text:
# Effects of non-thermal tails in Maxwellian electron distributions on synchrotron and Compton processes

## 1 Introduction

Power-law–like X-ray spectra extending up to soft $`\gamma `$-rays (hereafter X$`\gamma `$) are a common feature of the emission from accreting black holes in active galactic nuclei (AGNs) and black-hole binaries (BHBs). It is generally accepted that those spectra are produced by the Comptonization process taking place in optically thin accreting plasmas with typical electron energies of $`100`$ keV. Similar spectra are also observed from weakly-magnetized accreting neutron stars (e.g. Barret et al. 2000), where they can originate either in the accretion flow or in a boundary layer near the neutron star surface.

The X$`\gamma `$ spectra of Seyfert 1s and X-ray binaries in their hard (low) states can be modelled by Comptonization on thermal electrons, see e.g. Zdziarski (1999). However, both observations and theoretical considerations suggest that a small non-thermal component may be present in the electron distribution, giving rise to a weak non-thermal high-energy tail in the observed Comptonization spectrum. Such a tail above $`1`$ MeV has been detected in the spectrum of Cyg X-1 in the hard state by the COMPTEL detector aboard CGRO (McConnell et al. 1994, 2000a). In other BHBs and in Seyferts, the existing upper limits are compatible with the presence of weak non-thermal tails (e.g. Zdziarski et al. 1998; Gondek et al. 1996; Johnson et al. 1997).

A non-thermal tail in the electron distribution in compact X$`\gamma `$ sources can arise due to many different physical mechanisms. One such mechanism is acceleration in the process of dissipation of magnetic field (similarly to the case of the solar corona) in an optically-thin accretion flow or in active coronal regions above an optically thick disc. Then, the energy loss and subsequent thermalization of the accelerated electrons leads to a steady-state electron distribution consisting of a Maxwellian and a non-thermal tail. Apart from Coulomb interactions, an effective thermalizing process is synchrotron self-absorption (Ghisellini, Guilbert & Svensson 1988).

A strong support for the occurrence of non-thermal processes in black-hole accretion flows comes from observations of X$`\gamma `$ spectra of black-hole binaries in the soft state. In that case, Compton scattering by non-thermal electrons appears to be the dominant process (Poutanen & Coppi 1998; Gierliński et al. 1999). In the well-studied case of Cyg X-1, the power-law–like spectrum from this process becomes dominant at $`10`$ keV (Gierliński et al. 1999) and it extends to at least $`10`$ MeV (McConnell et al. 2000b).

An important issue concerns the source of soft photons undergoing Comptonization in the accreting plasma. A natural candidate is thermal emission from an optically thick accretion disc, either below the corona or outside the region of optically thin accretion flow. However, as the presence of magnetic fields of the order of that corresponding to equipartition with gas pressure (i.e. $`10^3`$ G in AGNs and $`10^7`$ G in BHBs) is expected, one has to take into account Comptonization of synchrotron photons (hereafter the CS process) when calculating the emergent spectra. In models of advection-dominated accretion flows, the CS process may be the main radiative process in the innermost part of the accretion flow, where most of the available gravitational energy is liberated (e.g. Narayan & Yi 1995).

Wardziński & Zdziarski (2000, hereafter WZ00) have investigated the efficiency of the thermal CS process under physical conditions relevant for accreting black holes and found that it can dominate only in low-luminosity AGNs and in luminous stellar X-ray sources with the hardest spectra. In the latter case, however, the presence of a correlation between the spectral index and the strength of X-ray reflection (Zdziarski, Lubiński & Smith 1999; Gilfanov, Revnivtsev & Churazov 1999; Lubiński & Zdziarski 2001) appears to imply that blackbody emission of cold media is the main source of seed photons for thermal Comptonization, and the CS process is negligible.

In this work, we extend our previous results to the case of a hybrid, quasi-thermal plasma and investigate the influence of a weak non-thermal component in the electron distribution on the efficiency of the CS process. In Section 2, we discuss the shape of the electron distribution function resulting from acceleration processes operating in accretion flows. The influence of non-thermal tails on the CS emission of the plasma is discussed in Section 3. We then investigate in detail synchrotron emission and Comptonization spectrum of hybrid electrons in Sections 4 and 5.1, respectively. This allows us to estimate the efficiency of the CS process in hybrid plasmas (Section 5.2), and those results are applied to astrophysical black-hole sources in Section 5.3. Finally, we summarize our results in Section 6.

## 2 Hybrid electron distributions

In accretion flows, hybrid electron distributions can be produced e.g. by stochastic acceleration of thermal electrons by interaction with plasma instabilities, e.g. with whistlers via gyroresonance (as applied to the case of Cyg X-1 by Li, Kusunose & Liang 1996) or with fast-mode Alfvén waves via transit-time damping (Li & Miller 1997). Such process can result from dissipation of magnetic field in an accretion disc corona (where magnetic fields are raised from the disc by buoyant forces, Galeev, Rosner & Vaiana 1979; Haardt, Maraschi & Ghisellini 1994) or in an optically thin accretion flow (Bisnovatyi-Kogan & Lovelace 1997), provided a significant fraction of the released energy generates turbulence.

Li & Miller (1997) considered parameters typical for a plasma in an accretion disc corona around a stellar black hole. They have shown, by solving the Fokker-Planck equation, that electrons can be effectively accelerated from the thermal background to relativistic energies, so that in the stationary state a non-thermal tail of a power-law–like shape develops above certain Lorentz gamma factor $`\gamma _{\mathrm{nth}}2`$ (which value is determined by the efficiency of thermalization processes, see Section 3,) and the power-law index, $`p`$, depends strongly on the plasma optical depth. The power-law component is cut off at a $`\gamma _\mathrm{f}10`$ due to radiative losses overcoming acceleration. Dermer, Miller & Li (1996) predict somewhat larger cut-off energies in the case of electron acceleration by whistlers and in a plasma with parameters relevant for AGNs.

Similar shape of steady-state electron distribution is also predicted from thermalization by synchrotron self-absorption and Coulomb interaction of an injected power-law electron distribution (produced by any acceleration mechanism), as shown by Ghisellini et al. (1988).

The exact shape of the non-thermal tail depends strongly on the plasma parameters and the (uncertain) mechanism of magnetic field dissipation and electron acceleration. We shall therefore consider a generic electron distribution function of the form qualitatively reproducing that predicted by detailed modelling of electron acceleration,

$`n_\mathrm{e}(\gamma )=\{\begin{array}{cc}n_\mathrm{e}^{\mathrm{th}}(\mathrm{\Theta },\gamma ),\hfill & \gamma \gamma _{\mathrm{nth}}\text{;}\hfill \\ n_\mathrm{e}^{\mathrm{th}}(\mathrm{\Theta },\gamma _{\mathrm{nth}})\left(\frac{\gamma }{\gamma _{\mathrm{nth}}}\right)^p\mathrm{exp}\left(\frac{\gamma \gamma _{\mathrm{nth}}}{\gamma _\mathrm{f}}\right),\hfill & \gamma >\gamma _{\mathrm{nth}}\text{.}\hfill \end{array}`$ (1)
Here $`\mathrm{\Theta }kT/m_\mathrm{e}c^2`$ is the dimensionless electron temperature and $`n_\mathrm{e}^{\mathrm{th}}(\mathrm{\Theta },\gamma )`$ is the relativistic Maxwellian normalized by the relation, $`n_\mathrm{e}=_1^{\mathrm{}}n_\mathrm{e}(\gamma )d\gamma `$, where $`n_\mathrm{e}`$ is the total electron density. This hybrid electron distribution consists thus of a thermal part below $`\gamma _{\mathrm{nth}}`$ and an e-folded power-law above it (see also Section 6).

We note that the cases $`p3`$ and $`p>3`$ are qualitatively different. In the former, most of the non-thermal synchrotron luminosity is produced by electrons with $`\gamma \gamma _\mathrm{f}`$ and the non-thermal synchrotron spectrum is hard. In the latter case, the luminosity is produced mostly by low-energy electrons, $`\gamma \gamma _{\mathrm{nth}}`$, and the spectrum is soft. Observations of Cyg X-1 (Gierliński et al. 1999; McConnell et al. 2000a) suggest $`p>3`$. In the following we assume $`p3`$.

The degree to which a hybrid distribution is non-thermal can be measured by the quantity, $`(\gamma _{\mathrm{nth}}1)/\mathrm{\Theta }`$, which shows how far in the Maxwellian tail the non-thermal power law starts. Another such quantity is the ratio of the energy densities in the non-thermal and thermal parts of the electron distribution, $`\delta `$, which, for mildly relativistic temperatures and $`\gamma _{\mathrm{nth}}1\mathrm{\Theta }`$ can be approximated as

$`\delta {\displaystyle \frac{_{\gamma _{\mathrm{nth}}}^{\mathrm{}}(\gamma 1)n_\mathrm{e}(\gamma )d\gamma }{\frac{6+15\mathrm{\Theta }}{4+5\mathrm{\Theta }}\mathrm{\Theta }n_\mathrm{e}}}`$ (2)
$`{\displaystyle \frac{4+5\mathrm{\Theta }}{6+15\mathrm{\Theta }}}{\displaystyle \frac{\left(\gamma _{\mathrm{nth}}^21\right)^{1/2}}{\mathrm{\Theta }^2K_2\left(\frac{1}{\mathrm{\Theta }}\right)}}\left({\displaystyle \frac{\gamma _{\mathrm{nth}}^3}{p2}}{\displaystyle \frac{\gamma _{\mathrm{nth}}^2}{p1}}\right)\mathrm{e}^{\frac{\gamma _{\mathrm{nth}}}{\mathrm{\Theta }}},`$
where $`K_2`$ is a modified Bessel function and we employed an approximation for the thermal energy density from Gammie & Popham (1998). In the above derivation (as well as in all analytical approximations below) we neglected, for simplicity, the exponential cut-off in the distribution (1). Typical relations between $`\delta `$ and $`(\gamma _{\mathrm{nth}}1)/\mathrm{\Theta }`$ are shown in Fig. 1.

Since the exponential term in equation (2) varies fastest, we have an approximate dependence,

$$\gamma _{\mathrm{nth}}C\mathrm{\Theta }\mathrm{ln}\delta ,$$
(3)
where $`C1`$ is slowly increasing with increasing $`\mathrm{\Theta }`$ and $`\gamma _{\mathrm{nth}}`$. We then see that $`(\gamma _{\mathrm{nth}}1)/\mathrm{\Theta }`$ only weakly depends on $`\mathrm{\Theta }`$ and the shape of the tail for a given $`\delta `$, see Fig. 1.

## 3 General form of emission from a hybrid electron distribution

In accreting black holes, synchrotron emission from thermal plasmas optically thin to scattering is usually strongly self-absorbed up to a turnover frequency, $`\nu _\mathrm{t}`$, below which the plasma is optically thick to absorption,

$$\alpha _{\nu _\mathrm{t}}R=1,$$
(4)
where $`\alpha _\nu `$ is the absorption coefficient and $`R`$ is the characteristic size of the source. Hereafter, we assume that $`\nu _\mathrm{t}`$ is determined only by synchrotron absorption, i.e. other radiative processes are negligible at $`\nu \nu _\mathrm{t}`$. Typically, accretion flow models predict (assuming equipartition magnetic fields) $`\nu _\mathrm{t}10^{15}`$ Hz for BHBs and $`\nu _\mathrm{t}10^{12}`$ Hz for AGNs, see e.g. WZ00. This corresponds to the value of $`\nu _\mathrm{t}/\nu _\mathrm{c}`$ from tens to hundreds, where $`\nu _\mathrm{c}=eB/2\pi m_\mathrm{e}c`$ is the cyclotron frequency for the magnetic field of the strength $`B`$ and $`m_\mathrm{e}`$ and $`c`$ are, respectively, the rest mass of the electron and the velocity of light. The ratio, $`\nu _\mathrm{t}/\nu _\mathrm{c}`$, is larger in AGNs as a result of their magnetic field weaker than that in BHBs. The Lorentz factor, $`\gamma _\mathrm{t}`$, of the electrons that emit most of their synchrotron radiation at $`\nu _\mathrm{t}`$ is of the order of $`\gamma _\mathrm{t}(\nu _\mathrm{t}/\nu _\mathrm{c})^{1/2}`$. The value of $`\gamma _\mathrm{t}`$ is of the order of a few, again larger in the case of AGNs than of BHBs.

First, let us consider a homogenous cloud of thermal plasma. The CS spectrum below the thermal turnover frequency, $`\nu _\mathrm{t}^{\mathrm{th}}`$, can be approximated as a Rayleigh-Jeans spectrum. Above $`\nu _\mathrm{t}^{\mathrm{th}}`$, the synchrotron emission drops rapidly and the spectrum is often dominated by a power-law component resulting from Comptonization, which is cut off above an energy of $`kT`$, where $`T`$ is the electron temperature and $`k`$ is the Boltzmann constant.

Now, let us assume that the electron distribution is given by equation (1). If $`\gamma _{\mathrm{nth}}<\gamma _\mathrm{t}`$, the number of electrons with $`\gamma \gamma _\mathrm{t}`$ will increase (compared to the thermal case) as well as the absorption coefficient and in consequence the new turnover frequency, $`\nu _\mathrm{t}^{\mathrm{nth}}`$, will be larger. Also the shape of the synchrotron spectrum is modified in this case. At low frequencies both emission and absorption are dominated by thermal electrons so that the optically thick spectrum is still the Rayleigh-Jeans one. However, at higher frequencies the emission of non-thermal electrons starts to dominate, while absorption is still dominated by thermal electrons \[which is a consequence of the thermal distribution falling off much more steeply than the non-thermal one, see equation (5)\], and the optically-thick spectrum quickly increases. This leads to the regime in which both the non-thermal emission and absorption dominate, in which case the spectrum is proportional to the source function of non-thermal electrons. Finally, the spectrum becomes optically thin ($`\nu ^{(p1)/2}`$), turning over above $`\nu _\mathrm{t}^{\mathrm{nth}}`$.

This behaviour is shown in Fig. 2 where the resulting synchrotron spectra are calculated for two different energy contents of the non-thermal tail and, as a reference, for the purely thermal case. Note that the non-thermal tail has to be sufficiently strong for the optically thick spectrum to reach the non-thermal source function before it turns over.

Both optically thick and optically thin parts of the synchrotron spectrum will undergo Comptonization. This process will thus be modified with respect to the thermal case due to two effects. Firstly, Compton scattering itself will involve both thermal electrons as well as non-thermal ones. The effect of non-thermal electrons on Comptonization will be most prominent at energies $`h\nu kT`$ (where $`h`$ is the Planck constant), where a high-energy tail (above the thermal cut-off) will develop. Secondly, the seed-photon flux will be higher, so luminosity in the thermal Comptonization spectrum will also increase. Since $`\gamma _\mathrm{t}`$ in optically-thin plasmas around black holes is expected to be well above the mean energy of the thermal electrons, even a non-thermal tail carrying a tiny fraction of the total electron energy, and thus only weakly modifying the shape of the Comptonization spectrum, can increase the turnover frequency and thus the luminosity produced by the CS process by a large factor.

Then a question arises whether the presence of strong synchrotron self-absorption, acting as a thermalizing mechanism for electrons of $`\gamma \gamma _\mathrm{t}`$ (i.e. in the regime of optically thick emission), will allow for an electron distribution deviating substantially from a pure Maxwellian for $`\gamma \gamma _\mathrm{t}`$. This issue was addressed by Ghisellini et al. (1988), who have shown that an initial power-law distribution of injected electrons will develop (after the thermalization time scale equal to a few synchrotron cooling timescales, $`t_\mathrm{c}^{\mathrm{syn}}`$, defined for the mean energy of the electrons radiating in the optically thick range of the spectrum) into a Maxwellian stationary distribution for $`\gamma <\gamma _\mathrm{t}`$, provided the synchrotron cooling dominates and $`t_\mathrm{c}^{\mathrm{syn}}t_{\mathrm{esc}}`$ (where $`t_{\mathrm{esc}}`$ is the escape time from the source). Assuming equipartition of magnetic field pressure with gas pressure, we find the latter criterion is always fulfilled in accretion flows onto black holes, independent of the form of the injected distribution.

However, when the inverse Compton process is dominant as a cooling mechanism, both $`\gamma _{\mathrm{nth}}`$ and $`kT`$ diminish and a power-law tail starting below $`\gamma _\mathrm{t}`$ is formed, see Ghisellini & Svensson (1990). A sufficient condition for the dominance of Compton cooling is $`\alpha <1`$, which is typically fulfilled in astrophysical accreting black holes. Even when $`\alpha >1`$, Comptonization of external soft photons (e.g. from a cold accretion disk) can dominate electron cooling (see, e.g. Ghisellini et al. 1998). Thus, $`\gamma _{\mathrm{nth}}<\gamma _\mathrm{t}`$ is a likely condition in those sources.

## 4 Synchrotron turnover frequency

Let $`\overline{\eta }_\nu (\gamma )`$ be the synchrotron emission coefficient of a single electron averaged both over the electron velocity direction and the direction of emission. Then, the absorption coefficient for an isotropic distribution of electrons in a chaotic magnetic field is (e.g. Ghisellini & Svensson 1991) is

$$\alpha _\nu =\frac{1}{2m_\mathrm{e}\nu ^2}\underset{1}{\overset{\mathrm{}}{}}\gamma \left(\gamma ^21\right)^{\frac{1}{2}}\overline{\eta }_\nu (\gamma )\frac{\mathrm{d}}{\mathrm{d}\gamma }\left[\frac{n_\mathrm{e}(\gamma )}{\gamma (\gamma ^21)^{\frac{1}{2}}}\right]d\gamma .$$
(5)
We assume a homogenous spherical source of the radius, $`R`$, and the Thomson optical depth, $`\tau _\mathrm{T}`$, where $`\tau _\mathrm{T}n_\mathrm{e}\sigma _\mathrm{T}R`$ and $`\sigma _\mathrm{T}`$ is the Thomson cross section. Then, we can numerically calculate the absorption coefficient with equation (5) for the distribution of equation (1) and then solve equation (4) for the turnover frequency, $`\nu _\mathrm{t}^{\mathrm{nth}}`$. This can be compared with the thermal turnover frequency, which for plasma parameters typical for accreting black holes can be approximated as (WZ00)

$$\nu _\mathrm{t}^{\mathrm{th}}7\times 10^3\mathrm{\Theta }^{0.95}\tau _\mathrm{T}^{0.05}\nu _\mathrm{c}^{0.91}\mathrm{Hz}.$$
(6)
On the other hand, we can obtain an analytical estimate, $`\nu _\mathrm{t}^{\mathrm{pl}}`$, of the turnover frequency for the hybrid electron distribution by considering a purely power-law electron distribution normalized to match the Maxwellian at $`\gamma _{\mathrm{nth}}`$. For the synchrotron absorption coefficient for power-law electrons (Rybicki & Lightman 1979),

$$\nu _\mathrm{t}^{\mathrm{pl}}=3^{\frac{1+p}{4+p}}2^{\frac{6}{4+p}}\pi ^{\frac{1}{4+p}}\nu _\mathrm{c}^{\frac{2+p}{4+p}}\left[G_1Rr_\mathrm{e}cn_\mathrm{e}(\gamma _{\mathrm{nth}})\gamma _{\mathrm{nth}}^p\right]^{\frac{2}{4+p}},$$
(7)
which, with equation (2), leads to

$`\nu _\mathrm{t}^{\mathrm{pl}}3^{\frac{1+p}{4+p}}2^{\frac{6}{4+p}}\pi ^{\frac{1}{4+p}}\nu _\mathrm{c}^{\frac{2+p}{4+p}}\left({\displaystyle \frac{\gamma _{\mathrm{nth}}}{p2}}{\displaystyle \frac{1}{p1}}\right)^{\frac{2}{4+p}}`$
$`\times \left({\displaystyle \frac{6+15\mathrm{\Theta }}{4+5\mathrm{\Theta }}}{\displaystyle \frac{G_1r_\mathrm{e}c\tau _\mathrm{T}\delta \mathrm{\Theta }\gamma _{\mathrm{nth}}^{p1}}{\sigma _\mathrm{T}}}\right)^{\frac{2}{4+p}},`$ (8)
where

$$G_1=\frac{\mathrm{\Gamma }(\frac{6+p}{4})\mathrm{\Gamma }(\frac{2+3p}{12})\mathrm{\Gamma }(\frac{22+3p}{12})}{\mathrm{\Gamma }(\frac{8+p}{4})}1,$$
(9)
$`r_\mathrm{e}`$ is the classical electron radius and $`\mathrm{\Gamma }`$ is Euler’s gamma function.

Provided $`p3`$ (so that radiation from electrons with $`\gamma >\gamma _\mathrm{f}`$ is negligible) and emission at $`\nu _\mathrm{t}^{\mathrm{nth}}`$ is dominated by non-thermal electrons (i.e. $`\gamma _\mathrm{f}>\gamma _\mathrm{t}>\gamma _{\mathrm{nth}}`$, which condition can be checked a posteriori by comparing the source function of power-law electrons to the Rayleigh-Jeans intensity), $`\nu _\mathrm{t}^{\mathrm{nth}}\nu _\mathrm{t}^{\mathrm{pl}}`$. On the other hand, $`\nu _\mathrm{t}^{\mathrm{pl}}<\nu _\mathrm{t}^{\mathrm{th}}`$ when emission from thermal electrons is significant at $`\nu _\mathrm{t}^{\mathrm{nth}}`$ and, in general,

$$\nu _\mathrm{t}^{\mathrm{nth}}\mathrm{max}(\nu _\mathrm{t}^{\mathrm{th}},\nu _\mathrm{t}^{\mathrm{pl}}).$$
(10)
The accuracy of this approximation is typically $`10`$ per cent (provided $`\gamma _\mathrm{t}<\gamma _\mathrm{f}`$ and $`p3`$). See WZ00 for a discussion of approximations of $`\nu _\mathrm{t}^{\mathrm{th}}`$.

We can obtain the dependences of $`\nu _\mathrm{t}^{\mathrm{nth}}/\nu _\mathrm{t}^{\mathrm{th}}`$ on the plasma parameters by considering the ratio $`\nu _\mathrm{t}^{\mathrm{pl}}/\nu _\mathrm{t}^{\mathrm{th}}`$ \[given by approximations (6) and (7)\]. First, at constant $`\delta `$, this ratio decreases with increasing $`\mathrm{\Theta }`$ since then $`\gamma _{\mathrm{nth}}`$ increases with $`\mathrm{\Theta }`$ \[see equation (3)\]. The dependence on $`\tau _\mathrm{T}`$ is stronger for $`\nu _\mathrm{t}^{\mathrm{pl}}`$ than for $`\nu _\mathrm{t}^{\mathrm{th}}`$ as a result of much steeper optically thin synchrotron spectrum in the latter case. This leads to $`\nu _\mathrm{t}^{\mathrm{pl}}/\nu _\mathrm{t}^{\mathrm{th}}\tau _\mathrm{T}^{2/(4+p)0.05}`$, i.e. the larger the optical depth, the larger the ratio $`\nu _\mathrm{t}^{\mathrm{pl}}/\nu _\mathrm{t}^{\mathrm{th}}`$. In a hot accretion flow, we generally expect $`\tau _\mathrm{T}`$ increasing and $`\mathrm{\Theta }`$ decreasing with increasing accretion rate, so the effect of a non-thermal tail (of a given form) on the value of the turnover frequency can be expected to be the largest for the most luminous sources.

Since $`\nu _\mathrm{t}^{\mathrm{pl}}/\nu _\mathrm{t}^{\mathrm{th}}\nu _\mathrm{c}^{(2+p)/(4+p)0.91}`$, the non-thermal increase of the turnover frequency grows with decreasing $`\nu _\mathrm{c}`$. This effect is due to the increase of $`\gamma _\mathrm{t}`$ with decreasing $`\nu _\mathrm{c}`$ \[as a result of increasing $`\nu _\mathrm{t}^{\mathrm{th}}/\nu _\mathrm{c}`$, see equation (6)\], which leads to a much larger ratio of the number of non-thermal to thermal electrons at $`\gamma _\mathrm{t}`$ for a given $`\gamma _{\mathrm{nth}}`$. Thus, the effect of a non-thermal tail with a given form will be much stronger in AGNs than in BHBs (due to much weaker magnetic field in the former objects).

Fig. 3 shows the values of $`\nu _\mathrm{t}^{\mathrm{nth}}/\nu _\mathrm{t}^{\mathrm{th}}`$ as a function of $`(\gamma _{\mathrm{nth}}1)/\mathrm{\Theta }`$ for three sets of thermal plasma parameters. In the first and second set, we assume parameters relevant for luminous BHBs, $`\tau _\mathrm{T}`$=1, $`R=3\times 10^7`$ cm, $`B=10^7`$ G, and $`\mathrm{\Theta }=0.1`$ and 0.2, respectively. In the third set, we assume parameters relevant to AGNs, $`\tau _\mathrm{T}`$=1, $`R=3\times 10^{14}`$ cm, $`B=10^3`$ G and $`\mathrm{\Theta }=0.1`$. For each set, we consider two different slopes of the non-thermal tail, $`p=3`$ and 4 and $`\gamma _\mathrm{f}=30`$. Additionally, we consider the cases of $`\gamma _\mathrm{f}=10`$ and 100 for $`p=3`$ (for $`p>3`$ the influence of the cut-off is much weaker). We see that even a weak non-thermal component, with $`(\gamma _{\mathrm{nth}}1)/\mathrm{\Theta }12`$ (which corresponds to only $`1`$ per cent of the total energy density of the electrons in the non-thermal tail) can lead to an increase of the turnover frequency by a factor of $`1.5`$–2 for BHBs and $`10`$ for AGNs.

The dependence of the relative increase of the turnover frequency on $`p`$ and $`\gamma _\mathrm{f}`$ is relatively weak. This reflects the fact that the turnover frequency depends mostly on the number of electrons at $`\gamma _\mathrm{t}`$ and since for small values of $`\delta `$, $`\gamma _\mathrm{t}`$ is only slightly larger than $`\gamma _{\mathrm{nth}}`$, $`n_\mathrm{e}(\gamma _\mathrm{t})`$ depends weakly on either $`p`$ or $`\gamma _\mathrm{f}`$. This suggests that the degree of the dependence on $`p`$ and $`\gamma _\mathrm{f}`$ should increase with increasing $`\nu _\mathrm{t}^{\mathrm{nth}}/\nu _\mathrm{t}^{\mathrm{th}}`$, in agreement with our results. It is important to note that weak dependence of $`\nu _\mathrm{t}^{\mathrm{nth}}/\nu _\mathrm{t}^{\mathrm{th}}`$ on the shape of the non-thermal tail makes our results weakly dependent on details of the acceleration mechanism.

## 5 Comptonization by a hybrid electron distribution

### 5.1 The general case

For Comptonization by a thermal plasma, we employ here a simple treatment of this process of Zdziarski (1985, 1986). In that approximation, the Comptonization spectrum above the frequency, $`\nu _{\mathrm{inj}}`$, at which soft seed photons are injected (which, in the case of the CS process, $`\nu _\mathrm{t}^{\mathrm{th}}`$), can be approximated as a sum of an e-folded power-law with energy index $`\alpha `$ and a Wien spectrum,

$$\frac{\mathrm{d}L_\mathrm{C}^{\mathrm{th}}}{\mathrm{d}x}\left[\left(\frac{x}{\mathrm{\Theta }}\right)^\alpha +\frac{\mathrm{\Gamma }(\alpha )P_{\mathrm{sc}}}{\mathrm{\Gamma }(2\alpha +3)}\left(\frac{x}{\mathrm{\Theta }}\right)^3\right]\mathrm{e}^{x/\mathrm{\Theta }},$$
(11)
where $`xh\nu /m_\mathrm{e}c^2`$ (hereafter all indices of $`x`$ have the same meaning as those of $`\nu `$) and $`P_{\mathrm{sc}}`$ is the volume-averaged scattering probability. Its integration yields the thermal-Compton luminosity, which, in the case of a spherical source, is

$$L_\mathrm{C}^{\mathrm{th}}(𝒞)4\pi R^2𝒞\mathrm{\Theta }\left[\mathrm{\Gamma }(1\alpha ,\frac{x_{\mathrm{inj}}}{\mathrm{\Theta }})+\frac{6\mathrm{\Gamma }(\alpha )P_{\mathrm{sc}}}{\mathrm{\Gamma }(2\alpha +3)}\right],$$
(12)
where $`𝒞`$ is a constant depending on the flux of the injected seed photons (we calculate $`𝒞`$ in Section 5.2 for the case of synchrotron seed photons).

The presence of non-thermal electrons modifies the Comptonization spectrum and since for the electron distributions we consider (i.e. $`\tau _\mathrm{T}1`$ and $`\delta 1`$) the Thomson optical depth for scattering off non-thermal electrons is $`1`$, the resulting spectrum can be approximated as a convolution of the thermal Comptonization spectrum with a spectrum resulting from single scattering of photons off non-thermal electrons. Therefore, for $`x\mathrm{\Theta }`$, the Comptonization spectrum becomes harder, while for $`x\mathrm{\Theta }`$, above the thermal cut-off, a power-law tail in the spectrum develops. As a result, the overall luminosity produced by Comptonization increases.

In calculations of the luminosity from Comptonization on non-thermal electrons, $`L_\mathrm{C}^{\mathrm{pl}}`$, where the Klein-Nishina effect has to be taken into account, we use the approximation for the rate of energy change, $`\mathrm{d}\gamma /\mathrm{d}t`$, of a single electron via Compton interaction with isotropic photons of energy density, $`U_{\mathrm{ph}}`$, and mean energy, $`x`$,

$$\frac{\mathrm{d}\gamma }{\mathrm{d}t}=\frac{4}{3}\frac{\sigma _\mathrm{T}\left(\gamma ^21\right)U_{\mathrm{ph}}}{m_ec}\left[1\frac{63}{10}\frac{\gamma x^2}{x}\right],$$
(13)
where the Klein-Nishina cross-section was approximated using the first-order correction to the Thomson-limit cross-section (Rybicki & Lightman 1979). We assume that the photons undergoing scattering off non-thermal electrons are those from the thermal Comptonization spectrum above $`x_{\mathrm{inj}}`$. Then

$$U_{\mathrm{ph}}(x)=\frac{3\left(\mathrm{d}L_\mathrm{C}^{\mathrm{th}}/\mathrm{d}x\right)\left(3/4+\tau _\mathrm{T}/5\right)}{4\pi cR^2},$$
(14)
where the factor $`\left(3/4+\tau _\mathrm{T}/5\right)`$ accounts for the change of the escape time due to scatterings in the source and is a matching formula between the optically thin case (where the escape time is $`3R/4c`$) and the optically thick one ($`\tau _\mathrm{T}R/5c`$, Sunyaev & Titarchuk 1980). We can further simplify the calculations assuming that the spectrum of photons undergoing non-thermal Comptonization is a pure power-law and then neglect scatterings in the Klein-Nishina limit (see below).

Now, for each $`\gamma `$ we calculate $`U_{\mathrm{ph}}`$, $`x`$ and $`x^2`$ as an integral over the power-law spectrum from $`x_{\mathrm{inj}}`$ up to the energy, $`x_{\mathrm{max}}`$, of the limit of the Klein-Nishina regime, which we obtain from the condition of $`163/(10\gamma x_{\mathrm{max}})=0`$, in which the numerical coefficient is chosen for consistency with equation (13).

Then the formula (13) is integrated over $`n_\mathrm{e}(\gamma )`$, from $`\gamma _{\mathrm{nth}}`$ to infinity, which leads to

$`L_\mathrm{C}^{\mathrm{pl}}(𝒞)={\displaystyle \frac{16\pi }{3}}\left({\displaystyle \frac{10}{63}}\right)^{1\alpha }{\displaystyle \frac{6+15\mathrm{\Theta }}{4+5\mathrm{\Theta }}}\left({\displaystyle \frac{3}{4}}+{\displaystyle \frac{\tau _\mathrm{T}}{5}}\right){\displaystyle \frac{R^2𝒞\delta \tau _\mathrm{T}\mathrm{\Theta }^{\alpha +1}}{(1\alpha )(2\alpha )}}`$
$`\times \left({\displaystyle \frac{\gamma _{\mathrm{nth}}}{p2}}{\displaystyle \frac{1}{p1}}\right)^1\left({\displaystyle \frac{\gamma _{\mathrm{nth}}^{\alpha +1}}{p\alpha 2}}{\displaystyle \frac{\gamma _{\mathrm{nth}}^{\alpha 1}}{p\alpha }}\right),`$ (15)
where equation (2) has been employed.

The total luminosity, $`L_\mathrm{C}^{\mathrm{nth}}=L_\mathrm{C}^{\mathrm{th}}+L_\mathrm{C}^{\mathrm{pl}}`$, as well as the ratio $`L_\mathrm{C}^{\mathrm{pl}}/L_\mathrm{C}^{\mathrm{th}}`$ can now be directly obtained from expressions (12) and (5.1). In the range of $`0.4\alpha 0.9`$, both the high-energy part of the power law component dominates the luminosity as well as the Wien component can still be neglected. Then, we obtain (see WZ00)

$`{\displaystyle \frac{L_\mathrm{C}^{\mathrm{pl}}}{L_\mathrm{C}^{\mathrm{th}}}}{\displaystyle \frac{4}{3}}\left({\displaystyle \frac{10}{63}}\right)^{1\alpha }{\displaystyle \frac{6+15\mathrm{\Theta }}{4+5\mathrm{\Theta }}}\left({\displaystyle \frac{3}{4}}+{\displaystyle \frac{\tau _\mathrm{T}}{5}}\right){\displaystyle \frac{\delta \tau _\mathrm{T}\mathrm{\Theta }^\alpha }{2\alpha }}`$
$`\times \left({\displaystyle \frac{\gamma _{\mathrm{nth}}}{p2}}{\displaystyle \frac{1}{p1}}\right)^1\left({\displaystyle \frac{\gamma _{\mathrm{nth}}^{\alpha +1}}{p\alpha 2}}{\displaystyle \frac{\gamma _{\mathrm{nth}}^{\alpha 1}}{p\alpha }}\right).`$ (16)
We find the accuracy of this approximation is typically $`30`$ per cent. Figure 4 shows that while the ratio $`L_\mathrm{C}^{\mathrm{pl}}/L_\mathrm{C}^{\mathrm{th}}`$ depends strongly on $`\delta `$, the dependence on the non-thermal tail shape is very weak.

In Fig. 5, we show as an example Comptonization spectra from a homogenous plasma cloud, produced with the Monte Carlo code of Gierliński (2000). We see that as long as $`\delta <0.1`$ the Comptonization spectrum at $`x\mathrm{\Theta }`$ is hardly modified by the presence of non-thermal electrons and the onset of the power-law tail is no higher than an order of magnitude below the peak in the $`xL(x)`$ spectrum.

From the observational point of view, an important feature measuring the deviation from purely thermal Comptonization spectrum is the fraction of the luminosity radiated in the high-energy tail beyond the thermal spectrum. The ratio $`L_\mathrm{C}^{\mathrm{pl}}/L_\mathrm{C}^{\mathrm{th}}`$ is larger than this fraction as it incorporates also the non-thermal luminosity at energies below the tail. This effect is stronger for smaller $`\delta `$, when a larger fraction of $`L_\mathrm{C}^{\mathrm{pl}}`$ is hidden below the thermal peak. For example, in the spectra in Fig. 5 corresponding to $`\delta =0.1`$ and $`\delta =0.01`$ the high-energy tail develops at $`300`$ keV and $`500`$ keV, respectively. Then $`L_\mathrm{C}^{\mathrm{pl}}/L_\mathrm{C}^{\mathrm{th}}`$ is larger than he fraction of the luminosity in the tail by factors 3 and 5, respectively. Thus, $`L_\mathrm{C}^{\mathrm{pl}}/L_\mathrm{C}^{\mathrm{th}}`$ provides only an upper limit to the fraction of the luminosity radiated in the high-energy tail.

### 5.2 Comptonization of synchrotron photons

As we have seen above, the presence of a non-thermal component in the thermal distribution of electrons in a plasma where synchrotron radiation is produced and then Comptonized leads to an increase of the luminosity available from this process in two ways. First, it increases the number of soft photons (both from optically thick emission up to $`x_\mathrm{t}^{\mathrm{nth}}`$ and from optically thin emission above $`x_\mathrm{t}^{\mathrm{nth}}`$) and second, it increases the efficiency of Comptonization, the former effect being significantly stronger.

In the case of purely thermal synchrotron Comptonization, the normalization of the Comptonization spectrum can be assumed (Zdziarski 1985, 1986) to be proportional to the Rayleigh-Jeans flux at the turnover frequency, so that the constant $`𝒞`$ in the expression (12), which we now denote as $`𝒞^{\mathrm{th}}`$, reads,

$$𝒞^{\mathrm{th}}=\frac{2\pi m_\mathrm{e}c^3}{\lambda _\mathrm{C}^3}\phi \mathrm{\Theta }^{1\alpha }\left(x_\mathrm{t}^{\mathrm{th}}\right)^{2+\alpha },$$
(17)
where

$$\phi (\mathrm{\Theta })\frac{1+(2\mathrm{\Theta })^2}{1+10(2\mathrm{\Theta })^2}$$
(18)
is a heuristic formula obtained by Zdziarski (1985) to match Monte Carlo simulation, and $`\lambda _\mathrm{C}`$ is the Compton wavelength. The total luminosity from the CS process is then

$$L_{\mathrm{CS}}^{\mathrm{th}}=L_\mathrm{S}^{\mathrm{th}}+L_\mathrm{C}^{\mathrm{th}}(𝒞^{\mathrm{th}}),$$
(19)
where $`L_\mathrm{S}^{\mathrm{th}}`$ is the luminosity in the thermal synchrotron spectrum.

If thermal and hybrid synchrotron photons experience the same average number of scatterings off thermal electrons before reaching $`x\mathrm{\Theta }`$ (i.e. when $`x_\mathrm{t}^{\mathrm{nth}}/x_\mathrm{t}^{\mathrm{th}}\mathrm{\Theta }/x_\mathrm{t}^{\mathrm{th}}`$ and the synchrotron spectrum is sufficiently narrow, which corresponds to $`p>3`$ and which we assume hereafter in this section), the ratio of the luminosity from Comptonization off thermal electrons in the hybrid case to that in the thermal case will equal the corresponding ratio of the synchrotron luminosities. Then,

$$\frac{𝒞^{\mathrm{nth}}}{𝒞^{\mathrm{th}}}=\frac{L_\mathrm{S}^{\mathrm{nth}}}{L_\mathrm{S}^{\mathrm{th}}},$$
(20)
where $`L_\mathrm{S}^{\mathrm{nth}}`$ is the hybrid synchrotron luminosity and $`𝒞^{\mathrm{nth}}`$ is the normalization constant for thermal Comptonization spectrum of hybrid synchrotron photons.

To calculate the ratio (20) we assume the optically thick hybrid synchrotron emission below $`x_\mathrm{t}^{\mathrm{nth}}`$ to be $`S^{\mathrm{pl}}(x)`$, the source function of power-law electrons,

$$S^{\mathrm{pl}}(x)=\frac{1}{\mathrm{2\hspace{0.17em}3}^{1/2}}\frac{G_2}{G_1}\frac{m_\mathrm{e}c^3}{\lambda _\mathrm{C}^3}x_\mathrm{c}^{\frac{1}{2}}x^{\frac{5}{2}},$$
(21)
where

$$G_2=\frac{\mathrm{\Gamma }(\frac{5+p}{4})\mathrm{\Gamma }(\frac{3p+19}{12})\mathrm{\Gamma }(\frac{3p1}{12})}{\mathrm{\Gamma }(\frac{7+p}{4})}1.$$
(22)
The emission above $`\nu _\mathrm{t}^{\mathrm{nth}}`$ is $`\nu ^{(p1)/2}`$ and then we have

$$\frac{𝒞^{\mathrm{nth}}}{𝒞^{\mathrm{th}}}=\frac{3^{\frac{1}{2}}(p+4)G_2\left(x_\mathrm{t}^{\mathrm{nth}}/x_\mathrm{c}\right)^{7/2}}{14(p3)G_1\mathrm{\Theta }\left(x_\mathrm{t}^{\mathrm{th}}/x_\mathrm{c}\right)^3}.$$
(23)
The above expressions were derived using the emission and absorption coefficients for power-law electrons (e.g. Rybicki & Lightman 1979). Since we neglected the cut-off in the electron distribution, the above relation holds for $`x_\mathrm{t}^{\mathrm{nth}}/x_\mathrm{c}\gamma _\mathrm{f}^2`$, i.e. when the contribution from electrons with $`\gamma >\gamma _\mathrm{f}`$ to the emission at $`x_\mathrm{t}^{\mathrm{nth}}`$ can be neglected.

The total luminosity of the CS process in the hybrid case is then

$$L_{\mathrm{CS}}^{\mathrm{nth}}=L_\mathrm{S}^{\mathrm{nth}}+L_\mathrm{C}^{\mathrm{th}}(𝒞^{\mathrm{nth}})+L_\mathrm{C}^{\mathrm{pl}}(𝒞^{\mathrm{nth}}).$$
(24)
When $`\alpha <1`$, we can neglect $`L_\mathrm{S}^{\mathrm{th}}`$ and $`L_\mathrm{S}^{\mathrm{nth}}`$ in equations (19) and (24) respectively, since most of the luminosity is then produced at much higher energies, $`x\mathrm{\Theta }`$. Then, the ratio $`L_{\mathrm{CS}}^{\mathrm{nth}}/L_{\mathrm{CS}}^{\mathrm{th}}`$ can be calculated as

$`{\displaystyle \frac{L_{\mathrm{CS}}^{\mathrm{nth}}}{L_{\mathrm{CS}}^{\mathrm{th}}}}={\displaystyle \frac{𝒞^{\mathrm{nth}}}{𝒞^{\mathrm{th}}}}\left[1+{\displaystyle \frac{L_\mathrm{C}^{\mathrm{pl}}\left(𝒞^{\mathrm{th}}\right)}{L_\mathrm{C}^{\mathrm{th}}\left(𝒞^{\mathrm{th}}\right)}}\right].`$ (25)
We can neglect the last term in the brackets if $`\delta `$ is sufficiently small, in which case $`L_\mathrm{C}^{\mathrm{pl}}\left(𝒞^{\mathrm{th}}\right)L_\mathrm{C}^{\mathrm{th}}\left(𝒞^{\mathrm{th}}\right)`$. We then have

$`{\displaystyle \frac{L_{\mathrm{CS}}^{\mathrm{nth}}}{L_{\mathrm{CS}}^{\mathrm{th}}}}{\displaystyle \frac{3^{\frac{1}{2}}(p+4)\left(x_\mathrm{t}^{\mathrm{nth}}/x_\mathrm{c}\right)^{7/2}}{14(p3)\mathrm{\Theta }\left(x_\mathrm{t}^{\mathrm{th}}/x_\mathrm{c}\right)^3}}.`$ (26)
Note that in calculations of $`𝒞^{\mathrm{nth}}`$ we used the source function for power-law electrons. Thus, this approximation does not have the correct thermal limit when $`\gamma _{\mathrm{nth}}\gamma _\mathrm{t}`$, i.e. when $`x_\mathrm{t}^{\mathrm{nth}}x_\mathrm{t}^{\mathrm{th}}`$ (see Fig. 2b, where a synchrotron spectrum for such a case is shown). Note also that for harder power-law electron distributions ($`p3`$), for which the approximation (20) does not hold, the amplification of CS luminosity will be even larger than that predicted above.

### 5.3 Applications to NGC 4151, GX 339–4 and Cyg X-1

We now calculate the ratio $`L_{\mathrm{CS}}^{\mathrm{nth}}/L_{\mathrm{CS}}^{\mathrm{th}}`$ as a function of the energy content of the electron power-law tail in three different cases. They correspond to the models of thermal CS emission from a hot accretion disc with magnetic field in equipartition with the gas internal energy (dominated by hot ions) of WZ00 for the Seyfert galaxy NGC 4151 and two BHBs in their hard states, GX 339–4 and Cyg X-1. The CS luminosities predicted by the thermal model (fig. 8a in WZ00) were orders of magnitude below the ones observed from NGC 4151 and GX 339–4, and less by a factor of $`3`$ for Cyg X-1.

For NGC 4151, we assume (see WZ00 for details and references) $`R=1.2\times 10^{14}`$ cm, $`B=3.8\times 10^3`$ G, $`\alpha =0.85`$, $`\mathrm{\Theta }=0.1`$, $`p=4`$, $`\gamma _\mathrm{f}=100`$. For GX 339–4, the same values of $`\mathrm{\Theta }`$, $`p`$ and $`\gamma _\mathrm{f}`$ are adopted while $`R=8.8\times 10^6`$ cm, $`B=1.5\times 10^7`$ G, $`\alpha =0.75`$. In the case of Cyg X-1, $`R=3\times 10^7`$ cm, $`B=6.3\times 10^6`$ G, $`\alpha =0.6`$, and we use the electron distribution obtained in McConnell et al. (2000a) by fitting the spectra from the COMPTEL and OSSE detectors aboard CGRO, i.e. $`\mathrm{\Theta }=0.17`$, $`\gamma _{\mathrm{nth}}=2.12`$, $`p=4.5`$, at assumed $`\gamma _\mathrm{f}=10^3`$, which corresponds to a rather large value of $`\delta =0.19`$.

The results are shown in Fig. 6. In the case of NGC 4151 and GX 339–4, a non-thermal tail in the electron distribution with $`\delta 10^3`$ amplifies the power produced by the CS process sufficiently to reproduce the luminosities of those sources. Therefore, unlike the results obtained by WZ00 from investigating the purely thermal CS process, in the case of electron distributions with relatively weak non-thermal tails the CS process appears to be capable of producing the observed X-ray spectra of accreting black holes, independent of the Eddington ratio or the black hole mass.

On the other hand, the estimated thermal CS luminosity in the case of Cyg X-1 in the hot disc model is comparable to that observed from the source. Then, we find that the non-thermal tail observed by McConnell et al. (2000a) would produce the CS luminosity $`500`$ times higher than that actually observed (see Fig. 2a, where thermal and hybrid synchrotron spectra are plotted for plasma parameters assumed in this model). The requirement that $`L_{\mathrm{CS}}^{\mathrm{nth}}`$ cannot exceed the source luminosity constrains the magnetic field strength in the model to be $`B2.3\times 10^5`$ G, which is $`30`$ times smaller than the equipartition value estimated by WZ00.

A similar constraint for Cyg X-1 can be obtained in a model of active coronal regions above the disc with magnetic field dissipation (see WZ00 for details of the corresponding model of CS emission). Assuming $`N=10`$ active regions of radius $`R_\mathrm{b}=3\times 10^6`$ cm and height $`H_\mathrm{b}=3\times 10^5`$ cm, we obtain $`B3.4\times 10^6`$ G, much less than the field strength predicted assuming dissipation of magnetic field at a fraction of the Alfvén speed, $`B_{\mathrm{diss}}=5\times 10^7`$ G. Though varying $`N`$, $`R_\mathrm{b}`$, $`H_\mathrm{b}`$ would change the upper limit on the field strength as well as the value of $`B_{\mathrm{diss}}`$, the relation between these quantities would remain similar. This suggests that the mechanism of dissipation of disc field in small active regions above the disc cannot be responsible for the energy release in the corona. This stems from the fact that $`B_{\mathrm{diss}}`$ is the minimum field strength necessary for the magnetic field to dissipate enough energy so as to produce the observed luminosity. Therefore, either the model of active regions above a cold disc is ruled out for the hard state of Cyg X-1 or another mechanism of energy release must operate. We point out that if the CS process were to be negligible (as suggested by observations, see Section 6), the magnetic field would have to be even weaker, which would strengthen our conclusion even more. On the other hand, a weak magnetic field does not present analogous difficulties in the hot accretion disc model as the disc heating is directly by gravity.

We then check how sensitive the above results are to details of the observed spectra of Cyg X-1. First, the spectral index of Cyg X-1 varies, and the value $`\alpha =0.6`$ assumed above is the hardest one and spectra with $`\alpha =0.7`$ have been also observed. Second, the relative normalization between the COMPTEL and OSSE spectra remains relatively uncertain, and it may be possible that the former is less by $`2`$ (McConnell et al. 2000b) than the one used in the calculations above. This would then diminish the energy content in the electron non-thermal tail twice, which corresponds to the new value of $`\gamma _{\mathrm{nth}}=2.31`$. Thus, we have repeated our calculations for $`\alpha =0.7`$ and $`\gamma _{\mathrm{nth}}=2.31`$. The resulting constraints on the magnetic field strength in the hot disc and active regions models are then $`B5.8\times 10^5`$ G and $`B6.8\times 10^6`$ G, respectively. These values are only about twice the previous upper limits, which does not affect our conclusions above.

We note that another constraint on the CS emission may be provided by extrapolating the X-ray power law spectrum to lower energies (but $`h\nu _\mathrm{t}`$), e.g. to the UV or optical ranges. That emission should not, in any case, exceed that observed. In the case of Cyg X-1, we have extrapolated the observed X-ray power-law spectrum of Cyg X-1 (Gierliński et al. 1997) to $`\nu 10^{15}`$ Hz and have found it is more than 2 orders of magnitude below the UV spectrum of the companion supergiant as observed by IUE (Treves et al. 1980). Thus, those data are compatible with the CS origin of the X-ray spectrum in Cyg X-1.

On the other hand, GX 339–4 is a low-mass X-ray binary with a low upper limit on the emission of the companion star (e.g. Zdziarski et al. 1998). Then, simultaneous (as the object is highly variable, e.g. Corbet et al. 1987) observations in X-rays and at longer wavelengths can put constraints on the CS model in this case (e.g. Fabian et al. 1982).

## 6 Conclusions and discussion

We have studied the influence of non-thermal tails in the electron distribution on the Comptonized synchrotron emission from plasmas of parameters typical for accreting black holes. In those plasmas, most of the self-absorbed thermal synchrotron emission is produced by electrons far in the tail of the Maxwellian distribution. As a result, even a weak non-thermal component far beyond the thermal peak can greatly enhance the synchrotron emission. This emission provides seed photons for Comptonization, which is also correspondingly enhanced. However, the shape the Comptonization spectrum remains only weakly modified, with the main signature being a relatively weak high-energy tail.

We have shown that this effect becomes stronger with the decreasing plasma temperature and increasing optical depth, which typically correspond to an increasing Eddington ratio. It is also much stronger in the case of AGNs than of accreting stellar black holes, as a result of weaker magnetic field in the former objects. The CS luminosity can be amplified by a factor $`10^5`$ and $`10^3`$, respectively, in the case of only $`1`$ per cent of the electron energy in the non-thermal tail.

Then, the Comptonized synchrotron emission becomes energetically capable of producing the observed luminosities of accreting black holes independently of their Eddington ratio or the black-hole mass. This conclusion is qualitatively different from that obtained in the case of purely thermal plasma.

It is then possible to obtain upper limits on the magnetic field strength in optically thin plasmas around black holes by constraining the parameters of the non-thermal component in the electron distribution from observations of the high-energy tail in Comptonization spectra. Based on the CGRO data for Cyg X-1 in the hard state, we have obtained the upper limits on the field strength at least an order of magnitude below both the value corresponding to equipartition (in the model of a two-temperature accretion disc) and the minimum value required by dissipation of magnetic fields (in the model of active coronal regions). Therefore, the latter model appears to be ruled out for the hard state of Cyg X-1.

We note that the presence of a strong correlation between reflection strength and the X-ray spectral index, observed in Cyg X-1 (Gilfanov et al. 1999), suggests that the CS emission is a negligible source of photons for Comptonization and thus the magnetic field strength must remain significantly below the upper limits derived above. The same correlation seen in other objects (Zdziarski et al. 1999; Gilfanov et al. 2000) implies that either the high-energy tails or the magnetic fields are weak enough for the CS process not to dominate the energy output.

In this work, we have assumed, for the sake of simplicity and compatibility with studies of, e.g., McConnell et al. (2000a), the electron distribution of the non-thermal electron tails to be $`\gamma ^p`$ at any value of $`\gamma `$. However, when the power-law tail extends down to non-relativistic energies, $`\gamma _{\mathrm{nth}}1`$, the distributions expected from acceleration processes are power laws in either the kinetic energy or the momentum, e.g. $`(\gamma 1)^p`$ for the former. Clearly, the smaller $`\gamma _{\mathrm{nth}}`$ (i.e. the larger $`\delta `$), the larger the differences between those distributions. Then, the actual shift of the turnover frequency at a given value of $`\gamma _{\mathrm{nth}}`$ would be somehow smaller. Also, the relation between $`\gamma _{\mathrm{nth}}`$ and $`\delta `$ would change. Nevertheless, the above effects would modify our results only quantitatively without affecting our conclusions.

The influence of a weak, non-thermal, component in the electron distribution on radiation spectra has also been independently considered by Özel, Psaltis & Narayan (2000). Their study has been devoted to the case of advection-dominated accretion flows (ADAF), and their presented spectra are integrated over all radii of the flow. In contrast to our results, they do not note any shift in the turnover frequency at the peak of the synchrotron spectrum. This appears to result from their study being constrained to the ADAF model at low accretion rates, in which case electrons in the innermost parts of the accretion flow reach rather high temperatures (at which even thermal electrons reach relatively high Lorentz factors). Those temperatures are significantly higher than those we consider based on observational data from luminous accreting black holes. On the other hand, Özel et al. (2000) find a significant excess non-thermal emission at frequencies well below the peak of the integrated synchrotron spectrum. This excess results from radiation emitted at large radii, where the electron temperature is low, and where the synchrotron emission (both optically thick and optically thin) around the local turnover frequency can be strongly amplified by the presence of an electron tail, as illustrated in Fig. 2 above.

## ACKNOWLEDGMENTS

This research has been supported in part by the Foundation for Polish Science and KBN grants 2P03D01619 and 2P03D00624. We thank Marek Gierliński for providing us with his Monte-Carlo Comptonization code and for help with its use. We are also grateful to Andrei Beloborodov, Joanna Mikołajewska and Juri Poutanen for valuable discussions.


format: text

text:
# Operator method for solution of the Schrödinger equation with the rational potential

## I Introduction

In this paper we generalize the operator method for solving the Schrödinger equation (SE) with a nonpolynomial structure of the Hamiltonian. As an example of such a problem we consider one-dimensional Hamiltonian with the rational potential

$$V_L(x)=\frac{x^2}{2}+\frac{\lambda x^{2L}}{1+gx^2}.$$
(1)
The problem chosen is of interest in different areas of physics, i.e. laser physics l15  and non-linear quantum field theory l16 . It can be also exploited for verifying new nonperturbative methods, which allow us to find approximate solution of the Schrödinger equation in the whole range of Hamiltonian parameters. Such methods can be used for treating many physical problems of current interest. The test problem chosen is very suitable for these purposes because the SE

$$\left[\frac{p^2}{2}+V_L(x)E\right]|\psi =0,$$
(2)
can be solved exactly if some algebraic relations between parameters $`g`$ and $`\lambda `$ are assumed l11 ; l17 ; l12 . The exact solutions allow to find out the accuracy of approximate methods.

The Schrödinger equation with the potentials $`V_1(x)`$ and $`V_2(x)`$ has been treated by different approaches based on finite-difference schemes, variational and asymptotic methodsl10 ; l3 ; l2 ; l7 ; l8 ; l9 ; l1 ; l4 . To the best of our knowledge there is only one work l4 , where the case $`L>2`$ was studied in the wide range of the parameter $`\alpha =g/\lambda `$ by very efficient method of finite-differences. The case $`L=+3`$ has been studied for $`g1,\lambda 1`$ by using the Padé-approximant and hypervirial methodsl20 . Three-dimensional generalization of Eq.(2) can be found in Ref. l19, .

In the present paper the SE Eq.(2) is solved by the operator method (OM), developed by Komarov and Feranchukl5 . The OM has been successfully applied to different problemsl18 ; thu ; l21 ; l14 , but the most problems considered had a polynomial structure of the Hamiltonian. The polynomial structure allows one to calculate matrix elements of the SE equation analytically in the Fock basisl1 . But this simple analytical scheme can not be directly applied to the nonpolynomial equation Eq.(2). We can still perform the OM calculationsl18 , but some extra numerical scheme must be exploited to calculate matrix elements. This problem can be resolved by the use of the nonunitary transformation $`C`$

$$|\psi =C|\phi ,C=(1+gx^2),$$
(3)
which leads to the following expression for Eq.(2)

$$\left[\frac{p^2}{2}+\frac{x^2}{2}+g\left(\frac{p^2x^2}{2}+\frac{x^4}{2}\right)+\lambda x^{2L}(1+gx^2)E\right]|\phi =0.$$
(4)
So we transformed the nonpolynomial operator equation Eq.(2) to the polynomial one Eq.(4), where the operator at the left side is nonhermitian with respect to the scalar product

$$\varphi |\psi =dx\varphi ^{}(x)\psi (x).$$
(5)
The idea of using nonunitary transformations was recently suggested in Ref. ander, . It was also developed in Ref. lee,  as nonhermitian technique of canonical transformations. The basic idea is that the commutation relation $`[x,p]=i`$ can be preserved by any unitary transformation, but the same is also true if we make use of the similarity transformation $`X=CxC^1`$, $`P=CpC^1`$. Below we show that the polynomial operator equation Eq.(4) is very easy to handle within the framework of the operator method in order to find both the exact numerical and approximate analytical solutions. The OM allows us to treat the problem in the wide range of Hamiltonian parameters $`g`$ and $`\lambda `$, where the SE Eq.(2) has a discrete energy spectrum, i.e. when $`\mathrm{}<L+1`$, $`|\lambda |<+\mathrm{}`$, $`g0`$; $`L>+2`$, $`\lambda 0`$, $`g0`$; $`L=+2`$, $`\lambda >g/2`$, $`g0`$. Below we present results of our nonperturbative calcutations for four cases $`1L4`$.

## II Iterative scheme of the operator method

In this section we apply the operator method for solving the equation Eq.(4) using the scheme developed in Ref. l5, . If we introduce the annihilation operator $`a(\omega )`$, the creation operator $`a^+(\omega )`$ and the excitation number operator $`n(\omega )=a^+(\omega )a(\omega )`$

$$a^+(\omega )=\frac{1}{\sqrt{2\omega }}(\omega xip),a(\omega )=\frac{1}{\sqrt{2\omega }}(\omega x+ip),$$
(6)
which satisfy the commutation relations

$$[a(\omega ),a^+(\omega )]=1,[a(\omega ),n(\omega )]=a(\omega ),[a^+(\omega ),n(\omega )]=a^+(\omega ),$$
(7)
one have for Eq.(4)

$`(^{}`$ $``$ $`E^{\prime \prime })|\phi =0,`$ (8)
$`^{}`$ $`=`$ $`\left({\displaystyle \frac{1}{4\omega }}+{\displaystyle \frac{\omega }{4}}\right)(2n+1)+{\displaystyle \frac{3g}{8\omega ^2}}(2n^2+2n+1)+{\displaystyle \frac{g}{8}}(2n^2+2n1)+{\displaystyle \frac{g}{8}}({\displaystyle \frac{1}{\omega ^2}}1)`$ (9)
$`\times `$ $`(a^{+4}+a^4)+\left({\displaystyle \frac{1}{4\omega }}{\displaystyle \frac{\omega }{4}}+{\displaystyle \frac{3g}{4\omega ^2}}+{\displaystyle \frac{g}{2\omega ^2}}n\right)(a^{+2}+a^2){\displaystyle \frac{g}{\omega ^2}}a^{+2}`$
$`+`$ $`{\displaystyle \frac{g}{2}}(a^{+2}a^2)+\lambda x^{2L},`$
$`^{\prime \prime }`$ $`=`$ $`1+{\displaystyle \frac{g}{2\omega }}(a^{+2}+a^2+2n+1),`$ (10)
where $`x^{2L}`$ is written below for the four particular cases $`1L4`$

$`x^2`$ $`=`$ $`{\displaystyle \frac{1}{2\omega }}(2n+1)+{\displaystyle \frac{1}{2\omega }}(a^{+2}+a^2),`$ (11)
$`x^4`$ $`=`$ $`{\displaystyle \frac{3}{4\omega ^2}}(1+2n+2n^2)+{\displaystyle \frac{1}{4\omega ^2}}(a^{+4}+a^4+(4n2)a^{+2}+(4n+6)a^2),`$ (12)
$`x^6`$ $`=`$ $`{\displaystyle \frac{5}{8\omega ^3}}(3+8n+6n^2+4n^3)+{\displaystyle \frac{1}{8\omega ^3}}(a^{+6}+a^6+(6n9)a^{+4}+(6n+15)a^4`$ (13)
$`+15(1n+n^2)a^{+2}+15(3+3n+n^2)a^2),`$
$`x^8`$ $`=`$ $`{\displaystyle \frac{1}{16\omega ^4}}(105+280n+350n^2+140n^3+70n^4)`$ (14)
$`+{\displaystyle \frac{1}{16\omega ^4}}(a^{+8}+a^8+(8n20)a^{+6}+(8n+28)a^6`$
$`+(9884n+28n^2)a^{+4}+(210+140n+28n^2)a^4`$
$`+(84+196n84n^2+56n^3)a^{+2}+(420+532n+252n^2+56n^3)a^2).`$
Any degree of $`x^2`$ can be derived by the use of a recursive relation between $`x^{2L}`$ and $`x^{2L2}`$. The normalized eigenvectors $`|n,\omega `$ of the excitation number operator $`n`$ are

$$|n,\omega =\frac{\left(a^+(\omega )\right)^n}{\sqrt{n!}}|0,\omega ,a(\omega )|0,\omega =0.$$
(15)
In contrast to the equation Eq.(2) we can calculate the matrix elements of the equation Eq.(8) exactly in the Fock basis Eq.(15). The eigenvector $`|\phi `$ can be represented in the Fock basis Eq.(15) as

$$|\phi =\underset{p=0}{\overset{\mathrm{}}{}}C_p|p,\omega ,$$
(16)
with the coefficient $`C_p`$ to be defined by the equation Eq.(8). The arguments of $`a^+(\omega )`$, $`a(\omega )`$, $`n(\omega )`$ and $`|n,\omega `$ will be omitted in order to simplify all expressions. Substitute Eq.(16) in Eq.(8) and find its projection on bra-vector $`k,\omega |`$

$$C_kk|^{}|kC_kk|^{\prime \prime }|kE+\underset{pk}{}C_pk|^{}|pE\underset{pk}{}C_pk|^{\prime \prime }|p=0.$$
(17)
The system of the linear equations Eq.(17) can be solved by the iterative methodl14

$`C_{k,kn}^{(s)}={\displaystyle \frac{_{p=km,pk}^{k+m}C_p^{(s1)}k|^{}|pE_n^{(s1)}_{p=km,pk}^{k+m}C_p^{(s1)}k|^{\prime \prime }|p}{k|^{}|kE_n^{(s1)}k|^{\prime \prime }|k}},`$
$`E_n^{(s)}={\displaystyle \frac{n|^{}|n+_{p=nm,pn}^{n+m}C_p^{(s1)}n|^{}|p}{n|^{\prime \prime }|n+_{p=nm,pn}^{n+m}C_p^{(s1)}n|^{\prime \prime }|p}},C_k^{(0)}=\delta _{n,k}`$ (18)
where $`C_k^{(s)}`$ is the coefficient of the eigenvector expansion and $`E_n^{(s)}`$ is the energy of $`n`$-level calculated in $`s`$-th iteration; $`m=2L`$; $`s=0,1,\mathrm{},s_{\mathrm{max}}`$. The results of numerical calculations of the energy eigenvalues for the ground ($`n=0`$) and first excited ($`n=1`$) states in the wide range of values of $`g`$, $`\lambda `$ are summarized in Tables 1-4. For convenience of the comparison with the results in Refs. l1, ; l20,  we present our results in the form $`2\lambda `$, $`2E`$. As is clear from Table 1 and Table 2, where some exact eigenvalues are given, the accuracy of the operator method is very high $`|E_n^{(s_{\mathrm{max}})}E_n|/E_n<10^{15}`$. But the case $`g1`$, $`\lambda 1`$ needs the increase of the iteration number $`s_{\mathrm{max}}`$. So it is important to make use of the parameter $`\omega `$, which allows us to speed up the convergence of the iterative scheme. Since the exact eigenvalues of the Hamiltonian do not depend on the choice of $`\omega `$, the following condition has to be satisfied for exact energy eigenvaluesl5

$$\frac{E_n}{\omega }=0.$$
(19)
Because $`E_n`$ is equal to $`E_n^{(\mathrm{})}(\omega )`$, a good accuracy can be achieved at the extremum point $`\omega _{extr}`$ of the function $`E_n^{(s_{\mathrm{max}})}(\omega )`$ (see Fig. 1). Notice that this statement is only valid in the case of large iteration number. As soon as the iteration number is relatively small $`E_n^{(s_{\mathrm{max}})}(\omega )`$ is described by an oscillating function crossing $`E=E_n(\mathrm{exact})`$ for some values of $`\omega `$ (see Fig. 1 and Ref. l21, ). In the case of large iteration number the magnitude of oscillations is very small and goes to zero quickly (see Fig. 1).

We emphasize that even for small number of iterations the points of extremum are near the exact energy eigenvalue and belong to the range of stable convergence of the iterative scheme. They can be chosen as starting values for $`\omega `$. As is clear from Table 2 the values of $`s_{\mathrm{max}}`$ and $`\omega `$ must be increased for excited states. Large values of $`g`$ and $`\lambda `$ require the increase of $`\omega `$ as well. Such behavior is typical for all potentials $`V_L(x)`$. Different criteria for choosing $`\omega `$ can be also found in Refs. l5, ; l18, ; l21, .

We’d like to notice that the normalized eigenfunctions $`\mathrm{\Psi }_n(x)`$ of the original SE Eq.(2) can be easily calculated

$`\mathrm{\Psi }_n(x)={\displaystyle \frac{\psi _n(x)}{\sqrt{\psi _n|\psi _n}}}{\displaystyle \frac{1}{N}}(1+gx^2){\displaystyle \underset{k=0}{\overset{s_{\mathrm{max}}}{}}}C_k^{(s_{\mathrm{max}})}x|k,\omega `$ (20)
$`N^2=\psi _n|\psi _n{\displaystyle _+\mathrm{}^{\mathrm{}}}dx(1+gx^2)^2\left[{\displaystyle \underset{k=0}{\overset{s_{\mathrm{max}}}{}}}C_k^{(s_{\mathrm{max}})}x|k,\omega \right]^2.`$ (21)
The polynomial structure of the operator equation Eq.(4) also allows us to obtain some analytical expressions for the energy eigenvalues with the explicit dependence on all parameters. Such expressions can be found by making a limited number of iterations with any system of computer algebra like Maple or Mathematica.

## III Conclusion

The Schrödinger equation with the rational potential was solved by the operator method. The energy of the ground state and the first excited state were calculated with very high accuracy without any limitations on the parameters of the Hamiltonian. The numerical solutions are in excellent agreement with the exact solutions known for some values of $`g`$ and $`\lambda `$. The criteria for choosing the convergence parameter $`\omega `$ was analyzed. So we have shown that the spectrum of the original nonpolynomial Schrödinger equation Eq.(2) is identical to the one of the polynomial operator equation Eq.(4). The results obtained demonstrate that the operator method is very efficiently for solving the problems with the nonpolynomial structure of the Hamiltonian if the nonunitary transformation is used. The application of the operator method to the three-dimensional nonpolynomial Hamiltonian is straightforward.

## Acknowledgments

The author is grateful to L. I. Komarov for useful discussions.


text:
# Self-trapping and stable localized modes in nonlinear photonic crystals

## Abstract

We predict the existence of stable nonlinear localized modes near the band edge of a two-dimensional reduced-symmetry photonic crystal with a Kerr nonlinearity. Employing the technique based on the Green function, we reveal a physical mechanism of the mode stabilization associated with the effective nonlinear dispersion and long-range interaction in the photonic crystals.

Photonic crystals are usually viewed as an optical analog of semiconductors that modify the properties of light similar to a microscopic atomic lattice that creates a semiconductor band-gap for electrons . It is therefore believed that by replacing relatively slow electrons with photons as the carriers of information, the speed and band-width of advanced communication systems will be dramatically increased, thus revolutionizing the telecommunication industry. To employ the high-tech potential of photonic crystals, it is crucially important to achieve a dynamical tunability of their band gap . This idea can be realized by changing the light intensity in the so-called nonlinear photonic crystals, having a periodic modulation of the nonlinear refractive index . Exploration of nonlinear properties of photonic band-gap (PBG) materials may open new applications of photonic crystals for all-optical signal processing and switching, suggesting an effective way to create tunable band-gap structures operating entirely with light.

One of the important ideas to control all-optical switching in the nonlinear regime is to explore the possibility of nonlinearity-induced self-trapping and nonlinear localized modes in photonic crystals. Existence of nonlinear localized modes for the frequencies in the forbidden gaps is usually associated with gap solitons, studied for one-dimensional  and even two-dimensional (2D) models , in the framework of the coupled-mode theory. Validity of the coupled-mode theory is usually restricted by a weak modulation of the refractive index (the so-called shallow-grating case), and therefore this theory is not directly applicable to the PBG crystals where modulation of the refractive index is of the order of its mean value. This observation calls for a systematic analysis of the (still open) problem of stable self-trapping and nonlinear localized modes in PBG materials, where the effects of discreteness  and long-range interaction  have been recently shown to be of crucial importance.

Nonlinear localized modes (also called intrinsic localized modes or discrete breathers) are associated with the energy localization that may occur in the absence of any disorder and solely due to nonlinearity . Such nonlinear modes can be easily identified as approximate (or sometimes exact) analytical solutions of coupled-oscillator nonlinear lattice models and also in numerical molecular-dynamics simulations, but only very recently the first observations of spatially localized nonlinear modes have been reported in the physical systems of a very different nature . The main purpose of this Letter is to predict the existence of nonlinear localized modes, analogous to gap solitons in the continuum limit, in 2D nonlinear photonic crystals, and to describe their unique properties including stability.

We study nonlinear properties of 2D photonic crystals, assuming that their symmetry is reduced by inserting the rods made from a Kerr-type nonlinear dielectric material characterized by the third-order nonlinear susceptibility $`\chi ^{(3)}`$. Specifically, we consider a periodic square lattice with the lattice spacing $`a`$ which consists of two types of infinitely long cylindrical rods: the rods of radius $`r_1`$ made from a linear material and placed at the corners of the lattice, and the rods of radius $`r_2`$ ($`r_2<r_1`$) made from a nonlinear material and placed at the center of each unit cell (see top right inset in Fig. 1). Linear properties of such photonic crystals are known .

We assume that the rods are parallel to the $`x_3`$ axis, so that the system is characterized by the dielectric constant $`ϵ(𝒙)=ϵ(x_1,x_2)`$. In this case the evolution of the $`E`$-polarized (with the electric field $`𝑬||𝒙_3`$) light propagating in the $`(x_1,x_2)`$-plane is governed by the scalar wave equation

$$^2E(𝒙,t)\frac{1}{c^2}_t^2\left[ϵ(𝒙)E\right]=0,$$
(1)
where $`^2_{x_1}^2+_{x_2}^2`$ and $`E`$ is the $`x_3`$ component of $`𝑬`$. Taking the electric field in the form $`E(𝒙,t)=e^{i\omega t}E(𝒙,t|\omega )`$, where $`E(𝒙,t|\omega )`$ is a slowly varying envelope, i.e. $`_t^2E(𝒙,t|\omega )\omega _tE(𝒙,t|\omega )`$, Eq. (1) can be reduced to

$$\left[^2+ϵ(𝒙)\left(\frac{\omega }{c}\right)^2\right]E(𝒙,t|\omega )2iϵ(𝒙)\frac{\omega }{c^2}\frac{E}{t}.$$
(2)
In the stationary case, i.e. when the r.h.s. vanishes, Eq. (2) reduces to the eigenvalue problem which can be solved in the linear limit , i.e. when $`ϵ(𝒙)`$ does not depend on the light intensity. In this case, the frequency spectrum has a band-gap structure shown in Fig. 1. It supports two band gaps, the lower of which extends from $`\omega =0.426\times 2\pi c/a`$ to $`\omega =0.453\times 2\pi c/a`$.

A low-intensity light cannot propagate through a photonic crystal if the light frequency falls into a band gap. However, it has been recently suggested  that in the case of a 2D periodic medium with a Kerr-type nonlinear material, high-intensity light with the frequency inside the photonic gap can propagate in the form of finite energy solitary waves2D gap solitons. Such solitary waves were analyzed in the framework of the coupled-mode continuum-limit equations valid for a weak modulation of the dielectric constant $`ϵ(𝒙)`$. However, in real photonic crystals the modulation of $`ϵ(𝒙)`$ is comparable to its average value, so that the existence and stability of localized modes in such structures is still an open problem.

More specifically, the coupled-mode equations are valid if and only if the band gap $`\mathrm{\Delta }\omega ^2`$ is vanishingly small, i. e. $`\mathrm{\Delta }\omega ^2A^2`$ where $`A`$ is an effective amplitude of the mode that is treated in the multi-scale asymptotic expansions  as a small parameter. If we apply this model to describe nonlinear modes in a wider gap (see, e.g., discussions in Ref. ), we obtain a 2D cubic nonlinear Schrödinger (NLS) equation known to possess no stable localized solutions. Moreover, 2D localized modes of the coupled-mode equations are expected to possess an oscillatory instability recently discovered for a broad class of 1D coupled-mode equations . Thus, it is clear that, if nonlinear localized modes do exist in realistic PBG materials, their existence and stability should be associated with different physical mechanisms not accounted for by simplified continuum coupled-mode models.

To study the nonlinear modes in such structures, we consider the nonlinear rods of small radius $`r_2`$ as “defects” embedded into the linear photonic crystal formed by a square lattice of the rods of larger radius $`r_1`$ (”diatomic crystal”). Then, writing the dielectric constant $`ϵ(𝒙)`$ as a sum of two periodic terms, $`ϵ(𝒙)=ϵ_1(𝒙)+ϵ_2(𝒙|E),`$ where $`ϵ_1(𝒙)`$ describes the linear photonic crystal and $`ϵ_2(𝒙|E)`$ corresponds to a lattice of nonlinear defect rods, one can present Eq. (2) in the form

$$E(𝒙,t|\omega )=d^2𝒚G(𝒙,𝒚|\omega )\widehat{}E(𝒚,t|\omega ),$$
(3)
where we introduce the Green function $`G(𝒙,𝒚|\omega )`$ of the linear photonic crystal (see, e.g., Ref.  for its properties) and the linear operator

$`\widehat{}=\left({\displaystyle \frac{\omega }{c}}\right)^2ϵ_2(𝒙|E)+2iϵ(𝒙){\displaystyle \frac{\omega }{c^2}}{\displaystyle \frac{}{t}}.`$ (4)
Now, using the indices $`n`$ and $`m`$ for numbering the nonlinear rods in the $`x_1`$ and $`x_2`$ directions, we can describe their positions by the vectors $`𝒙_{n,m}=n𝒂_1+m𝒂_2`$, where $`𝒂_1`$ and $`𝒂_2`$ are the primitive lattice vectors of the 2D photonic crystal, and write

$$ϵ_2(𝒙|E)=\{ϵ_2^{(0)}+\chi ^{(3)}|E(𝒙,t|\omega )|^2\}\underset{n,m}{}\theta (𝒙𝒙_{n,m}),$$

where $`\theta (𝒙)=1`$ for $`|𝒙|r_2`$ and $`\theta (𝒙)=0`$ otherwise. The parameter $`ϵ_2^{(0)}`$ is the dielectric constant of the defect rods in the linear limit, while the term $`\chi ^{(3)}|E|^2`$ takes into account a contribution due to the Kerr nonlinearity. If the radius $`r_2`$ of the defect rods is sufficiently small, the electric field $`E(𝒙,t|\omega )`$ inside them is almost constant, and Eq. (3) can be approximated  by the discrete nonlinear equation

$`i\sigma {\displaystyle \frac{}{t}}E_{n,m}E_{n,m}`$ (5)
$`+{\displaystyle \underset{k,l}{}}J_{nk,ml}(\omega )\left\{ϵ_2^{(0)}+\chi ^{(3)}|E_{k,l}|^2\right\}E_{k,l}=0,`$ (6)
for the amplitudes of the electric field $`E_{n,m}(t|\omega )E(𝒙_{n,m},t|\omega )`$ calculated at the defect rods. The parameter $`\sigma `$ and the coupling coefficients

$$J_{n,m}(\omega )=\left(\frac{\omega }{c}\right)^2_{r_2}d^2𝒚G(𝒙_{0,0},𝒙_{n,m}+𝒚|\omega )$$
(7)
are determined by the Green function $`G(𝒙,𝒚|\omega )`$ which we calculate numerically by the FDTD method  with the spatial step $`\mathrm{\Delta }x=0.01`$ to $`0.03`$ and time step $`\mathrm{\Delta }t=0.005`$, for the lattice $`24a\times 24a`$. The Green function and, therefore, the coupling coefficients $`J_{n,m}(\omega )`$ in 2D photonic crystals are usually long-ranged functions. For instance, for the case of Fig. 2 we obtain $`J_{n,0}0.012(1)^n\mathrm{exp}(0.66|n|)`$ for $`n2`$, so that one should take into account the interaction between at least 10 neighbors to achieve a good accuracy for the spectrum. By this means, Eq. (5) is a nontrivial long-range generalization of the 2D discrete NLS equations extensively studied during the last decade . We have checked the accuracy of the approximation provided by Eq. (5) by solving it in the linear limit. The low-frequency part of this dependence for $`ϵ_2^{(0)}=11.4`$ is depicted in Fig. 1 by a dashed line; it has a minimum at $`\omega =0.446\times 2\pi c/a`$ being in a good agreement with the band edge calculated directly from Eq. (2). This lends a support to the validity of Eq. (5) and allows us to use this discrete model for studying nonlinear properties.

Stationary nonlinear localized modes have been calculated numerically by the Newton-Raphson iteration scheme on the lattice $`170a\times 170a`$. Assuming that the linear rods are made from GaAs whereas the nonlinear rods are made from some nonlinear material (which we do not specify varying $`ϵ_2^{(0)}`$), we find a continuous family of nonlinear localized modes, and a typical example \[smoothed by continuous optimization for Eq. (3)\] of such a mode is shown in Fig. 3. At first glance, this mode can be regarded as a donor state created by a single defect rod with larger dielectric constant. However, in the nonlinear case the mode stability becomes a critical issue. It can be determined from the so-called Vakhitov-Kolokolov stability criterion extended to the 2D discrete NLS models in Ref. . According to this criterion, the nonlinear localized states with $`dQ/d\omega <0`$ are stable, and they are unstable otherwise. Here

$$Q(\omega )=\underset{n,m}{}|E_{n,m}|^2,$$
(8)
is the conserved mode power proportional to the energy of the electric field accumulated in the nonlinear localized mode.

As is well known (see, e.g., Refs. ) in the 2D discrete cubic NLS equation, only high-amplitude localized modes are stable, whereas no stable modes exist in the continuum limit. For our model, the high-amplitude modes are also stable (see inset in Fig. 4), but they are not accessible under realistic conditions: To excite such modes one should increase the refractive index at the mode center in more than 2 times. Thus, for realistic conditions and relatively small values of $`\chi ^{(3)}`$, only low-amplitude localized modes become a subject of much interest since they can be excited in experiment. However, such modes in unbounded 2D NLS models are always unstable and either collapse or spread out . Here we reveal that, in a sharp contrast to the 2D discrete NLS models discussed earlier in various applications, the low-amplitude localized modes of Eq. (5) can be stabilized due to nonlinear long-range dispersion inherent to the photonic crystals. It should be emphasized that such stabilization does not occur in the models with only linear long-range dispersion .

In order to gain a better insight into the stabilization mechanism, we have carried out the studies of Eq. (5) for the exponentially decaying coupling coefficients $`J_{n,m}`$. Our results show that the most important factor which determines stability of the low-amplitude localized modes is a ratio of the coefficients at the local nonlinearity ($`J_{0,0}`$) and the nonlinear dispersion ($`J_{0,1}`$). If the coupling coefficients $`J_{n,m}`$ decrease with the distances $`n`$ and $`m`$ rapidly, the low-amplitude modes of Eq. (5) with $`ϵ_2^{(0)}=11.4`$ are essentially stable for $`J_{0,0}/J_{0,1}13`$. However, this estimation is usually lowered because the stabilization is favored by the presence of long-range interactions.

It should be mentioned that the stabilization of low-amplitude 2D localized modes is not inherent to all types of nonlinear photonic crystals. On the contrary, the photonic crystals must be carefully designed to support stable low-amplitude nonlinear modes. For example, in the photonic crystal considered above such modes are stable at least for $`11<ϵ_2^{(0)}<12`$, however they become unstable for $`ϵ_2^{(0)}12`$ (see Fig. 4). The stability of these modes can also be controlled by varying $`r_2`$, $`r_1`$, or $`ϵ_1`$. Thus, experimental observation of the nonlinear localized modes would require not only the use of photonic materials with a relatively large nonlinear refractive index (such as GaAs/AlAs periodic structures  or polymer PBG crystals ), but also a fine adjustment of the parameters of the photonic crystal. The latter can be achieved, in principle, by employing the surface coupling technique  that is able to provide coupling to specific points of the dispersion curve, opening up a very straightforward way to access nonlinear effects.

In conclusion, we have developed a consistent theory of nonlinearity-induced self-trapping effects in 2D nonlinear photonic crystals and predicted the possibility of the energy localization near the bandgap edge in the form of stable 2D nonlinear localized modes.

Yuri Kivshar thanks V. Astratov, K. Busch, S. John, A. McGurn, M. Scalora, and C. Soukoulis for encouraging comments and discussions. Serge Mingaleev is indebted to Yu.B. Gaididei for useful remarks. The work has been partially supported by the Australian Research Council and the Performance and Planning Foundation grant of the Institute of Advanced Studies.