text:
# DYNAMICAL SUSCEPTIBILITIES IN STRONG COUPLING APPROACH: GENERAL SCHEME AND FALIKOV-KIMBALL MODEL

## I Introduction

The development of the Dynamical Mean Field Theory (DMFT), which is exact in the $`d=\mathrm{}`$ limit, clarified some problems connected with the simultaneous consideration of the electron hopping and strong local correlations and stimulated a large progress in the understanding of the strongly correlated electron systems . It was shown by Metzner and Vollhardt  that in the $`d=\mathrm{}`$ limit self-energies are single-site quantities (do not depend on wave vector) which leads to a significant simplification. But such a local self-energy only probes local properties in this limit and cannot detect instabilities associated with a specific wave-vector  that requires the calculation of susceptibilities.

An analytical schemes developed for the descriptions of the strongly correlated electrons can be divided into two types: (i) weak-coupling theories, that are based on the expansion over local many-electron interactions (Dyson approach) and are close to the standard Fermi-liquid theory, and (ii) strong coupling theories, that start from the expansion over electron hopping around the atomic limit (see Ref.  and references therein). In the weak-coupling approach susceptibilities are solutions of the Bethe-Salpeter equation and it is known that in the $`d\mathrm{}`$ limit the corresponding irreducible four-vertices are also local .

The aim of this article is to develop a general scheme to calculate dynamical susceptibilities within a strong-coupling DMFT approach for strongly correlated electron systems. A special case of the binary alloy \[Falicov-Kimball (FK)\] type models, that can be solved exactly in the limit of infinite dimensions , will be also considered. Static susceptibilities for the FK model have been already investigated by Brandt and Mielsch , who found an Ising-like phase transition to a chess-board charge-density-wave phase at half filling, and Freericks , who showed that system also displayed the incommensurate order and phase separation at other fillings.

The preliminary short version of this article was published in the Proceedings of M2S-HTSC-VI Conference . Here we present the details of the general scheme.

## II Two approaches in many-electron theory

In general, the Hamiltonian of the electronic system with strong local correlations can be written in the following form

$$H=H_t+H_{\text{loc}},$$
(1)
where

$$H_t=\underset{ij\sigma }{}t_{ij}a_{i\sigma }^{}a_{j\sigma }$$
(2)
describes an intersite electron hopping and $`H_{\text{loc}}`$ is a sum of the single site Hamiltonians

$$H_{\text{loc}}=\underset{i}{}H_i$$
(3)
which describe local electron correlations and/or interaction with other local excitations (lattice vibrations, (pseudo)spins, etc.).

### A Weak-coupling approach

As a rule, the first term in (1) is considered as an initial (zero-order) Hamiltonian and the second one is considered as a perturbation. Such approach is well known as the Dyson weak-coupling approach, where the single electron Green’s functions are determined by the Dyson equation

$$G_{\sigma 𝒌}(\omega _\nu )\text{}$$
(4)
or

$`G_{\sigma 𝒌}(\omega _\nu )`$ $`=`$ $`𝒢_{\sigma 𝒌}(\omega _\nu )+𝒢_{\sigma 𝒌}(\omega _\nu )\mathrm{\Sigma }_{\sigma 𝒌}(\omega _\nu )G_{\sigma 𝒌}(\omega _\nu )`$ (5)
$`=`$ $`{\displaystyle \frac{1}{i\omega _\nu +\mu t_𝒌\mathrm{\Sigma }_{\sigma 𝒌}(\omega _\nu )}}.`$ (6)
Here

$$𝒢_{\sigma 𝒌}(\omega _\nu )=\frac{1}{i\omega _\nu +\mu t_𝒌}$$
(7)
and the self-energy $`\mathrm{\Sigma }_{\sigma 𝒌}(\omega _\nu )`$ describes all scattering processes originated from the second term in (1) and cannot be divided into parts by cutting off one zero-order Green’s function $`𝒢_{\sigma 𝒌}(\omega _\nu )`$ line.

Also, it is well known that within such a weak-coupling approach the electron susceptibilities (charge, spin, etc.) can be presented in the form

$$\chi _𝒒^{AB}(\omega _m)=\text{},$$
(8)
where the type of susceptibility is determined by the ending parts. The full four vertex is a solution of the Bethe-Salpeter equation

(9)
with irreducible four vertex D that, in a same way as self-energy, describes scattering processes originated from local correlations and cannot be divided into parts by cutting off two zero-order Green’s function lines.

### B Strong-coupling approach

On the other hand, an alternative approach based on the expansion over electron hopping around the strong-coupling limit can be built . In this case, the single-electron Green’s functions are determined by the Larkin equation

$$G_{\sigma 𝒌}(\omega _\nu )\text{}$$
(10)
or

$`G_{\sigma 𝒌}(\omega _\nu )`$ $`=`$ $`\mathrm{\Xi }_{\sigma 𝒌}(\omega _\nu )+\mathrm{\Xi }_{\sigma 𝒌}(\omega _\nu )t_𝒌G_{\sigma 𝒌}(\omega _\nu )`$ (11)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Xi }_{\sigma 𝒌}^1(\omega _\nu )t_𝒌}},`$ (12)
where $`\mathrm{\Xi }_{\sigma 𝒌}(\omega _\nu )`$ is an irreducible according to Larkin part which cannot be divided into parts by cutting off one hopping (wavy) line.

Now, susceptibilities can be presented in the following form

$$\chi _𝒒^{AB}(i\omega _m)=\text{},$$
(13)
where

(14)
is a sum of the chains of hopping lines and full four vertex is a solution of the equation

$$\text{}.$$
(15)
In Eqs. (13) and (15) quantities

### C Connection between approaches

From Eqs. (5) and (11) one can get a connection between the self-energy and the irreducible vertex function according to the Larkin equation:

$$\mathrm{\Xi }_{\sigma 𝒌}^1(i\omega _\nu )=i\omega _\nu +\mu \mathrm{\Sigma }_{\sigma 𝒌}(i\omega _\nu ).$$
(16)
The connection between the four verticies in the weak and strong coupling approaches is more complicated. For the full four verticies one can get

$$\text{},$$
(17)
whereas the irreducible four verticies are connected by equation

$$\text{},$$
(18)
where $`D^{}`$ is a solution of equation

$$\text{}.$$
(19)
Here

(20)
In the case when the operators $`\widehat{A}`$ and $`\widehat{B}`$ are constructed only by the electronic operators, then for three and two verticies we find

(21)
and

(22)
which can be treated as some kind of sum rules.

## III Dynamical mean field theory

In the case of infinite dimensions $`d\mathrm{}`$ one should scale the hopping integral according to

$$t_{ij}\frac{t_{ij}}{\sqrt{d}}$$
(23)
in order to obtain a finite density-of-states and it was shown by Metzner  that in this limit the irreducible part becomes local

$$\mathrm{\Xi }_{ij\sigma }(\tau \tau ^{})=\delta _{ij}\mathrm{\Xi }_\sigma (\tau \tau ^{})$$
(24)
or

$$\mathrm{\Xi }_{\sigma 𝒌}(\omega _\nu )=\mathrm{\Xi }_\sigma (\omega _\nu ).$$
(25)
Such a site-diagonal function, as it was shown by Brandt and Mielsch , can be calculated by mapping the $`d\mathrm{}`$ lattice problem (1) with intersite hopping (23) onto the atomic model with single-site hopping

$$H_t_0^\beta 𝑑\tau _0^\beta 𝑑\tau ^{}\underset{\sigma }{}\zeta _\sigma (\tau \tau ^{})a_\sigma ^{}(\tau )a_\sigma (\tau ^{}).$$
(26)
Here $`\zeta _\sigma (\tau \tau ^{})`$ is the auxiliary Kadanoff-Baym field (dynamical mean field) which has to be self-consistently determined from the condition that the same function $`\mathrm{\Xi }_\sigma (\omega _\nu )`$ defines the Green’s function for the lattice and atomic model. The self-consistent set of equations for $`\mathrm{\Xi }_\sigma (\omega _\nu )`$ and $`\zeta _\sigma (\omega _\nu )`$ (e.g., see Ref.  and references therein) is the following:

$`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝒌}{}}{\displaystyle \frac{1}{\mathrm{\Xi }_\sigma ^1(\omega _\nu )t_𝒌}}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Xi }_\sigma ^1(\omega _\nu )\zeta _\sigma (\omega _\nu )}}`$ (27)
$`=`$ $`G_\sigma ^{(a)}(\omega _\nu ,\{\zeta _\sigma (\omega _\nu )\}),`$ (28)
where $`G_\sigma ^{(a)}(\omega _\nu ,\{\zeta _\sigma (\omega _\nu )\})`$ is the Green’s function for the atomic model (26).

In the same way as it was done by Metzner  for the irreducible part $`\mathrm{\Xi }_\sigma (\omega _\nu )`$, one can prove that in the $`d\mathrm{}`$ limit all irreducible verticies in the strong-coupling approach (Eqs. (13) and (15)) are also single site quantities. So, they can be also be calculated from the atomic model. On the other hand, using Eqs. (16)–(22) it is easy to show that the self-energy and irreducible verticies in the weak-coupling approach are also local quantities (in Ref.  that was proved using different approach).

In order to calculate these single-site irreducible verticies first of all one has to calculate an irreducible many particle Green’s functions for atomic model (26), i.e.

$``$ $`Ta_1^{}a_2a_3^{}a_4_a`$ (29)
$``$ $`Ta_1^{}a_2_aTa_3^{}a_4_a+Ta_1^{}a_4_aTa_3^{}a_2_a,`$ (30)
$$\text{}T\widehat{A}_0a_1^{}a_2_a\widehat{A}_0_aTa_1^{}a_2_a$$
(31)
and

$$\text{}T\widehat{A}_0\widehat{B}_1_a\widehat{A}_0_a\widehat{B}_1_a.$$
(32)
Then irreducible four vertex $`D`$ for the weak-coupling approach can be obtained from equation (9), where now arrows indicate single-electron Green’s functions for the atomic model (27) and full four vertex is determined from the many-particle Green’s function by

$$\text{}.$$
(33)
On the other hand, irreducible four vertex $`L`$ for the strong-coupling approach can be obtained from equation (15), where the full four vertex is determined from the many-particle Green’s function by

(34)
or

$$\text{}.$$
(35)
Now, the thin wavy lines are the dynamical mean field $`\zeta _\sigma (\omega _\nu )`$ and the thick ones, in contrast to (14), represent the sum of chains of wavy lines

$$\stackrel{~}{\zeta }_\sigma (\omega _\nu )=\frac{\zeta _\sigma (\omega _\nu )}{1\mathrm{\Xi }_\sigma (\omega _\nu )\zeta _\sigma (\omega _\nu )}=\zeta _\sigma (\omega _\nu )\frac{G_\sigma ^{(a)}(\omega _\nu )}{\mathrm{\Xi }_\sigma (\omega _\nu )}.$$
(36)
From the expression for the irreducible many particle Green’s function (31), (32) and(29) for atomic model one can find irreducible verticies $`|`$ and

(37)
and

$$\text{},$$
(38)
respectively, that complete the calculation of the irreducible verticies for the dynamical susceptibilities.

Finally, for the lattice dynamical susceptibility in the weak (8) and strong (13) coupling approaches we get

$$\text{},$$
(39)
and

$$\text{},$$
(40)
respectively. Here, $`(\mathrm{})^1`$ denote an inverse kernels of the corresponding integral equations,

$$\text{}=G_\sigma ^{(a)}(\omega _\nu )G_\sigma ^{(a)}(\omega _{\nu +m}),$$
(41)
and

$$\text{}=\frac{1}{N}\underset{𝒌}{}G_{\sigma 𝒌}(\omega _\nu )G_{\sigma 𝒌+𝒒}(\omega _{\nu +m}).$$
(42)
## IV Binary alloy type models

To test the possibilities of the above approach we consider a binary alloy type model

$$H_i=\frac{g}{2}\left(P_i^+P_i^{}\right)n_i\frac{h}{2}\left(P_i^+P_i^{}\right),$$
(43)
where $`P_i^\pm =\frac{1}{2}\pm S_i^z`$ for the $`U=0`$ pseudospin-electron (PE) model , $`P_i^+=c_i`$, $`P_i^{}=1c_i`$ for a binary alloy, and $`a_{i\sigma }d_i`$, $`n_i=d_i^{}d_i`$, $`P_i^+=f_i^{}f_i`$, $`P_i^{}=1f_i^{}f_i`$ for the Falicov-Kimball (FK) model . The main difference between these models is in the way that the averaging procedure is performed (a statistical one for the PE and FK models and a configurational one for the binary alloy) and in the thermodynamical regimes ($`h=\mathrm{const}`$ for the PE model, $`c=\mathrm{const}`$ for the binary alloy and $`n_d=\mathrm{const}`$ and/or $`n_f=\mathrm{const}`$ for the FK model).

The single-particle Green’s function for the effective atomic model (26) is a coherent sum of the Green’s functions for subspaces $`S_i^z=\pm \frac{1}{2}`$ and is equal to (see Refs. )

$`G_\sigma ^{(a)}(\omega _\nu )`$ $`=`$ $`{\displaystyle \frac{P^+}{i\omega _\nu +\mu \zeta _\sigma (\omega _\nu )\frac{g}{2}}}`$ (45)
$`+{\displaystyle \frac{P^{}}{i\omega _\nu +\mu \zeta _\sigma (\omega _\nu )+\frac{g}{2}}}`$
which allows us to find solutions of Eq. (27) for the irreducible part $`\mathrm{\Xi }_\sigma (\omega _\nu )`$ and dynamical mean field $`\zeta _\sigma (\omega _\nu )`$.

The Fourier transform of the two-electron Green’s function $`Ta_1^{}a_2a_3^{}a_4`$ is also a coherent sum of the two-electron Green’s functions for subspaces

$`G_{\sigma _1\sigma _1^{}}^{\sigma \sigma ^{}}\left({\displaystyle \genfrac{}{}{0pt}{}{\omega _\nu }{\omega _{\nu +m}}}{\displaystyle \genfrac{}{}{0pt}{}{\omega _\nu ^{}}{\omega _{\nu ^{}+m}}}\right)={\displaystyle \underset{\pm }{}}`$ (46)
$`\{{\displaystyle \frac{\delta _{m0}\delta _{\sigma \sigma _1}\delta _{\sigma ^{}\sigma _1^{}}P^\pm }{\left[i\omega _\nu +\mu \zeta _\sigma (\omega _\nu )\frac{g}{2}\right]\left[i\omega _\nu ^{}+\mu \zeta _\sigma ^{}(\omega _\nu ^{})\frac{g}{2}\right]}}`$ (47)
$`{\displaystyle \frac{\delta _{\nu \nu ^{}}\delta _{\sigma \sigma ^{}}\delta _{\sigma _1\sigma _1^{}}P^\pm }{\left[i\omega _\nu +\mu \zeta _\sigma (\omega _\nu )\frac{g}{2}\right]\left[i\omega _{\nu +m}+\mu \zeta _{\sigma _1}(\omega _{\nu +m})\frac{g}{2}\right]}}\}`$ (48)
that gives for the irreducible two-electron Green’s function (29) an expression which is a sum of the two separable functions

$`\stackrel{~}{G}_{\sigma _1\sigma _1^{}}^{\sigma \sigma ^{}}\left({\displaystyle \genfrac{}{}{0pt}{}{\omega _\nu }{\omega _{\nu +m}}}{\displaystyle \genfrac{}{}{0pt}{}{\omega _\nu ^{}}{\omega _{\nu ^{}+m}}}\right)=`$ $`\delta _{m0}\delta _{\sigma \sigma _1}\delta _{\sigma ^{}\sigma _1^{}}\mathrm{\Lambda }_\sigma (\omega _\nu )\mathrm{\Lambda }_\sigma ^{}(\omega _\nu ^{})`$ (49)
$`\delta _{\nu \nu ^{}}\delta _{\sigma \sigma ^{}}`$ $`\delta _{\sigma _1\sigma _1^{}}\mathrm{\Lambda }_\sigma (\omega _\nu )\mathrm{\Lambda }_{\sigma _1}(\omega _{\nu +m}),`$ (50)
where

$$\mathrm{\Lambda }_\sigma (\omega _\nu )=\frac{g\sqrt{P^+P^{}}}{\left[i\omega _\nu +\mu \zeta _\sigma (\omega _\nu )\right]^2\frac{g^2}{4}}.$$
(51)
Such a separable form allows us to calculate the inverse kernels of the integral equations and to find all quantities in an analytical form.

Within the weak-coupling approach, starting from the expression (50) with the use of Eq. (33), one can find solution of the Bethe-Salpeter equation for the irreducible four vertex in an analytical form:

$`\mathrm{\Gamma }_{\sigma _1\sigma _1^{}}^{\sigma \sigma ^{}}\left({\displaystyle \genfrac{}{}{0pt}{}{\omega _\nu }{\omega _{\nu +m}}}{\displaystyle \genfrac{}{}{0pt}{}{\omega _\nu ^{}}{\omega _{\nu ^{}+m}}}\right)={\displaystyle \frac{\mathrm{\Lambda }_\sigma (\omega _\nu )}{G_\sigma ^2(\omega _\nu )+\mathrm{\Lambda }_\sigma ^2(\omega _\nu )}}\times `$ (52)
$`{\displaystyle \frac{\delta _{m0}\delta _{\sigma \sigma _1}\delta _{\sigma ^{}\sigma _1^{}}}{1\underset{\nu ^{\prime \prime }\sigma ^{\prime \prime }}{}\frac{\mathrm{\Lambda }_{\sigma ^{\prime \prime }}^2(\omega _{\nu ^{\prime \prime }})}{G_{\sigma ^{\prime \prime }}^2(\omega _{\nu ^{\prime \prime }})+\mathrm{\Lambda }_{\sigma ^{\prime \prime }}^2(\omega _{\nu ^{\prime \prime }})}}}{\displaystyle \frac{\mathrm{\Lambda }_\sigma ^{}(\omega _\nu ^{})}{G_\sigma ^{}^2(\omega _\nu ^{})+\mathrm{\Lambda }_\sigma ^{}^2(\omega _\nu ^{})}}`$ (53)
$`{\displaystyle \frac{\delta _{\nu \nu ^{}}\delta _{\sigma \sigma ^{}}\delta _{\sigma _1\sigma _1^{}}}{G_\sigma (\omega _\nu )G_{\sigma _1}(\omega _{\nu +m})}}\times `$ (54)
$`{\displaystyle \frac{\mathrm{\Lambda }_\sigma (\omega _\nu )\mathrm{\Lambda }_{\sigma _1}(\omega _{\nu +m})}{G_\sigma (\omega _\nu )G_{\sigma _1}(\omega _{\nu +m})+\mathrm{\Lambda }_\sigma (\omega _\nu )\mathrm{\Lambda }_{\sigma _1}(\omega _{\nu +m})}},`$ (55)
which must be substituted into Eq. (9) for the full four vertex for lattice and allows to calculate electron susceptibilities, e.g. for charge susceptibility one can get

$$\chi _𝒒^{nn}(\omega _m)=\delta _{m0}\frac{\mathrm{\Delta }_n^2}{T\mathrm{\Theta }(T,𝒒)}+K_𝒒^{nn}(\omega _m),$$
(56)
$$K_𝒒^{nn}(\omega _m)=\frac{1}{\beta }\underset{\nu \sigma }{}\frac{1}{\chi _{\sigma 𝒒}^1(\omega _\nu ,\omega _m)\mathrm{\Gamma }_\sigma (\omega _\nu ,\omega _m)},$$
(57)
where

$`\mathrm{\Gamma }_\sigma (\omega _\nu ,\omega _m)`$ $`=`$ (59)
$`{\displaystyle \frac{\mathrm{\Lambda }_\sigma (\omega _\nu )\mathrm{\Lambda }_\sigma (\omega _{\nu +m})}{\stackrel{~}{\chi }_\sigma (\omega _\nu \omega _m)\left[\mathrm{\Lambda }_\sigma (\omega _\nu )\mathrm{\Lambda }_\sigma (\omega _{\nu +m})\stackrel{~}{\chi }_\sigma (\omega _\nu \omega _m)\right]}}`$
originates from the last term in (52) and

$`\mathrm{\Delta }_n`$ $`=`$ (61)
$`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\nu \sigma }{}}{\displaystyle \frac{\mathrm{\Lambda }_\sigma (\omega _\nu )\stackrel{~}{\chi }_\sigma (\omega _\nu 0)\chi _{\sigma 𝒒}(\omega _\nu 0)}{\stackrel{~}{\chi }_\sigma ^2(\omega _\nu 0)+\mathrm{\Lambda }_\sigma ^2(\omega _\nu )\left[\chi _{\sigma 𝒒}(\omega _\nu 0)\stackrel{~}{\chi }_\sigma (\omega _\nu 0)\right]}},`$
$`\mathrm{\Theta }(T,𝒒)`$ $`=`$ (63)
$`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\nu \sigma }{}}{\displaystyle \frac{\mathrm{\Lambda }_\sigma ^2(\omega _\nu )\left[\chi _{\sigma 𝒒}(\omega _\nu 0)\stackrel{~}{\chi }_\sigma (\omega _\nu 0)\right]}{\stackrel{~}{\chi }_\sigma ^2(\omega _\nu 0)+\mathrm{\Lambda }_\sigma ^2(\omega _\nu )\left[\chi _{\sigma 𝒒}(\omega _\nu 0)\stackrel{~}{\chi }_\sigma (\omega _\nu 0)\right]}},`$
$$\chi _{\sigma 𝒒}(\omega _\nu \omega _m)=\frac{1}{N}\underset{𝒌}{}G_{\sigma 𝒌}(\omega _\nu )G_{\sigma 𝒌+𝒒}(\omega _{\nu +m}),$$
(64)
$$\stackrel{~}{\chi }_\sigma (\omega _\nu \omega _m)=G_\sigma ^{(a)}(\omega _\nu )G_\sigma ^{(a)}(\omega _{\nu +m}).$$
(65)
But such a diagrammatic weak-coupling approach does not allow one to calculate pseudospin and mixed susceptibilities.

On the other hand, the strong-coupling approach allows us to calculate all susceptibilities. To do this, we calculate the irreducible many-particle Green’s functions

$`Tn_\sigma (\tau )a_\sigma ^{}^{}(\tau ^{})a_\sigma ^{}(\tau ^{\prime \prime })_{\text{irr}}\delta _{m0}{\displaystyle \underset{\nu ^{}}{}}\mathrm{\Lambda }_\sigma (\omega _\nu ^{})\mathrm{\Lambda }_\sigma ^{}(\omega _\nu )`$ (66)
$`\delta _{\sigma \sigma ^{}}\left[\mathrm{\Lambda }_\sigma (\omega _\nu )\mathrm{\Lambda }_\sigma (\omega _{\nu +m})+G_\sigma ^{(a)}(\omega _\nu )G_\sigma ^{(a)}(\omega _{\nu +m})\right],`$ (67)
$`TS^z(\tau )a_\sigma ^{}(\tau ^{})a_\sigma (\tau ^{\prime \prime })_{\text{irr}}`$ $``$ (69)
$`\delta _{m0}\sqrt{P^+P^{}}\mathrm{\Lambda }_\sigma (\omega _\nu ),`$
$$TS^z(\tau )S^z(\tau ^{})_{\text{irr}}\delta _{m0}P^+P^{},$$
(70)
$`Tn_\sigma (\tau )n_\sigma ^{}(\tau ^{})_{\text{irr}}\delta _{m0}{\displaystyle \underset{\nu }{}}\mathrm{\Lambda }_\sigma (\omega _\nu ){\displaystyle \underset{\nu ^{}}{}}\mathrm{\Lambda }_\sigma ^{}(\omega _\nu ^{})`$ (71)
$`\delta _{\sigma \sigma ^{}}{\displaystyle \underset{\nu }{}}\left[\mathrm{\Lambda }_\sigma (\omega _\nu )\mathrm{\Lambda }_\sigma (\omega _{\nu +m})+G_\sigma ^{(a)}(\omega _\nu )G_\sigma ^{(a)}(\omega _{\nu +m})\right]`$ (72)
$`TS^z(\tau )n_\sigma (\tau ^{})_{\text{irr}}`$ $``$ (74)
$`\delta _{m0}\sqrt{P^+P^{}}{\displaystyle \underset{\nu }{}}\mathrm{\Lambda }_\sigma (\omega _\nu ),`$
After substituting (67)–(74) into Eq. (40) we get for the charge susceptibility the same expression (56), whereas for the pseudospin and mixed susceptibilities

$`\chi _𝒒^{S^zS^z}(\omega _m)`$ $`=`$ $`\delta _{m0}{\displaystyle \frac{\mathrm{\Delta }_{S^z}^2}{T\mathrm{\Theta }(T,𝒒)}}`$ (75)
$`=`$ $`\delta _{m0}{\displaystyle \frac{P^+P^{}}{T\mathrm{\Theta }(T,𝒒)}},`$ (76)
$$\chi _𝒒^{nS^z}(\omega _m)=\chi _𝒒^{S^zn}(\omega _m)=\delta _{m0}\frac{\mathrm{\Delta }_{S^z}\mathrm{\Delta }_n}{T\mathrm{\Theta }(T,𝒒)},$$
(77)
where

$$\mathrm{\Delta }_{S^z}=\sqrt{P^+P^{}}.$$
(78)
Expression (63) for $`\mathrm{\Theta }(T,𝒒)`$ coincide with the one obtained by Freericks  from the equations for the static susceptibilities ($`\omega _m=0`$) derived by Brandt and Mielsch . Irreducible four vertex (59) was also derived by Freericks and Miller  within the Baym-Kadanoff formalism and used to find the exact solution for the nonresonant Raman scattering for the FK model .

## V Discussion

For the binary alloy (Falikov-Kimball) model (43) expressions (56), (75) and (77) define the so-called isothermal susceptibilities . Furthermore, the pseudospin (75) and mixed (77) susceptibilities are only static (with factor $`\delta _{m0}`$) because the pseudospin operator $`S_i^z`$ commutes with the Hamiltonian and is an integral of motion. It should be noted, that the binary alloy model (43) can be reduced to the Ising type model with an effective multisite retarded pseudospin interactions by taking trace of the statistical operator with the Hamiltonian (43) over electron (fermion) variables . These explains the Ising type expression obtained for the pseudospin susceptibility (75), but now the expression for the critical temperature $`\mathrm{\Theta }(T,𝒒)`$ is more complicated.

The expression for the charge susceptibility (56) contains two terms. The first one is static and can be written as

$$\chi ^{nS^z}\left[\chi ^{S^zS^z}\right]^1\chi ^{S^zn}.$$
(79)
This describes the contribution from the pseudospin subsystem to the charge susceptibility. It gives the main contribution to the static susceptibilities. The second term $`K_𝒒^{nn}(\omega _m)`$ in (56) gives the pure electron response and describes the so-called isolated (Kubo) susceptibility . In general, terms with the factor $`\delta _{m0}`$ give the difference between the isothermal and isolated susceptibilities (see Appendix in Ref. ).

## Acknowledgement

This work was partially supported by the Fundamental Researches State Fund of the Ministry of Ukraine for Science and Education (Project No 02.07/266).


format: text

text:
# Two-Band Superconductivity in MgB2

\[

## Abstract

The study of the anisotropic superconductor MgB2 using a combination of scanning tunneling microscopy and spectroscopy reveals two distinct energy gaps at $`\mathrm{\Delta }_1`$=2.3 meV and $`\mathrm{\Delta }_2`$=7.1 meV. Different spectral weights of the partial superconducting density of states (PDOS) are a reflection of different tunneling directions in this multi-band system. Our experimental observations are consistent with the existence of two-band superconductivity in the presence of interband superconducting pair interaction and quasiparticle scattering. Temperature evolution of the tunneling spectra follows the BCS scenario  with both gaps vanishing at the bulk Tc. Indeed, the study of tunneling junctions exhibiting only the small gap (c-axis tunneling) clearly and reproducibly show that this gap persists up to the bulk Tc. The data confirm the importance of Fermi-surface sheet dependent superconductivity in MgB2 proposed in the multigap model by Liu et al.  .

\]

The discovery of superconductivity in MgB2  at 39K sparked great interest in the fundamental physics and practical applications of this material. There has already been rapid progress in understanding the physical properties of this superconductor. Specific heat measurements  show that MgB2 is an s-wave superconductor and the presence of the isotope effect  points towards phonon-mediated pairing. Tunneling and photoemission spectroscopy directly measures the superconducting energy gap and can provide further understanding of the origin of the superconductivity in this material. Earlier tunneling spectroscopy measurements show a large spread in the gap values  each consistent with the BCS form. More recent experiments, including STM tunneling spectroscopy , point-contact spectroscopy , specific heat measurements , and Raman spectroscopy  point towards the existence of two distinct gaps. This scenario has been predicted theoretically by Liu et al. . First principle calculations show that the Fermi surface of MgB2 consists of 2D cylindrical sheets arising from $`\sigma `$ antibonding states of B $`p_{xy}`$ orbitals, and 3D tubular networks arising from $`\pi `$ bonding and antibonding states of B $`p_z`$ orbitals. In this theoretical framework  two different energy gaps exist, the smaller one being an induced gap associated with the 3D bands and the larger one associated with the superconducting 2D bands. Furthermore both superconducting gaps should vanish at the bulk critical temperature Tc. Due to this highly anisotropic band structure the superconducting gaps should be momentum-dependent reflecting the strength of the electron-phonon coupling of the carriers in the different bands. Up to now there has been no direct experimental evidence of the orientation dependence of the order parameter in this material. Moreover, the temperature dependence of the two gaps would give further insights into the nature of superconductivity in MgB2. Scanning tunneling spectroscopy is a unique technique that allows direct measure of the DOS near the Fermi energy with high spatial and energy resolution. Since the tunneling current couples only to bands with an appreciable Fermi velocity in the tunneling direction STM can distinguish between the 2D bands with little or no c-axis component of Fermi velocity and 3D bands with considerable c-axis component. Thus, in the clean limit($`l2\pi \xi `$) the direction of the tunneling current selects which band and which superconducting gap is probed. In this paper we report direct evidence of orientation-dependent double-gap structure in the quasi-particle energy spectra as determined from tunneling spectroscopy. Our data are consistent with the theoretical predictions of multigap superconductivity in the clean limit. Temperature dependent tunneling shows that the two distinct gaps vanish simultaneously near the bulk Tc.

Compact samples of MgB2 were synthesized from amorphous B powder (4N’s purity) and high purity Mg. The B powder was pressed into pellets under 6 kbar pressure. These pellets were reacted with Mg vapor at 850 C for 2 hr in a BN container under 50 bar of Ar. During the reaction the pellets broke up in pieces of several mm in size. The typical critical temperature of these gold-colored pellets is 39K.

The tunneling measurements were performed on these samples using different surface preparation methods. One procedure consisted of chemically etching the pellets for 50 seconds in bromine (Br 1% in pure ethanol), then rinsing in pure ethanol and drying in N2 gas. Alternative procedures consisted of cleaving or mechanical polishing of the sample in an inert atmosphere. After these treatments the samples were mounted on the STM stage in He exchange gas and cooled down to 4.2 K. These procedures all yielded equivalent tunneling results. The measurements were performed with a home-built STM described elsewhere .

The typical conductance spectra recorded on different grains within the sample are shown in Figure 1. The typical grain size is of the order of 50 nm and the spectra are remarkably reproducible within the same grain both with respect to location and tunneling resistance. All the spectra are normalized to the conductance value at –20mV. They reveal a double gap structure, flat background, and very low zero-bias conductance with very little broadening other than thermal smearing. Two peaks are present at V$`{}_{1}{}^{}`$ 3 mV and V$`{}_{2}{}^{}`$ 7.5 mV, symmetrically for both injection and emission of electrons and they stay within 10% of these values when changing between grains and/or samples. This gap size distribution is in line with the proposed distribution of the gap values over the Fermi surface by Choi et al. .

Figure 2(a) is the representative surface topography (150 nm x 150 nm) showing several different grains. The roughness within an individual grain is about 1nm and the grains are separated by steps that are 1-10nm high. Each of the four grains is characterized by a tunneling spectrum that is consistent within a particular grain. Within this scanning area we focus on the tunneling spectra acquired on grains 1 and 2 (Fig. 2a). The spectrum recorded in grain 2 shows a clear two-gap structure dominated by the high-energy gap while in grain 1 only the smaller gap is pronounced and the high-energy peak becomes a weak satellite feature. The main characteristics of the spectra, including the spectral weight of each peak, do not change with the tunneling resistance from 0.1 to 2 G$`\mathrm{\Omega }`$. Having this remarkable spatial reproducibility of the tunneling conductance we focus our attention on the origin of the different spectra in these grains. First, we point out that both spectra are in agreement with the multi-gap scenario proposed by Liu et al.  where in the clean limit the impurity and surface scattering do not significantly influence the spectra. Second, we find a striking similarity between the spectrum in grain 1 and the typical spectrum we recorded earlier on c-axis oriented films . In the c-axis tunneling on MgB2 epitaxial films the contribution from the 3D Fermi surface is expected to dominate the tunneling conductance. Thus we believe that the distinct spectra attributable to grain 1 and 2 are a reflection of different crystallographic grain orientation with respect to the tunneling direction. This evidence of momentum-dependent tunneling is further supported by the spatial and temperature-dependent analysis of the spectra presented below.

To elucidate the origin of the different spectra in the two grains we studied the spatial evolution of the spectra across the grain boundary. In Figure 2(b,c) we show the conductance maps at the superconducting peak values of V1=+3.3 mV and V2=+7.3 mV in the region indicated by the square spanning the two grains. In these maps the bright zones mark regions where the tunneling spectra show higher conductance at the energy corresponding to V1 (Fig. 2(b)) or V2 (Fig. 2(c)), respectively. The transition from one type of spectra to the other evolves over a length scale on the order of 5 nm. We associate these two different types of spectra with vacuum tunneling into grains with two distinct orientations. The smooth transition between the spectra across the grain boundary on the length scale of $`\xi `$ is further evidence of the sample being in the clean limit.

The distinct tunneling conductance spectra that we observe are consistent with the two-gap BCS model  taking into account interband impurity scattering (see, e.g., ). According to this model the partial densities of states (PDOS) of two bands $`N_{1(2)}(\omega )`$

$`N_{1(2)}(\omega )=N_{1(2)}(0)Re\left[{\displaystyle \frac{u_{1(2)}}{\sqrt{u_{1(2)}^21}}}\right],`$are determined by two dimensionless functions $`u_{1(2)}(\omega )`$ that obey equations

$`u_1\mathrm{\Delta }_1`$ $`=`$ $`\omega +i\mathrm{\Gamma }+i\mathrm{\Gamma }_{12}{\displaystyle \frac{u_2u_1}{\sqrt{u_2^21}}},`$
$`u_2\mathrm{\Delta }_2`$ $`=`$ $`\omega +i\mathrm{\Gamma }+i\mathrm{\Gamma }_{21}{\displaystyle \frac{u_1u_2}{\sqrt{u_1^21}}}.`$
Here $`\mathrm{\Delta }_{1,2}`$ are the gap parameters of two bands, $`\mathrm{\Gamma }_{12}`$ and $`\mathrm{\Gamma }_{21}`$ are the interband scattering rates, with $`\mathrm{\Gamma }_{12}/\mathrm{\Gamma }_{21}=N_2(0)/N_1(0)`$. It is important to note that strong interband scattering suppresses the superconducting transition temperature. In particular, at $`\mathrm{\Gamma }_{21}T_c`$ and $`\mathrm{\Delta }_1\mathrm{\Delta }_2`$, $`\delta T_c\pi \mathrm{\Gamma }_{21}/8`$ (see, e. g., ).

We use a linear combination of the two PDOS to obtain the total DOS

$`N_{tot}(\omega )=\alpha N_1(\omega )+\beta N_2(\omega )`$

that we compare with the experimental tunneling conductance spectra. The PDOS’ spectral weight represented by the parameters $`\alpha `$ and $`\beta `$ ($`\alpha +\beta `$=1) is associated with the relative orientation of the crystallographic axis of a grain with respect to the tunneling barrier. Thus in the fits shown below the parameters $`\alpha `$ and $`\beta `$ are kept constant for all spectra obtained on a particular grain at different temperatures.

The linear combination of the two PDOS shown above accounts for all major features in our tunneling spectra including the temperature dependent data. The evolution of the tunneling spectra in the temperature range between 4.2K and 42K are displayed in Figure 3. The two-gap signature in the spectra of grain 1 becomes hardly distinguishable with increasing the temperature (Fig. 3(a)). The small gap shows a rapid decrease in the peak height and the conductance peaks of the two gaps merge into a broad peak at T$``$10K. In order to separate the temperature evolution of each PDOS we used the temperature dependence of the spectrum on grain 2 (Fig. 3(b)) where only the small gap is present. Each curve was normalized to the conductance value at -20 mV and then compared with the theoretical curve at the same temperature. The theoretical curves were obtained by fixing the interband scattering rates $`\mathrm{\Gamma }_{12}`$, $`\mathrm{\Gamma }_{21}`$ and the smearing parameter $`\mathrm{\Gamma }`$ at values optimized for the 4.2 K conductance spectra. In all the fits we put $`\mathrm{\Gamma }_{12}`$= 0.73 $`\mathrm{\Gamma }_{21}`$ since their ratio is proportional to the ratio of the density of states at the Fermi level of the 3D and the 2D bands . The theoretical curves reproduce very well the shape inside the gap but not the peak heights. By adjusting the parameters $`\mathrm{\Gamma }_{12}`$ and $`\mathrm{\Gamma }`$ we could reproduce very well the peak heights but the zero bias value would be higher than the experimental one. Fits to the two set of parameters have been used to get an estimate of the error bars for the two gap values: $`\mathrm{\Delta }_1`$= 2.3$`\pm `$0.2 mV and $`\mathrm{\Delta }_2`$= 7.1$`\pm `$0.2. A very good fit of both peaks’ height and zero bias value can be obtained only for interband scattering values that are very large ($``$2.5 mV) and values of the 3D gap is extremely small. According to  such a large interband scattering value should lower the critical temperature by $``$9 K which contradicts our experimental findings. The set of temperature dependent tunneling spectra clearly shows that both superconducting gaps exist up to 37K.

Next, we turn our attention to the spectral weights of the two PDOS to the total DOS. The particular choice of the two grains in Fig. 2 was based on the peculiarity of the conductance spectra of these crystallites. They closely represent the two extreme limits of partial contribution by each PDOS that we have observed so far on a large number of pellets and thin films. The theoretical fits obtained by our model in these two limits (Fig. 3a,b) show that the contribution of the 3D PDOS varies from nearly 100% (Fig. 3b and ) down to 87% (Fig. 3a). Therefore it appears that momentum averaged STM tunneling conductance is dominated by the 3D band in this momentum dependent two-band tunneling process. This could be due to the specific geometry of the tunneling cone emanating from the STM tip and its projection on the highly anisotropic Fermi surface of the MgB2 crystallite. Obviously, the relative weights of the PDOS are temperature independent and the parameters $`\alpha `$ and $`\beta `$ remain fixed when the temperature dependent fits of the tunneling spectra are obtained.

The values of $`\mathrm{\Delta }_1`$(T) and $`\mathrm{\Delta }_2`$(T) are extracted from the theoretical curves and the results are reported in Figure 3(c). The high-energy gap follows a BCS behavior. The smaller gap remains constant up to 20K and follows the BCS predictions with a critical temperature far beyond the one expected for such a small gap (T$`{}_{c}{}^{}`$14K) with a ratio 2$`\mathrm{\Delta }`$/kT$`{}_{c}{}^{}`$1.4 similar to earlier findings from point spectroscopy . Such behavior clearly demonstrates that the gap in the 3D band is mainly induced by the interband pairing interaction rather than produced by pairing interaction inside this band, in agreement with theoretical expectations . This experimental evidence directly excludes the possibility that the grains or their surfaces are associated with different stoichiometry. The high-energy gap follows a BCS behavior with a ratio $`2\mathrm{\Delta }_2(0)/kT_c4.3`$, significantly larger than the BCS value of 3.52. There are two possible sources for this increase. Usually such enhancement is attributed to strong coupling effect. However, this ratio also increases due to the presence of additional superconducting bands even in the weak coupling limit. In the case of $`\mathrm{\Delta }_1\mathrm{\Delta }_2`$ one can derive a simple formula from the two-band BCS theory

$`{\displaystyle \frac{2\mathrm{\Delta }_2}{T_c}}3.52\left(1+{\displaystyle \frac{N_1\mathrm{\Delta }_1^2}{N_2\mathrm{\Delta }_2^2}}\mathrm{ln}{\displaystyle \frac{\mathrm{\Delta }_2}{\mathrm{\Delta }_1}}\right).`$Using $`N_2=0.73N_1`$ from  and experimental values of the gaps, we estimate from this formula $`2\mathrm{\Delta }_2/T_c4.1`$. Therefore, we conclude that the main part of the $`2\mathrm{\Delta }/T_c`$ increase originates from the multiband nature of the superconductivity in this compound.

Finally, we would like to point out that in the case of two-band superconductivity in the dirty limit ($`l\xi )`$ the conductance spectra are very broad and similar to those reported earlier . These spectra are mostly obtained on pellets without any surface treatment. Impurities near the surface may provide sufficiently strong interband scattering to smear out the two-gap features in this case.

In conclusion, the existence of two gaps has been observed by tunneling spectroscopy in compact pellets of MgB2. Tunneling spectra show different ratios between the two gap heights on different grains within the sample. This result can be interpreted in terms of different tunneling directions with respect to the crystallographic orientation of the grain, therefore supporting the two-gap scenario . Different tunneling direction should probe the two bands with different weights. Moreover, the two gaps both vanish at a temperature close to the bulk value. The observed temperature independence of the small gap up to T$``$ 20 K implies that no transition temperature corresponding to the BCS ratio 2$`\mathrm{\Delta }_1`$=3.52 kTc is present. This suggests that the coupling between the two bands is moderately strong and that the superconductivity is dominant in the 2D sheets.

This work was supported by US DOE Basic Energy Science - Material Science under contract No. W-31-109-ENG-38.


text:
# Transitions to Measure Synchronization in Coupled Hamiltonian Systems

## Abstract

Transitions to measure synchronization in two coupled $`\varphi ^4`$ lattices are investigated based on numerical simulations. The relationship between measure synchronization (MS), phase locking and system’s total energy is studied both for periodic and chaotic states. Two different scalings are discovered during the process to MS according to phase locking. Random walk like phase synchronization in chaotic measure synchronization is found, and phase locking interrupted by phase slips irregularly is also investigated. Meanwhile, related analysis is qualitative given to explain this phenomenon.
PACS numbers: 05.45.Xt, 05.45.Pq, 05.45. Ac, 05.45.-a

The phenomenon of synchronization has been studied since 17th , but the related research never stopped. Recently, in order to explore the collective behavior in complex systems, this topic has been paied more attention than ever before and became a hot focus in many science fields . Although synchronization is first mentioned in periodic systems, synchronization in chaotic systems attracted more attentions due to its rich dynamical properties and potential applications in practical systems . During the research process, the original concept has been largely extended, many new phenomena had been discovered, phase synchronization , lag synchronization , generalize synchronization  had been the familiar words in new related works. Till to now, most attentions had been concentrated on dissipative systems, the related research for conservative system, the hamiltonian system, has been neglected and only few works appeared . Hamiltonian system stands as the central topic of both classical and quantum mechanics, and a large variety of practical systems can be well approximated by Hamiltonian formalism even at weak dissipation. Thus, it is important to understand the behavior of Hamiltonian systems in their possible synchronization process.

The former works abut Hamiltonian systems are mainly about their collective behavior in globally coupled high dimension systems or studied their transitions to equipartition from non-equilibrium state based on statistical methods , few paper had concerned to investigate their detail dynamical behaviors during the system’s evolutions. The reason mainly because that Liouvill’s theory request the phase space be conservative, it prevents the full collapse of the orbits which related with complete synchronization and phase synchronization, this means the classical concept for synchronization developed in dissipative system is unsuitable for this kind of system. But this situation has been changed in 1999, a new phenomenon, measure synchronization (MS), for Hamiltonian systems is discovered by Hampton and Zenette . In their work, MS is defined as all orbits cover the same region in the phase space with identical invariant measures, and the transition to this state is characterized by a critical coupling strength. Although MS has been found and critical coupling is defined, that are still many questions still leaving unsolved. For example, the route from desynchronization to MS is unclear, and how about system’s variables behave and what’s the innate character for MS still need to be cleared. Also, the relation between chaos and MS need to be considered again, especially, the dynamical behavior for MS in chaotic state is also valuable to study carefully. In this paper, we will concerned mainly on these questions, and try to explore the way from non-MS to MS in coupled Hamiltonian systems.

Let’s consider two coupled classical $`\varphi ^4`$ lattices, the system’s Hamiltonian reads

$$H=\frac{p_1^2+p_2^2}{2}+\frac{q_1^4+q_2^4}{4}+\frac{\epsilon }{2}(q_1q_2)^2$$
(1)
the direct physical image is two coupled point masses each of them oscillated in an quartic potential. Due to the quartic potential, system could displays complex behavior and may express new property out of intuition. By changing coupling strength $`\epsilon `$ and initial values, both period and chaos can be realized. The corresponding canonical equations are

$`\underset{1}{\overset{.}{q}}`$ $`=`$ $`p_1,\text{ }\underset{1}{\overset{.}{p}}=q_1^3+\epsilon (q_2q_1)`$ (2)
$`\underset{2}{\overset{.}{q}}`$ $`=`$ $`p_2,\text{ }\underset{2}{\overset{.}{p}}=q_2^3+\epsilon (q_1q_2)`$ (3)
In this model, the system’s total energy is conserved, $`E(t)=\frac{p_1^2(t)+p_2^2(t)}{2}+\frac{q_1^4(t)+q_2^4(t)}{4}+\frac{\epsilon }{2}(q_1(t)q_2(t))^2=E(0)=`$Const, and initial values act as an independent additional parameters. If we set each lattice’s initial positions at zero, then total energy will be independend to coupling strength, this will be used in the following.

Once two lattices reached MS, then any macroscopic values will equal if only the average time long enough. One example is the energy of each lattice, defined as

$$e_i=\frac{1}{T}_0^T[p_i^2(t)/2+q_i^4(t)/4]𝑑t,\text{ }i=1,2$$
(4)
before MS, two orbits occupy different domains which never intersect, of course $`e_1e_2`$, this situation continued till MS reached, after the critical coupling, two lattice’s energy be equaled due to their identical invariant measures. In Fig. 1 we choose initial condition

$`q_1`$ $`=`$ $`0,\text{ }q_2=0,`$ (5)
$`p_1`$ $`=`$ $`0.1,\text{ }p_2=0.2`$ (6)
and plotted each lattice’s energy $`e_i`$ versus $`\epsilon `$. During our numerical simulations, third order bilateral symplectic algorithm  and small time step $`h=0.001`$ are used to keep the system’s conservative property. The critical coupling, $`\epsilon _c=0.0032,`$ which separate non-MS and MS is clearly shown in this figure. Meanwhile, we also plotted two lattice’s total energy $`e=e_1+e_2`$ versus $`\epsilon `$ in the same figure for comparison. It is shown that energy $`e`$ has its minimum on the critical point $`\epsilon _c`$, below and above this point $`e`$ decreased and increased monotonically. This made us reminisced phase synchronization which has been warmly studied recently , made it interesting to investigate the phase locking, or named phase synchronization, in this system. In Hampton’s work, the phase difference $`\mathrm{\Delta }\theta =\theta _1\theta _2,`$ with $`\theta _{1,2}`$ denotes each orbit’s angle and defined on the real line rather than on the circle, between two orbits is found to obey different rules. Before MS, $`\mathrm{\Delta }\theta `$ increases linearly with time goes on, and above MS, the difference oscillates with a characteristic frequency which around a constant value. But only these can’t explain the minimum, the deep reason need to be explored.

More information about MS and phase locking can be found in Fig. 2, in Fig. 2(a) and (b) we adopted $`\epsilon =0.002<\epsilon _c,`$ two orbits stay in different domain, and phase difference increases linearly with speed $`\eta =0.1`$. While $`\epsilon =0.0035>\epsilon _c,`$ MS reached and phase difference is oscillates within $`(\pi ,\pi )`$ with frequency $`f=2.12`$, as what is shown in Fig. 2(c) and (d). In order to explore the relation between phase difference and coupling strength, we plot $`\eta `$, which is defined as the speed of phase difference between two orbits before MS, and $`f,`$ which is defined as the oscillate frequency of phase difference after MS, versus $`\epsilon `$ in Fig. 2(e) and (f). It is observed that, before MS, $`\eta `$ decreases gradually as $`\epsilon `$ increases and trend to zero on $`\epsilon _c`$, while $`f`$ increases gradually from zero as $`\epsilon `$ increases from $`\epsilon _c`$, this is coincide to Ref. .

During our numerical simulations, we found that phase locking exists two different trends corresponding to non-MS and MS respectively, future study shows that its right the reason why the energy $`e`$ show its minimum at $`\epsilon _c.`$ Before the critical point, MS doesn’t exist and orbits have different main frequencies, but we found that: except $`\eta `$ decreases as $`\epsilon `$ increases, the distribution of phase difference $`\phi `$ (model $`\mathrm{\Delta }\theta `$ with $`2\pi `$) has a obviously trend. In Fig. 3(a) we plotted the distributions for $`\phi `$ under three different coupling which all below the critical point. From this figure, we can find that as $`\epsilon `$ increases, the probability for $`\phi =\pi `$ changes to higher, we remark here that $`\pi `$ is the biggest phase difference for two orbits after modeled $`\theta _1`$ and $`\theta _2`$ with $`2\pi .`$ In Fig. 3(b) we also plotted the probability distributions for distance difference $`\mathrm{\Delta }q=|q_1q_2|,`$ it is found that with $`\phi `$ trends to $`\pi `$ in distribution, the probability for large $`\mathrm{\Delta }q`$ increases as well. After $`ϵ>\epsilon _c,`$ as coupling increases, the oscillate frequency rise and the amplitude decrease, meanwhile, the distribution for $`\mathrm{\Delta }q`$ contracted to small scope.

Already had these ideas in mind, the reason for the minimum of energy $`e`$ on the critical point is straightforward. As the coupling increases, the relation between two phases become stronger, although they are unlocked, but they have the trend to reach phase locking, this can be found from $`\eta `$ decreases as $`\epsilon `$ increase. As we already known, for Hamiltonian system, the classical phase synchronization can’t be reached due to its conservative character, but as coupling increases, two orbits trend to reach it. On the way to phase synchronization, the first $`\mathrm{\Delta }\theta `$ should be $`\pi `$, because $`\pi `$ stands for the largest phase difference if we model phase by $`2\pi .`$ This directly made the probability of large $`\mathrm{\Delta }q`$ increased and thus induced energy $`e`$ decrease in this process. Advanced increase the coupling, $`\epsilon >\epsilon _c,`$ then two orbit become MS and be phase locked, under this situation, the coupling strength is strong enough to reach MS, but is weaker to bind the phase oscillate in a small amplitude, continue increase $`\epsilon `$, the amplitude will drops from large value $`\pi `$, where $`\epsilon `$ slightly larger than $`\epsilon _c,`$ to zero. During this process, phases are locked and the oscillate amplitude decreased, induced the distribution for $`\mathrm{\Delta }q`$ inclined to short distance, at last made energy $`e`$ tend upward. In order to give a more convincing expression, we plot the probability of $`\pi `$ (denote by $`\rho _\pi `$) versus $`\epsilon `$ in Fig. 4(a), $`\mathrm{log}(\rho _\pi )`$ for vertical axis and $`\mathrm{log}(\epsilon _c\epsilon )`$ for horizontal axis. A fine scaling is found

$$\rho _\pi =\frac{\chi }{(\epsilon _c\epsilon )^\nu }$$
(7)
with $`\nu =1.962.0`$ and $`\chi `$ denotes as constant. Also, in Fig. 4(b) we plotted amplitude $`A`$ versus $`\epsilon `$ in MS region, this time the scaling is found to be a exponential mode

$$A=A_0+Be^{\frac{\epsilon }{\alpha }}$$
(8)
with $`\alpha =1.0\times 10^3`$ and $`A_0,`$ $`B`$ constants.

Under the situation mentioned above, $`\epsilon _c`$ not only stands as the critical point for MS and the energy minimum, but also stands as the critical point for two different trends of phase locking, i.e., before MS, two orbits are trend to achieve the first possible locking $`\mathrm{\Delta }\theta =\pi ,`$ after MS, the amplitude of oscillation dropped from $`\pi `$ and trend to approach zero. But we should also noticed here that all what we have discussed is just within period, no chaos appeared, so a natural question is how about these phenomena within chaos, does different phase locking trend still exist, and what’s the new characters for MS when system is chaos. All these questions aroused our interests to investigate MS in coupled chaotic Hamiltonoian systems.

In order to get chaos, this time we choose initial condition

$`q_1`$ $`=`$ $`0,\text{ }q_2=0,`$ (9)
$`p_1`$ $`=`$ $`0.1,\text{ }p_2=0.4`$ (10)
Similar to Fig. 1, we plotted energy $`e_i`$ and total energy $`e`$ versus coupling $`\epsilon `$ in Fig. 5(a), different to the above result, energy $`e_i`$ behaves a disordered way in chaotic regions, made it difficult to get the critical point for MS. Through tracing the time evolution of their orbits, we found that the uneven behavior of $`e_i`$ is closely related to chaos, in other word, its just the conjoint action between chaos and MS that lead to this irregular behavior. On one hand, the definition of MS requests each lattice’s energy be equaled, just as what happened in periodic state, but on the other hand, the internal stochastic character keep their orbits moving in an irregular way. Thus, $`e_1=e_2`$ can be reached only under limit $`T\mathrm{}.`$ However, we found the behavior of summation $`e`$ is more regular than $`e_i`$, despite of some abrupt windows. Through calculate system’s largest Lyapunov exponent $`\lambda `$, it is found that these windows are corresponding to period which embed in chaos. In Fig. 5(b) we plotted $`e`$ and $`\lambda `$ versus $`\epsilon `$ in the same figure, compare with Fig. 5(a), it is shown that $`e`$ decreases more smoothly then $`e_i`$ in chaotic regions. This result is somewhat similar to the behavior of $`e`$ in region $`\epsilon <\epsilon _c`$ plotted in Fig. 1. As we have known, the trends of phase locking plays important role for the minimum of $`e`$ in periodic state, so we may ask does these the same to chaos, if not, what’s the reason made $`e`$ decrease in chaotic state and how about the phase locking now?

On important difference between this work and Hampton’s is, in our model, chaos is closely related with MS, and the largest Lyapunov exponent is a suitable order parameter for MS. We choose $`\epsilon =0.025>0.012=\epsilon _c`$, with initial condition Eqs. 7, system is chaos with $`\lambda =0.08`$. Phase difference $`\mathrm{\Delta }\theta `$ versus time $`t`$ under this coupling strength is plotted in Fig. 6. It is shown that $`\mathrm{\Delta }\theta `$ behaves a random walk process around zero, and there also exist some platforms with different wide. The inset figure is enlargement of the rectangle region in original figure, there the details about the platform is shown. In fact, each platform is in the state of phase locking, in these regions, although $`\mathrm{\Delta }\theta `$ oscillate all the time, but its amplitude is small than $`\pi `$. We also found that even in regions where phase are unlocked, the phase difference shows a $`2\pi `$ like slips. The same phenomenon also can be found in others chaotic states with different $`\epsilon >\epsilon _c`$, the only difference is, as $`\epsilon `$ increases, the random walk more near zero, and more time two orbits be phase locked.

The $`2\pi `$ slips in phase locking made us reminisced the phase slips in phase synchronization of dissipative systems  where various models were used and rich phenomena were discovered. But in our opinion, these results are unsuitable to Hamiltonian system due to its conservative property, here we only give a qualitative analysis for this phenomenon without proofs. When system is MS and within chaotic state, the coupling is stronger and trend to lock phases, this is exhibited by the platforms in Fig. 6(a), but due to the stochastic property of chaos, made the locking unstable, this is exhibited by the phase slips. Because the first possible locking is near $`\pi ,`$ and sometimes $`\pi `$ locking is unstable, made $`\phi `$ stays a long time within $`(0,\pi )`$ (phase locking), then broke suddenly because the stochastic property of chaos, after $`2\pi `$ increase, $`\phi `$ travels again near $`\pi ,`$ this time it has two choices, continue this slip if its still local unstable, or stay on locking state except local stable be changed. In our studies, we found that, in chaotic state, although phase locking doesn’t exist, but $`\phi `$ still owns the trend to be $`\pi `$ locked as $`\epsilon `$ increases, and the distribution for $`\mathrm{\Delta }q`$ also incline toward long distance, these properties made $`e`$ decreased in chaotic state with $`\epsilon `$ increases.

In conclusion, we have investigated the way from non-MS to MS in coupled $`\varphi ^4`$ lattices both for period and chaos. Due to Hamiltinian system’s conservative property, exactly phase locking, $`\mathrm{\Delta }\theta =0`$, is unreachable, but we found that, in periodic state, $`\phi `$ trend to $`\pi `$ as $`\epsilon `$ approach $`\epsilon _c`$, after $`\epsilon >\epsilon _c,`$ $`\mathrm{\Delta }\theta `$ is oscillate around zero and the related amplitude decreases with $`\epsilon `$ increases, on the critical point $`\epsilon _c`$, energy $`e`$ own its minimum, also, different scalings has been found to describe these trends. MS is found to be connected closely with chaos in this model, random walk like phase locking is discovered, after analysed system’s dynamical behaviors, we gave the qualitative explaining for this phenomenon. We argued that in some models, like $`\varphi ^4,`$ Lyapunov exponent can be a suitable parameter for MS. But there are still some works ahead waiting for future research. The first work is to explain the phase locking and slip in chaotic MS strictly, we noticed that a recent work about phase synchronization also shown a random walk like phase locking , its gave an example of nontrivial interaction between microscales and macroscales, so another work is to investigate the scalings for phase locking in chaotic states, also, to study these phenomena in large systems in also a interesting work. We wish this paper could deepen people’s understanding of MS, and the observed phenomenon could potentially be investigated in many other related research in coupled Hamiltonian systems, e.g., elliptical galaxies, Hamiltonian version Kuramoto’s model, transitions to equipartition in coupled large Hamiltonian systems, and we also hope this work could helpful in the research of quantum chaos.

This research was supported by the National Natural Science Foundation of China, the Nonlinear Science Project of China, and the Foundation of Doctoral training of Educational Bureau of China.

Captions of figures

Fig. 1 The time average energy for each lattice $`e_i=\underset{T\mathrm{}}{lim}\frac{1}{T}_0^T(\frac{p_i^2}{2}+\frac{q_i^2}{4})dt,`$ $`i=1,2`$ versus coupling strength $`\epsilon ,`$ critical coupling for MS is $`\epsilon _c=0.0032.`$ Total energy $`e=e_1+e_2`$ is plotted on the left vertical axis (denoted by hollow circle), on the critical point, $`e`$ has its minimum, before and above $`\epsilon _c`$, $`e`$ decreases and increases monotonically.

Fig. 2 With the initial condition Eqs. 4, when $`\epsilon =0.002<\epsilon _c`$, (a) plot the phase difference $`\mathrm{\Delta }\theta `$ versus $`t`$, the slope rate $`\eta =0.1,`$ (b) plot phase orbits for two lattices, it is shown that two orbits moving in different domains without overlap. (c) and (d) corresponding to $`\epsilon =0.0035>\epsilon _c`$, $`\mathrm{\Delta }\theta `$ oscillate with frequency $`f=2.12,`$ and two orbits be MS. (e) plote $`\eta `$ versus $`\epsilon `$ for $`\epsilon <\epsilon _c`$ and (f) plot $`f`$ versus $`\epsilon `$ for $`\epsilon >\epsilon _c`$.

Fig. 3 (a) Probability distributions of $`\phi `$ for three different $`\epsilon `$s, circle for $`\epsilon =0.001,`$ upward triangle for $`\epsilon =0.002`$ and downward triangle for $`\epsilon =0.025`$, all below $`\epsilon _c.`$ (b) Probability distributions for distance difference $`\mathrm{\Delta }q,`$ the couplings same to (a).

Fig. 4 (a) The $`\pi `$ distribution of $`\phi `$ versus coupling for $`\epsilon <\epsilon _c`$, after logarithm, all data lie on one line, the related scaling reads $`\rho _\pi =\frac{\chi }{(\epsilon _c\epsilon )^\nu }`$ with $`\gamma =1.96.`$ (b) Oscillation amplitude $`A`$ versus $`\epsilon ,`$ scaling relation reads $`A=A_0+Be^{\frac{\epsilon }{\alpha }}`$ with $`\alpha =1\times 10^3.`$

Fig. 5 Initial condition as Eqs. 7. (a) The summational energy $`e`$ and each lattice’s energy $`e_i`$ versus $`\epsilon `$, hollow circle represents the summation. (b) $`e`$ and the largest Lyapunov exponent versus $`\epsilon .`$ It can be found that in chaotic MS regions the summation decreases as $`\epsilon `$ increases.

Fig. 6 $`\epsilon =0.025,`$ the time evolution of phase locking in chaotic MS, $`\lambda =0.08.`$ The inset window is enlargement of the rectangle region in original figure, phase locking and slips can be found.


text:
# Super-Eddington Fluxes from Thin Accretion Disks?

## 1. Introduction

Magnetized, radiation pressure-dominated atmospheres tend to become highly inhomogeneous on scales much smaller than the radiation pressure scale height, due to the nonlinear development of “photon bubble instability” (Arons 1992; Gammie 1998). In this mechanism, gravity squeezes the radiation out of slightly overdense regions, while magnetic tension enhances the overdensity by preventing the gas from spreading sideways. Simultaneously, slightly underdense regions grow more tenuous as radiation force accelerates them; hence the density contrast increases. In a nonlinear analogue of this instability (Begelman 2001; hereafter, Paper I), gas cycles between narrow regions of very high density and much broader regions of very low density. Gravity and gas pressure dominate the force balance in the dense regions, which resemble gas pressure-supported exponential atmospheres, while the low density regions resemble radiation pressure-driven winds constrained to slide along magnetic field lines. To complete the cycle, low density matter rejoins the dense phase by passing through an isothermal, gas pressure-dominated shock, while gas moves from the high-density phase to the tenuous phase more gradually, via a slow magnetosonic critical point. Thus, the inhomogeneities represent nonlinear wave patterns rather than discrete, long-lived filaments.

Thanks to radiative diffusion, the presence of large density contrasts has little effect on the radiation pressure gradient. Under such conditions the Eddington limit can be circumvented. In the diffusion limit, the radiation flux is inversely proportional to the local density. Radiation flows readily through tenuous regions while avoiding regions of high density. If the low- and high-density regions are coupled dynamically, the upward radiation force in the tenuous zones can balance the downward gravitational force acting on the dense regions, and global dynamical equilibrium is maintained. In the steady-state, nonlinear plane-wave solutions presented in Paper I, a combination of ram pressure and magnetic tension forces maintains this coupling. If most of the mass resides in the dense phase while most of the volume is filled with tenuous gas, the total flux through the system can exceed the Eddington limit by a large factor (Shaviv 1998, 2000).

This effect should be relevant in luminous accretion disks, the inner regions of which are radiation pressure-dominated (Shakura & Sunyaev 1973). If there are no small-scale inhomogeneities, as assumed in the standard theory, the vertical scale height $`h`$ at a given radius is proportional to the local radiation flux. If the luminosity approaches the Eddington limit for the central mass, such a disk must inflate to a thickness of order its radius, probably triggering substantial mass loss. But in the presence of strong density inhomogeneities the disk becomes “leaky,” and more radiation can escape without inflating the disk. One can even imagine a scenario in which the luminosity exceeds the Eddington limit, while the disk remains geometrically thin.

In this Letter we analyze the possible effects of nonlinear photon bubble instability on thin accretion disks. In § 2 we summarize some generic features of the resulting inhomogeneities in radiation-dominated atmospheres. In § 3 we consider their effects on the structure and luminosity of accretion disks, and in § 4 we discuss the results and propose possible applications to astrophysical systems. In particular, we suggest that the observed “ultraluminous X-ray sources” are accretion disks radiating at above the Eddington limit for their central masses.

## 2. Inhomogenous Radiation-Dominated Atmospheres

Consider a plane-parallel radiation-dominated atmosphere with column density $`\mathrm{\Sigma }`$ and scale height $`h`$ in a gravitational field $`g(h)`$. The radiation flux passing through the atmosphere is $`F`$ and the Eddington flux is $`F_E(h)=cg(h)/\kappa `$, where $`\kappa `$ is the opacity. The mean mass density is $`\overline{\rho }=\mathrm{\Sigma }/2h`$. If most of the volume of the atmosphere is filled with tenuous gas with density $`\rho _{}\overline{\rho }`$, then the condition for overall dynamical equilibrium is (Shaviv 1998)

$$\frac{p_r}{h}\overline{\rho }g\frac{\rho _{}\kappa F}{c},$$
(1)
where $`p_r`$ is the radiation pressure. Since gas cycles between low and high densities via isothermal shocks, there is an additional, local pressure balance condition between the ram pressure of the tenuous phase, $`\rho _{}v_{}^2`$, and the postshock pressure, $`\rho _+c_g^2`$, in the dense phase, where $`c_g`$ is the thermal sound speed in the gas. Downstream of the shock, matter accelerates along each magnetic flux tube (which is assumed to be approximately straight on scales of order the gas pressure scale height, $`c_g^2/gh`$), and passes through a sonic point where $`v_s=c_g`$ (Paper I). This implies a constant mass flux $`\dot{m}_s=\rho _sc_g`$, where $`\rho _s`$ is the density at the sonic point. If radiative transfer obeys the diffusion law with $`F\rho ^1`$, then $`\rho _s=\overline{\rho }`$. We therefore have

$$\rho _+c_g^2\rho _{}v_{}^2=\frac{\dot{m}^2}{\rho _{}}=\frac{\overline{\rho }^2}{\rho _{}}c_g^2.$$
(2)
Thus, $`\overline{\rho }`$ is the geometric mean between $`\rho _{}`$ and $`\rho _+`$.

We solve the two pressure balance equations for $`\rho _+`$ in terms of $`F`$ and the other parameters. These results are verified by detailed calculations of steady-state, plane-parallel inhomogeneities (Paper I). This yields the “Eddington enhancement factor,” $`\mathrm{}F/F_E\rho _+/\overline{\rho }`$, giving the factor by which the radiation flux through the atmosphere exceeds that through a homogeneous atmosphere with the same column density and scale height.

The above analysis applies only if the standard radiative diffusion equation is valid locally throughout the flow, i.e., if the gas is optically thick across any density scale length along the flow, $`\rho ^2\kappa /|\rho |>1`$. Where this condition is violated the radiation field does not “see” the change in density and the flux does not adjust to local conditions. The local diffusion condition is first violated near the minimum density, $`\rho _{}`$ (Paper I), placing a lower limit on the critical mass flux, $`\dot{m}_sg/c_g\kappa `$. In this limit, the critical point occurs at a density larger than $`\overline{\rho }`$. The Eddington enhancement factor is therefore

$$\mathrm{}\frac{F}{F_E}\frac{\rho _+}{\overline{\rho }}\mathrm{min}(1,\frac{\overline{\rho }^2c_g^4\kappa ^2}{g^2}).$$
(3)
To complete the model we specify $`\rho _+`$ in terms of other atmospheric parameters. While the nonlinear theory in its current state does not predict the level of inhomogeneity, it clearly indicates an upper limit. If the maximum gas pressure $`\rho _+c_g^2`$ exceeds the magnetic pressure, the field will buckle and the crucial dynamical coupling will be lost. We conjecture that this limit also represents the natural saturation level of inhomogeneous structure. We therefore set $`\rho _+c_g^2\xi p_m`$, where $`p_m`$ is the magnetic pressure and $`\xi <1`$ is a parameter reflecting the “inefficiency” of the process by which inhomogeneities grow (e.g., due to the geometry and/or time-dependence of the magnetic field, etc.). The value of $`\xi `$ is not known, but we guess that it might lie in the range $`0.010.1`$. In anticipation of applying our model to accretion disks, we set $`p_m\alpha p_r`$ where $`\alpha `$ might also lie in the range $`0.010.1`$. Therefore, from eq. (3) we have

$$\mathrm{}\xi \alpha \frac{p_r}{\overline{p}_g}\mathrm{min}[1,\frac{\overline{p}_g^2\kappa ^2}{g^2}],$$
(4)
where $`\overline{p}_g=\overline{\rho }c_g^2p_r`$ is the mean gas pressure. The corresponding spacing between dense regions — the “wavelength” of the inhomogeneities — is given by $`\lambda \xi \alpha h\mathrm{min}(1,\overline{p}_g\kappa /g)`$, which is comfortably smaller than the radiation pressure scale height. Note that $`c_g`$ can be estimated by assuming that the gas is in LTE with the radiation field. The correct $`\alpha `$ to use here corresponds to the total magnetic pressure, rather than the $`B_rB_\varphi `$ stress associated with the transport of angular momentum. Simulations suggest that the latter might be $`10`$ times smaller than $`p_m`$, which is dominated by $`B_\varphi ^2/8\pi `$ (Brandenburg et al. 1995; Stone et al. 1996). In applying this model to accretion disks, we correct for the discrepancy by choosing a larger fiducial value of $`\xi `$, $`\xi 0.1`$.

## 3. Application to Accretion Disks

Now consider the atmosphere in the context of a thin, Keplerian accretion disk at radius $`r`$ with accretion rate $`\dot{M}`$. We will primarily consider a one-zone model in which the disk is characterized by a single mean density and scale height; later we will consider the effects of density stratification. We normalize $`r`$ to the gravitational radius of a central mass $`M`$, $`\dot{M}`$ to the Eddington accretion rate $`\dot{M}_E=4\pi GM/\kappa c`$, and $`h`$ to the radius:

$$x\frac{c^2r}{GM};\dot{m}\frac{\dot{M}}{\dot{M}_E};\delta \frac{h}{r}.$$
(5)
The local flux due to accretion is then

$$F=\frac{3}{8\pi }\frac{GM\dot{M}}{r^3}𝒟=\frac{3}{2}\frac{c^5}{GM\kappa }x^3\dot{m}𝒟,$$
(6)
where $`𝒟1(x_{\mathrm{in}}/x)^{1/2}`$ for the standard assumption of zero torque at the innermost stable orbit, $`x_{\mathrm{in}}`$ (Shakura & Sunyaev 1973). The Eddington enhancement factor is

$$\mathrm{}=\frac{3}{2}\frac{\dot{m}𝒟}{\delta x}.$$
(7)
In a standard radiation-dominated accretion disk $`\delta `$ adjusts so that $`\mathrm{}=1`$, but this will not be the case in the presence of small-scale inhomogeneities. For an $`\alpha `$model viscosity, the dissipated flux is given by

$$F\frac{3}{2}\alpha \mathrm{\Sigma }\mathrm{\Omega }^3h^2=\frac{3}{2}\frac{c^5}{GM}\alpha \mathrm{\Sigma }\delta ^2x^{5/2}.$$
(8)
Equating the accretion flux and the dissipated flux we obtain

$$\mathrm{\Sigma }\frac{\dot{m}𝒟}{\alpha \kappa }\delta ^2x^{1/2};\overline{\rho }=\frac{\mathrm{\Sigma }}{2h}\frac{c^2}{2GM\kappa }\frac{\dot{m}𝒟}{\alpha }\delta ^3x^{3/2}.$$
(9)
The radiation and gas pressure are calculated assuming vertical hydrostatic equilibrium and LTE, respectively. Setting $`\kappa =0.4`$ cm2 g-1 for electron scattering in ionized hydrogen, and letting $`mM/M_{}`$, we have

$`\overline{\rho }8\times 10^6{\displaystyle \frac{\dot{m}𝒟}{m\alpha }}\delta ^3x^{3/2}\mathrm{g}\mathrm{cm}^3;`$ (10)
$`p_r8\times 10^{15}{\displaystyle \frac{\dot{m}𝒟}{m\alpha }}\delta ^1x^{5/2}\mathrm{dyn}\mathrm{cm}^2;`$ (11)
$`\overline{p}_g5\times 10^{10}\left({\displaystyle \frac{\dot{m}𝒟}{m\alpha }}\right)^{5/4}\delta ^{13/4}x^{17/8}\mathrm{dyn}\mathrm{cm}^2.`$ (12)
Combining these formulae with equations (4) and (7), we obtain

$`\delta \mathrm{max}[0.3\xi _1^{4/13}\alpha _2^{5/13}(\dot{m}𝒟)^{5/13}m^{1/13}x^{5/26},`$
$`0.1\xi _1^{4/21}\alpha _2^{5/21}(\dot{m}𝒟)^{5/21}m^{1/21}x^{1/14}],`$ (13)
where we have normalized $`\xi `$ to 0.1 and $`\alpha `$ to 0.01. A similar equation for $`\mathrm{}`$ can be derived from eq. (3), using eq. (7). Inhomogeneities affect the disk structure only if $`\mathrm{}>1`$, implying

$`x<x_{\mathrm{inhom}}\mathrm{min}[7\xi _1^{8/21}\alpha _2^{10/21}\dot{m}^{16/21}m^{2/21},`$
$`9\xi _1^{8/45}\alpha _2^{2/9}\dot{m}^{32/45}m^{2/45}],`$ (14)
where we have set $`𝒟1`$ for $`x_{\mathrm{inhom}}x_{\mathrm{in}}`$.

Now consider the maximum possible mass flux that can be carried by the disk. Since inhomogeneities are ineffective outside $`x_{\mathrm{inhom}}`$, we have $`\dot{m}_{\mathrm{max}}(x_{\mathrm{inhom}})2x_{\mathrm{inhom}}/3`$, corresponding to a thick disk \[$`\delta (x_{\mathrm{inhom}})1`$\] with a local radiation flux equal to the Eddington limit. Any excess matter supplied at larger $`x`$ presumably will be driven away by radiation pressure before reaching $`x_{\mathrm{inhom}}`$. For $`\dot{m}`$ close to $`\dot{m}_{\mathrm{max}}`$ and $`\xi 0.1`$, local radiative diffusion should apply through most of the inhomogeneous region, and certainly in the inner regions where most energy is released. We therefore have $`\dot{m}_{\mathrm{max}}(x_{\mathrm{inhom}})580\xi _1^{8/5}\alpha _2^2m^{2/5}`$. By setting $`\delta (x)1`$ at $`x<x_{\mathrm{inhom}}`$, we find that the carrying capacity of the disk decreases with decreasing $`r`$, $`\dot{m}_{\mathrm{max}}x^{1/2}/𝒟`$, down to the radius $`x4x_{\mathrm{in}}`$ at which $`x^{1/2}/𝒟`$ has a local minimum. Thus, the mass flux reaching the central object may be much smaller than that entering the inhomogeneous region, $`\dot{m}_{\mathrm{in},\mathrm{max}}4(x_{\mathrm{in}}/x_{\mathrm{inhom}})^{1/2}\dot{m}_{\mathrm{max}}(x_{\mathrm{inhom}})`$. Nevertheless, the luminosity can exceed $`L_E`$.

According to the one-zone model, the maximum possible disk luminosity satisfies

$$\frac{L_{\mathrm{max}}}{L_E}>ϵ\dot{m}_{\mathrm{in},\mathrm{max}}30\xi _1^{4/5}\alpha _2\left(\frac{m}{10}\right)^{1/5}\left(\frac{ϵ}{0.1}\right)\left(\frac{x_{\mathrm{in}}}{6}\right)^{1/2},$$
(15)
where $`ϵ`$ is the accretion efficiency and the inequality reflects the fact that $`\dot{m}`$ increases with $`r`$. However, mass loss from the disk’s upper layers may limit the maximum attainable luminosity. Such a wind must arise, in a super-Eddington disk, from strata where the mean optical depth is so small that density inhomogeneities provide insufficient modulation of the flux. If the radiation becomes “trapped” in the accelerating outflow it will drive mass loss at a rate comparable to or exceeding the inflow rate, possibly regulating the accretion luminosity to a value $`L_E`$. But if the gas remains sufficiently inhomogeneous in the acceleration zone, radiative losses due to diffusion could reduce the energy supply available for unbinding the gas. The mass loss rate would then be lower (Shaviv 2001). Mass loss could be further inhibited if dynamo processes in the disk maintain a sufficiently dense array of closed magnetic flux loops in a magnetically dominated corona (Miller & Stone 2000). If $`\mathrm{}<p_m\kappa /g`$, then gas propelled upward into closed field lines anchored in the disk can become trapped by the magnetic tension. Once its optical depth reaches $`\mathrm{}`$, a trapped blob begins to be weighed down by gravity, and if $`\tau `$ reaches $`p_m\kappa /g`$ before the field opens up, the gas will overcome the tension and drop back onto the disk.

## 4. Discussion

We have argued that the radiation-dominated, inner regions of accretion disks may be able to approach or even exceed the Eddington limit, without becoming geometrically thick or suffering catastrophic mass loss. The disks become highly inhomogeneous on small scales, and correspondingly porous to radiation leakage, due to the nonlinear development of photon bubble instability (Gammie 1998; Paper I). The possible enhancement in radiation flux is limited by the strength of the magnetic field, which in turn is related to viscous dissipation and angular momentum transport via the magnetorotational instability (MRI, which operates concurrently with photon bubble instability \[Blaes & Socrates 2001; Turner, Stone, & Sano 2001\]). Greater porosity requires larger density constrasts, which place greater stresses on the magnetic field. Beyond a certain flux, the magnetic tension will be overwhelmed and the field will buckle. Simple one-zone models suggest that accretion disks around stellar-mass black holes could radiate as much as $`10`$ times the Eddington limit. A weak mass dependence implies that disks around supermassive black holes could exceed the Eddington limit by an even larger factor, perhaps $`100`$ or more.

Although the magnetic configurations in accretion disks are highly dynamical and do not resemble the quiescent, uniform fields assumed in our analytic models, this time-dependence need not impede the development of inhomogeneities. The linear growth time for photon bubble instability is $`(\lambda /g)^{1/2}`$, the free-fall time across the wavelength (Gammie 1998). This is shorter than the orbital time, which determines the evolution timescale of magnetic structures due to MRI (Blaes & Socrates 2001).

Catastrophic mass loss from the disk’s upper layers poses the most serious threat to the production of super-Eddington luminosities. If the luminosity is effectively trapped in the outflowing gas it could deplete the accretion flow, regulating the maximum luminosity to $`L_E`$. But if inhomogeneities persist into the acceleration zone the mass loss could be lower, leaving the accretion flow largely intact. A magnetically structured corona with a sufficient number of closed field lines could also inhibit mass loss by trapping and recycling some of the expelled gas. Both speculations will have to be checked by detailed calculations, including numerical simulations.

Fluxes $`10L_E`$ could explain the existence of “ultraluminous X-ray sources” (ULXs: Fabbiano 1989; Colbert & Mushotzky 1999; Makishima et al. 2000; Fabbiano, Zezas, & Murray 2001). The high luminosities of these sources have led to suggestions that they represent a new class of intermediate-mass black holes. However, the observed luminosities are also explainable in terms super-Eddington fluxes from disks around black holes of a few solar masses. This is consistent with anomalies reported by Makishima et al. (2000), who performed “multicolor blackbody disk” fits to ULX spectra and found that the best-fitting models were considerably hotter and more compact than expected for $`100M_{}`$ black holes. King et al. (2001) have also proposed that ULXs are sources producing super-Eddington fluxes, but they posit that the intense radiation is beamed within a small solid angle. In contrast to our proposal, which predicts a roughly isotropic flux, they require a large population of misdirected sources. However, we agree with them that these systems should not be very long-lived, since an object accreting at $`10`$ times the Eddington limit would double its mass in a few million years.

We make two other predictions based on the analysis presented in this paper. First, inhomogeneous accretion disks radiating above the Eddington limit should possess “cold” winds driven by continuum radiation pressure. In contrast to thermally driven outflows, these winds should be at or below the Compton temperature of the accretion disk radiation, and well below the virial temperature. They should have substantial optical depths and may therefore modify the underlying disk radiation. Moreover, if “hot” flares (energized by magnetic reconnection, for example) percolate through such a wind, Compton downscattering of the hard radiation could produce a prominent “soft excess”. It may be difficult to distinguish such winds from large-scale hydromagnetic flows (Blandford & Payne 1982; Königl & Kartje 1994), if the latter also possess large optical depths.

Second, there may be a population of “ultraluminous” AGNs, i.e., supermassive black holes radiating at well above their Eddington limits. Candidates could be identified by looking for discrepancies between the black hole mass as estimated from the properties of the host galaxy (Ferrarese & Merritt 2000; Gebhardt et al. 2000) or by reverberation mapping and line widths, and the lower bound on the mass according to the Eddington limit. An epoch of super-Eddington accretion at high $`z`$ could lead to the rapid growth of supermassive black holes by accretion.

Finally, we note that SS 433 may be an example of an object which is accreting at a rate that far exceeds even the upper limit attainable in a “leaky” disk (King, Taam, & Begelman 2000). Most of the matter initially in the accretion flow is therefore blown away before it reaches the vicinity of the central object. Nevertheless, our analysis suggests that the mass flux reaching the central object could exceed the Eddington limit by a substantial factor. Much of the energy released by the accreting matter would be available for driving an outflow and perhaps goes into powering the $`0.26c`$ jets, which are believed to carry a super-Eddington energy flux (Begelman et al. 1980). We presumably do not see a highly super-Eddington radiative flux because of the large optical depth of intervening matter.

## Acknowledgments

This work was supported in part by NSF grant AST-9876887. I am grateful to Neal Turner, Roger Blandford, and the referee, Nir Shaviv, for helpful comments.