text:
# Nonmonotonic critical temperature in superconductor/ferromagnet bilayers

## I Introduction

Superconductivity and ferromagnetism are two competing orders: while the former “prefers” an antiparallel spin orientation of electrons in Cooper pairs, the latter forces the spins to align in parallel. Therefore, their coexistence in one and the same material is possible only in a narrow interval of parameters; hence the interplay between superconductivity and ferromagnetism is most conveniently studied when the two interactions are spatially separated. In this case the coexistence of the two orders is due to the proximity effect. Recently, much attention has been paid to properties of hybrid proximity systems containing superconductors (S) and ferromagnets (F); new physical phenomena were observed and predicted in these systems.Ryazanov ; Kontos ; Radovic ; Tagirov\_PRL ; Buzdin\_DOS ; Nazarov\_DOS  One of the most striking effects in SF layered structures is highly nonmonotonic dependence of their critical temperature $`T_c`$ on the thickness $`d_f`$ of the ferromagnetic layers. Experiments exploring this nonmonotonic behavior were performed previously on SF multilayers such as Nb/Gd,Jiang  Nb/Fe,Muhge  V/V-Fe,Aarts  and Pb/Fe,Lazar  but the results (and, in particular, the comparison between the experiments and theories) were not conclusive.

To perform reliable experimental measurements of $`T_c(d_f)`$, it is essential to have $`d_f`$ large compared to the interatomic distance; this situation can be achieved only in the limit of weak ferromagnets. Active experimental investigations of SF bilayers and multilayers based on Cu-Ni dilute ferromagnetic alloys are carried out by several groups.Ryazanov\_new ; Aarts\_recent  In SF bilayers, they observed nonmonotonic dependence $`T_c(d_f)`$. While the reason for this effect in multilayers can be the $`0`$$`\pi `$ transition,Radovic  in a bilayer system with a single superconductor this mechanism is irrelevant, and the cause of the effect is interference of quasiparticle, specific to SF structures.

In the present paper, motivated by the experiments of Refs. Ryazanov\_new, ; Aarts\_recent,  we theoretically study the critical temperature of SF bilayers. Previous theoretical investigations of $`T_c`$ in SF structures were concentrated on systems with thin or thick layers (compared to the corresponding coherence lengths); with SF boundaries having very low or very high transparencies; the exchange energy was often assumed to be much larger than the critical temperature; in addition, the methods for solving the problem were usually approximate.Radovic ; Aarts ; Lazar ; Buzdin ; Demler ; Khusainov ; Tagirov ; Tagirov\_PRL  The parameters of the experiments of Refs. Ryazanov\_new, ; Aarts\_recent,  do not correspond to any of the above limiting cases. In the present paper we develop two approaches giving the opportunity to investigate not only the limiting cases of parameters but also the intermediate region. Using our methods, we find different types of nonmonotonic behavior of $`T_c`$ as a function of $`d_f`$, such as minimum of $`T_c`$ and even reentrant superconductivity. Comparison of our theoretical predictions with the experimental data shows good agreement.

A number of methods can be used for calculating $`T_c`$. When the critical temperature of the structure is close to the critical temperature $`T_{cs}`$ of the superconductor without the ferromagnetic layer, the Ginzburg–Landau (GL) theory applies. However, $`T_c`$ of SF bilayers may significantly deviate from $`T_{cs}`$, therefore we choose a more general theory valid at arbitrary temperature — the quasiclassical approach.Usadel ; LO ; RS  Near $`T_c`$ the quasiclassical equations become linear. In the literature the emerging problem is often treated with the help of the so-called “single-mode” approximation,Demler ; Khusainov ; Tagirov ; Tagirov\_PRL  which is argued to be qualitatively reasonable in a wide region of parameters. However, this method is justified only in a specific region of parameters which we find below. Moreover, below we show examples when this method fails even qualitatively. Thus there is need for an exact solution of the linearized quasiclassical equations. The limiting case of perfect boundaries and large exchange energy was treated by Radović et al.Radovic

Based on the progress achieved for calculation of $`T_c`$ in SN systems (where N denotes a nonmagnetic normal material),Golubov  we develop a generalization of the single-mode approximation — the multi-mode method. Although this method seems to be exact, it is subtle to justify it rigorously. Therefore we develop yet another approach (this time mathematically rigorous), which we call “the method of fundamental solution”. The models considered previouslyRadovic ; Aarts ; Lazar ; Buzdin ; Demler ; Khusainov ; Tagirov ; Tagirov\_PRL  correspond to limiting cases of our theory. A part of our results was briefly reported in Ref. JETPL, .

The paper is organized as follows. In Sec. II we formulate the Usadel equations and the corresponding boundary conditions. Section III is devoted to the exact multi-mode method for solving the general equations. An alternative exact method, the method of fundamental solution, is presented in Sec. IV. In Sec. V we describe results of our methods. In Sec. VI, a qualitative explanation of our results is presented, applicability of the results to multilayered structures is discussed, and the use of a complex diffusion constant is commented upon. Conclusions are presented in Sec. VII. Appendixes A,B contain analytical results for limiting cases. Finally, technical details of the calculations are given in Appendix C.

## II Model

We assume that the dirty-limit conditions are fulfilled, and calculate the critical temperature of the bilayer within the framework of the linearized Usadel equations for the S and F layers (the domain $`0<x<d_s`$ is occupied by the S metal , $`d_f<x<0`$ — by the F metal, see Fig. 1). Near $`T_c`$ the normal Green function is $`G=\mathrm{sgn}\omega _n`$, and the Usadel equations for the anomalous function $`F`$ take the form

$$\xi _s^2\pi T_{cs}\frac{d^2F_s}{dx^2}|\omega _n|F_s+\mathrm{\Delta }=0,0<x<d_s;$$
(1)
$$\xi _f^2\pi T_{cs}\frac{d^2F_f}{dx^2}(|\omega _n|+iE_{\mathrm{ex}}\mathrm{sgn}\omega _n)F_f=0,$$
(2)
$$d_f<x<0;$$

$$\mathrm{\Delta }\mathrm{ln}\frac{T_{cs}}{T}=\pi T\underset{\omega _n}{}\left(\frac{\mathrm{\Delta }}{|\omega _n|}F_s\right)$$
(3)
(the pairing potential $`\mathrm{\Delta }`$ is nonzero only in the S part). Here $`\xi _s=\sqrt{D_s/2\pi T_{cs}}`$, $`\xi _f=\sqrt{D_f/2\pi T_{cs}}`$ are the coherence lengths, while the diffusion constants can be expressed via the Fermi velocity and the mean free path: $`D=vl/3`$; $`\omega _n=\pi T(2n+1)`$ with $`n=0,\pm 1,\pm 2,\mathrm{}`$ are the Matsubara frequencies; $`E_{\mathrm{ex}}`$ is the exchange energy; and $`T_{cs}`$ is the critical temperature of the S material. $`F_{s(f)}`$ denotes the function $`F`$ in the S(F) region. We use the system of units in which Planck’s and Boltzmann’s constants equal unity, $`\mathrm{}=k_B=1`$.

Equations (1)–(3) must be supplemented with the boundary conditions at the outer surfaces of the bilayer:

$$\frac{dF_s(d_s)}{dx}=\frac{dF_f(d_f)}{dx}=0,$$
(4)
as well as at the SF boundary:KL

$`\xi _s{\displaystyle \frac{dF_s(0)}{dx}}=\gamma \xi _f{\displaystyle \frac{dF_f(0)}{dx}},`$ $`\gamma ={\displaystyle \frac{\rho _s\xi _s}{\rho _f\xi _f}},`$ (5)
$`\xi _f\gamma _b{\displaystyle \frac{dF_f(0)}{dx}}=F_s(0)F_f(0),`$ $`\gamma _b={\displaystyle \frac{R_b𝒜}{\rho _f\xi _f}}.`$ (6)
Here $`\rho _s`$, $`\rho _f`$ are the normal-state resistivities of the S and F metals, $`R_b`$ is the resistance of the SF boundary, and $`𝒜`$ is its area. The Usadel equation in the F layer is readily solved:

$$F_f=C(\omega _n)\mathrm{cosh}\left(k_f[x+d_f]\right),$$
(7)
$$k_f=\frac{1}{\xi _f}\sqrt{\frac{|\omega _n|+iE_{\mathrm{ex}}\mathrm{sgn}\omega _n}{\pi T_{cs}}},$$

and the boundary condition at $`x=0`$ can be written in closed form with respect to $`F_s`$:

$$\xi _s\frac{dF_s(0)}{dx}=\frac{\gamma }{\gamma _b+B_f(\omega _n)}F_s(0),$$
(8)
$$B_f=\left[k_f\xi _f\mathrm{tanh}(k_fd_f)\right]^1.$$

This boundary condition is complex. In order to rewrite it in a real form, we do the usual trick and go over to the functions

$$F^\pm =F(\omega _n)\pm F(\omega _n).$$
(9)
According to the Usadel equations (1)–(3), there is the symmetry $`F(\omega _n)=F^{}(\omega _n)`$ which implies that $`F^+`$ is real while $`F^{}`$ is a purely imaginary function.

The symmetric properties of $`F^+`$ and $`F^{}`$ with respect to $`\omega _n`$ are trivial, so we shall treat only positive $`\omega _n`$. The self-consistency equation is expressed only via the symmetric function $`F_s^+`$:

$$\mathrm{\Delta }\mathrm{ln}\frac{T_{cs}}{T}=\pi T\underset{\omega _n>0}{}\left(\frac{2\mathrm{\Delta }}{\omega _n}F_s^+\right),$$
(10)
and the problem of determining $`T_c`$ can be formulated in a closed form with respect to $`F_s^+`$ as follows. The Usadel equation for the antisymmetric function $`F_s^{}`$ does not contain $`\mathrm{\Delta }`$, hence it can be solved analytically. After that we exclude $`F_s^{}`$ from boundary condition (8) and arrive at the effective boundary conditions for $`F_s^+`$:

$$\xi _s\frac{dF_s^+(0)}{dx}=W(\omega _n)F_s^+(0),\frac{dF_s^+(d_s)}{dx}=0,$$
(11)
where

$$W(\omega _n)=\gamma \frac{A_s(\gamma _b+\mathrm{Re}B_f)+\gamma }{A_s|\gamma _b+B_f|^2+\gamma (\gamma _b+\mathrm{Re}B_f)},$$
(12)
$$A_s=k_s\xi _s\mathrm{tanh}(k_sd_s),k_s=\frac{1}{\xi _s}\sqrt{\frac{\omega _n}{\pi T_{cs}}}.$$

The self-consistency equation (10) and boundary conditions (11)–(12), together with the Usadel equation for $`F_s^+`$:

$$\xi _s^2\pi T_{cs}\frac{d^2F_s^+}{dx^2}\omega _nF_s^++2\mathrm{\Delta }=0$$
(13)
will be used below for finding the critical temperature of the bilayer.

The problem can be solved analytically only in limiting cases (see Appendix A). In the general case, one should use a numerical method, and below we propose two methods for solving the problem exactly.

## III Multi-mode method

### III.1 Starting point: the single-mode approximation and its applicability

In the single-mode approximation (SMA) one seeks the solution of the problem (10)–(13) in the form

$$F_s^+(x,\omega _n)=f(\omega _n)\mathrm{cos}\left(\mathrm{\Omega }\frac{xd_s}{\xi _s}\right),$$
(14)
$$\mathrm{\Delta }(x)=\delta \mathrm{cos}\left(\mathrm{\Omega }\frac{xd_s}{\xi _s}\right).$$
(15)
This anzatz automatically satisfies boundary condition (11) at $`x=d_s`$.

The Usadel equation (13) yields

$$f(\omega _n)=\frac{2\delta }{\omega _n+\mathrm{\Omega }^2\pi T_{cs}},$$
(16)
then the self-consistency Eq. (10) takes the form ($`\delta `$ and $`\mathrm{\Omega }`$ do not depend on $`\omega _n`$)

$$\mathrm{ln}\frac{T_{cs}}{T_c}=\psi \left(\frac{1}{2}+\frac{\mathrm{\Omega }^2}{2}\frac{T_{cs}}{T_c}\right)\psi \left(\frac{1}{2}\right),$$
(17)
where $`\psi `$ is the digamma function.

Boundary condition (11) at $`x=0`$ yields

$$\mathrm{\Omega }\mathrm{tan}\left(\mathrm{\Omega }\frac{d_s}{\xi _s}\right)=W(\omega _n).$$
(18)
The critical temperature $`T_c`$ is determined by Eqs. (17),(18).

Although this method is popular, it is often used without pointing out the limits of its applicability. We present the explicit formulation of the corresponding condition: the single-mode method is correct only if the parameters are such that $`W`$ can be considered $`\omega _n`$-independent \[because the left-hand side of Eq. (18) must be $`\omega _n`$-independent\].Buzdin

Appendix B demonstrates examples of the SMA validity and corresponding analytical results.

In one of experimentally relevant cases, $`E_{\mathrm{ex}}/\pi T_{cs}>1`$, $`d_f\xi _f`$, the SMA is applicable if $`\sqrt{E_{\mathrm{ex}}/\pi T_{cs}}1/\gamma _b`$ (see Appendix B for details).

### III.2 Inclusion of other modes

The single-mode approximation implies that one takes the (only) real root $`\mathrm{\Omega }`$ of Eq. (17). An exact (multi-mode) method for solving problem (10)–(13) is obtained if we also take imaginary roots into account — there is infinite number of these.Golubov

Thus we seek the solution in the form

$`F_s^+(x,\omega _n)`$ $`=`$ $`f_0(\omega _n)\mathrm{cos}\left(\mathrm{\Omega }_0{\displaystyle \frac{xd_s}{\xi _s}}\right)`$ (19)
$`+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}f_m(\omega _n){\displaystyle \frac{\mathrm{cosh}\left(\mathrm{\Omega }_m\frac{xd_s}{\xi _s}\right)}{\mathrm{cosh}\left(\mathrm{\Omega }_m\frac{d_s}{\xi _s}\right)}},`$
$`\mathrm{\Delta }(x)`$ $`=`$ $`\delta _0\mathrm{cos}\left(\mathrm{\Omega }_0{\displaystyle \frac{xd_s}{\xi _s}}\right)`$ (20)
$`+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\delta _m{\displaystyle \frac{\mathrm{cosh}\left(\mathrm{\Omega }_m\frac{xd_s}{\xi _s}\right)}{\mathrm{cosh}\left(\mathrm{\Omega }_m\frac{d_s}{\xi _s}\right)}}.`$
(The normalizing denominators in the $`\mathrm{cosh}`$-terms have been introduced in order to increase accuracy of numerical calculations.) This anzatz automatically satisfies boundary condition (11) at $`x=d_s`$.

Substituting the anzatz \[Eqs. (19)–(20)\] into the Usadel equation (13), we obtain

$`f_0(\omega _n)`$ $`={\displaystyle \frac{2\delta _0}{\omega _n+\mathrm{\Omega }_0^2\pi T_{cs}}},`$ (21)
$`f_m(\omega _n)`$ $`={\displaystyle \frac{2\delta _m}{\omega _n\mathrm{\Omega }_m^2\pi T_{cs}}},m=1,2,\mathrm{},`$
then the parameters $`\mathrm{\Omega }`$ are determined by the self-consistency equation (10) ($`\delta `$ and $`\mathrm{\Omega }`$ do not depend on $`\omega _n`$):

$`\mathrm{ln}{\displaystyle \frac{T_{cs}}{T_c}}`$ $`=\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\mathrm{\Omega }_0^2}{2}}{\displaystyle \frac{T_{cs}}{T_c}}\right)\psi \left({\displaystyle \frac{1}{2}}\right),`$ (22)
$`\mathrm{ln}{\displaystyle \frac{T_{cs}}{T_c}}`$ $`=\psi \left({\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{\Omega }_m^2}{2}}{\displaystyle \frac{T_{cs}}{T_c}}\right)\psi \left({\displaystyle \frac{1}{2}}\right),m=1,2,\mathrm{}`$
From Eqs. (22) and properties of the digamma functiondigamma  it follows that the parameters $`\mathrm{\Omega }`$ belong to the following intervals:

$$0<\mathrm{\Omega }_0^2<\frac{1}{2\gamma __E},$$
(23)
$$\frac{T_c}{T_{cs}}(2m1)<\mathrm{\Omega }_m^2<\frac{T_c}{T_{cs}}(2m+1),m=1,2,\mathrm{},$$

where $`\gamma __E1.78`$ is Euler’s constant.

Boundary condition (11) at $`x=0`$ yields the following equation for the amplitudes $`\delta `$:

$`\delta _0{\displaystyle \frac{W(\omega _n)\mathrm{cos}\left(\mathrm{\Omega }_0d_s/\xi _s\right)\mathrm{\Omega }_0\mathrm{sin}\left(\mathrm{\Omega }_0d_s/\xi _s\right)}{\omega _n+\mathrm{\Omega }_0^2\pi T_{cs}}}`$
$`+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\delta _m{\displaystyle \frac{W(\omega _n)+\mathrm{\Omega }_m\mathrm{tanh}\left(\mathrm{\Omega }_md_s/\xi _s\right)}{\omega _n\mathrm{\Omega }_m^2\pi T_{cs}}}=0.`$ (24)
The critical temperature $`T_c`$ is determined by Eqs. (22) and the condition that Eq. (24) has a nontrivial ($`\omega _n`$-independent) solution with respect to $`\delta `$.

Numerically, we take a finite number of modes: $`m=0,1,\mathrm{},M`$. To take account of $`\omega _n`$-independence of the solution, we write down Eq. (24) at the Matsubara frequencies up to the $`N`$th frequency: $`n=0,1,\mathrm{},N`$. Thus we arrive at the matrix equation $`K_{nm}\delta _m=0`$ with the following matrix $`\widehat{K}`$:

$$K_{n0}=\frac{W(\omega _n)\mathrm{cos}\left(\mathrm{\Omega }_0d_s/\xi _s\right)\mathrm{\Omega }_0\mathrm{sin}\left(\mathrm{\Omega }_0d_s/\xi _s\right)}{\omega _n/\pi T_{cs}+\mathrm{\Omega }_0^2},$$

$$K_{nm}=\frac{W(\omega _n)+\mathrm{\Omega }_m\mathrm{tanh}\left(\mathrm{\Omega }_md_s/\xi _s\right)}{\omega _n/\pi T_{cs}\mathrm{\Omega }_m^2},$$
(25)
$$n=0,1,\mathrm{},N,m=1,2,\mathrm{},M.$$

We take $`M=N`$, then the condition that Eq. (24) has a nontrivial solution takes the form

$$det\widehat{K}=0.$$
(26)
Thus the critical temperature $`T_c`$ is determined as the largest solution of Eqs. (22),(26).

## IV Method of fundamental solution

By definition, the fundamental solution $`G(x,y;\omega _n)`$ (which is also called the Green function) of problem (11)–(13) satisfies the same equations, but with the delta-functional “source”:Morse

$$\xi _s^2\pi T_{cs}\frac{d^2G(x,y)}{dx^2}\omega _nG(x,y)=\delta (xy),$$
(27)
$$\xi _s\frac{dG(0,y)}{dx}=W(\omega _n)G(0,y),\frac{dG(d_s,y)}{dx}=0.$$
(28)
The fundamental solution can be expressed via solutions $`v_1`$, $`v_2`$ of Eq. (27) without the delta-function, satisfying the boundary conditions at $`x=0`$ and $`x=d_s`$, respectively:

$`G(x,y;\omega _n)`$ $`=`$ $`{\displaystyle \frac{k_s/\omega _n}{\mathrm{sinh}(k_sd_s)+(W/k_s\xi _s)\mathrm{cosh}\left(k_sd_s\right)}}`$ (29)
$`\times \{\begin{array}{c}\hfill v_1(x)v_2(y),xy\\ \hfill v_2(x)v_1(y),yx\end{array},`$
where

$`v_1(x)`$ $`=\mathrm{cosh}(k_sx)+(W/k_s\xi _s)\mathrm{sinh}(k_sx),`$ (30a)
$`v_2(x)`$ $`=\mathrm{cosh}\left(k_s[xd_s]\right).`$ (30b)
Having found $`G(x,y;\omega _n)`$, we can write the solution of Eqs. (11)–(13) as

$$F_s^+(x;\omega _n)=2_0^{d_s}G(x,y;\omega _n)\mathrm{\Delta }(y)𝑑y.$$
(31)
Substituting this into the self-consistency equation (10), we obtain

$$\begin{array}{c}\mathrm{\Delta }(x)\mathrm{ln}\frac{T_{cs}}{T_c}\hfill \\ \hfill =2\pi T_c\underset{\omega _n>0}{}\left[\frac{\mathrm{\Delta }(x)}{\omega _n}_0^{d_s}G(x,y;\omega _n)\mathrm{\Delta }(y)𝑑y\right].\end{array}$$
(32)
This equation can be expressed in an operator form: $`\mathrm{\Delta }\mathrm{ln}(T_{cs}/T_c)=\widehat{L}\mathrm{\Delta }`$. Then the condition that Eq. (32) has a nontrivial solution with respect to $`\mathrm{\Delta }`$ is expressed by the equation

$$det\left(\widehat{L}\widehat{1}\mathrm{ln}\frac{T_{cs}}{T_c}\right)=0.$$
(33)
The critical temperature $`T_c`$ is determined as the largest solution of this equation.

Numerically, we put problem (32),(33) on a spatial grid, so that the linear operator $`\widehat{L}`$ becomes a finite matrix.

## V Numerical results

In Secs. III, IV we developed two methods for calculating the critical temperature of a SF bilayer. Specifying parameters of the bilayer we can find the critical temperature numerically. It can be checked that the multi-mode method and the method of fundamental solution yield equivalent results. However, at small temperatures $`T_cT_{cs}`$, the calculation time for the multi-mode method increases. Indeed, the size of the matrix $`\widehat{K}`$ \[Eq. (25)\] is determined by the number $`N`$ of the maximum Matsubara frequency $`\omega _N`$, which must be much larger than the characteristic energy $`\pi T_{cs}`$; hence $`NT_{cs}/T_c`$. Therefore, at low temperatures we use the method of fundamental solution.

### V.1 Comparison with experiment

Using our methods we fit the experimental data of Ref. Ryazanov\_new, ; the result is presented in Fig. 2. Estimating the parameters $`d_s=11`$ nm, $`T_{cs}=7`$ K, $`\rho _s=7.5`$$`\mu \mathrm{\Omega }`$ cm, $`\xi _s=8.9`$ nm, $`\rho _f=60`$$`\mu \mathrm{\Omega }`$ cm, $`\xi _f=7.6`$ nm, $`\gamma =0.15`$ from the experiment,private\_com  and fitting only $`E_{\mathrm{ex}}`$ and $`\gamma _b`$, we find good agreement between our theoretical predictions and the experimental data.

The fitting procedure was the following: first, we determine $`E_{\mathrm{ex}}130`$ K from the position of the minimum of $`T_c(d_f)`$; second, we find $`\gamma _b0.3`$ from fitting the vertical position of the curve.

The deviation of our curve from the experimental points is small; it is most pronounced in the region of small $`d_f`$ corresponding to the initial decrease of $`T_c`$. This is not unexpected because, when $`d_f`$ is of the order of a few nanometers, the thickness of the F film may vary significantly along the film (which is not taken into account in our theory), and the thinnest films can even be formed by an array of islands rather than by continuous material. At the same time, we emphasize that the minimum of $`T_c`$ takes place at $`d_f5`$ nm, when with good accuracy the F layer has uniform thickness.

### V.2 Various types of $`T_c(d_f)`$ behavior

The experimental results discussed above represent only one possible type of $`T_c(d_f)`$ behavior. Now we address the general case; we obtain different kinds of $`T_c(d_f)`$ curves depending on parameters of the bilayer.

To illustrate, in Fig. 3 we plot several curves for various values of $`\gamma _b`$ \[we recall that $`\gamma _bR_b`$, where $`R_b`$ is the resistance of the SF interface in the normal state — see Eq. (6)\]. The exchange energy is $`E_{\mathrm{ex}}=150`$ K; the other parameters are the same as in Fig. 2.

We observe three characteristic types of $`T_c(d_f)`$ behavior: 1) at large enough interface resistance, $`T_c`$ decays nonmonotonically to a finite value exhibiting a minimum at a particular $`d_f`$, 2) at moderate interface resistance, $`T_c`$ demonstrates the reentrant behavior: it vanishes in a certain interval of $`d_f`$, and is finite otherwise, 3) at low enough interface resistance, $`T_c`$ decays monotonically vanishing at finite $`d_f`$. A similar succession of $`T_c(d_f)`$ curves as in Fig. 3 can be obtained by tuning other parameters, e.g., the exchange energy $`E_{\mathrm{ex}}`$ or the normal resistances of the layers (the parameter $`\gamma `$).

A common feature seen from Fig. 3 is saturation of $`T_c`$ at large $`d_f\lambda _{\mathrm{ex}}`$. This fact has a simple physical explanation: the suppression of superconductivity by a dirty ferromagnet is only due to the effective F layer with thickness on the order of $`\lambda _{\mathrm{ex}}`$, adjacent to the interface (this is the layer explored and “felt” by quasiparticles entering from the S side due to the proximity effect).

It was shown by Radović et al.Radovic  that the order of the phase transition may change in short-periodic SF superlattices, becoming the first order. We also observe this feature in the curves of types 2) and 3) mentioned above. This phenomenon manifests itself as discontinuity of $`T_c(d_f)`$: the critical temperature jumps to zero abruptly without taking intermediate values (see Figs. 3,4). Formally, $`T_c`$ becomes a double-valued function, but the smaller solution is physically unstable (dotted curve in Fig. 4).

An interesting problem is determination of the tricritical point where the order of the phase transition changes. The corresponding result for homogeneous bulk superconductors with internal exchange field was obtained a long time ago in the framework of the Ginzburg–Landau theory.Sarma  However, the generalization to the case when the GL theory is not valid is a subtle problem which has not yet been solved. We note that the equations used in Refs. Radovic, ; Khusainov,  were applied beyond their applicability range because they are GL results valid only when $`T_c`$ is close to $`T_{cs}`$.

### V.3 Comparison between single- and multi-mode methods

A popular method widely used in the literature for calculating the critical temperature of SF bi- and multi-layers is the single-mode approximation. The condition of its validity was formulated in Sec. III.1. However, this approximation is often used for arbitrary system’s parameters. Using the methods developed in Secs. III,IV, we can check the actual accuracy of the single-mode approximation. The results are presented in Fig. 5.

We conclude that although at some parameters the results of the single-mode and multi-mode (exact) methods are close (Figs. 5 a,f), in the general case they are quantitatively and even qualitatively different \[Figs. 5 b,c,d,e — these cases correspond to the most nontrivial $`T_c(d_f)`$ behavior\]. Thus to obtain reliable results one should use one of the exact (multi-mode or fundamental-solution) techniques.

### V.4 Spatial dependence of the order parameter

The proximity effect in the SF bilayer is characterized by the spatial behavior or the order parameter, which can be chosen as

$$F(x,\tau =0)=T\underset{\omega _n}{}F(x,\omega _n),$$
(34)
where $`\tau `$ denotes the imaginary time \[in the S metal $`F(x,\tau =0)\mathrm{\Delta }(x)`$\]. This function is real due to the symmetry relation $`F(\omega _n)=F^{}(\omega _n)`$.

We illustrate this dependence in Fig. 6, which shows two cases differing by the thickness of the F layer $`d_f`$ (and by the corresponding $`T_c`$). Although the critical temperatures differ by more than the order of magnitude, the normalized order parameters are very close to each other, which means that the value of $`T_c`$ has almost no effect on the shape of $`F(x,\tau =0)`$. Details of the calculation are presented in Appendix C.

Another feature seen from Fig. 6 is that the order parameter in the F layer changes its sign when the thickness of the F layer increases (this feature can be seen for the dotted curve, although negative values of the order parameter have very small amplitudes). We discuss this oscillating behavior in the next section.

## VI Discussion

### VI.1 Qualitative explanation of the nonmonotonic $`T_c(d_f)`$ behavior

The thickness of the F layer at which the minimum of $`T_c(d_f)`$ occurs, can be estimated from qualitative arguments based on the interference of quasiparticles in the ferromagnet.

Let us consider a point $`x`$ inside the F layer. According to Feynman’s interpretation of quantum mechanics,Feynman  the quasiparticle wave function may be represented as a sum of wave amplitudes over all classical trajectories; the wave amplitude for a given trajectory is equal to $`\mathrm{exp}(iS)`$, where $`S`$ is the classical action along this trajectory. We are interested in an anomalous wave function of correlated quasiparticles, which characterizes superconductivity; this function is equivalent to the anomalous Green function $`F(x)`$. To obtain this wave function we must sum over trajectories that (i) start and end at the point $`x`$, (ii) change the type of the quasiparticle (i.e., convert an electron into a hole, or vice versa). There are four kinds of trajectories that should be taken into account (see Fig. 7).

Two of them (denoted 1 and 2) start in the direction toward the SF interface (as an electron and as a hole), experience the Andreev reflection, and return to the point $`x`$. The other two trajectories (denoted 3 and 4) start in the direction away from the interface, experience normal reflection at the outer surface of the F layer, move toward the SF interface, experience the Andreev reflection there, and finally return to the point $`x`$. The main contribution is given by the trajectories normal to the interface. The corresponding actions are

$$S_1=Qx\alpha ,$$
(35)
$$S_2=Qx\alpha ,$$
(36)
$$S_3=Q(2d_f+x)\alpha ,$$
(37)
$$S_4=Q(2d_f+x)\alpha $$
(38)
(note that $`x<0`$), where $`Q`$ is the difference between the wave numbers of the electron and the hole, and $`\alpha =\mathrm{arccos}(E/\mathrm{\Delta })`$ is the phase of the Andreev reflection. To make our arguments more clear, we assume that the ferromagnet is strong, the SF interface is ideal, and consider the clean limit first: then $`Q=k_ek_h=\sqrt{2m(E+E_{\mathrm{ex}}+\mu )}\sqrt{2m(EE_{\mathrm{ex}}+\mu )}2E_{\mathrm{ex}}/v`$, where $`E`$ is the quasiparticle energy, $`\mu `$ is the Fermi energy, and $`v`$ is the Fermi velocity. Thus the anomalous wave function of the quasiparticles is

$$F(x)\underset{n=1}{\overset{4}{}}\mathrm{exp}(iS_n)\mathrm{cos}(Qd_f)\mathrm{cos}\left(Q[d_f+x]\right).$$
(39)
The suppression of $`T_c`$ by the ferromagnet is determined by the value of the wave function at the SF interface: $`F(0)\mathrm{cos}^2(Qd_f)`$. The minimum of $`T_c`$ corresponds to the minimum value of $`F(0)`$ which is achieved at $`d_f=\pi /2Q`$. In the dirty limit the above expression for $`Q`$ is replaced by

$$Q=\sqrt{\frac{E_{\mathrm{ex}}}{D_f}}\frac{2\pi }{\lambda _{\mathrm{ex}}}$$
(40)
(here we have defined the wavelength of the oscillations $`\lambda _{\mathrm{ex}}`$); hence the minimum of $`T_c(d_f)`$ takes place at

$$d_f^{(\mathrm{min})}=\frac{\pi }{2}\sqrt{\frac{D_f}{E_{\mathrm{ex}}}}=\frac{\lambda _{\mathrm{ex}}}{4}.$$
(41)
For the bilayer of Ref. Ryazanov\_new,  we obtain $`d_f^{(\mathrm{min})}7`$ nm, whereas the experimental value is $`5`$ nm (Fig. 2); thus our qualitative estimate is reasonable.

The arguments given above seem to yield not only the minimum but rather a succession of minima and maxima. However, numerically we obtain either a single minimum or a minimum followed by a weak maximum (Fig. 3). The reason for this is that actually the anomalous wave function not only oscillates in the ferromagnetic layer but also decays exponentially, which makes the amplitude of the subsequent oscillations almost invisible.

Finally, we note that our arguments concerning oscillations of $`F(x)`$ also apply to a half-infinite ferromagnet, where we should take into account only the trajectories 1 and 2 (see Fig. 7). This yields $`F(x)\mathrm{cos}(Qx)`$ (another qualitative explanation of this result can be found, for example, in Ref. Demler, ).

### VI.2 Multilayered structures

The methods developed and the results obtained in this paper apply directly to more complicated symmetric multilayered structures in the $`0`$-state such as SFS and FSF trilayers, SFIFS and FSISF systems (I denotes an arbitrary potential barrier), and SF superlattices. In such systems an SF bilayer can be considered as a unit cell, and joining together the solutions of the Usadel equations in each bilayer we obtain the solution for the whole system (for more details see Sec. VIII of Ref. FF, ).

Our methods can be generalized to take account of possible superconductive and/or magnetic $`\pi `$-states (when $`\mathrm{\Delta }`$ and/or $`E_{\mathrm{ex}}`$ may change their signs from layer to layer). In this case the system cannot be equivalently separated into a set of bilayers. Mathematically, this means that the solutions of the Usadel equations lose their purely cosine form \[see Eqs. (7), (14), (15), (19), (20), (30b)\] acquiring a sine part as well.

### VI.3 Complex diffusion constant?

Finally, we comment on Refs. Tagirov, ; Tagirov\_PRL, ; Khusainov, ; Buzdin\_new, , where the authors considered (in the vicinity of $`T_c`$) diffusion equations with a complex diffusion constant $`D_f`$ for the F part of the structure. This implies small complex corrections to $`D_f`$ over $`E_{\mathrm{ex}}\tau 1`$ in the Usadel equations ($`\tau `$ is the time of the mean free path). However, we disagree with this method for the following reason: although the complex $`D_f`$ can indeed be formally obtained in the course of the standard derivation of the Usadel equationsUsadel  from the Eilenberger onesEilenberger  by expanding over the spherical harmonics, one can check that the higher harmonics neglected in the derivation have the same order of magnitude as the retained complex correction to $`D_f`$. Hence the complexity of $`D_f`$ in the context of the Usadel equations is the excess of accuracy. Below we present our arguments.

We give a brief derivation of the Usadel equations showing how the complex diffusion constant can be obtained and why this result cannot be trusted. In the “quasi-one-dimensional” geometry (which means that the parameters vary only as a function of $`x`$) the linearized Eilenberger equation in the presence of disorder and the exchange field has the form

$$\frac{v\mathrm{cos}\theta }{2}\frac{d}{dx}F+\left(\omega _n+\frac{1}{2\tau }+iE_{\mathrm{ex}}\right)F=\mathrm{\Delta }+\frac{F}{2\tau },$$
(42)
where, for simplicity, we assume a positive Matsubara frequency $`\omega _n>0`$, and $`\theta `$ is the angle between the $`x`$ axis and the direction of the Fermi velocity $`𝐯`$, while $`\mathrm{}`$ denotes angular averaging over the spherical angles. The disorder is characterized by the time of the mean free path $`\tau `$ and the mean free path $`l`$ (to be used below). In the dirty limit the anomalous Green function $`F`$ is nearly isotropic. However, to obtain the Usadel equation for the isotropic part of $`F`$, we must also take into account the next term from the full Legendre polynomial expansion:

$`F(x,\omega _n,\theta )`$ $`={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}F_k(x,\omega _n)P_k(\mathrm{cos}\theta )`$
$`F_0(x,\omega _n)+F_1(x,\omega _n)\mathrm{cos}\theta .`$ (43)
Here we have neglected the harmonics with $`k2`$ assuming them small; we shall check this assumption later.

Averaging Eq. (42) over the spherical angles first directly and second — after being multiplied by $`\mathrm{cos}\theta `$, we arrive at

$$\frac{v}{6}\frac{d}{dx}F_1+\left(\omega _n+iE_{\mathrm{ex}}\right)F_0=\mathrm{\Delta },$$
(44)
$$\frac{v}{2}\frac{d}{dx}F_0+\left(\omega _n+\frac{1}{2\tau }+iE_{\mathrm{ex}}\right)F_1=0.$$
(45)
Equation (45) yields

$$F_1=\left(\frac{l}{1+2\omega _n\tau +2iE_{\mathrm{ex}}\tau }\right)\frac{d}{dx}F_0,$$
(46)
then Eq. (44) leads to

$$\frac{D}{2}\frac{d^2}{dx^2}F_0\left(\omega _n+iE_{\mathrm{ex}}\right)F_0+\mathrm{\Delta }=0,$$
(47)
$$D=\frac{vl/3}{1+2\omega _n\tau +2iE_{\mathrm{ex}}\tau }.$$

Now we must check that the assumption $`|F_1/F_0|1`$, $`|F_2/F_1|1`$, etc. that we used is indeed satisfied. From Eq. (46) we obtain

$$\left|\frac{F_1}{F_0}\right|\frac{l/L}{\mathrm{max}(1,2\omega _n\tau ,2E_{\mathrm{ex}}\tau )},$$
(48)
where $`L`$ is the characteristic space scale on which $`F_0`$ varies. According to the Usadel equation (47), it is given by

$$L\frac{l}{\sqrt{\mathrm{max}(1,2\omega _n\tau ,2E_{\mathrm{ex}}\tau )\mathrm{max}(2\omega _n\tau ,2E_{\mathrm{ex}}\tau )}},$$
(49)
and the condition of the Usadel equation’s validity is written as

$$\left|\frac{F_1}{F_0}\right|\sqrt{\frac{\mathrm{max}(2\omega _n\tau ,2E_{\mathrm{ex}}\tau )}{\mathrm{max}(1,2\omega _n\tau ,2E_{\mathrm{ex}}\tau )}}1$$
(50)
\[similarly, we can also keep the term with $`k=2`$ in series (43), which yields $`|F_2/F_1||F_1/F_0|`$, etc.\].

Finally, condition (50) takes the form

$$2\pi T_{cs}\tau 1,2E_{\mathrm{ex}}\tau 1$$
(51)
(we have taken into account that the characteristic energy is $`\omega _n\pi T_{cs}`$).

Now we can analyze our results. If condition (51) is satisfied and the Usadel equation is valid, the neglected angular harmonics have the relative order of magnitude $`|F_2/F_0|\mathrm{max}(2\pi T_{cs}\tau ,2E_{\mathrm{ex}}\tau )`$; hence we cannot retain the terms of the same order in the diffusion constant \[see Eq. (47)\], and we should use the standard expression $`D=vl/3`$.

## VII Conclusions

In the present paper we have developed two methods for calculating the critical temperature of a SF bilayer as a function of its parameters (the thicknesses and material parameters of the layers, the quality of the interface). The multi-mode method is a generalization of the corresponding approach developed in Ref. Golubov,  for SN systems. However, the rigorous justification of this method is not clear. Therefore, we propose yet another approach — the method of fundamental solution, which is mathematically rigorous. The results demonstrate that the two methods are equivalent; however, at low temperatures (compared to $`T_{cs}`$) the accuracy requirements are stricter for the multi-mode method, and the method of fundamental solution is preferable. Comparing our method with experiment we obtain good agreement.

In the general case, we observe three characteristic types of $`T_c(d_f)`$ behavior: 1) nonmonotonic decay of $`T_c`$ to a finite value exhibiting a minimum at particular $`d_f`$, 2) reentrant behavior, characterized by vanishing of $`T_c`$ in a certain interval of $`d_f`$ and finite values otherwise, 3) monotonic decay of $`T_c`$ and vanishing at finite $`d_f`$. Qualitatively, the nonmonotonic behavior of $`T_c(d_f)`$ is explained by interference of quasiparticles in the F layer, which can be either constructive or destructive depending on the value of $`d_f`$.

Using the developed methods we have checked the accuracy of the widely used single-mode approximation. We conclude that although at some parameters the results of the single-mode and exact methods are close, in the general case they are quantitatively and even qualitatively different. Thus, to obtain reliable results one should use one of the exact (multi-mode or fundamental-solution) techniques.

The spatial dependence of the order parameter (at the transition point) is shown to be almost insensitive to the value of $`T_c`$.

The methods developed and the results obtained in this paper apply directly to more complicated symmetric multilayered structures in the $`0`$-state such as SFS and FSF trilayers, SFIFS and FSISF systems, and SF superlattices. Our methods can be generalized to take account of possible superconductive and/or magnetic $`\pi `$-states (when $`\mathrm{\Delta }`$ and/or $`E_{\mathrm{ex}}`$ may change their signs from layer to layer).

We argue that the use of the complex diffusion constant in the Usadel equation is the excess of accuracy.

In several limiting cases, $`T_c`$ is considered analytically.

###### Acknowledgements.

We thank V. V. Ryazanov and M. V. Feigel’man for stimulating discussions. We are especially indebted to V. V. Ryazanov for communicating the experimental result of his group to us prior to the detailed publication. We are also grateful to J. Aarts, A. I. Buzdin, M. Yu. Kupriyanov, Yu. Oreg, and L. R. Tagirov for useful comments. Ya.V.F. acknowledges financial support from the Russian Foundation for Basic Research (RFBR) under project No. 01-02-17759, and from Forschungszentrum Jülich (Landau Scholarship). The research of N.M.C. was supported by the RFBR (projects Nos. 01-02-06230 and 00-02-16617), by Forschungszentrum Jülich (Landau Scholarship), by the Netherlands Organization for Scientific Research (NWO), by the Einstein Center, and by the Swiss National Foundation.

## Appendix A Analytical results for a thin S layer

(i) When $`d_s\xi _s`$ and $`E_{\mathrm{ex}}\pi T_{cs}`$, problem (10)–(13) can be solved analytically. The first of the above conditions implies that $`\mathrm{\Delta }`$ can be considered constant, and $`F^+`$ weakly depends on the spatial coordinate; so $`F^+(x,\omega _n)=2\mathrm{\Delta }/\omega _n+A(\omega _n)\mathrm{cosh}(k_s[xd_s])`$. The boundary conditions determine the coefficient $`A`$; as a result

$$F^+(\omega _n)F^+(x=0,\omega _n)=\frac{2\mathrm{\Delta }}{\omega _n}\left[\frac{A_s(\omega _n)}{A_s(\omega _n)+W(\omega _n)}\right],$$
(52)
where $`k_s`$, $`A_s`$, and $`W`$ are defined in Eq. (12). Finally, the self-consistency equation for $`T_c`$ takes the form

$$\mathrm{ln}\frac{T_{cs}}{T_c}=\mathrm{Re}\psi \left(\frac{1}{2}+\frac{\gamma }{2}\frac{\xi _s}{d_s}\frac{1}{\gamma _b+B_f}\frac{T_{cs}}{T_c}\right)\psi \left(\frac{1}{2}\right),$$
(53)
where $`B_f`$ does not depend on $`\omega _n`$ due to the condition $`E_{\mathrm{ex}}\pi T_{cs}`$:

$$B_f=\left[k_f\xi _f\mathrm{tanh}(k_fd_f)\right]^1,k_f=\frac{1}{\xi _f}\sqrt{\frac{iE_{\mathrm{ex}}}{\pi T_{cs}}}.$$
(54)
(ii) If the F layer is also thin, $`d_f\sqrt{D_f/2E_{\mathrm{ex}}}`$, Eq. (53) is further simplified:

$$\mathrm{ln}\frac{T_{cs}}{T_c}=\mathrm{Re}\psi \left(\frac{1}{2}+\frac{\tau _f}{\tau _s}\left[\frac{1}{i+\tau _fE_{\mathrm{ex}}}\right]\frac{E_{\mathrm{ex}}}{2\pi T_c}\right)\psi \left(\frac{1}{2}\right),$$
(55)
where $`\tau _s`$, $`\tau _f`$ are defined similarly to Ref. FF, :

$$\tau _s=\frac{2d_sR_b𝒜}{\rho _sD_s},\tau _f=\frac{2d_fR_b𝒜}{\rho _fD_f},$$
(56)
and have the physical meaning of the escape time from the corresponding layer. They are related to the quantities $`\gamma `$, $`\gamma _b`$ used in the body of the paper as

$$\tau _s=\frac{\gamma _b}{\gamma }\frac{1}{\pi T_{cs}}\frac{d_s}{\xi _s},\tau _f=\gamma _b\frac{1}{\pi T_{cs}}\frac{d_f}{\xi _f}.$$
(57)
(iii) If the S layer is thin, $`d_s\xi _s`$, and the SF interface is opaque, $`\gamma _b\mathrm{}`$, the critical temperature of the bilayer only slightly deviates from $`T_{cs}`$. In this limit Eq. (52) applies with $`W=\gamma /\gamma _b1`$, and we finally obtain:

$$T_c=T_{cs}\frac{\pi }{4\tau _s}.$$
(58)
Interestingly, characteristics of the F layer ($`d_f`$, $`E_{\mathrm{ex}}`$, etc.) do not enter the formula. In particular, this formula is valid for an SN bilayerMcMillan ; T\_c\_McMillan  (where N is a nonmagnetic normal material, $`E_{\mathrm{ex}}=0`$) because Eq. (58) was obtained without any assumptions about the value of the exchange energy.

### A.1 Transparent interface

When both layers are very thin \[$`d_s\sqrt{D_s/2\omega __D}`$, $`d_f\mathrm{min}(\sqrt{D_f/2\omega __D},\sqrt{D_f/2E_{\mathrm{ex}}})`$, with $`\omega __D`$ the Debye energy of the S metal\] and the interface is transparent, the bilayer is equivalent to a homogeneous superconducting layer with internal exchange field. This layer is described by effective parameters: the pairing potential $`\mathrm{\Delta }^{(\mathrm{eff})}`$, the exchange field $`E_{\mathrm{ex}}^{(\mathrm{eff})}`$, and the pairing constant $`\lambda ^{(\mathrm{eff})}`$. In this subsection we develop the ideas of Ref. Bergeret-Efetov, , demonstrate a simple derivation of this description, and find the limits of its applicability.

The Usadel equations (1),(2) for the two layers can be written as a single equation:

$$\begin{array}{c}\frac{D_f\theta (x)+D_s\theta (x)}{2}\frac{d^2F}{dx^2}\hfill \\ \hfill |\omega _n|FiE_{\mathrm{ex}}\mathrm{sgn}(\omega _n)\theta (x)F+\mathrm{\Delta }\theta (x)=0,\end{array}$$
(59)
where $`\theta `$ is the Heaviside function \[$`\theta (x>0)=1`$, $`\theta (x<0)=0`$\]. The self-consistency equation (3) can be rewritten as

$$\mathrm{\Delta }(x)=\lambda \theta (x)\pi T\underset{\omega _n}{}F(x,\omega _n),$$
(60)
where $`\lambda `$ is the pairing constant.

First, we consider the ideal SF interface: $`\gamma _b=0`$ \[see Eq. (6)\], then $`F(x)`$ is continuous at the interface and nearly constant across the whole bilayer, i.e., $`F_s(x)F_f(x)=F`$. Applying the integral operator to Eq. (59):

$$\frac{\nu _f}{\nu _sd_s+\nu _fd_f}_{d_f}^0𝑑x+\frac{\nu _s}{\nu _sd_s+\nu _fd_f}_0^{d_s}𝑑x$$
(61)
(here $`\nu `$ is the normal-metal density of states), and cancelling gradient terms due to the boundary condition (5), we obtain the equations describing a homogeneous layer:

$$|\omega _n|F(\omega _n)iE_{\mathrm{ex}}^{(\mathrm{eff})}\mathrm{sgn}(\omega _n)F(\omega _n)+\mathrm{\Delta }^{(\mathrm{eff})}=0,$$
(62)
$$\mathrm{\Delta }^{(\mathrm{eff})}=\lambda ^{(\mathrm{eff})}\pi T\underset{\omega _n}{}F(\omega _n),$$
(63)
with the effective parameters (see also Ref. Bergeret-Efetov, ):

$$E_{\mathrm{ex}}^{(\mathrm{eff})}=\frac{\tau _f}{\tau _s+\tau _f}E_{\mathrm{ex}},\mathrm{\Delta }^{(\mathrm{eff})}=\frac{\tau _s}{\tau _s+\tau _f}\mathrm{\Delta },$$
(64)
$$\lambda ^{(\mathrm{eff})}=\frac{\tau _s}{\tau _s+\tau _f}\lambda ,T_{cs}^{(\mathrm{eff})}=\frac{\gamma __E}{\pi }2\omega __D\mathrm{exp}\left(\frac{1}{\lambda ^{(\mathrm{eff})}}\right),$$

where $`\gamma __E`$ is Euler’s constant and $`T_{cs}^{(\mathrm{eff})}`$ is the critical temperature of the layer in the absence of ferromagnetism (i.e., at $`E_{\mathrm{ex}}^{(\mathrm{eff})}=0`$). The critical temperature is determined by the equation

$$\mathrm{ln}\frac{T_{cs}^{(\mathrm{eff})}}{T_c}=\mathrm{Re}\psi \left(\frac{1}{2}+i\frac{E_{\mathrm{ex}}^{(\mathrm{eff})}}{2\pi T_c}\right)\psi \left(\frac{1}{2}\right).$$
(65)
Actually, the description in terms of effective parameters (64) is applicable at an arbitrary temperature (i.e., when the Usadel equations are nonlinear) and has a clear physical interpretation: the superconducting ($`\mathrm{\Delta }`$, $`\lambda `$) and ferromagnetic ($`E_{\mathrm{ex}}`$) parameters are renormalized according to the part of time spent by quasiparticles in the corresponding layer. This physical picture is based on interpretation of $`\tau `$ as escape times, which we present in the next subsection.

Now we discuss the applicability of the above description for a nonideal interface ($`\gamma _b0`$). In this case $`F`$ is nearly constant in each layer, but these constants are different: $`F_s(x)F_s+C_s(xd_s)^2`$, $`F_f(x)F_f+C_f(x+d_f)^2`$, where $`|F_s||C_s|d_s^2`$ and $`|F_f||C_f|d_f^2`$. Using the Usadel equation (59) and the boundary conditions (5),(6), we find the difference $`\delta FF_sF_f`$:

$$\delta F=\frac{\mathrm{\Delta }}{{\displaystyle \frac{1}{\tau _s}}+|\omega _n|\left[1+{\displaystyle \frac{1}{\tau _f\left(\left|\omega _n\right|+iE_{\mathrm{ex}}\mathrm{sgn}\omega _n\right)}}\right]}.$$
(66)
Finally, the homogeneous description is valid when $`|\delta F/F|1`$ \[with $`F`$ determined by Eq. (62)\], which yields:

$$\mathrm{max}(E_{\mathrm{ex}},\omega __D)\mathrm{max}(\tau _s,\tau _f)1$$
(67)
(here $`\omega _n\omega __D`$ has been taken as the largest characteristic energy scale in the quasi-homogeneous bilayer).

### A.2 Interpretation of $`\tau `$ as escape times

The quantities $`\tau _s`$, $`\tau _f`$ introduced in Eq. (56) may be interpreted as escape times from the corresponding layers. The arguments go as follows.

If the layers are thin, then the diffusion inside the layers is “fast” and the escape time from a layer is determined by the interface resistance. The time of penetration through a layer or the interface is determined by the corresponding resistance: $`R_{s(f)}`$ or $`R_b`$, hence the diffusion is “fast” if $`R_{s(f)}R_b`$.

Let us use the detailed balance approach, and consider an interval of energy $`dE`$. In the S layer, the charge in this interval is $`Q_s=e\nu _sdE𝒜d_s`$. Let us define the escape time from the S layer $`t_s`$, so that the current from S to F is equal to $`Q_s/t_s`$. On the other hand, this current can be written as $`dE/eR_b`$, hence

$$\frac{Q_s}{t_s}=\frac{dE}{eR_b},$$
(68)
and we immediately obtain

$$t_s=\frac{d_sR_b𝒜}{\rho _sD_s}.$$
(69)
Similarly, we obtain the expression for the escape time from the F layer $`t_f`$. As a result, the relations between the quantities $`\tau `$ defined in Eq. (56) and the escape times $`t`$ are simply

$$\tau _s=2t_s,\tau _f=2t_f.$$
(70)
Microscopic expressions for the escape times may be obtained using the Sharvin formula for the interface resistance. Assuming, for definiteness, that the Fermi velocity is smaller in the S metal, $`v_s<v_f`$, we obtain

$$R_b=\frac{\pi r_b}{e^2\nu _sv_s𝒜},$$
(71)
and consequently

$$t_s=\pi \frac{d_s}{v_s}r_b,t_f=\pi \frac{v_fd_f}{v_s^2}r_b,$$
(72)
where $`r_b`$ is the inverse transparency of one channel. The asymmetry in these expressions stems from our assumption $`v_s<v_f`$. In the opposite case the indices $`s`$ and $`f`$ in Eqs. (71),(72) should be interchanged.

## Appendix B Applicability of the single-mode approximation

As pointed out in Sec. III.1, the single-mode approximation (SMA) is applicable only if the parameters are such that $`W`$ \[see Eq. (12)\] can be considered $`\omega _n`$-independent. An example is the case when $`\gamma _b|B_f|`$, hence $`W=\gamma /\gamma _b`$.

The condition $`\gamma _b|B_f|`$ can be written in a simpler form; to this end we should estimate $`|B_f|`$. We introduce the real and imaginary parts of $`k_f`$: $`k_f=k_f^{}+ik_f^{\prime \prime }`$, and note that $`k_f^{}>k_f^{\prime \prime }`$. Then using the properties of the trigonometric functions and the estimate $`\mathrm{tanh}x\mathrm{min}(1,x)`$ we obtain

$$|B_f|\left[k_f^{}\xi _f\mathrm{tanh}(k_f^{}d_f)\right]^1,$$
(73)
and finally cast the condition $`\gamma _b|B_f|`$ into the form

$$\frac{1}{\gamma _b}\mathrm{min}\{\sqrt{\mathrm{max}(\frac{T_c}{T_{cs}},\frac{E_{\mathrm{ex}}}{\pi T_{cs}})};\frac{d_f}{\xi _f}\mathrm{max}(\frac{T_c}{T_{cs}},\frac{E_{\mathrm{ex}}}{\pi T_{cs}})\},$$
(74)
where the ratio $`T_c/T_{cs}`$ originates from $`\omega _n/\pi T_{cs}`$ with $`\omega _n\pi T_c`$ as the characteristic energy scale in the bilayer.

If condition (74) is satisfied, then the SMA is valid and $`T_c`$ is determined by the equations

$$\mathrm{ln}\frac{T_{cs}}{T_c}=\psi \left(\frac{1}{2}+\frac{\mathrm{\Omega }^2}{2}\frac{T_{cs}}{T_c}\right)\psi \left(\frac{1}{2}\right),$$
(75)
$$\mathrm{\Omega }\mathrm{tan}\left(\mathrm{\Omega }\frac{d_s}{\xi _s}\right)=\frac{\gamma }{\gamma _b}.$$
(76)
These equations can be further simplified in two limiting cases which we consider below.

(1) $`{\displaystyle \frac{\gamma }{\gamma _b}}{\displaystyle \frac{d_s}{\xi _s}}1`$:

in this case Eq. (76) yields $`\mathrm{\Omega }^2=\frac{\gamma }{\gamma _b}\frac{\xi _s}{d_s}`$, and Eq. (75) takes the form

$$\mathrm{ln}\frac{T_{cs}}{T_c}=\psi \left(\frac{1}{2}+\frac{1}{2}\frac{\gamma }{\gamma _b}\frac{\xi _s}{d_s}\frac{T_{cs}}{T_c}\right)\psi \left(\frac{1}{2}\right),$$
(77)
which reproduces the $`\gamma _b|B_f|`$ limit of Eq. (53).

(2) $`{\displaystyle \frac{\gamma }{\gamma _b}}{\displaystyle \frac{d_s}{\xi _s}}1`$:

in this case Eq. (76) yields $`\mathrm{\Omega }\frac{d_s}{\xi _s}=\frac{\pi }{2}`$, and Eq. (75) takes the form

$$\mathrm{ln}\frac{T_{cs}}{T_c}=\psi \left(\frac{1}{2}+\frac{\pi ^2}{8}\left[\frac{\xi _s}{d_s}\right]^2\frac{T_{cs}}{T_c}\right)\psi \left(\frac{1}{2}\right).$$
(78)
Equations (75)–(78) can be used for calculating the critical temperature $`T_c`$ and the critical thickness of the S layer $`d_s^{(\mathrm{cr})}`$ below which the superconductivity in the SF bilayer vanishes (i.e., $`T_c=0`$).

### B.1 Results for the critical temperature

In the limit when $`T_c`$ is close to $`T_{cs}`$, Eqs. (77),(78) yield

$$T_c=T_{cs}\left(1\frac{\pi ^2}{4}\frac{\gamma }{\gamma _b}\frac{\xi _s}{d_s}\right)\text{if }\frac{\gamma }{\gamma _b}\mathrm{min}(\frac{d_s}{\xi _s},\frac{\xi _s}{d_s}),$$
(79)
and

$$T_c=T_{cs}\left[1\left(\frac{\pi ^2}{4}\frac{\xi _s}{d_s}\right)^2\right]\text{if }\frac{d_s}{\xi _s}\mathrm{max}(1,\frac{\gamma _b}{\gamma }).$$
(80)
Using relations (57) one can check that result (79) is equivalent to Eq. (58).

### B.2 Results for the critical thickness

The critical thickness of the S layer $`d_s^{(\mathrm{cr})}`$ is defined as the thickness below which there is no superconductivity in the SF bilayer: $`T_c(d_s^{(\mathrm{cr})})=0`$. When $`T_c0`$, Eq. (75) yields $`\mathrm{\Omega }=1/\sqrt{2\gamma __E}`$ (where $`\gamma __E1.78`$ is Euler’s constant), and Eq. (76) takes the form

$$\frac{1}{\sqrt{2\gamma __E}}\mathrm{tan}\left(\frac{1}{\sqrt{2\gamma __E}}\frac{d_s^{(\mathrm{cr})}}{\xi _s}\right)=\frac{\gamma }{\gamma _b}.$$
(81)
Explicit results for $`d_s^{(\mathrm{cr})}`$ can be obtained in limiting cases:

$$\frac{d_s^{(\mathrm{cr})}}{\xi _s}=2\gamma __E\frac{\gamma }{\gamma _b}\text{if }\frac{\gamma }{\gamma _b}\frac{d_s}{\xi _s}1,$$
(82)
and

$$\frac{d_s^{(\mathrm{cr})}}{\xi _s}=\pi \sqrt{\frac{\gamma __E}{2}}\text{if }\frac{\gamma }{\gamma _b}\frac{d_s}{\xi _s}1.$$
(83)
## Appendix C Spatial dependence of the order parameter

According to the self-consistency equation, in the S layer the order parameter $`F(x,\tau =0)`$ is proportional to $`\mathrm{\Delta }(x)`$:

$$F_s(x,\tau =0)=\frac{\mathrm{\Delta }(x)}{\pi \lambda },$$
(84)
where $`\lambda `$ is the pairing constant which can be expressed via the Debye energy:

$$\lambda ^1=\mathrm{ln}\left(\frac{2\gamma __E\omega __D}{\pi T_{cs}}\right).$$
(85)
The pairing potential $`\mathrm{\Delta }(x)`$ can be found as the eigenvector of the matrix $`\widehat{L}\widehat{1}\mathrm{ln}(T_{cs}/T_c)`$ \[see Eq. (33)\], corresponding to the zero eigenvalue.

After that we can express $`F(x,\tau =0)`$ in the F layer via $`\mathrm{\Delta }(x)`$ in the superconductor. The Green function $`F_f(x,\omega _n)`$ in the F layer is given by Eq. (7), with $`C(\omega _n)`$ found from the boundary conditions:

$$C(\omega _n)=\left(\frac{B_f}{\gamma _b+B_f}\right)\frac{F_s(0,\omega _n)}{\mathrm{cosh}(k_fd_f)}.$$
(86)
The Green function at the S side of the SF interface is

$$F_s(0,\omega _n)=\frac{F_s^+(0,\omega _n)+F_s^{}(0,\omega _n)}{2}.$$
(87)
The symmetric part $`F_s^+`$ is given by Eq. (31). The antisymmetric part is

$$F_s^{}=C^{}(\omega _n)\mathrm{cosh}\left(k_s[xd_s]\right),$$
(88)
with $`C^{}(\omega _n)`$ found from the boundary conditions:

$$C^{}(\omega _n)=\left[\frac{i\gamma \mathrm{Im}B_f}{A_s|\gamma _b+B_f|^2+\gamma (\gamma _b+\mathrm{Re}B_f)}\right]\frac{F_s^+(0,\omega _n)}{\mathrm{cosh}(k_sd_s)}.$$
(89)
Finally, the order parameter in the F layer is the Fourier transform \[see Eq. (34)\] of

$$F_f(x,\omega _n)=\left[1+\frac{i\gamma \mathrm{Im}B_f}{A_s|\gamma _b+B_f|^2+\gamma (\gamma _b+\mathrm{Re}B_f)}\right]$$

$$\times \left(\frac{B_f}{\gamma _b+B_f}\right)\frac{\mathrm{cosh}\left(k_f[x+d_f]\right)}{\mathrm{cosh}(k_fd_f)}_0^{d_s}G(0,y;\omega _n)\mathrm{\Delta }(y)𝑑y.$$
(90)


format: text

text:
# Evidence for Circumburst Extinction of Gamma-Ray Bursts with Dark Optical Afterglows and Evidence for a Molecular Cloud Origin of Gamma-Ray Bursts

## 1 Evidence for Circumburst Extinction of Bursts with Dark Optical Afterglows

Optical afterglows have been detected for about 1/3 of the rapidly-, well-localized gamma-ray bursts (e.g., ). This data-rich subsample of the rapidly-, well-localized bursts has naturally been the focus of the vast majority of the field’s attention and resources over the past five years, but in the end, it is a biased sample. The nature of the dark optical afterglows (DOAs) has only recently become a subject of greater interest, and as might be expected given that these are by definition data-poor events, contradictory initial findings:  find that $`>3/4`$ of the limiting magnitudes for the DOAs are brighter than the detected, but faint, optical afterglow of GRB 000630, and conclude that the DOAs are consistent with a failure to image deeply enough quickly enough. However, , using a smaller sample of DOAs, find that the distribution of limiting magnitudes for the DOAs, even if treated as detections, is significantly fainter than the distribution of magnitudes of the detected optical afterglows, and conclude that the DOAs cannot be explained by a failure to image deeply enough quickly enough.

In this section, we first resolve this issue of whether the DOAs can be explained by a failure to image deeply enough quickly enough by applying Bayesian inference, a statistical formalism in which limits can be treated as limits, instead of detections (e.g., ). Next, we argue that extinction by circumburst111By circumburst, we mean within the circumburst cloud, which is probably parsecs to tens of parsecs across (; see also ). dust is the most likely explanation for the vast majority of the DOAs. Finally, we consider the effects of dust sublimation and fragmentation , and place constraints on the sizes and masses of the circumburst clouds of DOAs.

### 1.1 Evidence that DOAs are Fainter than Detected Optical Afterglows

The samples of  and  differ in the following ways: The sample of  is nearly complete – it contains limiting magnitudes for 95% of the DOAs preceding GRB 000630 – whereas the sample of  includes limiting magnitudes for BeppoSAX bursts only. Also, the sample of  relies less on GCN Reports, and more on published results. As BeppoSAX bursts appear to have been deeply imaged more often than bursts detected by other satellites, the different findings of these two papers can be attributed to sample differences, and the practice of comparing limits to detections, which can turn such sample differences into sample biases, as we demonstrate below.

But first, we combine the samples of  and , keeping only the bursts through GRB 000630 to maintain the completeness of the  sample. As is done in these papers, we scale the data to a common time: We use a temporal index of –1.4, the median temporal index in Table 1 of , and we scale the data to 18 hours after the burst, the median observation time in our combined sample. In the event of multiple limiting magnitudes for a single burst, we adopt the most constraining. We plot the combined sample in Figure 1: The hashed histogram is a binning of the R-band magnitude distribution of the detected optical afterglows. The unhashed histogram, which is added to the hashed histogram, is a binning of the limiting R-band magnitude distribution of the DOAs.

Although these distributions appear to be similar, since one is a distribution of limits, they can be consistent only if the vast majority of the DOAs could have been detected had they only been imaged a little more deeply. Being this unlucky, consistently, is of course very improbable. We quantify this improbability with Bayesian inference: Let $`n(\mathrm{R})`$ be the normalized R-band magnitude distribution of all of the bursts. The likelihood function is then given by

$$=\underset{i=1}{\overset{N}{}}_i,$$
(1)
where $`N`$ is the number of bursts in our combined sample, and

$$_i=\{\begin{array}{cc}n(\mathrm{R}_i)\hfill & \text{(}\mathrm{R}_i=\text{ detection)}\hfill \\ _{\mathrm{R}_i}^{\mathrm{}}n(\mathrm{R})𝑑\mathrm{R}\hfill & \text{(}\mathrm{R}_i=\text{ limit)}\hfill \end{array}$$
(2)
(e.g., ). We now consider a wide variety of two-parameter models for $`n(\mathrm{R})`$: a Gaussian, a boxcar, an increasing power law, and a decreasing power law. We have fitted these models to the data in Figure 1, we plot the best-fit models in Figure 1 (dotted curves), and we list the best-fit parameter values and 68% credible intervals in Table 1 of . Also in Table 1 of , we list the fractions $`f_{ALL}`$ of all bursts, and $`f_{DOA}`$ of DOAs, that have afterglows fainter than R $`=24`$ mag 18 hours after the burst, the magnitude scaled to this time of the faintest detected optical afterglow in our sample. The results are fairly independent of the assumed shape of the brightness distribution: On average, we find that $`57_{11}^{+13}`$% of all bursts, or $`82_{17}^{+22}`$% of DOAs, have afterglows that are fainter than R $`=24`$ mag 18 hours after the burst. Consequently, the DOAs cannot be explained by a failure to image deeply enough quickly enough: As a whole, they are fainter than the detected optical afterglows.

### 1.2 Evidence for Circumburst Extinction and Against Alternative Explanations

 find that the far-infrared luminosity-to-mass ratio for isolated and weakly interacting spiral galaxies is consistent with the average for Galactic molecular clouds, which suggests that most, if not all, of the star formation in such galaxies occurs in molecular clouds. Since bursts likely occur where massive stars form (see, e.g.,  for a review of some of the evidence), and since the optical depths through such clouds are typically many magnitudes (e.g., ), circumburst extinction is a likely explanation for the vast majority of the DOAs.

 put circumburst extinction to a simple test: Considering only (1) strongly collimated bursts that burn completely through the optical depth of their circumburst clouds (e.g., ), and (2) weakly collimated bursts that do not,  find that circumburst extinction can explain snapshot optical vs. X-ray and optical vs. radio brightness distributions of all of the afterglows, DOAs and dark radio afterglows included.

We now address a number of alternative explanations:

Galactic extinction:  show that the distribution of Galactic column densities for the DOAs is consistent with the distribution of Galactic column densities for the bursts with detected optical afterglows, ruling this out as an explanation for the vast majority of the DOAs.

Host galaxy extinction unrelated to the circumburst medium: If burst host galaxies are like our galaxy, this cannot explain the vast majority of the DOAs for the very same reason: Most lines of sight through our galaxy do not result in a great deal of extinction.

However, the burst host galaxies might not be like our galaxy:  argue that the vast majority of the DOAs might be occurring in the nuclear regions of ultraluminous infrared galaxies (ULIRGs), which have column densities of $`N_H10^{23}10^{24}`$ cm-2. We now test this idea. In Figure 2, we plot the source-frame column density probability distributions of the four DOAs in Table 1 of  (see below). Redshifts are not known for three of these bursts, so we adopt $`z=1`$ (top panel), 2 (middle panel), and 3 (bottom panel) for these three bursts. Given the redshifts that have been measured for bursts to date, the probabilities that all three of these bursts have $`z>1`$, 2, and 3 are 0.27, 0.013, and 0.0016, respectively. Even in the very unlikely event that all three of these bursts have $`z3`$, these four probability distributions are not consistent with the expected column density distribution for bursts that occur in the nuclear regions of ULIRGs, but are consistent with the expected column density distribution for bursts that occur in molecular clouds .

Another possibility is that burst host galaxies are very dusty. However, if this were the case, one would expect the X-ray bright afterglows to be DOAs as often as the X-ray dim afterglows, which is not the case (Figure 3). However, Figure 3 is naturally explained by circumburst extinction: The X-ray bright afterglows are probably strongly collimated, and consequently the burst and optical flash are more likely to burn through the optical depth of the circumburst cloud; the X-ray dim afterglows are probably weakly collimated, and consequently the burst and optical flash are less likely to do so .

High redshift effects: Lyman limit absorption in the source frame, absorption by the Ly$`\alpha `$ forest, absorption by excited molecular hydrogen in the circumburst medium , and source-frame extinction by the FUV component of the extinction curve (e.g., ) could all result in DOAs if at sufficiently high redshifts. However,  show that unless the burst history of the universe differs dramatically from the star-formation history of the universe, the redshift distribution of the bursts should be fairly narrowly peaked around $`z2`$, primarily because there is very little volume in the universe at low and high redshifts.  model the detection efficiency functions of BeppoSAX and IPN, and show that these satellites should detect even fewer bursts at high redshifts, pushing the expected typical redshift for bursts detected by these satellites down to the observed value of about one. However, based on their variability redshift estimates of 220 BATSE bursts,  (see also ) find that the burst history of the universe might differ dramatically from the star-formation history of the universe, with very many more bursts at high redshifts. However, using similar variability redshift estimates  for 907 BATSE bursts,  find that only $`15`$% of bursts above BATSE’s detection threshold have $`z>5`$ (and if the luminosity function is evolving, far fewer bursts below BATSE’s detection threshold have $`z>5`$). Consequently, the above high redshift effects probably affect $`<10`$% of the bursts in our BeppoSAX- and IPN-dominated sample. Furthermore,  find that the variability redshift estimates for all of the DOAs for which high resolution BATSE light curves are available have $`z<5`$.

A number of ad hoc hypotheses are also addressed in .

### 1.3 Evidence for a Large, Massive Cloud Origin of DOAs

Neglecting for a moment extinction exterior to the circumburst cloud, we now consider the optical depth to a burst that is embedded a distance $`r`$ within its circumburst cloud, along the line of sight. Prior to sublimation and fragmentation , the column to the burst is optically thin, which we take to mean $`\tau <0.3`$, if the hydrogen column density $`N_H<5\times 10^{20}`$ cm-2, and optically thick, which we take to mean $`\tau >3`$, if $`N_H>5\times 10^{21}`$ cm-2. We mark column densities with dotted lines in Figure 4. Since sublimation and fragmentation burn through $`10L_{49}^{1/2}`$ pc of optical depth , where $`L=10^{49}L_{49}`$ erg s-1 is the 1 – 7.5 eV isotropic-equivalent peak luminosity of the optical flash, the post-sublimation/fragmentation column to the burst is optically thin if $`r`$ is less than this distance or $`N_H<5\times 10^{20}`$ cm-2, and optically thick if $`r`$ is greater than this distance and $`N_H>5\times 10^{21}`$ cm-2. We show this in Figure 4 by plotting $`\tau =0.3`$ and 3 for $`L_{49}=0.1`$ (thin curves, left), 1 (thick curves), and 10 (thin curves, right; see  for details).

Consequently, modulo the value of $`L`$, bursts that occur to the lower left of the $`\tau =0.3`$ curve have either relatively unextinguished optical afterglows, or optical afterglows that are extinguished by dust elsewhere in the host galaxy, or in our galaxy. Bursts that occur to the upper right of the $`\tau =3`$ curve have highly extinguished afterglows. Since circumburst extinction appears to be responsible for most of the DOAs, and Galactic extinction and host galaxy extinction unrelated to the circumburst medium account for no more than perhaps a few of the DOAs detected to date, we find that most DOAs have $`r>10L_{49}^{1/2}`$ pc and $`N_H>5\times 10^{21}`$ cm-2.

We now estimate the sizes and masses of the circumburst clouds of the DOAs. Taking the clouds to be spherical, and taking the bursts to be located at the cloud centers, both of which are reasonable approximations on average, we find that most DOAs occur in clouds of radius $`R>10L_{49}^{1/2}`$ pc and mass $`M>3\times 10^5L_{49}`$ M (Figure 4 of ). Clouds of this size and mass are typical of giant molecular clouds (e.g., ), and are active regions of star formation.

## 2 Evidence for a Molecular Cloud Origin of Bursts

In this section, we first consider the equilibrium properties of clouds, and show that clouds of these sizes and masses, and hence the circumburst clouds of DOAs, cannot be diffuse, but rather are probably molecular, in which case they are active regions of star formation. Next, we show that Galactic-like molecular clouds have sufficiently high column densities to make bursts dark at optical wavelengths, and that sublimation of dust by the optical flash and fragmentation of dust by the burst and afterglow can leave a comparable number of bursts – bursts that are sufficiently near the earth-facing side of the cloud and/or in sufficiently small clouds – relatively unextinguished at optical wavelengths, which is in agreement with what is observed. Finally, we model and constrain the column density distribution of the bursts with detected optical afterglows, and show that these bursts are also consistent with a molecular cloud origin, and are not consistent with a diffuse cloud origin. Consequently, we find that all but perhaps a few bursts, dark or otherwise, probably occur in molecular clouds.

### 2.1 Evidence for a Molecular Cloud Origin of DOAs: Cloud Equilibrium, Gravitational Collapse, and Fragmentation-Driven Bursts of Star Formation

We first consider the equilibrium properties of diffuse clouds: i.e., gas that is bound by pressure equilibrium with a warm or hot phase of the interstellar medium. Specifically,  considers the equilibrium properties of a uniformly magnetized, isothermal, non-rotating diffuse cloud, and introduces factors $`c_1=0.53`$ and $`c_2=0.60`$ that correct for the fact that in equilibrium the cloud will not be of uniform hydrogen density $`n_H`$, but will be centrally condensed, and will not be spherical with radius $`R`$, but will be moderately flattened along the direction of the magnetic field.  shows that clouds with $`R>R_m`$ are unstable to gravitational collapse, where

$$\frac{GM^2}{R_m^4}=\frac{25p_m}{1(M_c/M)^{2/3}},$$
(3)
$`M=4\pi \rho R_m^3/3`$ is the mass of the cloud, $`\rho =\mu _Hm_Hn_H`$, $`\mu _H=1.87`$ is the mean molecular weight per hydrogen atom,

$$p_m=\frac{3.15c_2(kT/\mu )^4}{G^3M^2[1(M_c/M)^{2/3}]^3},$$
(4)
$`T`$ is the temperature of the cloud, $`\mu =1.44`$ is the mean molecular weight,

$$M_c=\frac{0.0236(c_1B)^3}{G^{3/2}\rho ^2},$$
(5)
and $`B`$ is the strength of the magnetic field within the cloud. We plot $`R_m(n_H)`$ for $`T=10`$ and 100 K and $`B=3`$ $`\mu `$G in Figure 5, and for $`T=30`$ K and $`B=1`$ and 10 $`\mu `$G in Figure 2 of : Most Galactic clouds have $`10<T<100`$ K, and $`B=3`$ $`\mu `$G is typical of the Galactic interstellar medium.

Also in Figures 5 and 2 of , we plot the curves of constant post-sublimation/fragmentation optical depth $`\tau `$ from Figure 4. Above, we argue that most of the DOAs occur to the upper right of the $`\tau =3`$ curve, which implies relatively large sizes and high densities for their circumburst clouds. We show in Figures 5 and 2 of  that diffuse clouds of these sizes and densities are unlikely to be in equilibrium, or even near equilibrium, for typical temperatures and magnetic field strengths. Such clouds would be unstable to gravitational collapse, resulting in the collapse and fragmentation of the cloud until star formation re-establishes pressure equilibrium within the fragments, and the fragments are bound by self-gravity: This is precisely the structure of molecular clouds (e.g., ). Furthermore, molecular clouds regularly have sizes $`R>10`$ pc and masses $`M>3\times 10^5`$ M, and are active regions of star formation. Consequently, DOAs probably occur in collapsing diffuse clouds/forming molecular clouds, and/or molecular clouds that re-established equilibrium some time ago. We confirm that molecular clouds have sufficiently high column densities to make a large fraction of bursts dark at optical wavelengths below.

 find the mass scale of the molecular cloud fragments to be on the order of a few solar masses, which supports the idea that they are supported against further collapse and fragmentation by star formation. Since the lifetimes of massive stars are comparable to the free fall time of the progenitor cloud for these cloud masses, the first generation of massive star formation and supernovae should occur on this short timescale: $`10^710^8`$ yr. Consequently, DOAs are probably a byproduct of this burst of star formation if the molecular cloud formed recently, and/or the result of lingering or latter generation star formation if the molecular cloud formed some time ago.

### 2.2 Consistency Check: The Column Density and Optical Depth Distributions of Bursts in Galactic-Like Molecular Clouds

We now consider the column density and optical depth properties of bursts in Galactic-like molecular clouds. First, we compute the mean column densities of the 273 molecular clouds in the sample of , which was constructed using data from the Massachusetts-Stony Brook CO Galactic Plane Survey. The mean column density of a cloud is given by $`N_H=M/(\mu _Hm_HA)`$, where $`M`$ is the virial mass ( show that the clouds are in or near virial equilibrium), $`A=11.56(D\mathrm{tan}\sqrt{\sigma _l\sigma _b})^2`$ is the effective cross section of the cloud , $`D`$ is the distance to the cloud, and $`\sigma _l`$ and $`\sigma _b`$ are the angular sizes of the cloud. We plot the distribution (dotted histogram) of mean column densities in Figure 6, and confirm the finding of  that $`\mu _HN_H`$ is narrowly peaked around $`170`$ M pc-2, independent of cloud mass and size.

Next, we correct this distribution for a number of effects. First, since bursts more likely trace cloud mass than cloud number, we weight these mean column densities by cloud mass, and renormalize the distribution (dashed histogram). Since molecular cloud column densities appear to be fairly independent of cloud mass, the distribution is not significantly changed. Also, we weight the $`M<7\times 10^4`$ M clouds by an additional factor of $`(M/7\times 10^4`$ M$`{}_{}{}^{})^{3/2}`$ to correct for an undercounting of low mass clouds on the far side of the Galaxy . However, this affects the mass-weighted mean column density distribution negligibly.

Finally, we correct for geometrical effects: (1) bursts occur in their circumburst clouds, not behind them; and (2) molecular clouds are centrally condensed. First, we adopt the same density distribution for the clouds that  adopted to compute the virial masses: $`\rho r^1`$ within the effective radius of the cloud, and $`\rho =0`$ beyond this radius. This density distribution results in a surface brightness profile that is similar to what is observed . Next, we trace this density distribution with $`10^4`$ randomly-placed bursts for each cloud, and compute for each burst the distance and column density to the earth-facing side of the cloud. The effect of placing bursts in the clouds, as opposed to behind the clouds, is to lower the mean column densities on average by a factor of two, and sometimes (when bursts occur near the earth-facing side of the cloud) considerably more. However, this is partially offset by the central condensation of the clouds, which places more bursts in above average density (and column density) environments. We plot the final distribution – the expected column density distribution for bursts in Galactic-like molecular clouds – also in Figure 6 (solid histogram).

Lastly, we confirm that Galactic-like molecular clouds, the column densities of which appear to be fairly universal, at least within the Galaxy, have sufficiently high column densities to make a large fraction of bursts dark at optical wavelengths. To this end, we plot in Figure 7 the expected distribution (1 and 2 $`\sigma `$ dashed contours) of column densities and distances to the earth-facing side of the cloud for bursts in Galactic-like molecular clouds. Also in Figure 7, we plot the curves of constant post-sublimation/fragmentation optical depth from Figures 4 and 5, and lines of constant column density. It appears that for low values of the isotropic-equivalent peak luminosity of the optical flash ($`L_{49}<0.1`$), Galactic-like molecular clouds can make a large fraction of bursts dark at optical wavelengths. However, even for very low values of $`L`$, a reasonable fraction of bursts would be left relatively unextinguished at optical wavelengths. Since comparable fractions of DOAs and detected optical afterglows are observed (about 2/3 and 1/3, respectively; e.g., ), this suggests that the bursts with detected optical afterglows might also occur in molecular clouds. We confirm this by modeling and constraining the distribution of column densities, measured from absorption of the X-ray afterglow, of the bursts with detected optical afterglows below.

### 2.3 Evidence for a Molecular Cloud Origin for Bursts with Detected Optical Afterglows: The Column Density Distribution

In this section, we model and constrain the column density distribution of the bursts with detected optical afterglows, and compare this distribution to the expected distribution for bursts in Galactic-like molecular clouds (solid histogram of Figure 6). To this end, we have reviewed the literature, and list in Table 1 of  observer-frame column densities and uncertainties that have been measured from the X-ray afterglows of 15 bursts, 11 of which have detected optical afterglows (this information is not yet available for the vast majority of the DOAs). Before modeling these data, we convert each measurement $`N_{H,i}`$ and upper and lower 1-$`\sigma `$ uncertainty $`\sigma _{u,i}`$ and $`\sigma _{l,i}`$ into probability distributions, first in the observer frame, and then in the source frame by correcting the observer-frame probability distribution for the Galactic column density along the line of sight, and the redshift of the burst if known. If $`\sigma _{u,i}=\sigma _{l,i}`$ (or $`\sigma _{l,i}=N_{H,i}`$), we take the observer-frame probability distribution to be a Gaussian of mean $`N_{H,i}`$ and standard deviation $`\sigma _{u,i}`$: $`p_i(N_{H,obs})=G(N_{H,obs},N_{H,i},\sigma _{u,i})`$. Otherwise, we take the observer-frame probability distribution to be a skewed Gaussian: $`p_i(N_{H,obs})=G(N_{H,obs}^s,N_{H,i}^s,\sigma _i^s)dN_{H,obs}^s/dN_{H,obs}`$, where $`s`$ is the skewness parameter, and $`s`$ and $`\sigma _i`$ are given by solving

$`\sigma _i^s`$ $`=`$ $`(N_{H,i}+\sigma _{u,i})^sN_{H,i}^s`$ (6)
$`=`$ $`N_{H,i}^s(N_{H,i}\sigma _{l,i})^s,`$ (7)
i.e., we take the observer-frame probability distribution to be a Gaussian in $`N_{H,obs}^s`$, which reduces to a Gaussian in $`N_{H,obs}`$ when $`\sigma _{u,i}=\sigma _{l,i}`$ ($`s=1`$). The source-frame probability distribution $`p_i(N_H)`$ is then given by substituting $`N_{H,obs}=N_{H,MW}+N_H(1+z)^{2.6}`$, where $`N_{H,MW}`$ is the Galactic column density along the line of sight (interpolated222We use the $`n_H`$ FTOOL at http://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl. from the maps of ), $`N_H`$ is the source-frame column density, and $`(1+z)^{2.6}`$ scales $`N_H`$ to the observer frame (e.g., ). If multiple column density measurements are available for a burst, we average their source-frame probability distributions if the measurements were made from the same data, and we take the product of these distributions if the measurements were made from independent data (GRB 000926). We plot the source-frame probability distributions, peak-normalized to one, for the bursts with detected optical afterglows in the top panel of Figure 8. The dotted curves mark bursts of unknown redshift, in which case we fix $`z=0`$: These distributions can be thought of as fuzzy lower limits (see below).

We model the column density distribution of the bursts with detected optical afterglows with a broken power law. We use a three-parameter function for flexibility, and a simple function (as opposed to, e.g., the more elegant skewed Gaussian) so the inner integral of Equation (9) can be evaluated analytically (otherwise, the error analysis would be computationally taxing). We note that the broken power law is log-normalized \[$`_{\mathrm{}}^{\mathrm{}}n(\mathrm{log}N_H)d\mathrm{log}N_H=1`$\] to facilitate comparison of the fitted distribution to the solid histogram of Figure 6 (see below).

We fit the model $`n(\mathrm{log}N_H)`$ to the data $`p_i(N_H)`$ using Bayesian inference (e.g., ). The likelihood function is given by Equation (1), where $`N=11`$ is the number of bursts with detected optical afterglows in our sample, and

$$_i=\{\begin{array}{cc}_{\mathrm{}}^{\mathrm{}}p_i(N_H)n(N_H)𝑑N_H\hfill & \text{(}z_i\text{ known)}\hfill \\ _{\mathrm{}}^{\mathrm{}}p_i(N_H^{})\left[_{N_H^{}}^{\mathrm{}}n(N_H)𝑑N_H\right]𝑑N_H^{}\hfill & \text{(}z_i\text{ unknown)}\hfill \end{array}$$
(8)
(e.g., ), or equivalently

$$_i=\{\begin{array}{cc}_{\mathrm{}}^{\mathrm{}}p_i(N_H)n(\mathrm{log}N_H)d\mathrm{log}N_H\hfill & \\ _0^{\mathrm{}}p_i(N_H^{})\left[_{\mathrm{log}N_H^{}}^{\mathrm{}}n(\mathrm{log}N_H)d\mathrm{log}N_H\right]𝑑N_H^{}\hfill & \end{array},$$
(9)
since $`n(\mathrm{log}N_H)=0`$ for $`N_H<0`$. We find that $`\mathrm{log}(a/\mathrm{cm}^2)=22.16_{0.33}^{+0.20}`$, arctan$`b=0.92_{0.35}^{+0.38}`$ ($`b1.31_{0.68}^{+2.25}`$), and arctan$`c=1.46_{0.11}^{+0.53}`$ \[$`c<1.33`$ (1 $`\sigma `$)\]. We plot the likelihood-weighted average model distribution (solid curve), which is also log-normalized, and the 1-$`\sigma `$ uncertainty in this distribution (dotted curves) in the bottom panel of Figure 8. We find that the column density distribution of the bursts with detected optical afterglows peaks around $`10^{22}`$ cm-2, and spans a factor of a few to an order of magnitude or so. These findings verify and improve upon the earlier finding of  that these column densities span $`10^{22}10^{23}`$ cm-2.

Lastly, we compare the fitted column density distribution of the bursts with detected optical afterglows to the expected column density distribution for bursts in Galactic-like molecular clouds (solid histogram, replotted from Figure 6), and find these distributions to be consistent within the 1-$`\sigma `$ uncertainties. Furthermore, the fitted distribution is not consistent with the expectation for bursts that occur in diffuse clouds, the canonical column density of which is a few times $`10^{20}`$ cm-2. Consequently, we find that all but perhaps a few bursts, dark and otherwise, probably occur in molecular clouds.