text:
# Enhancement of the upper critical field by nonmagnetic impurities in dirty two-gap superconductors

## I Introduction

Recent discovery of $`MgB_2`$ has renewed interest in superconductors with multicomponent order parameters. In particular, extensive ab-initio calculations tg1 ; tg2 , along with ample experimental evidences based on STM tg3 ; tg4 , point contact tg5 , and Raman tg6  spectroscopy, heat capacity tg7  and rf response tg8  measurements strongly indicate that $`MgB_2`$ has two distinct s-wave superconducting gaps $`\mathrm{\Delta }_\sigma (0)7.2`$mV and $`\mathrm{\Delta }_\pi (0)2.3`$mV. These gaps reside on different disconnected sheets of the Fermi surface (FS), which comprises nearly cylindrical 2D parts formed by in-plane $`\sigma `$ antibonding $`p_{xy}`$ orbitals of B, and a more isotropic 3D tubular network formed by out-of-plane $`\pi `$ bonding and antibonding $`p_z`$ orbitals of B. So far all attempts to increase the critical temperature $`T_c`$ of $`MgB_2`$ by doping have been unsuccessful rev . On the other hand, the significant potential of $`MgB_2`$ for applicationsappl  is limited by comparatively low upper critical fields $`H_{c2}^{}(0)35T`$ and $`H_{c2}^{||}(0)1519T`$ of hexagonal $`MgB_2`$ single crystals a1 ; a2 ; a3 ; a4 ; a5 ; a6 ; a7 ; a8 ; a9 ; a10 ; a11 , where the indices $``$ and $`||`$ correspond to the magnetic field H perpendicular and parallel to the ab plane, respectively. As far as $`H_{c2}`$ is concerned, it can be increased by adding nonmagnetic impurities whh ; a15 , which are especially effective in the dirty limit, $`2\pi k_BT_c<\mathrm{}/\tau `$, where $`\tau `$ is the elastic scattering time, $`\mathrm{}`$ is the Planck constant, and $`k_B`$ is the Boltzmann constant. In this case, there is a simple universal relation between the zero-temperature value $`H_{c2}(0)`$, the slope $`H_{c2}^{}=dH_{c2}/dT`$ at $`T_c`$, and the normal state residual resistivity $`\rho _n`$ whh :

$$H_{c2}(0)=0.69T_cH_{c2}^{},H_{c2}^{}=\frac{4ek_B}{\pi \mathrm{}}N_F\rho _n,$$
(1)
where $`N_F`$ is the density of states at the FS and $`e`$ is the electron charge. The fact that adding nonmagnetic impurities increases $`H_{c2}`$ has been used to significantly improve $`H_{c2}`$ of $`NbTi`$ and $`Nb_3Sn`$ compounds a15 . The same approach has been also taken to improve the high-field performance of $`MgB_2`$ in which the additional scattering was introduced by proton perk  and neutron buket  irradiation, and atomic substitutions on both B and Mg sites in thin films sat , bulk samples rev  and wires jc1 ; jc2 . For instance, in highly resistive c-axis oriented $`MgB_2`$ films sat , $`\rho _n`$ was increased from the nominal value $`1\mu \mathrm{\Omega }`$cm fin  to more than $`700\mu \mathrm{\Omega }`$cm, which in turn increased the slopes $`H_{c2}^{}`$ to rather high values $`H_{c2}^{}1`$T/K and $`H_{c2}^{}1.8`$T/K, while reducing $`T_c`$ down to $`30`$K. Based on these numbers, the extrapolation (1) would give the upper critical field, $`H_{c2}^{}(0)20`$T, which is still below $`H_{c2}(0)30`$T of $`Nb_3Sn`$a15 . However, Eq. (1) is a result of the one-gap theory whh  which does not take into account the multiple scattering channels in two-gap superconductors. In this paper I show that $`H_{c2}(0)`$ in two gap superconductors can be significantly higher than what follows from Eq. (1), thus the upper critical fields of $`MgB_2`$ alloys can in fact exceed $`H_{c2}(0)`$ of $`Nb_3Sn`$ even for the values $`H_{c2}^{}1`$T/K which have already been achievedsat ; gur .

The theory of the two-band superconductivity has been developed long ago suhl , and more recently generalized to include the effects of impuritiesimp1 ; imp2 ; imp3 ; lanis1 ; lanis2 ; goluba  and strong electron-phonon coupling tg2 ; lanis2 ; goluba ; shulga . The Fermi surface of $`MgB_2`$ provides three different impurity scattering channels: intraband scattering within each $`\sigma `$ and $`\pi `$ FS sheets and interband scattering between them. Strong impurity scattering in the dirty limit causes intraband electron diffusion over the respective FS sheets, which reduces the intrinsic anisotropy of $`\mathrm{\Delta }_\sigma `$ and $`\mathrm{\Delta }_\pi `$, but does not affect $`T_c`$, in accordance with the Anderson theorem. By contrast, the pairbreaking interband scattering reduces $`T_c`$ imp1 ; imp2 , but its effect in $`MgB_2`$ seems to be not very pronounced due to orthogonality of $`\sigma `$ and $`\pi `$ orbitalsimp3 . It is the multiple scattering channels, which provide the essential flexibility to increase the $`H_{c2}`$ of $`MgB_2`$ to a much greater extent than in one-gap superconductors not only by the usual increase of $`\rho _n`$, but also by optimizing the relative weight of the $`\sigma `$ and $`\pi `$ scattering rates. The latter could be done by selective atomic substitution on B sites (with C, O, etc. ) which mostly affect the in-plane $`\sigma `$ scattering, and Mg site (Mg vacancies, Al, etc) rev , which mostly affects out-of-plane $`\pi `$ scattering imp3 .

In this work I use a weak coupling multiband BCS model suhl  to develop a theory of magnetic properties of two-gap superconductors in the dirty limit, taking $`MgB_2`$ as an example. This approach generalizes the well-known theory developed for dirty one-gap superconductorswhh . The paper is organized as follows. In section II the quasiclassic Uzadel equations are obtained for an anisotropic two-gap superconductor with the account of both interband and intraband scattering by nonmagnetic impurities. In section III these equations are used to derive the Ginzburg-Landau equations and the critical temperature $`T_c`$ in the dirty limit. In section IV the linearized Uzadel equations are solved exactly to obtain a general equation for the upper critical field parallel to the c-axis. It is shown that the temperature dependence of $`H_{c2}(T)`$ can be very different from the one-band theory whh , and $`H_{c2}(0)`$ can be considerably higher than what follows from Eq. (1). In section V I calculate the angular dependence of $`H_{c2}(\theta )`$ and propose a scaling rule, which shows that the anisotropy ratio $`H_{c2}^{}/H_{c2}^{}`$ can both increase and decrease with temperature, depending on the relation between intraband diffusivities. In section VI the results of this work are used to address the observed temperature and orientational dependences of $`H_{c2}`$ of $`MgB_2`$.

## II Uzadel equations

A powerful tool for investigating inhomogeneous states of superconductors is the quasiclassic Eilenberger equations eilen ; book  for the Green functions $`f(𝐤,𝐫,\omega )`$, $`f^+(𝐤,𝐫,\omega )`$, $`g(𝐤,𝐫,\omega )`$ and the order parameters $`\mathrm{\Delta }(𝐤,𝐫)`$ which depend on the coordinates $`𝐫`$, the Matsubara frequency $`\omega =\pi T(2n+1)`$, and the wave vector $`𝐤`$ on the FS (hereafter the units with $`k_B=\mathrm{}=1`$ are used). The essential dependence of $`f(𝐤,𝐫,\omega )`$, $`g(𝐤,𝐫,\omega )`$ and $`\mathrm{\Delta }(𝐤,𝐫)`$ on the direction of $`𝐤`$, makes them sensitive to the shape of the FS, which greatly complicates solving the nonlinear Eilenberger equations. Thus, the Eilenberger equations are not so useful in the case of complicated FS, such as the tubular network of disconnected FS sheets of $`MgB_2`$ for which the numerical solutions of the Eliashberg equation have shown the significance of the FS anisotropy tg2 . However the situation simplifies in the dirty limit for which the impurity scattering within each FS sheet averages out the angular dependences of $`f(𝐤,𝐫,\omega )`$, $`f^+(𝐤,𝐫,\omega )`$, $`g(𝐤,𝐫,\omega )`$, making them independent of $`𝐤`$. In this case the Eilenberger equations reduce to much simpler Uzadel equations uzadel  in which all microscopic details are hidden in the electronic diffusivity tensors $`D_m^{\alpha \beta }`$ for each m-th FS sheet and the interband scattering rates $`\gamma _{mm^{}}`$ which reflect the underlying FS symmetry allen .

In this paper I consider a dirty two-gap anisotropic superconductor in the simplest case of two disconnected sheets 1 and 2 of the FS. This model corresponds to the $`\sigma `$ and $`\pi `$ bands in $`MgB_2`$ tg1 ; tg2 , for which the superconducting gaps take the constant values $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$ on the sheets 1 and 2, respectively. The Uzadel equations derived in Appendix A have the form

$`2\omega f_1D_1^{\alpha \beta }[g_1\mathrm{\Pi }_\alpha \mathrm{\Pi }_\beta f_1f_1_\alpha _\beta g_1]=`$
$`2\mathrm{\Delta }_1g_1+\gamma _{12}(g_1f_2g_2f_1)`$ (2)
$`2\omega f_2D_2^{\alpha \beta }[g_2\mathrm{\Pi }_\alpha \mathrm{\Pi }_\beta f_2f_2_\alpha _\beta g_2]=`$
$`2\mathrm{\Delta }_2g_2+\gamma _{21}(g_2f_1g_1f_2),`$ (3)
Eqs. (2) and (3) are supplemented by the self-consistency equations for the order parameters $`\mathrm{\Delta }_m=|\mathrm{\Delta }_m|\mathrm{exp}(i\phi _m)`$:

$$\mathrm{\Delta }_m=2\pi T\underset{\omega >0}{\overset{\omega _D}{}}\underset{m}{}\lambda _{mm^{}}f_m^{}(𝐫,\omega ),$$
(4)
normalization condition

$$|f_m|^2+g_m^2=1,$$
(5)
and the expression for the supercurrent density

$$J^\alpha =4\pi ieTIm\underset{\omega }{}\underset{m}{}N_mD_m^{\alpha \beta }f_m^+\mathrm{\Pi }_\beta f_m.$$
(6)
Here the band index m runs from 1 and 2, the functions $`f_m(𝐫,\omega )`$ and $`g_m(𝐫,\omega )`$ in the m-th band depend on $`𝐫`$ and $`\omega `$ but not on $`𝐤`$, $`D_m^{\alpha \beta }`$ are the intraband diffusivity tensors due to nonmagnetic impurity scattering, $`\gamma _{mm^{}}`$ are the interband scattering rates, $`N_m`$ is the partial density of states, $`𝚷=+2\pi i𝐀/\varphi _0`$, $`𝐀`$ is the vector potential, $`\varphi _0`$ is the flux quantum, and the summation over the Greek Cartesian indices is implied.

Eqs. (4) contains the matrix of the BCS superconducting coupling constants $`\lambda _{mm^{}}=\lambda _{mm^{}}^{(ep)}\mu _{mm^{}}`$, where $`\lambda _{mm^{}}^{(ep)}`$ are electron-phonon constants, and $`\mu _{mm^{}}`$ is the matrix of the Coulomb pseudopotential. Here the diagonal terms $`\lambda _{11}`$ and $`\lambda _{22}`$ quantify the intraband superconducting coupling, and off-diagonal terms $`\lambda _{12}`$ and $`\lambda _{21}`$ describe the interband coupling. The eigenvalues of $`\lambda _{mm^{}}`$ are assumed positive, the band 1 having the highest coupling constant $`\lambda _{11}`$. The indices 1 and 2 thus correspond to the $`\sigma `$ and $`\pi `$ bands of $`MgB_2`$ for which ab-initio calculations yield $`\lambda _{\sigma \sigma }0.81`$, $`\lambda _{\pi \pi }0.285`$, $`\lambda _{\sigma \pi }0.119`$, and $`\lambda _{\pi \sigma }0.09`$goluba . The mixed components $`\lambda _{12}`$ and $`\lambda _{21}`$ satisfy the symmetry relation suhl  (see also Eq. (81)):

$$N_1\lambda _{12}=N_2\lambda _{21}$$
(7)
where $`N_1`$ and $`N_2`$ are partial densities of states in the bands 1 and 2 ($`N_\pi 1.3N_\sigma `$ for $`MgB_2`$). A similar approach based on the Uzadel equations was recently proposed to describe vortices in two-gap superconductors koshgol .

Formulas for the interband scattering rates $`\gamma _{mm^{}}`$ and the intraband diffusivity tensors $`D_m^{\alpha \beta }`$ expressed in terms of microscopic material parameters are given in Appendix A. They will be used in the next sections as input parameters, which can be either calculated from first principles or extracted from the observed temperature dependences of $`H_{c2}(T)`$ and resistivity $`\rho `$, as discussed below. Here we just emphasize the fact that, although both $`D_1^{\alpha \beta }`$ and $`D_2^{\alpha \beta }`$ reflect the underlying symmetry of the FS, the features of atomic orbitals which form the bands 1 and 2 can manifest themselves in very different properties of $`D_1^{\alpha \beta }`$ and $`D_2^{\alpha \beta }`$. For instance, the principal value $`D_\sigma ^{(c)}`$ along the c-axis in $`MgB_2`$ is much smaller than two principal values $`D_\sigma ^{(a)}`$ and $`D_\sigma ^{(b)}`$ in the ab plane, because of the nearly 2D nature of the $`\sigma `$-band formed by in-plane bonding and antibonding $`p_{xy}`$ orbitals of B, as was shown by STM tg3 . By contrast, the difference in the principal values $`D_\pi ^{(c)}`$, $`D_\pi ^{(a)}`$ and $`D_\pi ^{(b)}`$ is less pronounced for the 3D $`\pi `$-band formed by the out-of-plane $`p_z`$ orbitals of B. The resulting relation $`D_\sigma ^{(c)}/D_\sigma ^{(a)}D_\pi ^{(c)}/D_\pi ^{(a)}`$ can manifest itself in anomalous temperature dependence of the anisotropy of $`H_{c2}`$, as shown in the next sections.

## III Critical temperature and Ginzburg-Landau equations

The necessity to satisfy the self-consistency condition (4) complicates solving the nonlinear Eqs. (2)-(3) for inhomogeneous $`\mathrm{\Delta }_m(𝐫)`$. The situation simplifies near $`T_c`$ where Eqs. (2)-(4) reduce to the Ginzburg-Landau (GL) equations for the order parameters $`\mathrm{\Delta }_m`$. We derive the GL equations neglecting the interband scattering terms $`\gamma _{12}\pi T_c`$ which are usually small in $`MgB_2`$ due to orthogonality of the $`\sigma `$ and $`\pi `$ orbitals imp3 . Expanding the solution of Eqs. (2) and (3) in powers of $`\mathrm{\Delta }_m`$ and its spatial derivatives, and using $`g_m1|f_m|^2/2`$, we obtain

$$f_m=\frac{\mathrm{\Delta }_m}{\omega }+\frac{D_m^{\alpha \beta }}{2\omega ^2}\mathrm{\Pi }_\alpha \mathrm{\Pi }_\beta \mathrm{\Delta }_m\frac{\mathrm{\Delta }_m|\mathrm{\Delta }_m|^2}{2\omega ^3}$$
(8)
Inserting Eqs. (8) into Eq. (4) and summing up over $`\omega `$, results in the following equations

$`\mathrm{\Delta }_1=\lambda _{11}R_1+\lambda _{12}R_2,`$ (9)
$`\mathrm{\Delta }_2=\lambda _{21}R_1+\lambda _{22}R_2,`$ (10)
$`R_m=\mathrm{\Delta }_ml+{\displaystyle \frac{\pi D_m^{\alpha \beta }}{8T}}\mathrm{\Pi }_\alpha \mathrm{\Pi }_\beta \mathrm{\Delta }_m{\displaystyle \frac{7\zeta (3)}{8\pi ^2T^2}}|\mathrm{\Delta }_m|^2\mathrm{\Delta }_m,`$ (11)
where $`R_m=_\omega f_m`$, $`l=\mathrm{ln}(2\gamma \omega _D/\pi T)`$, $`\mathrm{ln}\gamma =0.577`$ is the Euler constant, and $`\zeta (3)=1.202`$. The supercurrent density at $`TT_cT_c`$ is obtained by inserting the first term in the right hand side of Eq. (44) into Eq. (6) and summing up over $`\omega `$:

$$J^\alpha =\frac{\pi e}{4T}Im\underset{m}{}N_mD_m^{\alpha \beta }\mathrm{\Delta }_m^{}\mathrm{\Pi }_\beta \mathrm{\Delta }_m,$$
(12)
where the asterisk means complex conjugation. As follows from Eq. (12), the current density near $`T_c`$ is a sum of independent intraband contributions.

Equilibrium equations (9)-(10) for the order parameters $`\mathrm{\Delta }_m`$ are not yet the GL equations which are obtained by varying the free energy functional, $`\delta F/\delta \mathrm{\Delta }^{}=0`$. The reason is that the off-diagonal terms with respect to the band index m in Eqs. (9)-(10) violate the necessary symmetry conditions which the variational equations $`\delta F/\delta \mathrm{\Delta }^{}=0`$ must satisfy. To find the proper linear combinations of Eqs. (9) and (10) which satisfy all necessary symmetry conditions, we first solve Eqs. (9)-(11) for the quantities $`N_1R_1=(\lambda _{22}\mathrm{\Delta }_1\lambda _{12}\mathrm{\Delta }_2)N_1/w`$ and $`N_2R_2=(\lambda _{11}\mathrm{\Delta }_2\lambda _{21}\mathrm{\Delta }_1)N_2/w`$, where $`w=\lambda _{11}\lambda _{22}\lambda _{12}\lambda _{21}`$ is the determinant of the matrix $`\lambda _{mm^{}}`$. Expressing then $`R_m`$ via $`\mathrm{\Delta }_m`$ using Eq. (11), we obtain the GL equations for an anisotropic two-band superconductor:

$`a_1\mathrm{\Delta }_1+b_1\mathrm{\Delta }_1|\mathrm{\Delta }_1|^2+c_1^{\alpha \beta }\mathrm{\Pi }_\alpha \mathrm{\Pi }_\beta \mathrm{\Delta }_1+g\mathrm{\Delta }_2=0,`$ (13)
$`a_2\mathrm{\Delta }_2+b_2\mathrm{\Delta }_2|\mathrm{\Delta }_2|^2+c_2^{\alpha \beta }\mathrm{\Pi }_\alpha \mathrm{\Pi }_\beta \mathrm{\Delta }_2+g\mathrm{\Delta }_2=0.`$ (14)
Here the GL expansion coefficients are given by

$`a_1=\left(\mathrm{ln}{\displaystyle \frac{2\gamma \omega _D}{\pi T}}{\displaystyle \frac{\lambda _{22}}{w}}\right)N_1,b_1={\displaystyle \frac{7\zeta (3)N_1}{8\pi ^2T^2}},`$ (15)
$`a_2=\left(\mathrm{ln}{\displaystyle \frac{2\gamma \omega _D}{\pi T}}{\displaystyle \frac{\lambda _{11}}{w}}\right)N_2,b_2={\displaystyle \frac{7\zeta (3)N_2}{8\pi ^2T^2}},`$ (16)
$`c_1^{\alpha \beta }=\pi N_1D_1^{\alpha \beta }/8T,c_2^{\alpha \beta }=\pi N_2D_2^{\alpha \beta }/8T,`$ (17)
$`g=\lambda _{21}N_2/w=\lambda _{12}N_1/w.`$ (18)
The equality in Eq. (18) results from Eq. (7) for the interband coupling constants $`\lambda _{12}=V_iN_2`$, and $`\lambda _{21}=V_iN_1`$. The GL equations (13) and (14) can now be obtained by varying the free energy $`F=[F_1+F_2+F_i]𝑑V`$, which contains the usual GL intraband contributions

$$F_m=\frac{a_m}{2}|\mathrm{\Delta }_m|^2+\frac{b_m}{4}|\mathrm{\Delta }_m|^4+\frac{c_m^{\alpha \beta }}{2}\mathrm{\Pi }_\alpha \mathrm{\Delta }_m\mathrm{\Pi }_\beta ^{}\mathrm{\Delta }_m^{},$$
(19)
and the interband interaction term

$$F_i=\frac{g}{4}(\mathrm{\Delta }_1\mathrm{\Delta }_2^{}+\mathrm{\Delta }_2\mathrm{\Delta }_1^{}).$$
(20)
Notice that Eq. (12) for the current density obtained from the Uzadel equations coincides with the general expression $`𝐉=\delta F/\delta 𝐀`$, where $`F`$ is defined by Eqs. (19) and (20). Static GL equations (13)-(20) and their time-dependent generalization gurvin  have been suggested phenomenologically to describe various distributions of the order parameter and vortex properties in two-gap superconductorstanaka ; gurvin .

Eqs. (13) and (14) give the following equation $`a_1a_2g^2=0`$ for the critical temperature $`T_c`$, which results in the quadratic equation for $`l_c=\mathrm{ln}(2\gamma \omega _D/\pi T_c)`$:

$$1\lambda _+l_c+wl_c^2=0$$
(21)
where $`\lambda _\pm =\lambda _{11}\pm \lambda _{22}`$, and $`w=\lambda _{11}\lambda _{22}\lambda _{12}\lambda _{21}`$. Solution of Eq. (21) reproduces the well-known result by Suhl, Matthias and Walker suhl :

$$T_{c0}=1.14\omega _D\mathrm{exp}[(\lambda _+\lambda _0)/2w],$$
(22)
where $`\lambda _0=(\lambda _{}^2+4\lambda _{12}\lambda _{21})^{1/2}`$. As follows from Eq. (22), the interband coupling always increases $`T_c`$ as compared to noninteracting bands $`(\lambda _{12}=\lambda _{21}=0)`$, while the intraband impurity scattering does not affect $`T_c`$, in accordance with the Anderson theorem book .

## IV Upper critical field for $`𝐇c`$

Now we turn to the calculation of the upper critical field $`H_{c2}`$ applied along the c-axis of a hexagonal crystal. Here $`H_{c2}`$ is the maximum eigenvalue of the linearized Eqs. (2) and (3) for $`f1`$ and $`g1`$. In Appendix B, these equations are solved exactly and an equation for $`H_{c2}`$ is obtained for arbitrary relation between $`D_1`$, $`D_2`$ and $`\gamma _{mm^{}}`$. This general equation for $`H_{c2}`$ is very cumbersome, so to make the essential effects of two-gap superconductivity more transparent, I consider here a much simpler case of negligible $`\gamma _{mm^{}}`$ for which the nonmagnetic impurity scattering does not affect $`T_c`$. If $`\gamma _{12}=\gamma _{21}=0`$, the linearized Eqs. (2) and (3) take the form

$`2\omega f_1D_1^{\alpha \beta }\mathrm{\Pi }_\alpha \mathrm{\Pi }_\beta f_1=2\mathrm{\Delta }_1,`$ (23)
$`2\omega f_2D_2^{\alpha \beta }\mathrm{\Pi }_\alpha \mathrm{\Pi }_\beta f_2=2\mathrm{\Delta }_2.`$ (24)
For $`𝐇||c`$, the relevant solutions of Eqs. (23) and (24) depends only on isotropic in-plane diffusivities, $`D_m^{\alpha \beta }=D_m\delta _{\alpha \beta }`$. In the gauge $`A_y=Hx`$ the solutions of Eqs. (4), (23), and (24) are

$`f_m(x,\omega )=\mathrm{\Delta }_m(x)/(\omega +\pi HD_m/\varphi _0),`$ (25)
$`\mathrm{\Delta }_m(x)=\stackrel{~}{\mathrm{\Delta }}_m\mathrm{exp}(\pi Hx^2/\varphi _0).`$ (26)
Inserting Eqs. (25) and (26) into the gap equation (4) yields two linear equations for $`\stackrel{~}{\mathrm{\Delta }}_1`$ and $`\stackrel{~}{\mathrm{\Delta }}_2`$ in which the summation over $`\omega `$ is performed using the identity

$`2\pi T{\displaystyle \underset{\omega >0}{\overset{\omega _D}{}}}{\displaystyle \frac{1}{\omega +\mathrm{\Omega }}}=\mathrm{ln}{\displaystyle \frac{2\gamma \omega _D}{\pi T}}U\left({\displaystyle \frac{\mathrm{\Omega }}{2\pi T}}\right),`$ (27)
$`U(x)=\psi (1/2+x)+\psi (1/2),`$ (28)
where $`\psi (x)`$ is the di-gamma function. The equations for $`\stackrel{~}{\mathrm{\Delta }}_1`$ and $`\stackrel{~}{\mathrm{\Delta }}_2`$ become

$`\stackrel{~}{\mathrm{\Delta }}_1=\lambda _{11}[lU(h)]\stackrel{~}{\mathrm{\Delta }}_1+\lambda _{12}[lU(\eta h)]\stackrel{~}{\mathrm{\Delta }}_2,`$ (29)
$`\stackrel{~}{\mathrm{\Delta }}_2=\lambda _{22}[lU(\eta h)]\stackrel{~}{\mathrm{\Delta }}_2+\lambda _{21}[lU(h)]\stackrel{~}{\mathrm{\Delta }}_1,`$ (30)
where $`l=\mathrm{ln}(2\gamma \omega _D/\pi T)`$, $`h=H_{c2}D_1/2\varphi _0T`$, and $`\eta =D_2/D_1`$. The solvability condition of Eqs. (29) and (30) gives an equation for $`H_{c2}`$, in which it is convenient to express $`\omega _D`$ via $`T_c`$ using Eqs. (21) and (22). As a result, the equation for $`H_{c2}`$ takes the form

$`2a_0[\mathrm{ln}t+U(h)][\mathrm{ln}t+U(\eta h)]+`$
$`a_2[\mathrm{ln}t+U(\eta h)]+a_1[\mathrm{ln}t+U(h)]=0.`$ (31)
where $`a_1=1+\lambda _{}/\lambda _0`$, $`a_2=1\lambda _{}/\lambda _0`$, and $`a_0=w/\lambda _0`$. For equal diffusivities, $`\eta =1`$, Eq. (31) reduces to the equation $`\mathrm{ln}t+U(h)=0`$ for $`H_{c2}`$ in one-gap dirty superconductors whh . A somewhat different form of Eq. (31) was recently obtained from the Uzadel equations koshgol .

Shown in Fig. 1 are $`H_{c2}(T)`$ curves calculated from Eq. (31) for different diffusivity ratios $`\eta =D_2/D_1`$. One can clearly see the evolution of the temperature dependence of $`H_{c2}(T)`$ from the dirty-limit one-band behavior at $`D_1=D_2`$, to rather different $`H_{c2}(T)`$ curves which have portions with both upward and downward curvatures for either $`D_1D_2`$ and $`D_1D_2`$. In the latter two cases, the zero-temperature value $`H_{c0}(0)`$ can be significantly higher than the one-gap extrapolation (1). To show this, we obtain expressions for $`H_{c2}(T)`$ near $`T_c`$ and $`T=0`$. For $`TT_c`$, Eq. (31) can be expanded in small terms $`h`$, using $`U(h)\pi ^2h/2`$, and then solving for $`H_{c2}`$:

$$H_{c2}=\frac{8\varphi _0(T_cT)}{\pi ^2(a_1D_1+a_2D_2)},$$
(32)
For $`D_1=D_2`$, Eq. (32) yields the result of the one-band theory, $`H_{c2}=4(T_cT)\varphi _0/\pi ^2D`$ whh . Furthermore, if $`\lambda _{12}=\lambda _{21}=0`$, then $`a_1=2`$ and $`a_2=0`$, thus $`H_{c2}(T)4(T_cT)\varphi _0/\pi ^2D_1`$ is determined by $`D_1`$ for the band with the highest coupling constant $`\lambda _{11}`$. However, if all constants $`\lambda _{mm^{}}`$ are finite, but $`D_1`$ and $`D_2`$ are very different, then $`H_{c2}`$ is determined by the maximum diffusivity, for example, by $`D_2D_1`$ if the scattering in the ”weak” band 2 dominates.

The zero-temperature $`H_{c2}(0)`$ can be obtained using the asymptotic behavior of $`U(x)\mathrm{ln}(4\gamma h)`$ for $`h\mathrm{}`$. Then all terms proportional to $`\mathrm{ln}t`$ cancel out, and Eq. (31) reduces to a quadratic equation for $`\mathrm{ln}H_{c2}`$, whence

$`H_{c2}(0)={\displaystyle \frac{\varphi _0T_c}{2\gamma \sqrt{D_1D_2}}}\mathrm{exp}({\displaystyle \frac{g}{2}}),`$ (33)
$`g=\left({\displaystyle \frac{\lambda _0^2}{w^2}}+\mathrm{ln}^2{\displaystyle \frac{D_2}{D_1}}+{\displaystyle \frac{2\lambda _{}}{w}}\mathrm{ln}{\displaystyle \frac{D_2}{D_1}}\right)^{1/2}{\displaystyle \frac{\lambda _0}{w}}.`$ (34)
For $`D_1=D_2`$, Eqs. (33)-(34) again yield the result $`H_{c2}(0)=\varphi _0T_c/2\gamma D`$ of the one-band dirty limit whh . However for $`D_1D_2`$, Eqs. (31)-(34) predict a significant enhancement of $`H_{c2}(0)`$ as compared to $`H_{c2}(0)=0.69H_{c2}^{}T_c`$ for the symmetric case $`D_1=D_2`$. Indeed, in the limit of very different diffusivities, we obtain

$`H_{c2}(0)={\displaystyle \frac{\varphi _0T_c}{2\gamma D_2}}e^{\frac{\lambda _{}\lambda _0}{2w}},D_2D_1e^{\frac{\lambda _0}{w}},`$ (35)
$`H_{c2}(0)={\displaystyle \frac{\varphi _0T_c}{2\gamma D_1}}e^{\frac{\lambda _{}\lambda _0}{2w}},D_1D_2e^{\frac{\lambda _0}{w}}.`$ (36)
Unlike $`H_{c2}(T)`$ at $`TT_c`$, the limiting value of $`H_{c2}(0)`$ is determined by the minimum diffusivity. It is the feature of a two-band superconductor, which causes both the upward curvature and the enhancement of $`H_{c2}`$ at low T in Fig. 1. In the limit $`D_2D_1`$, the diffusivity $`D_2`$ does not affect $`H_{c2}`$ practically for all temperatures $`T<T_c`$, except a narrow region $`TT_c`$, but for $`D_20`$, the upper critical field $`H_{c2}(T)`$ diverges at $`T=0`$.

## V Anisotropy of $`H_{c2}`$

.

### V.1 General equations

To address the anisotropy of $`H_{c2}`$, observed in $`MgB_2`$, we consider a uniaxial crystal in a field $`𝐇`$ inclined by the angle $`\theta `$ with respect to the c-axis. Then the tensors $`D_{\alpha \beta }`$ in Eq. (2)-(3) have two equal principal values $`D_m^{(a)}`$ in the ab plane and a different value $`D_m^c`$ along the c-axis. If the z-axis is taken along $`𝐇`$ and the a-axis coincides with the y-axis, $`D_m^{\alpha \beta }`$ has the following nonzero components

$`D_m^{yy}`$ $`=`$ $`D_m^{(a)},`$ (37)
$`D_m^{xx}`$ $`=`$ $`D_m^{(a)}\mathrm{cos}^2\theta +D_m^{(c)}\mathrm{sin}^2\theta ,`$ (38)
$`D_m^{zz}`$ $`=`$ $`D_m^{(c)}\mathrm{cos}^2\theta +D_m^{(a)}\mathrm{sin}^2\theta ,`$ (39)
$`D_m^{xz}`$ $`=`$ $`D_m^{(zx)}=(D_m^{(c)}D_m^{(a)})\mathrm{sin}\theta \mathrm{cos}\theta `$ (40)
In the gauge $`A_y=Hx`$, the Uzadel equations (23) and (24) take the form

$`D_m^{xx}_{xx}f_mD_m^{yy}(_y+2\pi iHx/\varphi _0)^2f_m`$
$`D_m^{zz}_{zz}f_m2D_m^{xz}_{xz}f_m+2\omega f_m=2\mathrm{\Delta }_m`$ (41)
The solution of this equation has the form $`f_m(𝐫)=\mathrm{exp}(ik_zzik_zxD_m^{xz}/D_m^{xx}+ik_yy)f_m(x)`$, where $`f_m(x)`$ is determined by the following equation

$`D_m^{xx}f_m^{\prime \prime }+D_m^{yy}(2\pi Hx/\varphi _0)^2f_m+2\omega f_m`$
$`+k_z^2(D_m^{zz}D_m^{xz2}/D_m^{xx})f_m=2\mathrm{\Delta }_m,`$ (42)
the prime denotes differentiation over x, and $`k_y`$ was absorbed by the shift of x ll . The upper critical field $`H_{c2}`$ is the maximum eigenvalue of Eq. (42). Because of the stability condition $`D_m^{zz}D_m^{xx}>D_m^{xz2}`$, the term quadratic in $`k_z`$ always decreases the eigenvalues of Eq. (42), so $`H_{c2}`$ corresponds to $`k_z=0`$. In this case the lowest eigenfunction $`\phi _0`$ of Eq. (42) is

$$\phi _0^{(m)}e^{q_m^2x^2/2},q_m=\left[\frac{2\pi H}{\varphi _0}\right]^{1/2}\left[\frac{D_m^{yy}}{D_m^{xx}}\right]^{1/4}$$
(43)
An interesting situation occurs if the ratio $`D_m^{yy}/D_m^{xx}`$ is different for different bands, as characteristic of the $`\sigma `$ and $`\pi `$ bands of $`MgB_2`$. In this case the solution $`f_m(x)\mathrm{\Delta }_m(x)\phi _0(x)`$ satisfies Eq. (43), but does not satisfy the self-consistency condition (4), because $`q_1q_2`$. This feature of two-gap superconductors can essentially complicate the calculation of $`H_{c2}`$ for inclined fields, as compared to the case $`𝐇c`$.

$`H_{c2}(\theta )`$ for an arbitrary field orientation can be obtained by expanding the solution of the inhomogeneous Eq. (43) in the series of orthogonal normalized eigenfunctions $`\phi _n(q_mx)`$ for the different Landau levels n:

$`f_m={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\phi _n(q_mx)}{\omega +(2n+1)\mathrm{\Omega }_m}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\mathrm{\Delta }_m(u)\phi _n(q_mu)𝑑u,`$ (44)
$`\phi _n(q_mx)=\left({\displaystyle \frac{q_m}{\sqrt{\pi }2^nn!}}\right)^{1/2}e^{q_m^2x^2/2}H_n(q_mx),`$ (45)
$`\mathrm{\Omega }_m=(D_m^{xx}D_m^{yy})^{1/2}\pi H/\varphi _0,`$ (46)
where $`H_n(qx)`$ is the Hermitian polynomials ll . Inserting Eq. (45) into Eq. (4) and summing up over $`\omega `$, gives $`H_{c2}`$ as the maximum eigenvalue of a matrix $`M_{nn^{}}`$, for which the determinant $`M=0`$. The matrix $`M_{nn^{}}`$ is given by (see Appendix C):

$`M_{nn^{}}=`$
$`(1\lambda _{11}F_n^{(1)})\delta _{nn^{}}\lambda _{12}\lambda _{21}F_n^{}^{(1)}{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{V_{ns}V_{n^{}s}F_s^{(2)}}{1\lambda _{22}F_s^{(2)}}},`$ (47)
$`F_n^{(m)}=\mathrm{ln}{\displaystyle \frac{2\gamma \omega _D}{\pi T}}\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\mathrm{\Omega }_m(2n+1)}{2\pi T}}\right)+\psi \left({\displaystyle \frac{1}{2}}\right),`$ (48)
$`V_{ns}={\displaystyle \phi _n(q_1x)\phi _s(q_2x)𝑑x}`$ (49)
where $`\delta _{nn^{}}`$ is the Kronecker symbol. For $`𝐇||c`$ or for equal anisotropy ratios $`D_1^{(a)}/D_1^{(c)}=D_2^{(a)}/D_2^{(c)}`$, we have $`q_1=q_2`$, so $`V_{nn^{}}=\delta _{nn^{}}`$, and all off-diagonal terms of the matrix M vanish. In this case the equation for the maximum eigenvalue $`H_{c2}`$ is simply $`M_{00}=0`$, which reduces to Eq. (31) of the previous paragraph. For arbitrary field orientation and general relations between the coupling constants and the anisotropy ratios $`D_m^{(a)}/D_m^{(c)}`$ the equation $`M=0`$ for $`H_{c2}`$ can be solved numerically. In this case the matrix $`M_{nn^{}}`$ should be first diagonalized to $`\stackrel{~}{M}_n\delta _{nn^{}}`$, and then $`H_{c2}`$ can be found as a root of $`\stackrel{~}{M}_0=0`$.

The equation for $`H_{c2}`$ greatly simplifies for the moderate anisotropy characteristic of $`MgB_2`$. To quantify the degree of anisotropy, we introduce the asymmetry parameter $`\zeta =[(q_2^2q_1^2)/(q_2^2+q_1^2)]^2`$, that is,

$$\zeta =\left[\frac{\sqrt{D_2^{xx}/D_2^{yy}}\sqrt{D_1^{xx}/D_1^{yy}}}{\sqrt{D_2^{xx}/D_2^{yy}}+\sqrt{D_1^{xx}/D_1^{yy}}}\right]^2$$
(50)
As shown in the next section, the observed anisotropy of $`H_{c2}`$ of $`MgB_2`$ near $`T_c`$ indicates that $`D_2^{(a)}D_2^{(c)}`$ and $`D_1^{(a)}4D_1^{(c)}`$, $`\zeta 1/9`$. In this case the matrix elements $`V_{0,2s}\zeta ^s`$ rapidly decreases with $`s`$ (see Appendix C). By contrast, the factor $`F_s^{(2)}/(1\lambda _{22}F_s^{(2)})`$, which varies much slower with s, can be replaced by its value at $`s=0`$ and taken out of the sum. Indeed, near $`T_c`$ where $`\mathrm{\Omega }_m0`$, the function $`F_s\mathrm{ln}(2\gamma \omega _D/\pi T)`$ is independent of $`s`$ over a very wide range of s, so the above procedure becomes exact. Even for $`T=0`$, the function $`F\mathrm{ln}[(2s+1)\mathrm{\Omega }_m/T_{c0}]`$ varies much weaker as compared to the exponential decay of $`V_{0s}^2`$, so the replacement of $`F_s^{(2)}/(1\lambda _{22}F_s^{(2)})`$ by a constant is still a very good approximation. The remaining summation over s in Eq. (47) can then be performed exactly, using the sum rule due to the orthogonality of the eigenfunctions $`\phi _n`$ (see Appendix C):

$$\underset{s=0}{\overset{\mathrm{}}{}}V_{ns}V_{n^{}s}=\delta _{nn^{}}$$
(51)
Because of the $`\delta _{nn^{}}`$ in Eq. (51), all off-diagonal elements of the matrix M vanish, thus the equation $`\stackrel{~}{M}_{00}=0`$ for $`H_{c2}`$ takes the form

$$(1\lambda _{11}F_0^{(1)})(1\lambda _{22}F_0^{(2)})=\lambda _{12}\lambda _{21}F_0^{(1)}F_0^{(2)}$$
(52)
In this equation all $`q_m(\theta )`$ cancel out, and the anisotropy manifests itself only via $`\mathrm{\Omega }_m(\theta )`$ in the functions $`F_0^{(m)}`$ defined by Eqs. (46) and (48). After subtraction Eq. (21) for $`l_c=\mathrm{ln}(2\gamma \omega _D/\pi T_{c0})`$, from Eq. (52), the latter reduces to Eq. (31) if the following rescaling $`D_m=(D_m^{xx}D_m^{yy})^{1/2}`$ is made. Therefore, all formulas of the previous section can be generalized to the anisotropic case of inclined field by replacing $`D_1`$ and $`D_2`$ with the angular-dependent diffusivities $`D_1(\theta )`$ and $`D_2(\theta )`$ for both bands:

$$D_m(\theta )=[D_m^{(a)2}\mathrm{cos}^2\theta +D_m^{(a)}D_m^{(c)}\mathrm{sin}^2\theta ]^{1/2}$$
(53)
For example, Eqs. (31) and (53) give the following angular dependence of $`H_{c2}(\theta )`$ near $`T_c`$:

$$H_{c2}(\theta )=\frac{8\varphi _0(T_cT)}{\pi ^2[a_1D_1(\theta )+a_2D_2(\theta )]},$$
(54)
For $`D_1=D_2`$, Eq. (54) reduces to the well-known result of the anisotropic GL theory, $`H_{c2}=4\varphi _0(T_cT)/\pi ^2D(\theta )`$ for a one-gap superconductor in the dirty limit gork .

### V.2 Temperature dependence of the anisotropy of $`H_{c2}`$

Eqs. (31) and (53) determine both the temperature and angular dependences of $`H_{c2}`$. Here the anisotropy of $`H_{c2}(\theta )`$ essentially depends on both T and the in-plane diffusivity ratio $`D_2^{(a)}/D_1^{(a)}`$. For instance, let us consider a simpler case of isotropic diffusivity $`D_2^{(a)}=D_2^{(c)}`$ for the band 2 and anisotropic diffusivity with $`D_1^{(a)}D_2^{(c)}`$ for the band 1, which qualitatively models the 3D $`\pi `$ band and the 2D $`\sigma `$ band of $`MgB_2`$. Then for $`D_1^{(a)}D_2^{(ab)}`$, the angular dependence of $`H_{c2}`$ is most pronounced at $`T_c`$ were it is determined by the large anisotropic $`D_1(\theta )`$ (see Eq. (54), while at lower T the angular dependence of $`H_{c2}`$ decreases, as it is mostly affected by the nearly isotropic small $`D_2`$ (see Eqs. (33)-(36)). However for $`D_1^{(a)}D_2^{(a)}`$, the situation reverses: the upper critical field is nearly isotropic at $`T_c`$, becoming more anisotropic at lower $`T`$. These two different types of angular dependencies of $`H_{c2}(\theta )`$ are shown in Fig. 2.

The above arguments indicate that the temperature dependence of the anisotropy parameter $`H_{c2}^{||}/H_{c2}^{}`$ can be rather sensitive to the ratio $`D_1/D_2`$, which in turn is determined by the type of intraband scatterers. Indeed, if the intraband scattering is dominated by impurities (like C or N) which mostly reduce $`D_1`$ in the main $`\sigma `$-band, ($`D_1D_1`$), then, according to Eqs. (32) and (33), the ratio $`H_{c2}^{||}/H_{c2}^{}`$ is minimum at $`T_c`$ and increases as T decreases. This behavior has indeed been observed on many, mostly clean, $`MgB_2`$ samples a1 ; a2 ; a3 ; a4 ; a5 ; a6 ; a7 ; a8 ; a9 ; a10 ; a11 . By contrast, if impurities (like Al, Mg vacancies) cause the strongest intraband scattering in the ”weak” $`\pi `$-band, $`(D_2D_1)`$, then the anisotropy parameter $`H_{c2}^{||}/H_{c2}^{}`$ would decrease as T decreases, The temperature dependence of $`H_{c2}^{||}/H_{c2}^{}`$ for these limiting cases is shown in Fig. 3.

## VI Discussion

The results of this work show that the magnetic behavior of two-gap superconductors can essentially depend on the intraband diffusivity ratio $`\eta =D_2^{(a)}/D_1^{(a)}`$. This fact can have important consequences for $`MgB_2`$ for which $`H_{c2}`$ can be increased to a much greater extent than in one-band superconductors not just by the usual increase of $`\rho _n`$, but also by tuning the ratio of intraband scattering rates via selective atomic substitutions on both Mg and B sites. Actually, this possibility is naturally provided by the electron structure of $`MgB_2`$ for which substitutions on Mg site mostly disturb the out-of-plane $`p_z`$ orbitals of B thus increasing the $`\pi `$ scattering. Likewise, atomic substitutions on the B sites mostly increase scattering in the 2D $`\sigma `$ band. As shown above, the upward curvature of $`H_{c2}(T)`$ for $`D_2^{(a)}D_1^{(a)}`$ can make the zero-temperature $`H_{c2}(0)`$ significantly higher than what the one-gap estimate (1) suggests. For example, if $`D_\pi =0.05D_\sigma `$ as in Fig. 1a, the field $`H_{c2}(0)2.3T_cH_{c2}^{}`$ could approach the paramagnetic limit of $`MgB_2`$, $`70`$T for the values $`H_{c2}^{}1T/K`$ and $`T_c30`$K observed on dirty c-axis oriented $`MgB_2`$ filmssat . In fact, recent high-field measurements on this film gave $`H_{c2}^{}33.5T`$ and $`H_{c2}^{}48T`$ which corresponds $`D_\pi 0.12D_\sigma `$ gur . Thus, dirty $`MgB_2`$ can exhibit upper critical fields significantly higher than $`H_{c2}(0)30`$T of $`Nb_2Sn`$ a15 , which can be very important for applications appl .

The intraband diffusivities $`D_1`$ and $`D_2`$ can be either calculated from first principles using Eqs. (80), or extracted from measurements of $`H_{c2}`$ and $`\rho _n`$. For negligible interband scattering, $`D_1`$ and $`D_2`$ can be obtained using Eq. (32) for $`H_{c2}^{}`$ and the relation between intraband diffusivities and the residual conductivity $`\sigma =1/\rho _n`$:

$$\sigma _{\alpha \beta }=e^2\underset{m}{}N_mD_m^{\alpha \beta }$$
(55)
For $`𝐇c`$, solving two linear equations (32) and (55) expresses $`D_1`$ and $`D_2`$ in terms of observed $`H_{c2}^{}`$ and $`\rho _n`$:

$`D_1={\displaystyle \frac{1}{a_1N_2a_2N_1}}\left({\displaystyle \frac{8k_B\varphi _0N_2}{\pi ^2\mathrm{}H_{c2}^{}}}{\displaystyle \frac{a_2}{e^2\rho _n}}\right),`$ (56)
$`D_2={\displaystyle \frac{1}{a_1N_2a_2N_1}}\left({\displaystyle \frac{a_1}{e^2\rho _n}}{\displaystyle \frac{8k_B\varphi _0N_1}{\pi ^2\mathrm{}H_{c2}^{}}}\right),`$ (57)
where the normal units are restored.

Using the results of ab-initio calculationsimp3 ; goluba  for $`MgB_2`$, $`\lambda _{\sigma \sigma }0.81`$, $`\lambda _{\pi \pi }0.285`$, $`\lambda _{\sigma \pi }0.119`$, $`\lambda _{\pi \sigma }0.09`$, $`N_\sigma 0.3`$ states/cell eV, and $`N_\pi 0.41`$ states/cell eV, we obtain $`\lambda _{}=\lambda _{\sigma \sigma }\lambda _{\pi \pi }=0.525`$, $`\lambda _0=(\lambda _{}^2+4\lambda _{\pi \sigma }\lambda _{\sigma \pi })^{1/2}=0.564`$, thus $`a_1=1+\lambda _{}/\lambda _0=1.93`$, $`a_2=1\lambda _{}/\lambda _0=0.07`$, and $`a_1N_2a_2N_1`$. Since both $`D_1`$ and $`D_2`$ must be positive, the slope $`B_{c2}^{}`$ for a given $`\rho _n`$ in confined within the following limits

$$\frac{8ek_B}{\pi a_1}N_1\rho _n<B_{c2}^{}<\frac{8ek_B}{\pi a_2}N_2\rho _n.$$
(58)
For negligible interband coupling $`a_20`$, Eq. (56) reduces to the second Eq. (1), of the one-gap theory, whereas the stability condition (58) requires that $`B_{c2}^{}`$ cannot be smaller than $`4ek_BN_1\rho _n/\pi `$. This condition gives $`B_{c2}^{}>0.6T/K`$ for $`\rho _n=1m\mathrm{\Omega }`$cm, and $`N_1=0.3`$ states/cell eV imp3 , thus $`B_{c2}(0)`$ for $`T_c=30K`$ sat  cannot be smaller than $`0.7B_{c2}^{}T_c=12.6`$T even for the symmetric case $`D_1=D_2`$ which gives the minimum $`H_{c2}(0)`$ for a given $`H_{c2}^{}`$. For asymmetric case, $`D_\pi D_\sigma `$ characteristic of dirty c-axis oriented $`MgB_2`$ films sat , $`H_{c2}(0)`$ can considerably exceed this minimum gur .

Another property for which the diffusivity ratio $`D_2^{(a)}/D_1^{(a)}`$ is essential is the temperature dependence of the anisotropy of $`H_{c2}`$. As shown in the previous section, the anisotropy parameter $`H_{c2}^{||}/H_{c2}^{}`$ increases as T decreases for $`D_\sigma D_\pi `$, but decreases as T decreases for $`D_\sigma D_\pi `$ (see Fig. 3). Both cases are in stark contrast with the one-gap theory in which the anisotropy ratio for both the upper critical field $`H_{c2}^{||}/H_{c2}^{}`$ and the lower critical field $`H_{c1}^{||}/H_{c1}^{}`$ is determined by the ratio of the effective masses $`m_c/m_{ab}`$ and is temperature independent. By contrast, the two-gap $`MgB_2`$ exhibits a strong temperature dependence of $`H_{c2}^{||}/H_{c2}^{}`$ which is rather different from the temperature dependence of $`H_{c1}^{||}/H_{c1}^{}`$, as was shown both experimentally and theoreticallya2 ; lanis1 ; lanis2 . For the dirty c-axis oriented films, the observed $`H_{c2}^{||}/H_{c2}^{}`$ is about 2 sat , and $`D_\pi 0.1D_\sigma `$gur . In this case the anisotropy of $`H_{c2}`$ near $`T_c`$ is entirely determined by the anisotropy of $`D_\sigma `$, thus from Eq. (54) it follows that $`D_\sigma ^{(a)}4D_\sigma ^{(c)}`$.

## VII Acknowlegments

This work was supported by the NSF MRSEC (DMR 9214707), and by AFOSR MURI (49620-01-1-0464). I am grateful to D. C. Larbalestier for illuminating discussions.

## Appendix A Derivation of the Uzadel equations

We derive the Uzadel equations for an anisotropic two-gap superconductor from the general Eilenberger equations for the quasiclassical Green functions $`f(𝐯,𝐫,\omega )`$, $`f^+(𝐯,𝐫,\omega )`$ and $`g(𝐯,𝐫,\omega )`$ eilen

$`(2\omega +𝐯_𝐤𝚷)f(𝐤)`$ $``$ $`2\mathrm{\Delta }g(𝐤)=`$
$`n_i{\displaystyle }{\displaystyle \frac{dA_𝐪}{v_𝐪}}|u_{\mathrm{𝐤𝐪}}|^2[g(𝐤)f(𝐪)`$ $``$ $`f(𝐤)g(𝐪)],`$ (59)
$`(2\omega 𝐯_𝐤𝚷^{})f^+(𝐤)`$ $``$ $`2\mathrm{\Delta }^{}g(𝐤)=`$
$`n_i{\displaystyle }{\displaystyle \frac{dA_𝐪}{v_𝐪}}|u_{\mathrm{𝐤𝐪}}|^2[g(𝐤)f^+(𝐪)`$ $``$ $`f^+(𝐤)g(𝐪)].`$ (60)
Here $`dA_𝐤/v_𝐤`$ means integration over the FS with the local density of states $`1/v_𝐤`$, the integral terms account for nonmagnetic impurity scattering, $`u_{𝐤,𝐪}`$ is the scattering amplitude, $`n_i`$ is the density of impurities, $`v_𝐤`$ is the normal group velocity on an anisotropic FS, the wave vectors $`𝐤`$ and $`𝐪`$ lie on the FS, and the asterisk means complex conjugate. Eqs. (59) and (60) are supplemented by the normalization condition

$$ff^++|g|^2=1,$$
(61)
the equation for the superconducting gap

$$\mathrm{\Delta }(𝐤,𝐫)=2\pi T\underset{\omega >0}{\overset{\omega _D}{}}\frac{dA_𝐪}{v_𝐪}V(𝐤,𝐪)f(𝐪,𝐫,\omega ),$$
(62)
and the expression for the current density

$$𝐉=4\pi TeIm\underset{\omega >0}{\overset{\omega _D}{}}\frac{dA_𝐪}{v_𝐪}𝐯_𝐪g(𝐪,𝐫,\omega ),$$
(63)
where $`V(𝐤,𝐪)`$ is the pairing potential.

Now we define the Green functions $`f_1`$ and $`f_2`$ on two separate sheets of the FS, and write the collision integrals in the right hand side of Eqs. (59) and (60) in the form

$`St_1=n_i{\displaystyle \frac{dA_{𝐪_\mathrm{𝟏}}}{v_{𝐪_\mathrm{𝟏}}}|u_{𝐤_\mathrm{𝟏}𝐪_\mathrm{𝟏}}|^2[g_1(𝐤_\mathrm{𝟏})f_1(𝐪_\mathrm{𝟏})f_1(𝐤_\mathrm{𝟏})g_1(𝐪_\mathrm{𝟏})]}`$
$`+n_i{\displaystyle \frac{dA_{𝐪_\mathrm{𝟐}}}{v_{𝐪_\mathrm{𝟐}}}|u_{𝐤_\mathrm{𝟏}𝐪_\mathrm{𝟐}}|^2[g_1(𝐤_\mathrm{𝟏})f_2(𝐪_\mathrm{𝟐})f_1(𝐤_\mathrm{𝟏})g_2(𝐪_\mathrm{𝟐})]}`$ (64)
$`St_2=n_i{\displaystyle \frac{dA_{𝐪_\mathrm{𝟐}}}{v_{𝐪_\mathrm{𝟐}}}|u_{𝐤_\mathrm{𝟐}𝐪_\mathrm{𝟐}}|^2[g_2(𝐤_\mathrm{𝟐})f_2(𝐪_\mathrm{𝟐})f_2(𝐤_\mathrm{𝟐})g_2(𝐪_\mathrm{𝟐})]}`$
$`+n_i{\displaystyle \frac{dA_{𝐪_\mathrm{𝟏}}}{v_{𝐪_\mathrm{𝟏}}}|u_{𝐤_\mathrm{𝟐}𝐪_\mathrm{𝟏}}|^2[g_2(𝐤_\mathrm{𝟐})f_1(𝐪_\mathrm{𝟏})f_2(𝐤_\mathrm{𝟐})g_1(𝐪_\mathrm{𝟏})]}`$ (65)
The first integral term in Eqs. (64) and (65) describe the intraband scattering for which the wave vectors $`𝐤`$ and $`𝐪`$ lie on the same sheet (1 or 2) of the FS, while the second integral term describes the interband scattering for which $`𝐤`$ and $`𝐪`$ lie on the different FS sheets.

To derive the equations for $`f_m(𝐫,\omega )`$, $`f_m^+(𝐤,\omega )`$, $`g_m(𝐤,𝐫,\omega )`$, and $`\mathrm{\Delta }_m(𝐫)`$ in the dirty limit, we follow the procedure developed for one-band superconductors, expanding $`f(𝐤,𝐫,\omega )`$ and $`g(𝐤,𝐫,\omega )`$, in spherical harmonics, and keeping only the first dipole term uzadel

$`f_m(𝐤,𝐫,\omega )=f_m(𝐫,\omega )+(𝐯_m\delta 𝐟_m(𝐫,\omega )),`$ (66)
$`g_m(𝐤,𝐫,\omega )=g_m(𝐫,\omega )+(𝐯_m\delta 𝐠_m(𝐫,\omega )),`$ (67)
where the band index $`m`$ runs from 1 to 2, and the last terms in the right hand side describe first order anisotropic corrections to $`f(𝐤,𝐫,\omega )`$ and $`g(𝐤,𝐫,\omega )`$. These corrections are small if $`6\pi T_cv_m^2/D_m`$, where $`D_m`$ are intraband diffusivities (see below). The vector corrections $`\delta 𝐠_m`$ and $`\delta 𝐟_m`$ are related through the normalization condition (61):

$$2g_m\delta 𝐠_m=f_m^+\delta 𝐟_mf_m\delta 𝐟_m^+$$
(68)
Now we insert Eqs. (66) and (67) into Eqs. (59)-(60) and (64)-(65) and integrate over the respective FS sheets: $`𝑑A_𝐤/v_m(𝐤)`$, using $`𝐯_m(𝐤)𝑑A_𝐤=0`$. Since all $`f_m`$ and $`\delta 𝐟_m`$ are independent of $`𝐯`$, this yields

$`2\omega f_12\mathrm{\Delta }_1g_1+r_1^{\alpha \beta }\mathrm{\Pi }_\alpha \delta f_1^\beta =(g_1f_2g_2f_1)\gamma _{12},`$ (69)
$`2\omega f_22\mathrm{\Delta }_2g_2+r_2^{\alpha \beta }\mathrm{\Pi }_\alpha \delta f_2^\beta =(g_2f_1g_1f_2)\gamma _{21},`$ (70)
$`r_m^{\alpha \beta }={\displaystyle v_m^\alpha (𝐤)v_m^\beta (𝐤)\frac{dA_𝐤}{v_m(𝐤)N_m}}`$ (71)
where $`N_m=𝑑A_𝐤/v_m(𝐤)`$ are partial densities of states, and $`\gamma _{mm^{}}`$ are the interband scattering rates allen

$$\gamma _{mm^{}}=\frac{n_i}{N_m}\frac{dA_𝐪}{v_m(𝐪)}\frac{dA_𝐤}{v_m^{}(𝐤)}|u_{𝐤,𝐪}|^2$$
(72)
Then we substitute Eqs. (66) and (67) into Eqs. (59)-(60) and (64)-(65), multiply them by $`𝐯_m(𝐤)`$ and integrate over the respective FS sheets, using the fact that $`𝐯_1(𝐤)𝐯_2(𝐤)𝑑A_𝐤=0`$. In the dirty limit, this yields the following relations

$`\mathrm{\Gamma }_m^{\alpha \beta }(f_m\delta g_m^\beta g_m\delta f_m^\beta )=r_m^{\alpha \beta }\mathrm{\Pi }_\beta f_m,`$ (73)
$`\mathrm{\Gamma }_m^{\alpha \beta }(f_m^+\delta g_m^\beta g_m\delta f_m^{+\beta })=r_m^{\alpha \beta }\mathrm{\Pi }_\beta ^{}f_m^+,`$ (74)
Here the tensor $`\mathrm{\Gamma }_m^{\alpha \beta }`$ is defined by

$$\mathrm{\Gamma }_m^{\alpha \beta }=\frac{n_i}{N_m}\frac{|u_{𝐪,𝐤}|^2v_m^\alpha (𝐤)[v_m^\beta (𝐤)v_m^\beta (𝐪)]}{v_m(𝐤)v_m(𝐪)}𝑑A_𝐤𝑑A_𝐪$$
(75)
Eqs. (73) and (74) can be rewritten as

$`f_m\delta g_m^\alpha g_m\delta f_m^\alpha =F_m^{\alpha \beta }\mathrm{\Pi }_\beta f_m,`$ (76)
$`f_m^+\delta g_m^\alpha g_m\delta f_m^{+\alpha }=F_m^{\alpha \beta }\mathrm{\Pi }_\beta ^{}f_m^+,`$ (77)
where $`F_m^{\alpha \beta }=(\mathrm{\Gamma }_m^{\alpha \gamma })^1r_m^{\gamma \beta }`$. Multiplying Eq. (76) by $`f_m^+`$ and Eq. (77) by $`f_m`$, adding, and then using Eq. (68), we obtain

$$2\delta g_m^\alpha =F_m^{\alpha \beta }(f_m^+\mathrm{\Pi }_\beta f_mf_m\mathrm{\Pi }_\beta ^{}f_m^+)$$
(78)
Using Eq. (73), (74), (78), and the normalization condition $`g_m^2+|f_m|^2=1`$, we obtain

$$\delta f_m^\alpha =F_m^{\alpha \beta }(g_m\mathrm{\Pi }_\beta f_mf_m_\beta g_m).$$
(79)
Inserting Eq. (79) and its conjugate into Eqs. (69) and (70), we arrive at Eqs. (2) and (3) in which the intraband diffusivity tensors are given by

$$D_m^{\alpha \beta }=r_m^{\alpha \gamma }(\mathrm{\Gamma }_m^{\gamma \mu })^1r_m^{\mu \beta }.$$
(80)
Averaging the gap equation (62) over the FS, yields Eq. (4) in which the coupling constants $`\lambda _{mm^{}}`$ are given by

$$\lambda _{mm^{}}=\frac{1}{N_m}\frac{dA_𝐤dA_𝐪}{v_m(𝐤)v_m^{}(𝐪)}V(𝐤,𝐪),$$
(81)
whence $`N_1\lambda _{12}=N_2\lambda _{21}`$.

## Appendix B Equations for $`H_{c2}`$ and $`T_c`$ with the account of interband scattering

For $`H||c`$, the linearized Eqs. (2) and (3) take the form

$`\omega f_1{\displaystyle \frac{D_1}{2}}\mathrm{\Pi }^2f_1=\mathrm{\Delta }_1+(f_2f_1)\gamma _{12},`$ (82)
$`\omega f_2{\displaystyle \frac{D_2}{2}}\mathrm{\Pi }^2f_2=\mathrm{\Delta }_2+(f_1f_2)\gamma _{21},`$ (83)
In the gauge $`A_y=Hx`$, the solution of Eqs. (82) and (83) is

$$f_m(x)=\stackrel{~}{f}_me^{q^2x^2/2},\mathrm{\Delta }_m(x)=\stackrel{~}{\mathrm{\Delta }}_me^{q^2x^2/2}.$$
(84)
Here $`q^2=2\pi H/\varphi _0`$, and the amplitudes $`\stackrel{~}{f}_m`$ can be expressed via $`\stackrel{~}{\mathrm{\Delta }}_1`$ and $`\stackrel{~}{\mathrm{\Delta }}_2`$ from Eqs. (82)-(84):

$`\stackrel{~}{f}_1={\displaystyle \frac{g_+\stackrel{~}{\mathrm{\Delta }}_1g_{12}\stackrel{~}{\mathrm{\Delta }}_2}{\omega +\mathrm{\Omega }_+}}+{\displaystyle \frac{g_{}\stackrel{~}{\mathrm{\Delta }}_1+g_{12}\stackrel{~}{\mathrm{\Delta }}_2}{\omega +\mathrm{\Omega }_{}}},`$ (85)
$`\stackrel{~}{f}_2={\displaystyle \frac{g_{}\stackrel{~}{\mathrm{\Delta }}_2g_{21}\stackrel{~}{\mathrm{\Delta }}_1}{\omega +\mathrm{\Omega }_+}}+{\displaystyle \frac{g_+\stackrel{~}{\mathrm{\Delta }}_2+g_{21}\stackrel{~}{\mathrm{\Delta }}_1}{\omega +\mathrm{\Omega }_{}}},`$ (86)
where

$`2\mathrm{\Omega }_+=\omega _++\gamma _++\mathrm{\Omega }_0,2\mathrm{\Omega }_{}=\omega _++\gamma _+\mathrm{\Omega }_0,`$ (87)
$`\mathrm{\Omega }_0=(\omega _{}^2+\gamma _+^2+2\gamma _{}\omega _{})^{1/2},`$ (88)
$`\gamma _\pm =\gamma _{12}\pm \gamma _{21},\omega _\pm =(D_1\pm D_2)\pi H/\varphi _0,`$ (89)
$`2g_\pm =1\pm {\displaystyle \frac{\omega _{}+\gamma _{}}{\mathrm{\Omega }_0}},g_{12}={\displaystyle \frac{\gamma _{12}}{\mathrm{\Omega }_0}},g_{21}={\displaystyle \frac{\gamma _{21}}{\mathrm{\Omega }_0}}.`$ (90)
Inserting Eqs. (85) and (86) into Eq. (4) we first sum up over $`\omega `$, and express the result in terms of the functions

$`F_\pm =2\pi T{\displaystyle \underset{\omega }{}}{\displaystyle \frac{1}{\omega +\mathrm{\Omega }_\pm }}=\mathrm{ln}{\displaystyle \frac{2\gamma \omega _D}{\pi T}}`$
$`\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\mathrm{\Omega }_\pm }{2\pi T}}\right)+\psi \left({\displaystyle \frac{1}{2}}\right).`$ (91)
Then the self-consistency equations (4) reduce to

$`G_{11}\stackrel{~}{\mathrm{\Delta }}_1+G_{12}\stackrel{~}{\mathrm{\Delta }}_2=0,`$ (92)
$`G_{21}\stackrel{~}{\mathrm{\Delta }}_1+G_{21}\stackrel{~}{\mathrm{\Delta }}_2=0,`$ (93)
where the matrix $`G_{mm^{}}`$ is given by

$`G_{11}=1\lambda _{11}(g_+F_++g_{}F_{})+\lambda _{12}g_{21}(F_+F_{}),`$ (94)
$`G_{12}=\lambda _{11}g_{12}(F_+F_{})\lambda _{12}(g_{}F_++g_+F_{}),`$ (95)
$`G_{22}=1\lambda _{22}(g_+F_++g_{}F_{})+\lambda _{21}g_{12}(F_+F_{}),`$ (96)
$`G_{21}=\lambda _{22}g_{21}(F_+F_{})\lambda _{21}(g_+F_++g_{}F_{}).`$ (97)
Using Eqs. (94)-(97), the final equation $`G_{11}G_{22}G_{12}G_{21}=0`$ for $`H_{c2}`$ can be presented in the form

$`1\lambda _{11}(g_+F_++g_{}F_{})\lambda _{22}(g_+F_++g_{}F_{})`$
$`+w(g_+F_++g_{}F_{})(g_+F_++g_{}F_{})`$
$`+(\lambda _{12}g_{21}+\lambda _{21}g_{12})(F_+F_{})=0,`$ (98)
where $`H_{c2}`$ is the maximum root of Eq. (98). Setting $`H_{c2}=0`$ in Eq. (98), gives the equation for $`T_c`$ with the account of interband scattering for which $`\omega _\pm =\mathrm{\Omega }_{}=0`$, $`\mathrm{\Omega }_+=\gamma _{12}+\gamma _{21}`$, $`g_+=2\gamma _{12}/(\gamma _{12}+\gamma _{21})`$, $`g_{}=2\gamma _{21}/(\gamma _{12}+\gamma _{21})`$. Subtracting Eq. (21) for $`T_{c0}`$ from Eq. (98), the equation for $`T_c`$ can be reduced to

$`U\left({\displaystyle \frac{\gamma _{12}+\gamma _{21}}{2\pi T}}\right)={\displaystyle \frac{(\lambda _0w\mathrm{ln}t)\mathrm{ln}t}{pw\mathrm{ln}t}},`$ (99)
$`p={\displaystyle \frac{\gamma _{12}(\lambda _12\lambda _{21})+\gamma _{21}(\lambda _22\lambda _{12})}{2(\gamma _{12}+\gamma _{21})}},`$ (100)
where $`t=T_c/T_{c0}`$, $`\lambda _1=\lambda _0+\lambda _{}`$, and $`\lambda _2=\lambda _0\lambda _{}`$. The effect of interband scattering on $`T_c`$ suppression was considered beforeimp1 ; imp2 ; lanis1

## Appendix C Matrix M

Inserting Eq. (44) into the gap equation (4) yields

$$\mathrm{\Delta }_m=2\pi T\underset{\omega }{}\underset{m^{}}{}\lambda _{mm^{}}\underset{n=0}{\overset{\mathrm{}}{}}\frac{\phi _n(q_m^{}x)C_n^{(m^{})}}{\omega +(2n+1)\mathrm{\Omega }_m^{}},$$
(101)
where $`C_n^{(m)}=_{\mathrm{}}^{\mathrm{}}\phi _n(q_mx)\mathrm{\Delta }_m(x)𝑑x`$. Now we multiply Eq. (101) by $`\phi _s(q_mx)`$, integrate over x using the orthogonality condition $`\phi _n(q_mx)\phi _n^{}(q_mx)𝑑x=\delta _{nn^{}}`$, and then sum up over $`\omega `$. This gives two sets of linear equations

$`(1\lambda _{11}F_s^{(1)})C_s^{(1)}=\lambda _{12}{\displaystyle \underset{n}{}}V_{sn}F_n^{(2)}C_n^{(2)},`$ (102)
$`(1\lambda _{22}F_s^{(2)})C_s^{(2)}=\lambda _{21}{\displaystyle \underset{n}{}}V_{ns}F_n^{(1)}C_n^{(1)},`$ (103)
where $`V_{sn}`$ and $`F_n^{(m)}`$ are defined as follows

$`V_{sn}={\displaystyle _{\mathrm{}}^{\mathrm{}}}\phi _s(q_1x)\phi _n(q_2x)𝑑x`$ (104)
$`F_n^{(m)}={\displaystyle \underset{\omega >0}{\overset{\omega _D}{}}}{\displaystyle \frac{1}{\omega +(2n+1)\mathrm{\Omega }_m}}`$ (105)
Eq. (105) can be expressed in terms of the $`\psi `$ function as in Eq. (48). Inserting $`C_n^{(2)}`$ from Eq. (103) into Eq. (102) yields the closed set of linear equations for $`C_n^{(1)}`$

$$\underset{n^{}=0}{\overset{\mathrm{}}{}}M_{nn^{}}C_n^{}^{(1)}=0,$$
(106)
where the matrix $`M_{nn^{}}`$ is given by Eq. (47). The solvability condition of Eq. (106) gives the equation for $`H_{c2}`$, which is the maximum eigenvalue of the matrix M.

A general formula for $`V_{nn^{}}`$ is rather cumbersome, so to illustrate general behavior of $`V_{nn^{}}`$ we consider the case of $`n^{}=0`$ for which the integral (104) equals ryzik

$$V_{n0}=\left(\frac{2q_1q_2(2n)!}{2^{2n}(n!)^2(q_1^2+q_2^2)}\right)^{1/2}\left(\frac{q_2^2q_1^2}{q_2^2+q_1^2}\right)^n,$$
(107)
if $`n=2s`$, and $`V_{n0}=0`$ for $`n=2s+1`$, where s is any positive integer. As follows from Eq. (107), $`V_{2s,0}\zeta ^s`$ exponentially decays with s. The matrix elements $`V_{nn^{}}`$ obey the useful sum rule $`_{s=0}^{\mathrm{}}V_{ns}V_{n^{}s}=\delta _{nn^{}}`$, which can be proven by inserting Eq. (104) under the sum:

$$\underset{s=0}{\overset{\mathrm{}}{}}_{\mathrm{}}^{\mathrm{}}𝑑x𝑑x^{}\phi _n(q_1x)\phi _s(q_2x)\phi _n^{}(q_1x^{})\phi _s(q_2x^{})$$
(108)
To sum up over s, we use the identity ll

$$\underset{s=0}{\overset{\mathrm{}}{}}\phi _s(q_2x)\phi _s(q_2x^{})=\delta (xx^{})$$
(109)
After inserting Eq. (109) into Eq. (108), the latter reduces to $`\phi _n(q_1x)\phi _n^{}(q_1x)𝑑x`$, giving the sum rule (51).

Another equation for $`H_{c2}`$ valid for weak interband coupling $`\lambda _{12}\lambda _{21}\lambda _{11}\lambda _{22}`$ and any anisotropy can be obtained as follows. Since off-diagonal terms of M in Eq. (47) are proportional to the small factor $`\lambda _{12}\lambda _{21}`$, the matrix element $`\stackrel{~}{M}_0`$ can be obtained in the second order perturbation theory:

$$\stackrel{~}{M}_0=M_{00}\underset{s=1}{\overset{\mathrm{}}{}}\frac{M_{0,2s}M_{2s,0}}{M_{2s,2s}}$$
(110)
Here $`M_{0,2s}`$ is proportional to the small parameter $`\lambda _{12}\lambda _{21}\zeta `$, so the sum in Eq. (110) is by the factor $`\lambda _{12}\lambda _{21}\zeta ^2`$ smaller than the sum in the right hand side of Eq. (47) for $`M_{00}`$. Thus, the contribution of the off-diagonal terms to $`\stackrel{~}{M}_0`$ can be neglected, and the equation for $`H_{c2}`$ is again $`M_{00}=0`$, that is,

$$1\lambda _{11}F_0^{(1)}=\lambda _{12}\lambda _{21}F_0^{(1)}\underset{s=0}{\overset{\mathrm{}}{}}\frac{V_{0s}^2F_s^{(2)}}{1\lambda _{22}F_s^{(2)}}$$
(111)
For moderate anisotropy $`(\zeta 1)`$, Eq. (111) reduces to Eq. (52), as the convergence of the sum in Eq. (111) is due to the rapidly decreasing factor $`V_{0s}^2\zeta ^s`$. However, for very strong anisotropy $`\zeta 1`$ at low T when the dependence of $`F_s^{(2)}`$ on s is essential, Eq. (111) may be simpler for numerical solution than the full Eq. (106).


format: text

text:
# A Thickness of Stellar Disks of Edge-on Galaxies and Position of Their Truncation Radii

## 1 Introduction

Stellar disks in spiral galaxies consist mainly of old stars and usually include the main part of their stellar content. Being collisionless, they nevertheless may experience the evolution of mass distribution - both along the radial and vertical coordinates due to internal dynamical processes and outer interactions with neighbour systems. What determines the relative disk thickness and disk radial extension, remains an open question. Are these parameters a result of a long lasting evolution or a product of the very early period of the stellar disks formation? The analysis of connection between these parameters and the kinematic properties of galaxies helps to find an answer.

## 2 A thickness of stellar disks of edge-on galaxies

A flatness of stellar disks differ strongly from one galaxy to another. As it was found by some authors, photometric axes ratioes for flat edge-on galaxies correlate with their morphological type: disks of late-type galaxies (Sc–Sd) are on average ”thinner” than that of early-type objects (Karachentsev et al.; de Grijs ; Ma et al., ). At the same time the data of the Revised Flat Galaxy Catalog (Karachentsev et al., , RFGC hereafter) does not reveal any correlation of axes ratioes with the rotation velocity or the luminosity. To illustrate it, in Fig. 1 we compare the $`B`$-band axes ratio $`a/b`$ according to the catalog RFGC with the known HI line width ($`W_{50}`$) and with the absolute magnitude $`M_B`$. Both these parameters were adopted from the LEDA catalog. The relationship between the thickness of the disks and their kinematics, if exists, has to be more sophisticated.

A thickness of a self-gravitating disk depends on its density and the stellar vertical velocity dispersion. It is worth noting that the dynamical heating of a collisionless stellar disk is the one-way process – after the increasing of the stellar velocity dispersion the disk never becomes cooler or thinner again. There are several heating mechanisms known which might increase the vertical scale of the disk, such as a scattering of stars either by massive gas clouds or by other density irregularities, a merging of small satellites, or a gravitational interaction with neighbours. One may easily find the discussion of these mechanisms in the literature. However although the efficiency of any of them should depend on the distance from the galactic center, the observed thickness of disks of edge-on galaxies is usually nearly constant within a wide interval of radial distances (although some exceptions may also exist).

Let’s imagine for instance that we have a disk where the initial velocity dispersion of stars in z-direction $`C_z`$ is very low, say 3 km/s instead of more realistic value 30 km/s. Let also we have neither gas clouds nor satellites or close neighbours. Will a disk stay razor-thin in this case? The answer is definitely NO, and numerical experiments demonstrate it very clearly: the disk always heats up within a short time interval (a couple of rotation periods) as a result of the bending instability which inevitably develops in a thin disk. The thickness ceases to grow when some critical value of $`C_z`$ is reached. This instability determines the minimum possible thickness for collisionless disks. The other heating processes, if involved, may only increase the vertical scalelength of the disk.

As it is known, the bending instability is in some sense an addition to the classical Jeans instability in the disk plane. Indeed, the increasing of radial velocity dispersion $`C_r`$ makes the disk more stable to Jeans perturbations and less stable to the bending ones. Hence a growth of $`C_r`$ leads to the parallel increasing of the vertical dispersion $`C_z`$. As a result, one can expect the proportionality between $`C_r`$ and $`C_z`$ in the final stable state (although, as N-body simulations show, their ratio may depends on the relative masses of spherical and flat components, see Mikhailova et al., 2001). The simplified model of the disk leads to the ratio $`C_z/C_r0.4`$ (Polyachenko & Shukhman ), whereas the observations of real galaxies inferred the values 0.5 - 0.8 with rather large uncertainties. On the other hand, the minimal value of $`C_r`$ is determined by the local gravitational stability of the disk. In general case the critical value of $`C_r`$ is proportional to the surface density over the epicyclic frequency ratio $`\mathrm{\Sigma }_{disk}/\kappa `$. The latter in turn is a function of circular velocity $`V(R)`$, which makes the minimum possible thickness of the disk connected with the kinematic properties of a galaxy.

Let us use the simple proportions written for some fixed dimensionless radius $`R/h`$ of a galaxy where $`h`$ is its radial scalelength:

– the vertical velocity dispersion $`C_zC_r`$;

– the vertical scalelength $`Z_0C_z^2/\mathrm{\Sigma }_{disk}`$;

– the minimal radial velocity dispersion $`C_r\mathrm{\Sigma }_{disk}/\kappa `$;

– the total disk mass of a galaxy $`M_d\mathrm{\Sigma }_{disk}h^2`$;

– the total disk mass of a galaxy $`M_tV^2h`$, and

– the epicyclic frequency $`\kappa V/h`$.

Combining these expressions, one can find that the minimal relative thickness $`Z_0/h`$ is expected to be proportional to the relative mass of the disk: $`Z_0/hM_d/M_t`$. This conclusion agrees well with the results of the 3D $`N`$-body numerical simulations of collisionless initially thin disks embedded into a rigid halo (see Zasov & Morozov , Zasov et al. , Mikhailova et al. , Zasov et al. ). Fig. 2 shows the diagram taken from Mikhailova et al. (), where the mass ratio of spherical and disk components $`M_s/M_d`$ is compared with the vertical over radial scalelengths ratio. In accordance with the expectations, the model galaxies with the most massive halo own the thinniest stable disks. The question is whether the real galaxies obey this rule as well.

To verify it, we composed and analysed two samples of edge- on galaxies with available rotation curves and a surface photometry (see Zasov et al.  for details). All galaxies were taken from RFGC and thought to be disk-dominated objects with a small, if any, bulge contribution to the integrated luminosity. The nearby galaxies, the Virgo cluster members, and galaxies with asymmetric or distorted isophotes were excluded from the consideration.

The first sample includes the objects for which the R-band surface CCD photometry was performed with the 6m telesope at the Special Astrophysical Observatory of RAS (Karachentsev et al. ). We estimated the radial scale length $`h`$ by fitting the photometric major-axis profiles to the analytical expression obtained for a model edge-on exponential disk (van der Kruit&Searle \[1981a\], \[1981b\]). Once the radial scalelength is known, the vertical scale $`Z_0`$ can be found from the observed axes ratio. We will call this sample as the ”BTA-sample” hereafter. The second (”2MASS”) sample consists of galaxies whose infrared images had been taken from the 2MASS survey in the $`K_s`$ band. The radial and the vertical scales for the galaxies were restored from their photometric cuts made with the 2MASS $`K_s`$-band images (avoiding major axes of the disks to reduce the influence of the dust) by Bizyaev & Mitronova  (all references may be found in this paper). The comparison of two samples has revealed that the galaxies in the $`K_s`$ -band look systematically thinner than those in the bluer one (in accordance with the results obtained by de Grjis, ), and hence, the samples should be considered separately.

The starting point for our comparison of the disk thickness with the relative mass of halo in galaxies is that the presence of a dark halo always increases the total mass-to-luminosity ratio $`M/L`$ of a galaxy. For definiteness we will estimate all masses within four disk scalelengths $`4h`$, which encompasses practically all stellar mass and luminosity of the galaxy. In the first approximation the total mass of the galaxy is $`M(4h)=W_H^2h/G`$. The luminosity was estimated from the photometric parameters of the disks found after the data processing. The diagrams $`M/L`$ (that is $`W_H^2h/GL`$) versus $`h/Z_0`$ are shown in Fig. 3a,b.

One can see that, although the dispersion is rather high, the stellar disks indeed tend to be thicker when $`M/L`$ is lower. It is better seen for the 2MASS sample, evidently due to the more reliable photometric models and low internal absorption in the infrared band. The existence of this dependence is what we expect if the bending instability is the major factor determining the disk thickness.

Fig. 4 is the same as Fig. 3b with three curves superposed. These curves are what we expect from the N-body numerical simulations for the initially unstable collisionless galactic disks. The prototype of the curves is the model curve shown in Fig. 2 above. It was recalculated to jump from $`M_s/M_d`$ to the total ratio $`M/L_{Ks}=W_H^2h/GL_{Ks}`$ using the values $`(M/L_{Ks})_{disk}=`$ 1, 2 and 3 solar units for a stellar disk.

One can see that the sample galaxies are actually situated in the domain between $`(M/L_{Ks})_{disk}`$ = 1 and 3 which is quite reasonable for old stellar systems. It indicates that the numerical models of marginally stable disks agree well with observations. In turn, it means that for most of the galaxies the mechanisms of additional disk heating along Z-axis (other than the bending instability) may have a little effect on the vertical disk structure. In this case the approximately constant disk thickness along the radius is just a result of two opposite tendencies: the radial decreasing of both the surface disk density and the vertical velocity dispersion. The influence of these factors on the vertical scale length cancels each other almost exactly within the extended range of radial distances.

Fig. 5 compares the relative disk thickness and the face-on central surface brightness (the latter was reduced to face-on orientation using the model parameters found from the $`R`$\- and the $`K_s`$-band surface brightness distributions). The correlation between these parameters is even more conspicuous than that between $`Z_0/h`$ and $`M/L`$ in Fig 3a,b. Note that the existence of such a relation is nothing unexpected if to take into account that the ”normal” and the low surface brightness galaxies were already shown to exhibit the close relation between $`S_0`$ and the integrated $`M/L`$, which characterizes the dark halo mass fraction (see discussion in MacGaugh & de Block ). The lower is the central surface brightness (and, consequently, the surface density), the higher is its flatness and hence the dark halo mass fraction in a galaxy (the latter dependence between the brightness and dark halo mass was also discussed by MacGaugh & de Block ).

## 3 Where do the old stellar disks end?

The connection between the radial extent of old stellar disks and their kinematic properties is not so evident as in the case of their thickness. It is well known that many late type galaxies experience the significant steepening of the photometric profiles at some radius (truncation radius) $`R_{cut}`$ (usually between 2.5 and 5 of photometric scalelengths $`h`$). As in the case of the disk flattening, the truncation of the disk can be studied most accurately in edge-on galaxies. The nature of the disk truncation is not well known. Here we meet the well known paradox - the less we know about some event, the easier we find its possible interpretations. Some authors consider the disk truncation as a result of a gas density treshold, beyond which the thermal instability, being responsible for the formation of gas cold phase, is absent. The others tie the observed edge of a disk with condition of the local gravitational instability of gaseous layer (see the references in Pohlen et al., Bizyaev & Zasov ). The scarcity of the data makes it difficult to verify the different scenarios. Both approaches mentioned above may explain more or less successfully the radius where the observed present day star formation rate drops down in nearby galaxies. Nevertheless, this radius do not coincide with the truncation radius of an old stellar disk.

It is rather difficult to verify the role of different mechanisms regulating the star formation rate in the past when the presently observed old stellar disk has formed: the velocity dispersion and the density of a gaseous proto-disk, not even saying about the ultraviolet ionizing field, should strongly differ from those we have today. It is possible however to check the fulfilment of the gravitational stability condition for the proto-disks near the truncation radius. For this purpose we used the fact that the critical density of the primeval gaseous disk depends on the same parameters as the observed thickness of the old stellar disk. In other words, we admit that the current stellar vertical velocity dispersion at the truncation radius is equal to (or, in general case, is not higher than) the velocity dispersion of the parent gas. Indeed, the critical surface density of a thin gaseous disk is $`\mathrm{\Sigma }_{crit}C_{gas}\kappa /Q_T\pi G`$ and $`\mathrm{\Sigma }_{disk}C_z^2/\pi GZ_0F(\rho _{halo})`$, where $`Q_T`$ is the Toomre stability parameter, and $`F(\rho _{halo})`$ is the coefficient taking into account the influence of a halo on the disk thickness. To estimate $`\mathrm{\Sigma }_{crit}`$, we used Polyachenko et al.() criterium which gives $`Q_T\sqrt{3}`$ for a flat rotation curve.

We analysed available data for 16 galaxies with the truncation radii and vertical scalelength values taken from van der Kruit & Searle (\[1981a\], \[1981b\],  (NGC 891, 4013, 4217, 4244, 4565, 5907, 5907), and from Pohlen et al. (NGC 2424, ESO187-008, 269-015, 319-026, 321-010, 446-018, 446-044, 528-017 and 581-006), and whose rotation curves were found in the literature. The rotation curves were used to evaluate $`\mathrm{\Sigma }_{disk}`$, $`F(\rho _{halo})`$ and $`\kappa `$ at $`R=R_{cut}`$ for the galaxies chosen. It is well known that the rotation curve decomposing, which we used to find $`\mathrm{\Sigma }_{disk}`$, is ambiguous. Indeed, a comparison of the best fit model with the maximum disk and the minimum disk models (the latter is chosen as the model where ($`M/L_{Ks})_{disk}`$ of the disk is not less than 0.5) shows that the typical uncertainties for the disk mass estimates are about 30%.

The calculations of $`\mathrm{\Sigma }_{disk}`$ and $`\mathrm{\Sigma }_{crit}`$ showed that their ratio decreases slowly toward the periphery for all galaxies considered, and is close to unit (within a factor of about 1.5) at $`R=R_{cut}`$. Their values at $`R=R_{cut}`$ are shown in Fig 6. The filled symbols are related to the galaxies studied by van der Kruit &Searle (\[1981a\], \[1981b\], ), whereas the open ones mark the more distant galaxies taken from Pohlen et al. (). The cross marks the position of our Galaxy.

Curiously, the mean velocity dispersion of stars at the truncation radius for the galaxies we considered appears to be close to the turbulent gas velocity in the present interstellar medium for star forming galaxies: 10- 15 km/s, although it varies from one galaxy to another.

A comparison of the values $`\mathrm{\Sigma }_{disk}`$ and $`\mathrm{\Sigma }_{crit}`$ enables us to conclude that the observational data well agree with the assumption that the observed disk truncation in late type galaxies considered here is connected with the disk kinematics. It supports an idea that the sharp decrease of the surface brightness at $`R>R_{cut}`$ might be a result of rather inefficient star formation beyond this radius due to the gravitational stability of outer parts of a primeval gaseous disk at the time of intense formation of the first generation stars.

## 4 Conclusions

The thickness of stellar disks and their radial extension in most of galaxies are rather conservative parameters reflecting the conditions of the initial formation of galactic disks rather than their long-lasting evolution. In this case the observed dependence of the disk thickness on the relative mass of the dark halo may be naturally explained. The low surface brigtness galaxies are to be among those which possesses both the thinniest disks and the highest dark-to-luminous mass ratioes. Concerning the galactic truncation radii $`R_{cut}`$, although we cannot prove their direct relation to the local Jeans instability of the primeval gas disk at $`R<R_{cut}`$, however this assumption demonstrates a good compatibility with the observed properties of the edge-on galaxies.

This work was supported by the Russian Foundation of Basic Research (grant 01-02-17597).


text:
# The Quantum Query Complexity of 0-1 Knapsack and Associated Claw Problems

## 1 Introduction

The goal of the present paper is to explore the possibility of developing quantum algorithms for different NP-hard problems that are faster than what we get by a direct application of Grover’s search algorithm . This can be seen as part of an active research theme in quantum computing: discovering tasks for which there are quantum algorithms that are significantly faster than classical algorithms (c.f. ). The difficulty in the area is that there are as yet only a few known techniques for designing quantum algorithms. The two central methods are essentially from Shor’s factoring algorithm  and Grover’s search algorithm . Other known algorithms are essentially based on these methods.

In this paper, we consider quantum algorithms for certain NP-hard problems which have a special divide-and-conquer structure that can be exploited to design significantly faster quantum algorithms than what a direct application of Grover’s search algorithm would yield. For example, we give an $`\stackrel{~}{O}(2^{n/3})`$ quantum algorithm for the 0-1 Knapsack problem with $`n`$ variables. This is a consequence of a more general result for 0-1 Integer Linear Programs with $`n`$ variables and $`d`$ inequalities for which we give an $`\stackrel{~}{O}(2^{n/3}n^d)`$ quantum algorithm. For $`d=o(n/\mathrm{log}n)`$ this running time is bounded by $`\stackrel{~}{O}(2^{n(1/3+ϵ)})`$ for every $`ϵ>0`$ and in particular it is significantly better than the $`\stackrel{~}{O}(2^{n/2})`$ time bound if we directly apply Grover’s search algorithm. Our algorithms are, of course, based on Grover’s search algorithm and the powerful method of amplitude amplification .

The next question we address is whether faster quantum algorithms for the above NP-hard problems are possible. To study this, we formulate a new claw-like problem (see  for claw problems) which we term as the *symmetric* claw problem and study its quantum query complexity. The symmetric claw problem essentially captures the structure of 0-1 Knapsack and we study its complexity in the quantum query model . For the symmetric claw problem we are able to prove a lower bound of $`\stackrel{~}{O}(2^{n/4})`$ in the quantum query model using Ambainis’ method . The problem also has an $`\stackrel{~}{O}(2^{n/3})`$ upper bound given by essentially the same quantum algorithm that works for 0-1 Knapsack.

Finally, we consider CNF satisfiability for the case of CNF formulas $`F`$ with number of clauses linearly bounded in the number of variables. Note that we do not make any assumptions on the size of clauses in $`F`$. More precisely, let $`\text{CNF}_c`$ denote CNF formulas $`F`$ such that the number of clauses in $`F`$ is bounded by $`cn`$ for the given constant $`c`$, where $`n`$ is the number of variables. We give a $`2^{(1\alpha )n/2}`$ quantum algorithm for satisfiability of inputs from $`\text{CNF}_c`$, where $`\alpha `$ is a constant depending on $`c`$.

### 1.1 Preliminaries

The set $`\{0,1\}^n`$ denotes the set of binary strings of length $`n`$, and for a positive integer $`N`$ we denote the set $`\{1,2,\mathrm{},N\}`$ by $`[N]`$.

We design new quantum algorithms for the following NP-hard problems.

0-1 Knapsack

Input: A list of positive integers $`c_1,c_2,\mathrm{},c_n`$ and a positive integer $`K`$.

Problem: Is there a subset $`S[n]`$ such that $`_{iS}c_i=K`$.

Actually, we consider is 0-1 Integer Linear Programs in general.

0-1 ILP

Input: Integers $`a_{ij}`$ and $`b_i`$, $`1jn`$ and $`1id`$.

Problem: Is the following set of linear inequalities feasible?

$$\underset{j=1}{\overset{n}{}}a_{ij}x_jb_i,x_j\{0,1\},1jn1id.$$

We use standard definitions from quantum computing from . In particular, we use the quantum query model (c.f. ). In this model the input is a function $`f:[N][M]`$ and the values of $`f`$ is accessed by oracle queries. The complexity of computing some property of $`f`$ is measured as the number of oracle queries. A quantum computation with $`T`$ queries can be seen as a sequence of unitary transformations:

$$U_0OU_1\mathrm{}OU_T,$$

where the unitary transform $`O`$ implements the oracle access to $`f`$ and the $`U_i`$ are arbitrary unitary transformations which do not depend on the input $`f`$.

A central idea from quantum computing, which is a generalized form of Grover’s search algorithm that we use throughout the paper is amplitude amplification : Essentially, if we have a quantum algorithm $`𝒜`$ with success probability $`p`$ then the success probability can be amplified to a constant by $`O(\sqrt{1/p})`$ calls to $`𝒜`$ and $`𝒜^1`$.

A *range tree* (see  for details) is a data structure for storing a set $`S`$ of $`N`$ elements $`\{x_1,x_2,\mathrm{},x_N\}`$, each of which is a $`d`$-tuple of integers $`x_i=(x_{i1},x_{i2},\mathrm{},x_{id})`$. The data structure can be be built in time $`N(log^{d1}N)`$ and needs space $`N(log^{d1}N)`$. The specific property of interest that we require is that we can process *range queries* on a range tree in time $`\mathrm{log}^dN`$ to retrieve one element, if it exists, of the set $`S`$ that satisfies the query bounds. More precisely, given a $`d`$-tuple of real numbers $`(a_1,a_2,\mathrm{},a_d)`$, in time $`O(\mathrm{log}^dN)`$ we can search for an element $`x_i=(x_{i1},x_{i2},\mathrm{},x_{id})`$ in the range tree such that $`x_{ij}a_j`$, $`1id`$.

## 2 New Quantum Algorithms for 0-1 Integer Linear Programs

###### Theorem 1

There is an $`\stackrel{~}{O}(2^{n/3}n^d)`$ quantum algorithm with constant error probability that solves 0-1 integer linear programs with $`n`$ variables and $`d`$ inequalities.

Proof.  Let the input instance be

$$\underset{j=1}{\overset{n}{}}a_{ij}x_jb_i,1id,x_j\{0,1\},1jn.$$

The goal is to find a feasible solution. We give a stepwise description of the quantum algorithm.

1. Partition the variables into two sets $`A=\{x_1,\mathrm{},x_{n/3}\}`$ and $`B=\{x_{n/3+1},\mathrm{},x_n\}`$, where note that $`|A|=n/3`$ and $`|B|=2n/3`$.
2. For each of the $`2^{n/3}`$ 0-1 assignments $`I`$ to the variables in $`A`$ define the $`d`$-tuple $`y_I=(y_I(1).\mathrm{},y_I(d))`$ of integers, where

$$y_I(i)=\underset{j=1}{\overset{n/3}{}}a_{ij}I_j,$$

where $`I_j`$ is the value of variable $`x_j`$ in assignment $`I`$.
3. Let $`N=2^{n/3}`$. The set $`X=\{y_II`$ is a 0-1 assignment to $`A\}`$ is of size $`N`$. In time $`\stackrel{~}{O}(N\mathrm{log}^{d1}N)`$ build a range tree (c.f. \[9, Theorem 2.11\]) to store the set $`X`$. The range tree has size $`\stackrel{~}{O}(N\mathrm{log}^{d1}N)`$ and range queries can be processed in time $`\mathrm{log}^dN`$. This completes the preprocessing phase.
4. Using the Hadamard transform prepare the uniform superposition over the set $`T`$ of 0-1 assignments to the $`2n/3`$ variables in $`B`$:

$$|\psi =\frac{1}{2^{n/3}}\underset{aT}{}|a.$$
5. Define a unitary transform $`U`$ as the standard reversible implementation of the following classical subroutine $`f`$ given $`uT`$ as input:

* Compute $`z_j=_{j=n/3+1}^na_{ij}u_j`$, $`1id`$.
    * Let $`\widehat{z_i}=b_iz_i`$ for $`1id`$, giving a $`d`$-tuple of integers $`(\widehat{z_1},\mathrm{},\widehat{z_d})`$.
    * Search in the range tree for a $`y_I`$ such that $`y_I(1)\widehat{z_1},\mathrm{},y_I(d)\widehat{z_d}`$. Notice that this is a range query and it can be processed in time $`\mathrm{log}^dN`$ in the range tree date structure \[9, Theorem 2.11\].
    * If such a tuple $`y_I`$ is in $`X`$ then $`f(u)=1`$. Otherwise, $`f(u)=0`$.
6. With the initial state as the uniform superposition $`|\psi `$ and using the unitary transform $`U`$, apply Grover’s search algorithm to search for $`uT`$ such that $`f(u)=1`$.

It follows from the well-known analysis of Grover’s algorithm and amplitude amplification  that the above algorithm has running time $`\stackrel{~}{O}(2^{n/3}\mathrm{log}^dN)=\stackrel{~}{O}(2^{n/3}n^d)`$, with constant success probability.

###### Corollary 2

There is an $`\stackrel{~}{O}(2^{n/3})`$ time quantum algorithms with constant success probability for 0-1 Knapsack. Additionally, any NP-hard optimization problem that takes the form of a 0-1 Integer Linear Program with constant number of constraints has an $`\stackrel{~}{O}(2^{n/3})`$ time quantum algorithms with constant success probability.

Proof.  As an instance of 0-1 Knapsack consists of one equation, it can be expressed as a 0-1 integer linear program with two inequalities. Hence the algorithm of Theorem 1 yields the claimed quantum algorithm.

Consider any optimization problem that has a linear optimality function and the constraints can be expressed as a 0-1 integer linear program with $`d`$ constraints. Using binary search we can reduce the optimization problem to feasibility of a 0-1 integer linear program with $`d+1`$ constraints, which can be solved using the algorithm of Theorem 1 with running time $`\stackrel{~}{O}(2^{n/3})`$ for constant $`d`$. Thus, the corresponding optimization problem can also be solved in time $`\stackrel{~}{O}(2^{n/3})`$ with constant success probability.

* Interestingly, if the 0-1 ILP is of the following form:

minimize/maximize $`{\displaystyle \underset{i}{}}c_ix_i`$
s.t. $`{\displaystyle \underset{j=1}{\overset{n}{}}}a_{ij}x_j=b_i,`$ $`1id,x_j\{0,1\},1jn,`$
then it is easy to see that essentially the same quantum algorithm presented in Theorem 1 solves this optimization problem in time $`\stackrel{~}{O}(2^{n/3})`$, independent of the number of equations $`d`$. The reason is that we do not have to maintain a range tree data structure as we do not have to process range queries. Just a sorted list would suffice as we only need to make exact queries in this case.

The above kind of 0-1 ILP is referred to as the 0-1 Group Problem in Nemhauser and Wolsey’s book . An interesting instance of this problem is the following NP-hard problem: given a CNF formula $`F`$ (no restrictions on clause size), search for a satisfying assignment that satisfies *exactly* one variable in each clause. This is clearly an instance of the above 0-1 ILP. Thus we have the following corollary.

###### Corollary 3

There is an $`\stackrel{~}{O}(2^{n/3})`$ quantum algorithm that takes a CNF formula $`F`$ with $`n`$ variables as input (no restrictions on clause size or number of clauses) and searches, with constant success probability, for a satisfying assignment that satisfies exactly one variable in each clause.

## 3 The Symmetric Claw Problem

In order to study how far we can exploit this idea of dividing the input, we examine in this section a black-box version of the 0-1 Knapsack problem which we term as the Symmetric Claw problem.

For a pair of strings $`x,y\{0,1\}^n`$ let $`x|_y`$ denote the substring of $`x`$ obtained by projecting it on the positions where $`y`$ is 1. Likewise, let $`x|_{\overline{y}}`$ denote the substring of $`x`$ obtained by projecting it on the positions where $`y`$ is 0.

###### Definition 4

The input to the *symmetric claw problem* is a function $`P:\{0,1\}^n\times \{0,1\}^n\{0,1\}^m\times \{0,1\}^m`$ where we write $`P(x,y)=(P_1(x,y),P_2(x,y))`$ for every $`x,y\{0,1\}^n`$. Additionally, the input function $`P`$ is such that it fulfills the following promise symmetry conditions:

1. For $`x\{0,1\}^n`$, if $`P_1(x,y)=P_2(x,y)`$ for some $`y\{0,1\}^n`$ then $`P_1(x,y)=P_2(x,y)`$ for every $`y\{0,1\}^n`$.
2. For any $`x,x^{},y^{}\{0,1\}^n`$ if $`x|_y=x^{}|_y`$ then $`P_1(x,y)=P_1(x^{},y)`$. Likewise, if $`x|_{\overline{y}}=x^{}|_{\overline{y}}`$ then $`P_2(x,y)=P_2(x^{},y)`$.

Given such an input $`P`$ the symmetric claw problem is to find an $`x\{0,1\}^n`$, if it exists, such that $`P_1(x,y)=P_2(x,y)`$ for some $`y\{0,1\}^n`$, by querying $`P`$ on $`(x,y)\{0,1\}^n\times \{0,1\}^n`$.

Notice that the 0-1 Knapsack problem is an instance of the symmetric claw problem.

* Consider a simpler version of the 0-1 Knapsack problem, where we seek a 0-1 assignment to variables $`x_1,\mathrm{},x_n`$ such that $`_{i=1}^na_ix_i=0`$, given integers $`a_1,\mathrm{},a_n`$ as input.

For each $`y\{0,1\}^n`$ define a partition of the variable set $`A=\{i[n]y_i=1\}`$ and $`B=\{i[n]y_i=0\}`$. Define $`P_1(x,y)=_{iA}a_ix_i`$ and $`P_2(x,y)=_{iB}a_ix_i`$.

It is easy to check that $`P_1`$ and $`P_2`$ fulfill conditions of Definition 4, making this version of 0-1 Knapsack an instance of the symmetric claw problem.

We first establish the upper bound result for symmetric claw problem. It is the same algorithm we described in the previous section.

###### Theorem 5

There is an $`\stackrel{~}{O}(2^{n/3})`$ quantum algorithm with constant error probability for the symmetric claw problem in the quantum query model.

Proof.  The algorithm is along exactly the same lines as the algorithm in Theorem 1 for 0-1 ILP. In fact it is an easier version of that algorithm as we do not require range trees here.

We fix $`y\{0,1\}^n`$ to be such that $`y_i=1`$ for $`1in/3`$ and $`y_i=0`$ for $`n/3+1in`$.

Let $`X_A=\{x\{0,1\}^nx_i=0`$ for $`n/3+1in\}`$ and $`X_B=\{x\{0,1\}^nx_i=0`$ for $`1in/3\}`$.

Compute and sort the list $`\{P_1(x,y)xX_A\}`$ of size $`2^{n/3}`$ using a classical sorting algorithm in time $`\stackrel{~}{O}(2^{n/3})`$. This is the preprocessing phase.

Prepare a superposition of $`xX_B`$ and search for $`P_2(x,y)`$ in the sorted list using a classical binary search implemented by a unitary transform $`U`$ defined as: For $`xX_B`$, $`U|x=1`$ if $`P_2(x,y)`$ occurs in the sorted list and $`U|x=0`$ otherwise.

Now, using Grover’s search we can find $`xX_B`$ such that $`P_1(x^{},y)=P_2(x,y)`$ for some $`x^{}X_A`$. At this point, using the properties of $`P_1`$ and $`P_2`$ we can put together a complete 0-1 assignment $`x^{\prime \prime }`$ from $`x`$ and $`x^{}`$ such that $`P_1(x^{\prime \prime },y)=P_2(x^{\prime \prime },y)`$. We can obtain $`x^{\prime \prime }`$ by concatenating the first $`n/3`$ bits from $`x^{}`$ and the last $`2n/3`$ bits from $`x`$.

We now turn to the question of lower bounds. Using the technique of Ambainis  we establish a lower bound of $`\mathrm{\Omega }(2^{n/4})`$ for the quantum query complexity of the symmetric claw problem. We first state the result in a form that we need for our setting.

###### Theorem 6

 Let $``$ be the set of input functions $`f`$, where $`f:[M][M]`$, and let $`\varphi :\times `$ be the function which we wish to compute. Let $`X,Y`$ be two subsets of $``$ such that $`\varphi (f)\varphi (g)`$ for all $`fX`$ and $`gY`$. Let $`RX\times Y`$ such that:

1. For every $`fX`$ there are at least $`m`$ different $`gY`$ such that $`(f,g)R`$.
2. For every $`gY`$ there are at least $`m^{}`$ different $`fX`$ such that $`(f,g)R`$.
3. For every $`fX`$ and $`x[N]`$, there are at most $`l`$ different $`gY`$ such that $`(f,g)R`$ such that $`f(x)g(x)`$.
4. For every $`gY`$ and $`x[N]`$, there are at most $`l^{}`$ different $`fX`$ such that $`(f,g)R`$ such that $`f(x)g(x)`$.

Any quantum algorithm that evaluates $`\varphi `$ with constant success probability must make $`\mathrm{\Omega }(\sqrt{\frac{mm^{}}{ll^{}}})`$ queries to the input function in the quantum query model.

Now we are ready to prove our lower bound result of this section.

###### Theorem 7

Any quantum algorithm that solves the symmetric claw problem with constant success probability needs $`\mathrm{\Omega }(2^{n/4})`$ queries in the quantum query model.

Proof.  We will prove the claimed lower bound by applying Theorem 6. The inputs consist of the class $``$ of functions $`P:\{0,1\}^n\times \{0,1\}^n\{0,1\}^n\times \{0,1\}^n`$, where we write $`P(x,y)=(P_1(x,y),P_2(x,y))`$, and every $`P`$ in $``$ fulfills the conditions of the Definition 4. The goal of a quantum algorithm $`\varphi `$ for the symmetric claw problem is to find $`x\{0,1\}^n`$ such that $`P_1(x,y)=P_2(x,y)`$ for some (and therefore, by the promise of Definition 4, every $`y`$). Let NO denote the subset of “no” instances of $``$ and YES the subset of “yes” instances. Define $`X\text{NO}`$ as:

$`X=\{PP_1(x,y)P_2(x,y)x,y\{0,1\}^n,P_1(x,y)P_1(x^{},y)`$ if $`x|_yx^{}|_y`$ and $`P_2(x,y)P_2(x^{},y)`$ if $`x|_{\overline{y}}x^{}|_{\overline{y}}\}`$.

We will define the subset $`Y`$ of YES by first defining a relation $`RX\times \text{YES}`$ and then setting $`Y=\{P^{}\text{YES}(P,P^{})R`$ for some $`PX\}`$. We now define the relation $`R`$.

Given $`PX`$, for every $`x\{0,1\}^n`$ we define a distinct $`P^{}\text{YES}`$ and include $`(P,P^{})`$ in $`R`$. Given $`PX`$ and $`x\{0,1\}^n`$, $`P^{}`$ is obtained from $`P`$ as follows:

* If $`_{i=1}^ny_in/2`$ then define $`P_1^{}(x,y)=P_1(x,y)`$, and $`P_2^{}(x,y)=P_1(x,y)`$. Furthermore, for every $`x^{}\{0.1\}^n`$ such that $`x^{}|_{\overline{y}}=x|_{\overline{y}}`$ define $`P_2^{}(x^{},y)=P_1(x,y)`$. For such $`x^{}`$, notice that by definition of $`X`$ we have $`P_1^{}(x^{},y)P_1(x,y)`$.
* If $`_{i=1}^ny_i<n/2`$ then define $`P_2^{}(x,y)=P_2(x,y)`$, and $`P_1^{}(x,y)=P_2(x,y)`$. Furthermore, for every $`x^{}\{0.1\}^n`$ such that $`x^{}|_y=x|_y`$ define $`P_1^{}(x^{},y)=P_2(x,y)`$. For such $`x^{}`$, notice that by definition of $`X`$ we have $`P_2^{}(x^{},y)P_2(x,y)`$.
* For all other pairs $`(u,v)\{0,1\}^n\times \{0,1\}^n`$ define $`P_1^{}(u,v)=P_1(u,v)`$ and $`P_2^{}(u,v)=P_2(u,v)`$.

Let $`P^{}`$ be defined as above for a given $`PX`$ and $`x\{0,1\}^n`$. Notice that by construction $`P^{}`$ is in YES and, furthermore, $`x`$ is the unique element of $`\{0,1\}^n`$ such that $`P_1^{}(x,y)=P_2^{}(x,y)`$ for $`y\{0,1\}^n`$. Thus we define $`2^n`$ distinct elements $`P^{}\text{YES}`$, one for every $`x\{0,1\}^n`$, such that $`(P,P^{})R`$. Therefore, in the terminology of Theorem 6, $`m=2^n`$. For the remainder of the proof, we denote by $`P_x^{}`$ the element $`P^{}\text{YES}`$ defined as above from an $`x\{0,1\}^n`$ and a given $`PX`$.

We next show that for any $`PX`$ and $`(x,y)\{0,1\}^n\times \{0,1\}^n`$ there are at most $`2^{n/2}`$ elements $`P^{}Y`$ such that $`P(x,y)P^{}(x,y)`$ and $`(P,P^{})R`$. In terms of Theorem 6 this implies $`l2^{n/2}`$. We argue this in two cases:

1. Suppose $`_{i=1}^ny_in/2`$. Then notice from the definition of $`R`$ that $`P_z^{}`$ differs from $`P`$ at $`(x,y)`$ exactly when $`x|_y=z|_y`$. Thus, there are at most $`2^{n/2}`$ such elements $`P_z^{}`$.
2. Suppose $`_{i=1}^ny_i<n/2`$. Again, notice from definition of $`R`$ that $`P_z^{}`$ differs from $`P`$ at $`(x,y)`$ exactly when $`x|_{\overline{y}}=z|_{\overline{y}}`$. Thus, again, there can be at most $`2^{n/2}`$ such elements $`P_z^{}`$.

Finally, notice that $`l^{}m^{}`$ always holds.

We can now apply the lower bound result Theorem 6. Since

$$\frac{mm^{}}{ll^{}}\frac{m}{l}\frac{2^n}{2^{n/2}},$$

The lower bound of $`\mathrm{\Omega }(\sqrt{\frac{mm^{}}{ll^{}}})`$ given by Theorem 6 yields $`\mathrm{\Omega }(2^{n/4})`$. This completes the proof.

## 4 Quantum Algorithm for CNF-SAT

In this section we use ideas from Section 2 to give a simple quantum search algorithm for CNF-SAT when the number of clauses is linearly bounded in the number of variables. More precisely, for any constant $`c>0`$ we define the following class of CNF formulas: $`\text{CNF}_c=\{FF`$ is in conjunctive normal form with $`n`$ variables and $`F=C_1C_2\mathrm{}C_m`$ where $`mcn\}`$.

###### Theorem 8

For any constant $`c>0`$ there is an $`\stackrel{~}{O}(2^{(1\alpha )n/2})`$ quantum algorithm for satisfiability for inputs $`F\text{CNF}_c`$, where $`\alpha `$ is a constant such that $`\alpha <1/6`$ and $`H(\alpha )\frac{1}{4c}`$. Here, $`H(\alpha )`$ is the entropy function defined as $`\alpha \mathrm{log}_e\alpha (1\alpha )\mathrm{log}_e(1\alpha )`$.

Proof.  Let $`\alpha <1/6`$ be a constant to be fixed later, as claimed in the statement. Let the input instance from $`\text{CNF}_c`$ be $`F=C_1C_2\mathrm{}C_m`$ from $`\text{CNF}_c`$ with variables $`x_1,x_2,\mathrm{},x_n`$. Let $`\frac{1}{\alpha }=k`$, where we assume for simplicity of analysis that $`k`$ is an integer. Partition the variable set $`\{x_1,x_2,\mathrm{},x_n\}`$ into $`k`$ equal size sets $`A_1,A_2,\mathrm{},A_k`$, where $`|A_i|=\alpha n`$, $`1ik`$. Suppose $`F`$ is satisfiable and $`a^{}`$ is some fixed satisfying assignment for $`F`$. Denote by $`\overline{A_i}`$ the set of variables $`\{x_1,\mathrm{},x_n\}A_i`$ for each $`i`$. Denote by $`b_i^{}`$ the partial assignment given by $`a^{}`$ when restricted to variables in $`A_i`$. Similarly, let $`c_i^{}`$ denote the partial assignment given by $`a^{}`$ when restricted to variables in $`\overline{A_i}`$.

Claim. There is an $`i:1ik`$ such that the partial assignment $`c_i^{}`$ satisfies at least $`(1\alpha )m`$ clauses of $`F`$.

Suppose the claim is false. Then $`\overline{A_i}`$ satisfies at most $`(1\alpha )m`$ clauses of $`F`$ for each $`i`$. Since $`a^{}`$ is a satisfying assignment, for each $`i`$ there is a set $`S_i`$ of *more than* $`\alpha m`$ clauses satisfied by the partial assignment $`b_i^{}`$ such that $`S_iS_j=\mathrm{}`$ for $`ij`$. This is impossible since there are only $`m`$ clauses in $`F`$. The claim follows.

We now give a step-wise description of the algorithm. The algorithm is a loop with index $`i:1ik`$, where in the $`i`$th iteration it considers the partition of the variable set into $`A_i`$ and $`\overline{A_i}`$. The algorithm succeeds when it considers a partition $`A_i`$ and $`\overline{A_i}`$ such that $`c_i^{}`$ satisfies at least $`(1\alpha )m`$ clauses of the input $`F`$ (such an index exists by the claim).

1. Partition the variable set into $`A_i`$ and $`\overline{A_i}`$.
2. For $`1l\alpha m`$ build a sorted table $`T_l`$ consisting of the following set of pairs $`(u,b)`$: $`u\{0,1\}^m`$ is am $`m`$-bit string with exactly $`l`$ 1’s and $`b`$ is a truth assignment to variables in $`A_i`$ such that $`b`$ satisfies each clause $`C_j`$ for which $`u_j=1`$, where $`1jm`$. All other pairs $`(u,b)`$ are discarded. Then, using a classical sorting algorithm and treating $`u`$ in $`(u,b)`$ as an $`m`$-bit integer key, sort the pairs in $`T_l`$ in increasing order. This entire step can be done in time $`\stackrel{~}{O}(\left(\genfrac{}{}{0pt}{}{m}{\alpha n}\right)2^{\alpha n})`$.
3. Using the Hadamard transform prepare a uniform superposition of the truth assignments $`v`$ to variables in $`\overline{A_i}`$. Let $`S`$ be the set of all such truth assignments. Notice that $`|S|=2^{(1\alpha )n}`$. The superposition is

$$|\psi =\sqrt{\frac{1}{2^{(1\alpha )n}}}\underset{vS}{}|v.$$
4. Define a unitary transform $`U`$ as the standard reversible implementation of the following classical subroutine $`f`$:

* Compute the vector $`u\{0,1\}^m`$ where $`u_j=1`$ if and only if $`C_j`$ is not satisfied by the partial assignment $`v`$.
    * If $`u`$ has more than $`\alpha m`$ 1’s then $`f(v)`$ returns value $`0`$ and stop.
    * Otherwise, if $`u`$ has $`l`$ 1’s in it, $`1l\alpha m`$, then do a binary search for $`u`$ in the table $`T_l`$. If $`(u,b)`$ is found in $`T_l`$ for some $`b`$ then $`f(v)=1`$ else $`f(v)=0`$. Notice that if $`f(v)=1`$ then $`v`$ and $`b`$ together give a satisfying assignment for $`F`$.
5. Now, starting with initial state as $`|\psi `$ and using the unitary transform $`U`$ to implement $`f`$, apply Grover search to find $`vS`$ such that $`f(v)=1`$.
6. Output the satisfying assignment $`(b,v)`$.

To argue correctness it suffices to notice by our earlier claim that for some $`i`$, $`A_i`$ and $`\overline{A_i}`$ is a partition such that in Step 3, $`f(v)=1`$ for at least the partial assignment $`v=c_i^{}`$ obtained from $`a^{}`$.

Thus, if $`F`$ is satisfiable the algorithm will find a satisfying assignment with constant success probability in time $`\stackrel{~}{O}(\left(\genfrac{}{}{0pt}{}{m}{\alpha m}\right)2^{\alpha n})+\stackrel{~}{O}(2^{(1\alpha )n/2})`$, where the first term is the preprocessing time in Step 2, and the second term is the time for the Grover search.

Since $`mcn`$, an easy analysis yields $`\left(\genfrac{}{}{0pt}{}{m}{\alpha m}\right)2^{H(\alpha )cn}`$ for large $`n`$. Thus, if we choose $`\alpha <1/6`$ such that $`(1\alpha )/2H(\alpha )c+\alpha `$, then the overall time taken is $`\stackrel{~}{O}(2^{(1\alpha )n/2})`$. It suffices to choose $`\alpha `$ such that $`H(\alpha )<\frac{1}{4c}`$.

## 5 Discussion

The question that arises is whether we can close the gap between the upper and lower bounds for the symmetric claw problem. Also, for the 0-1 ILP problem with $`n`$ variables and $`d`$ inequalities we do not have a better upper bound than $`\stackrel{~}{O}(2^{n/3})n^d`$. To study this issue in the quantum query model we define the simultaneous claw problem.

###### Definition 9

Let $`(f_i,g_i),1id`$ be functions where $`f_i:[N]X`$ and $`g_i:[N]X`$, where $`X`$ is some set. The *simultaneous claw* problem is to find an $`x[N]`$ such that $`f_i(x)=g_i(x)`$ for each $`i=1,2,\mathrm{},d`$.

We can also define a symmetric version of the simultaneous claw problem similar to Definition 4, but for simplicity we focus on this definition. It is known  that for $`d=1`$ there is an $`\stackrel{~}{O}(N^{3/4})`$ upper bound and the recent new techniques of Aaronson  followed by Shi’s sharpened results  imply an $`\mathrm{\Omega }(N^{2/3})`$ lower bound for the problem.

We show here that the simultaneous claw problem has a quantum upper bound of $`\stackrel{~}{O}(N^{3/4}\mathrm{log}^dN)`$. On the other hand, we are not able to strengthen the lower bound of $`\mathrm{\Omega }(N^{2/3})`$ which already holds for a single claw pair $`(f,g)`$. The following theorem is easy to prove.

###### Theorem 10

There is a quantum algorithm that takes as input a collection $`(f_i,g_i),1id`$ of functions where $`f_i:[N]X`$ and $`g_i:[N]X`$, where $`X`$ is some set, which makes $`\stackrel{~}{O}(N^{3/4}\mathrm{log}^dN)`$ queries in the quantum query model and outputs a simultaneous claw $`x[N]`$, if it exists, such that $`f_i(x)=g_i(x)`$ for each $`i=1,2,\mathrm{},d`$.

We recall the $`r`$-to-1 collision problem: given as input a function $`f:[N]X`$ with the promise that $`f`$ is either 1-1 or $`r`$-to-1, the problem is to find a pair $`xy[N]`$ such that $`f(x)=f(y)`$. The recent results of Aaronson  followed by Shi  imply an $`\mathrm{\Omega }(N^{1/3})`$ lower bound for the 2-to-1 collision problem.

Analogous to the simultaneous claw problem we can define the simultaneous 2-to-1 collision problem which is given as input functions $`f_i:[N]X`$, $`1id`$, with the promise that either all the $`f_i`$ are 1-1 or all are 2-to-1, find a pair $`xy[N]`$ if it exists such that $`f_i(x)=f_i(y)`$ for each $`i`$. The $`\mathrm{\Omega }(N^{1/3})`$ lower bound of Shi  clearly holds but we have nothing better than that. On the other hand, we have a straightforward algorithm that gives an upper bound for the problem.

###### Theorem 11

There is a quantum algorithm that takes as input a collection $`f_i,1id`$ of functions $`f_i:[N]X`$ where $`X`$ is some set, makes $`\stackrel{~}{O}(N^{3/4}\mathrm{log}^dN)`$ queries in the quantum query model and outputs a simultaneous collision $`xy[N]`$ such that $`f_i(x)=f_i(y)`$ for each $`i=1,2,\mathrm{},d`$, if it exists.