Orbital upper critical field and its anisotropy of clean one- and two-band superconductors

The seminal work of Helfand and Werthamer (HW) [1] on the temperature dependence of the upper critical field H_c2(T) is routinely applied to analyze data on new superconductors despite the fact that HW considered the isotropic s-wave case whereas in the majority of new materials both Fermi surfaces and order parameters are strongly anisotropic. The problem of H_c2(T) has been studied theoretically for anisotropic situations as well: for layered systems [5], for a few cases of hexagonal anisotropy of the Fermi surface and of the order parameter [6], for the two-band MgB₂ [2, 3, 7–9], for d-wave cuprates, see e.g., [10] and references therein, and in a comprehensive work in [11], to name a few.

A ubiquitous feature of H_c2 in many new superconductors is the temperature dependent anisotropy parameter γ = H_c2,ab/H_c2,c (for uniaxial materials), the property absent in conventional one-band isotropic s-wave materials. For example, γ(T) of MgB₂ decreases on warming [12], whereas for many Fe-based materials γ(T) increases with increasing T [13–15]. Up to date, the T dependence of the anisotropy parameter γ is considered by many to be caused by a multi-band character of the materials in question with the commonly given reference to the example of MgB₂. To the best of our knowledge, no explanation was offered for the 'unusual' γ(T) of pnictides.

In this work, we develop a method to evaluate both H_c2,c(T) and γ(T) that can be applied with minor modifications to various situations of different order parameter symmetries and Fermi surfaces, two bands included. Having in mind possible applications for data analysis, we provide only the necessary minimum of analytic developments resorting to numerical methods as soon as possible.

The upper critical field is affected by many factors: magnetic structures that may coexist or interfere with superconductivity, the paramagnetic limit, scattering, etc. In this work, we have in mind to establish a qualitative picture of how the general features of anisotropic Fermi surfaces and order parameters affect H_c2 and its anisotropy, in particular.

Since the paramagnetic effects and the possibility of Fulde–Ferrel–Larkin–Ovchinnikov phase are not included in our scheme (one can find a comprehensive discussion of these questions in [8]), applications to materials with high H_c2(0) should be carried out with care; our results can prove useful for interpretation of data at temperatures where H_c2(T) does not exceed the paramagnetic limit.

As far as applications to two-band materials are concerned, we note that our formalism applies only to superconductors with a single critical temperature T_c. More exotic possibilities of two component condensates with two distinct T_cs are out of the scope of this paper; these are considered on general group-theoretical symmetry grounds, e.g., in [16, 17].

Fine details of Fermi surfaces are unlikely to strongly influence H_c2(T) because, as is well known and shown explicitly below, only the integrals over the whole Fermi surface enter equations for H_c2(T). Circumstantial evidence for a weak connection between fine particularities of the Fermi surface and H_c2(T) is provided by the very fact that the HW isotropic model works so well for many one-band s-wave materials, although their Fermi surfaces hardly resemble a sphere, take for example Nb. Another example is given by the calculations of [2, 3] for MgB₂ based on different model Fermi surfaces, but giving similar results reasonably close to the measured H_c2(T). We, therefore, model actual Fermi surfaces of uniaxial materials of interest here by rotational ellipsoids (spheroids) choosing them so as to have averaged squared Fermi velocities equal to the measured values or to those calculated for realistic Fermi surfaces. This idea, in fact, has been employed by Miranovic et al for MgB₂ [2]. Also, we tested our method on a rotational hyperboloid as an example of open Fermi surfaces, appendix D. This work is still in progress. We intend to study more about effects of open Fermi surfaces (or two-band combinations of closed and open surfaces) on H_c2(T).

We focus in this work on the clean limit for two major reasons. First, commonly after discovery of a new superconductor, an effort is made to obtain as clean single crystals as possible since those provide a better chance to study the underlying physics. Second, almost as a rule, new materials are multi-band so that the characterization of scattering leads to a multitude of scattering parameters which cannot be easily controlled or separately measured. In addition, in general, scattering suppresses the anisotropy of H_c2, the central quantity of interest in this work. We refer readers interested in effects of scattering to a number of successful dirty-limit microscopic models, e.g. to [7] and references therein.

In the next section, we begin with the general discussion of the H_c2 problem for arbitrary Fermi surfaces and order parameters. We show that in the isotopic limit our approach is reduced to that of HW. The basic HW approach is then applied to anisotropic situations. The derivation involves rescaling of the coordinates and, therefore, necessitates co- and contravariant vector representations. In our view, disregarding this necessity may lead to incorrect conclusions. Another formal feature of our approach should be mentioned: we avoid the minimization relative to the actual coordinate dependence of the order parameter in the mixed state often employed to find H_c2(T) [2, 3, 9, 10, 18]. In this sense, our method is close to the original HW work that stresses that H_c2(T) is actually an eigenvalue problem.

Next, we formulate the problem for two s-wave bands and show that along with H_c2(T) and its anisotropy one can find the ratio of order parameters on two bands, a quantity that up to date has been studied only in zero field. We then formulate the problem for two bands with order parameters of different symmetries. To show how the method actually works, we consider in detail one or two bands with Fermi surfaces as rotational ellipsoids. The method is demonstrated on the well-studied example of MgB₂.

The anisotropy parameter γ is shown to depend on temperature even for the one-band case for other than s-wave order parameters. This dependence is weak in the d-wave materials with closed Fermi surfaces, is stronger for open Fermi shapes and is stronger yet for order parameters of the form Δ = Δ₀(1 + η cos k_za), one of the candidates for Fe-based materials [4]. Moreover, γ(T) increases or decreases on warming depending on the sign of the coefficient η, in other words, on whether Δ is maximum or minimum at k_z = 0. These results challenge the common belief that temperature dependence of γ is always related to the multi-band topology of Fermi surfaces.

For two bands, after checking the method on MgB₂, we focus on situations with dominant inter-band coupling, which is relevant for Fe-based materials. While in most cases we have considered, H_c2,c(T) is qualitatively similar to the HW curve, the anisotropy parameter γ(T) may show a non-monotonic T dependence even for s-wave order parameters depending on the Fermi surface shapes and densities of states (DOS). Most interesting are the order parameters Δ = Δ₀(1 + η cos k_za) which yield nearly linear increase in γ(T) on warming, a ubiquitous feature seen in many Fe-based superconductors.

Our approach is basically that of HW, although formally the equations we solve for H_c2(T) are different and can be applied to anisotropic and multi-band situations. To establish the link to HW, we start the discussion with the one-band case. The problem of the second-order phase transition at H_c2 is addressed here on the basis of linearized (the order parameter Δ → 0) quasiclassic Eilenberger equations [19]. For clean materials we have

$$\begin{equation} (2\omega +{\bit v}{\boldsymbol {\cdot}} \boldsymbol{\Pi})\,f=2\Delta /\hbar, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\Psi}{2\pi T} \ln \frac{T_{\rm c}}{T}= \sum_{\omega >0}^{\infty}\bigg(\frac{\Psi}{\hbar\omega}-\langle \Omega \, f \rangle \bigg). \end{equation} \tag{ 2 }$$

Here, v is the Fermi velocity, Π = ∇ + 2πiA/φ₀, A is the vector potential and φ₀ is the flux quantum. Δ(r, k_F) is the order parameter that in general depends on the position k_F at the Fermi surface (or on v). The function f(r, v, ω) originates from Gor'kov's Green's function integrated over energy near the Fermi surface. Further, N(0) is the total DOS at the Fermi level per spin; the Matsubara frequencies are ω = πT(2n + 1)/ℏ with an integer n. The averages over the Fermi surface are shown as 〈...〉. The Eilenberger function $g=\sqrt{1-ff^+}=1$ at H_c2. The temperature T is in energy units, i.e. k_B = 1.

The self-consistency equation (2) is written for the general case of anisotropic gaps: Δ = Ψ(r, T) Ω(k_F). The function Ω(k_F) which determines the k_F dependence of Δ is normalized so that

$$\begin{equation} \langle\Omega^2\rangle=1, \end{equation} \tag{ 3 }$$

for details see, e.g., [20]. Equation (2) corresponds to the factorizable coupling potential, V(k, k') = V₀Ω(k)Ω(k'). This popular approximation [21] works well for one band materials with anisotropic coupling. We show in sections 5 and 6 how this convenient shortcut can be generalized to a multi-band case.

We now recast the self-consistency equation (2) by writing the solution of equation (1) in the form

$$\begin{equation} f={2\over \hbar}\int_0^{\infty}{\rm d}\rho\,{\rm e}^{-\rho (2\omega +{\bit v}{\boldsymbol {\cdot}} \boldsymbol{\Pi})} \Delta, \end{equation} \tag{ 4 }$$

using the identity

$$\begin{equation} {1\over \hbar \omega}={2\over \hbar}\int_0^{\infty}{\rm d}\rho\,{\rm e}^{-2\omega\rho}, \end{equation} \tag{ 5 }$$

and by summing up over ω:

$$\begin{equation} - \Psi \ln t= \int_0^{\infty} \frac{{\rm d}u}{\sinh\,u} ( \Psi - \langle\Omega^2 {\rm e}^{-\rho {\bit v}{\boldsymbol {\cdot}} \boldsymbol{ \Pi}}\Psi\rangle ), \end{equation} \tag{ 6 }$$

where u = 2πTρ/ℏ and t = T/T_c. Hence, we got rid of the summation over ω, a convenient feature for further analysis.

The self-consistency equation (6) can be further rewritten in the form free of the divergent factor 1/ sinh u in the integrand. To this end, we integrate by parts the right-hand side (rhs) of equation (6) using du/ sinh u = d ln tanh(u/2). The first term on the rhs diverges:

$$\begin{equation} \int_0^{\infty} \frac{{\rm d}u}{\sinh u} \Psi =-\Psi \ln\tanh\frac{\pi T\rho}{\hbar}\bigg|_{\rho\to 0}. \end{equation} \tag{ 7 }$$

The second term transforms to

$$\begin{eqnarray} &&- \bigg\langle\Omega^2 \bigg[\ln\tanh\frac{\pi T\rho}{\hbar}\,{\rm e}^{-\rho {\bit v}\boldsymbol{ \Pi}}\bigg]_{\rho= 0}^{\infty} \Psi\bigg\rangle\nonumber\\ &&- \bigg\langle\Omega^2 {\bit v}\boldsymbol{ \Pi}\int_0^{\infty} {\rm d}\rho\, \ln\tanh\frac{\pi T\rho}{\hbar} {\rm e}^{-\rho {\bit v}\boldsymbol{ \Pi}}\Psi\bigg\rangle. \end{eqnarray} \tag{ 8 }$$

The first term here and the contribution (7) cancel, and we obtain instead of equation (6):

$$\begin{equation} \Psi \ln t= \int_0^{\infty}\rmd\rho\, \ln\tanh\frac{\pi T\rho}{\hbar}\, \langle\Omega^2 {\bit v}\boldsymbol{ \Pi}\,\rme^{-\rho {\bit v}{\boldsymbol {\cdot}} \boldsymbol{ \Pi}}\Psi\rangle. \end{equation} \tag{ 9 }$$

Here, the singularity at ρ → 0 is integrable.

2.1. H_c2 near T_c

In this domain, the gradients Π ∼ ξ⁻¹ → 0, and one can expand exp(−ρvΠ) in the integrand (9) and keep only the linear term:

$$\begin{equation} -\Psi \delta t= \frac{7\zeta (3)\hbar^2}{16\pi^2T_{\rm c}^2}\, \langle\Omega^2 ({\bit v}{\boldsymbol {\cdot}} \boldsymbol{\Pi})^2\Psi\rangle, \end{equation} \tag{ 10 }$$

where $\int_0^\infty {\rm d}x\, x\ln \tanh x=-7\zeta(3)/16$ with ζ(3) = 1.202. This is, in fact, the anisotropic version of the linearized Ginzburg–Landau (GL) equation

$$\begin{equation} - \xi^2_{ik} \Pi_i\Pi_k \Psi =\Psi, \end{equation} \tag{ 11 }$$

with

$$\begin{equation} \xi^2_{ik} = \frac{7\zeta (3)\hbar^2}{16\pi^2T_{\rm c}^2 \tau}\, \langle\Omega^2 v_iv_k \rangle,\tqs \tau=1-t, \end{equation} \tag{ 12 }$$

the result of Gor'kov and Melik-Barkhudarov [22]. Solving the eigenvalue problem for equation (11), which is similar to one for a charged particle in uniform magnetic field, see e.g. work by Tilley [23], one finds the critical fields in two principal directions of uniaxial materials:

$$\begin{eqnarray} &&H_{{\rm c2},c}=\frac{8\pi\phi_0T_{\rm c}^2\tau} {7\zeta(3)\hbar^2\langle\Omega^2 v_a^2\rangle},\\ &&H_{{\rm c2},a}=\frac{8\pi\phi_0T_{\rm c}^2\tau}{7\zeta(3)\hbar^2\sqrt{\langle\Omega^2 v_a^2\rangle\langle\Omega^2 v_{\rm c}^2\rangle}}, \end{eqnarray} \tag{ 13 }$$

so that

$$\begin{equation} \gamma ^2(T_{\rm c}) =\left(\frac{H_{{\rm c2},a}}{H_{{\rm c2},c}}\right)^2= \frac{\xi_{aa}^2}{\xi_{cc}^2}= \frac{\langle \Omega^2 v_a^2\rangle}{\langle \Omega^2 v_c^2\rangle}. \end{equation} \tag{ 14 }$$

The angular dependence

$$\begin{equation} H_{\rm c2}(\theta)=\frac{H_{{\rm c2},a}}{\sqrt{\sin^2\theta+\gamma^2\cos^2\theta}} \end{equation} \tag{ 15 }$$

is a direct consequence of equation (11) in which $\xi^2_{ik}$ is a second rank tensor (θ is the angle between the applied field and the c-axis). We argue below that equation (11) holds at H_c2, in fact, at all temperatures, and so should the angular dependence (15), the common practice to call it 'GL' notwithstanding. We show that this angular dependence holds for any order parameter symmetry and for any Fermi surface shape including two-band situations. These conclusions are, in fact, confirmed experimentally [24–30] and by calculations of [2].

This problem has been solved by HW for the whole curve H_c2(T) [1]. It is instructive and useful for the following generalization to the anisotropic case to reproduce their result within the quasiclassic scheme [31].

It was established in [1] that at H_c2(T) at any temperature, the order parameter satisfies a linear equation

$$\begin{equation} - \xi^2 \Pi^2\Delta=\Delta \end{equation} \tag{ 16 }$$

in which ξ(T) should be determined so as to satisfy the self-consistency equation. One can see that this equation is equivalent to the Schrödinger equation for a charge moving in uniform magnetic field and that the maximum field in which non-trivial solutions Δ exist is H_c2 = φ₀/2πξ². For the field along z we choose the gauge A_y = Hx. One readily verifies that in terms of operators

$$\begin{eqnarray} &&\fl \Pi^{\pm}=\Pi_x\pm {\rm i}\Pi_y,\tqs \Pi_x=\partial_x,\tqs \Pi_y=\partial_y +{\rm i}q^2x,\nonumber\\ &&\fl q^2= 2\pi H/\phi_0, \end{eqnarray} \tag{ 17 }$$

Equation (16) reads Π⁺Π⁻Δ = 0 provided q² = 1/ξ². Therefore, we obtain a useful property of Δ at H_c2(T):

$$\begin{equation} \Pi^-\Delta=\Pi^-\Psi=0. \end{equation} \tag{ 18 }$$

We now introduce v^± = v_x ± iv_y so that vΠ = (v⁻Π⁺ + v⁺Π⁻)/2 and evaluate the average 〈vΠ e^−ρvΠ Ψ〉 needed in the self-consistency equation (9). To this end, we use the known property of exponential operators:

$$\begin{equation} {\rm e}^{-\rho {\bit v}\boldsymbol{ \Pi}}\,\Psi={\rm e}^{P+Q}\,\Psi={\rm e}^P{\rm e}^Q{\rm e}^{[Q,P]/2}\,\Psi. \end{equation} \tag{ 19 }$$

Here, P = −ρv⁻Π⁺/2, Q = −ρv⁺Π⁻/2, the commutator $[Q,P] /2=- \rho^2v_{\perp}^2/4\xi^2$ and $v_{\perp}^2= v_x ^2+v_y ^2$.

Since Π⁻Ψ = 0 and e^QΨ = Ψ, we have

$$\begin{eqnarray} &&\fl {\bit v} \boldsymbol{ \Pi}{\rm e}^{-\rho {\bit v}\boldsymbol{ \Pi}}\Psi = {{\rm e}^{-\eta}\over 2}(v^-\Pi^++v^+\Pi^-) \sum_{n=0}^{\infty} \frac{(-\rho v^-\Pi^+)^n}{2^nn!}\Psi,\nonumber\\ &&\fl \eta=\frac{\rho^2v_{\perp}^2}{4\xi^2}. \end{eqnarray} \tag{ 20 }$$

After averaging over the Fermi sphere, only 〈v⁺v⁻〉 survives (use v^± = v_⊥e^±iφ):

$$\begin{equation} \big\langle {\bit v}\boldsymbol{ \Pi}\,{\rm e}^{-\rho {\bit v}\boldsymbol{ \Pi}}\,\Psi\big\rangle = -{\rho \over 2}\big\langle {\rm e}^{-\eta}v^+v^-\big\rangle\Pi^- \Pi^+\Psi. \end{equation} \tag{ 21 }$$

Now, with the help of equation (16), the self-consistency equation (9) yields

$$\begin{equation} \xi^2 \ln t= {1\over 2}\int_0^{\infty}{\rm d}\rho\,\rho\, \ln\tanh\frac{\pi T\rho}{\hbar}\, \left \langle {\rm e}^{-\eta}v_{\perp}^2\right\rangle. \end{equation} \tag{ 22 }$$

Introducing the dimensionless field

$$\begin{equation} h= \frac{\hbar^2v_{\rm F}^2}{4\pi^2T_{\rm c}^2\xi^2}. \end{equation} \tag{ 23 }$$

and a variable s = πT_cρ/ℏ, we rewrite equation (22) for ξ(t) as an equation for h(t):

$$\begin{equation} \ln t= 2h\int_0^{\infty}{\rm d}s\,s\, \ln\tanh(st) \langle \mu {\rm e}^{-\mu h s^2}\rangle, \end{equation} \tag{ 24 }$$

where $\mu=v_{\perp}^2/v_{\rm F}^2$ and v_F is the Fermi velocity. The field h is the upper critical field in units of $\phi_0/ (\hbar^2v_{\rm F}^2/2\pi T_{\rm c}^2)$. Thus, solving equation (24) with respect to h(t) we have the HW solution of the problem for the isotropic one-band case. In the following we refer to equation (24) as the HW result although in the original work they obtained a different form of the equation for H_c2(T).

At arbitrary T, equation (24) can be solved numerically; the exceptions are T → 0 and T → T_c. Since h(0) is finite, the integral over s is truncated at $s\sim 1/\sqrt{\mu h}$; therefore, for low enough t we have ln tanh(st) ≈ ln t + ln s. The integration over s then yields

$$\begin{equation} {{\bit C}\over 2}= \bigg\langle \ln \frac{2\pi T_{\rm c}\xi(0)}{\hbar v_{\perp}}\, \bigg\rangle, \end{equation} \tag{ 25 }$$

where C = 0.577 is the Euler constant. The averaging over the Fermi sphere with v_⊥ = v_F sin θ is readily performed:

$$\begin{equation} \xi^2(0)=\frac{\hbar^2 v_{\rm F}^2}{2\pi^2 T_{\rm c}^2}\,{\rm e}^{{\bit C}-2}. \end{equation} \tag{ 26 }$$

Thus, we obtain

$$\begin{equation} H_{\rm c2}(0)=\frac{\phi_0}{2\pi\xi^2(0)}=\frac{\phi_0\pi T_{\rm c}^2}{2\hbar^2v_{\rm F}^2}\,{\rm e}^{2-{\bit C}}, \end{equation} \tag{ 27 }$$

the value obtained variationally by Gor'kov [22] and proven to be exact by HW.

Close to T_c, H_c2 is obtained as isotropic limits of equations (13). It is instructive, however, to see how this can be deduced directly from equation (24). In this domain, ln t ≈ t − 1 = −τ and h ∝ τ. Then, the integral in equation (24) should be evaluated in zero order in τ, in other words, h in exp(−μhs²) can be set zero. Integration over s gives −(7ζ(3)/16)〈μ〉, whereas $\langle \mu \rangle=\langle v^2_\perp/v_{\rm F}^2 \rangle=2/3$:

$$\begin{equation} h= 12\,\tau/7\zeta(3). \end{equation} \tag{ 28 }$$

The reduced HW field

$$\begin{eqnarray} &&\fl h^*(0)={H_{\rm c2}(0)\over |H_{\rm c2}^{\prime}(T_{\rm c})|T_{\rm c}}=\frac{h(0)}{|h^\prime(1)|}= \frac {7\zeta(3) e^2}{48 \, {\rm e}^{{\bit C}}} \approx 0.727, \end{eqnarray} \tag{ 29 }$$

the correct HW value.

The central point of the HW paper is the proof that the linearized GL equation (16) holds not only near T_c but along the whole curve H_c2(T). For anisotropic materials equation (16) should be replaced with equation (11) where all components of the tensor (ξ²)_ik should be determined from the self-consistency equation. We consider here uniaxial materials for which the symmetric tensor $\xi_{ik}^2$ has two independent eigenvalues. One has for the field along one of the principal crystal directions, which we call z:

$$\begin{equation} -(\xi_{xx}^2 \Pi_x^2 + \xi_{yy}^2\Pi_y^2) \Delta=\Delta. \end{equation} \tag{ 30 }$$

We denote $\xi_{xx} = \xi/\sqrt{m_x}$, $\xi_{yy} = \xi/\sqrt{m_y}$, $\xi_{zz} = \xi/\sqrt{m_z}$ with dimensionless constants m_{x, y, z}. Since the three quantities, ξ_xx, ξ_yy and ξ_zz, are replaced with four, ξ and m_{x, y, z}, we can impose an extra condition: m_xm_ym_z = 1. For uniaxial materials of interest here with m_c ≠ m_a = m_b, we introduce the anisotropy parameter $\gamma = \sqrt{m_c/m_a}$, so that all 'masses' are expressed in terms of γ: m_a = γ^−2/3, m_c = γ^4/3. It is worth noting that ms here do not necessarily have the meaning of the band theory effective masses; rather, they are parameters describing the anisotropy of H_c2; near T_c they can be expressed in terms of Fermi velocities and the gap anisotropy, equation (14) [22]. Hence, if the ansatz (11) is correct, the self-consistency equation must provide equations to determine ξ(T) and γ(T).

To make use of the property (18) in anisotropic situations we rescale coordinates in equation (30):

$$\begin{equation} x'=x \sqrt{m_x},\tqs y'=y\sqrt{m_y}, \end{equation} \tag{ 31 }$$

$$\begin{equation} \frac{\partial^2\Delta}{\partial x^{\prime\, 2}}+ \bigg(\frac{\partial}{\partial y'} +\rmi q^{\prime\, 2}x'\bigg)^2\Delta=-\frac{\Delta}{\xi^2}, \end{equation} \tag{ 32 }$$

$$\begin{equation} q^{\prime\, 2}=\frac{2\pi H}{\phi_0\sqrt{m_xm_y}}. \end{equation} \tag{ 33 }$$

Formally, equation (32) is equivalent to the isotropic equation (16). The upper critical field is then determined by q^{' 2} = 1/ξ²:

$$\begin{equation} H_{\rm c2}=\frac{\phi_0\sqrt{m_xm_y}}{2\pi\xi^2}\,. \end{equation} \tag{ 34 }$$

Therefore, we have in the uniaxial case:

$$\begin{eqnarray} &&\fl H_{{\rm c2},{c}}=\frac{\phi_0}{2\pi\xi^2}\,m_{a},\tqs H_{{\rm c2},{b}}=\frac{\phi_0}{2\pi\xi^2}\sqrt{m_am_c}\, \end{eqnarray} \tag{ 35 }$$

and

$$\begin{equation} \frac{H_{{\rm c2},{b}}}{H_{{\rm c2},{c}}}=\sqrt{\frac{m_c}{m_a}}=\gamma. \end{equation} \tag{ 36 }$$

Thus, the formally introduced 'masses' are related to the measurable ratio of H_c2s.

Any coordinate transformation results in transformation of vectors (and tensors). The scaling transformation (31) necessitates the co- and contravariant representations for vectors, see, e.g., [33] or [34]. The covariant gradients π_x = ∂/∂x' and π_y = ∂/∂y' + iq^{' 2}x' have the same properties as their isotropic counterparts Π_{x, y}. Equation (32) acquires the 'isotropic' form:

$$\begin{equation} -\xi^2(\pi_x^2+\pi_y^2)\Delta =\Delta, \end{equation} \tag{ 37 }$$

$$\begin{equation} \pi^-\pi^+\Delta=-2\Delta/\xi^2,\tqs \pi^-\Delta=0. \end{equation} \tag{ 38 }$$

The following manipulation then is similar to what was carried out in the isotropic case; one, however, should keep track of differences between co- and contravariant components of vectors.

Since the scalar products are invariant, we have vΠ = νⁱπ_i = (ν⁻π⁺ + ν⁺π⁻)/2, ν^± = ν^x ± iν^y where

$$\begin{equation} \nu^x=v_x\sqrt{m_x},\tqs \nu^y=v_y\sqrt{m_y} \end{equation} \tag{ 39 }$$

are the contravariant components of the Fermi velocity that transform as coordinates, equation (31). Further, we use the property (19) of exponential operators with

$$\begin{eqnarray} &&\fl P=-{\rho\over 2}\nu^-\pi^+,\tqs Q=-{\rho\over 2}\nu^+\pi^-,\nonumber\\ &&{1\over 2}[Q,P] =-\frac{\rho^2\nu_{\perp}^2}{4\xi^2} \end{eqnarray} \tag{ 40 }$$

and

$$\begin{equation} \nu_{\perp}^2=(\nu^x)^2+(\nu^y)^2=v_x^2m_x+v_y^2m_y. \end{equation} \tag{ 41 }$$

Since π⁻Ψ = 0 and e^QΨ = Ψ, we have

$$\begin{eqnarray} &&\fl {\bit v} \boldsymbol{\Pi}{\rm e}^{-\rho {\bit v}\boldsymbol{ \Pi}}\Psi ={{\rm e}^{-\eta}\over 2}(\nu^-\pi^++\nu^+\pi^-) \sum_{n=0}^{\infty} \frac{(-\rho\nu^-\pi^+)^n}{2^nn!}\Psi,\nonumber\\ \end{eqnarray} \tag{ 42 }$$

where $\eta=\rho^2\nu_{\perp}^2/4\xi^2$.

The next step in the 'isotropic derivation' was to use the fact that $\langle v^+v^-\rangle=2v_{\rm F}^2/3$ whereas all other averages (such as 〈v⁺v⁻v⁻〉) vanish because v^± = v_⊥e^±iφ where φ is the azimuthal angle on a sphere. To prove rigorously that this is valid for a general Fermi surface is difficult, although it is probably true for uniaxial materials of interest here. As mentioned in the introduction, to make progress in evaluation of H_c2 and its anisotropy we resort to modeling actual Fermi surfaces with spheroids. The rescaling employed above, in fact, transforms spheroids to 'spheres' in rescaled coordinates so that we can still claim that the only surviving product after averaging of equation (42) is 〈ν⁺ν⁻〉 and we obtain

$$\begin{equation} \langle \Omega^2\, {\bit v}\boldsymbol{ \Pi}\,{\rm e}^{-\rho {\bit v}\boldsymbol{ \Pi}}\,\Psi\rangle = -{\rho \over 2}\langle \Omega^2\, {\rm e}^{-\eta}\nu^+ \nu^-\rangle\pi^- \pi^+ \Psi. \end{equation} \tag{ 43 }$$

The self-consistency equation (9) now yields with the help of equation (38):

$$\begin{equation} \xi^2 \ln t= {1\over 2}\int_0^{\infty}{\rm d}\rho\,\rho\, \ln\tanh\frac{\pi T\rho}{\hbar}\, \langle\Omega^2 {\rm e}^{-\eta}\nu_{\perp}^2\rangle. \end{equation} \tag{ 44 }$$

This equation (written for a certain field orientation) contains two unknown functions, ξ(T) and γ(T). Therefore, one has to write it for two field orientations: (a) for H along c with $\eta=\rho^2\nu_{\perp}^2/4\xi^2$ and

$$\begin{equation} \nu_{\perp}^2=m_a(v_x^2+v_y^2)=\gamma^{-2/3} (v_x^2+v_y^2), \end{equation} \tag{ 45 }$$

and (b) for H along b with

$$\begin{equation} \nu_{\perp}^2= m_av_x^2+m_cv_z^2 =\gamma^{-2/3} (v_x^2+\gamma^2v_z^2). \end{equation} \tag{ 46 }$$

In principle, by solving the system of these two equations one can determine both ξ(T) and γ(T), thus proving the correctness of the ansatz (11).

Since the Fermi velocity is not a constant at anisotropic Fermi surfaces, we normalize velocities on some value v₀ for which we choose

$$\begin{equation} v_{0}^3= \frac{2E_{\rm F}^2}{\pi^2\hbar^3N(0)} , \end{equation} \tag{ 47 }$$

where E_F is the Fermi energy and N(0) is the total DOS at the Fermi level per spin. One easily verifies that v₀ = v_F for the isotropic case.

We now write the self-consistency equation (44) for H || c, i.e. with ν_⊥ of equation (45), in dimensionless form. To this end, we go to the integration variable s = πT_c/ℏ, divide both parts by ξ², and multiply and divide the integrand by $v_0^2$:

$$\begin{eqnarray} &&\fl \ln t= 2 h_c \int_0^{\infty}s\, \ln\tanh (st) \, \langle\Omega^2 \mu_c {\rm e}^{-\mu_c h_c s^2}\rangle {\rm d}s, \end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray} &&\fl h_c=\frac{\hbar^2 v_0^2\gamma^{-2/3}}{4\pi^2T_{\rm c}^2\xi^2}, \tqs \mu_c=\frac{v_x^2+v_y^2}{v_0^2}. \end{eqnarray} \tag{ 49 }$$

One can see that h_c is in fact H_c2,c in units of $\phi_0/ (\hbar^2v_0^2/2\pi T_{\rm c}^2)$. An important feature of equation (48) should be noted: it does not contain the anisotropy γ explicitly so that it can be solved for h_c(t).

Writing equation (44) for H ⊥ c with ν_⊥ of equation (46) we obtain

$$\begin{eqnarray} &&\fl \ln t = 2 h_c \int_0^{\infty}s\, \ln\tanh (st) \langle\Omega^2 \mu_{\rm b} {\rm e}^{-\mu_b h_c s^2} \rangle {\rm d}s, \end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray} &&\fl \mu_b = \frac {v_x^2+\gamma ^2v_z^2}{v_0^2} . \end{eqnarray} \tag{ 51 }$$

For the isotropic s-wave case, γ = Ω = 1 and equations (48) and (50) coincide with each other and with equation (24) of HW. We note that the integrals on the rhs of equations (48), (50) differ only in µs; for brevity we denote

$$\begin{equation} {\cal I} =\int_0^{\infty}s\, \ln\tanh (st) \langle\Omega^2 \mu {\rm e}^{-\mu h_c s^2}\rangle {\rm d}s. \end{equation} \tag{ 52 }$$

Thus, the general scheme for solving for H_c2(T) and its anisotropy consists of (a) solving equation (48) for h_c(t) and (b) for the now known h_c(t), solving equation (50) for γ(t).

4.1. T → T_c and T → 0

Analysis of equations (48) and (50) for arbitrary temperatures is difficult because h_c(t) and t enter the integrals ${\cal I}$ in a nonlinear fashion. The situation is simpler near T_c where ln t ≈ −τ = t − 1 and h_c ∝ τ. Therefore, $\cal I$ can be evaluated in zero order in τ, in other words, h_c in the $\exp(-\mu^2_ch_cx^2)$ can be set zero and t = 1 in ln tanh(xt). We obtain after integration:

$$\begin{equation} {\cal I} =-\frac{7\zeta(3)}{16} \langle\Omega^2 \mu \rangle, \end{equation} \tag{ 53 }$$

where for H || c one takes μ = μ_c whereas μ = μ_b for H ⊥ c. Equation (48) now yields

$$\begin{eqnarray} &&\fl h_c= \frac{8\tau}{7\zeta(3) \langle \Omega^2\mu_c \rangle},\tqs h_c^\prime(1)= -\frac{8}{7\zeta(3) \langle \Omega^2\mu_c \rangle}. \end{eqnarray} \tag{ 54 }$$

It is readily shown that equation (50) for γ(T_c) reduces to

$$\begin{equation} \langle\Omega^2 \mu_b \rangle = \langle\Omega^2 \mu_c \rangle. \end{equation} \tag{ 55 }$$

Using μ_b and μ_c of equations (51) and (49) one reproduces the general result (14).

As t → 0, h_c → const, and the exponential factor in $\cal I$ truncates the integrand at a finite $x\sim 1/\sqrt{\mu h_c}$. Hence, at small enough t, ln tanh(xt) ≈ ln(xt) = ln t + ln x. One can now integrate over x:

$$\begin{equation} {\cal I} (\mu )= \frac{\ln t}{2h_c} - \frac{{{\bit C}}+ \langle\Omega^2 \ln (h_c\mu ) \rangle}{4h_c} . \end{equation} \tag{ 56 }$$

Substituting this in equation (48) one obtains

$$\begin{equation} h_c(0)=\exp(-{\bit C} - \langle\Omega^2 \ln \mu_c \rangle). \end{equation} \tag{ 57 }$$

The HW ratio for a general anisotropy reads

$$\begin{eqnarray} &&\fl h^*_c(0)= \frac{h_c(0)}{h_c^\prime(1)}=\frac{7\zeta(3)}{8{\rm e}^{{\bit C}}} \left\langle \Omega^2\mu_c \right\rangle \exp \left(- \big\langle\Omega^2 \ln \mu_c \big\rangle\right). \end{eqnarray} \tag{ 58 }$$

Thus, the HW number h^*(0) = 0.727 is corrected by both Fermi surface shape and by the order parameter symmetry.

Equation (50) for γ(0) along with the expression for h_c(0) readily gives at T = 0

$$\begin{equation} \langle\Omega^2 \ln \mu_c \rangle = \langle\Omega^2 \ln \mu_b \rangle. \end{equation} \tag{ 59 }$$

The general self-consistency equation for two bands reads

$$\begin{equation} \Delta_\alpha({\bit r},{\bit k})= 2\pi T \sum_{\beta,\omega} N_\beta \langle V_{\alpha\beta}({\bit k},{\bit k}^\prime) f_\beta ({\bit r},{\bit k}^\prime,\omega) \rangle_{{\bit k}^\prime}. \end{equation} \tag{ 60 }$$

Here, α, β = 1, 2 are band indices and N_β are the bands DOS. We consider elements of V_αβ as constants, so that our model is a weak-coupling two-band theory in which the s-wave (i.e. k independent) order parameters Δ_1,2(r, T, H) should be calculated self-consistently for a given coupling matrix V_αβ.

As is commonly done, it is convenient to rewrite this equation in the form containing the measured critical temperature which is related to the effective coupling V₀ via the standard BCS formula

$$\begin{equation} \Delta(0)=\pi {\rm e}^{-{\bit C}}T_{\rm c}=2\hbar {\omega_{\rm D} \rme}^{-1/N(0)V_0}, \end{equation} \tag{ 61 }$$

where ℏω_D is the energy scale of a 'glue boson'. To this end, we introduce the normalized coupling matrix λ_αβ = V_αβ/V₀ and use the relation identical to equation (61) for T_c:

$$\begin{equation} \frac{1}{N(0) V_0}= \ln \frac{T}{T_{\rm c}}+2\pi T\sum_{\omega >0}^{\omega_{\rm D}} \frac{1}{\hbar \omega}. \end{equation} \tag{ 62 }$$

We then obtain

$$\begin{equation} -\Delta_\alpha \ln t= 2\pi T \sum_{\omega} \left(\frac{\Delta_\alpha}{\hbar\omega}- \sum_\beta n_\beta\lambda_{\alpha\beta} \langle f_\beta \rangle\right) \end{equation} \tag{ 63 }$$

with n_β = N_β/N(0).

Solving the self-consistency equation in zero field and Δ → 0 one obtains a relation for T_c (or for V₀) in terms of couplings V_αβ and DOS n_α (appendix B):

$$\begin{equation} V_0^{-2} n_1n_2d - V_0^{-1} (n_1 V_{11} + n_2 V_{22})+1=0, \end{equation} \tag{ 47 }$$

$$\begin{equation} d = V_{11} V_{22}-V_{12}^2, \end{equation} \tag{ 65 }$$

Therefore, the normalized λ_αβ = V_αβ/V₀ obey

$$\begin{equation} n_1n_2 \delta- n_1 \lambda_{11} - n_2 \lambda_{22} +1=0, \end{equation} \tag{ 66 }$$

$$\begin{equation} \delta = \lambda_{11} \lambda_{22}-\lambda_{12}^2. \end{equation} \tag{ 65 }$$

This property, in fact, means that normalized couplings λ_αβ for a given T_c have only two independent components, which are chosen in the following as λ₁₁ and λ₂₂.

To avoid misunderstanding, we stress that our notation for the normalized λ_αβ = V_αβ/V₀ differs from $\lambda_{\alpha\beta} ^{\rm lit}$ used in the literature: the latter are $\lambda_{11} ^{\rm lit}=N_1V _{11}$, $\lambda_{22} ^{\rm lit}=N_2V _{22}$ and $\lambda_{12} ^{\rm lit}=N_1V _{12}$. It should also be noted that equation (63) is not valid for the unlikely situation of zero inter-band coupling, V₁₂ = 0, because two decoupled condensates in general have two different critical temperatures.

As was done in section 2, we can transform the self-consistency equation (63) to

$$\begin{eqnarray} &&\fl \Delta_\alpha \ln t=\sum_\beta n_\beta\lambda_{\alpha\beta} \int_0^\infty \rmd\rho\ln\tanh\frac{{\boldsymbol \pi} T\rho}{\hbar} \langle {\bit v}{\boldsymbol {\cdot}} {\Pi} \rme^{-\rho{\bit v} {\boldsymbol \Pi}}\Delta \rangle_\beta,\nonumber\\ \end{eqnarray} \tag{ 68 }$$

the averaging in the last term here is only over β-band.

We have generalized the HW isotropic one-band approach by showing that linearized GL equation −(ξ²)_ikΠ_iΠ_kΔ = Δ holds everywhere along the H_c2 line. Clearly, the tensor (ξ²)_ik gives the length scales of spatial variations of Δ at H_c2. Before solving the self-consistency equation for two-band systems (68) it is instructive to recall the situation in the GL domain of a two-band material, which has recently been discussed in some detail for the case of two isotropic bands [35–37].

The system of two GL equations for two order parameters written in terms of coefficients of the GL energy expansion looks–at first sight—as if it contains two different coherence lengths, i.e. each order parameter varies in space with its own length scale different from the other. It has been shown, however, that at T = T_c these two length scales coincide, provided of course that the system has a single T_c [35]. In fact, two GL equations can be written as one with a single coherence length ξ which is related in a non-trival manner to coefficients of the GL energy functional. Thus, for a material with two isotropic bands, the linearized GL equation is the same for both bands:

$$\begin{equation} - \xi^2 \Pi^2 \Delta_\alpha =\Delta_\alpha,\tqs\alpha=1,2. \end{equation} \tag{ 69 }$$

When the two bands are anisotropic, we can look for solutions of the self-consistency system (68) which satisfies at H_c2 the linear equation

$$\begin{equation} - (\xi^2)_{ik} \Pi_i\Pi_k \Delta_\alpha=\Delta_\alpha,\tqs \alpha=1,2. \end{equation} \tag{ 70 }$$

All components of the tensor (ξ²)_ik are to be determined from the self-consistency equations. One can consider equation (70) as an ansatz which should be substituted in the self-consistency relations. If one succeeds in finding such a (ξ²)_ik so that the latter are satisfied, the ansatz (70) is proven correct.

Repeating the derivation of section 4 we obtain

$$\begin{equation} \langle {\bit v}{\boldsymbol {\cdot}} \boldsymbol{ \Pi} \rme^{-\rho{\bit v}{\boldsymbol {\cdot}} {\boldsymbol \Pi}}\Delta \rangle_\beta= \frac{\rho\Delta_\beta}{2\xi^2}\langle \nu_\perp^2 \rme^{-\eta} \rangle_\beta, \end{equation} \tag{ 71 }$$

where ξ is the average coherence length related to the eigenvalues of $\xi_{ik}^2$: $\xi_{aa}^2=\xi^2\gamma^{2/3}$ and $\xi_{cc}^2=\xi^2\gamma^{-4/3}$. Further, $\eta=\rho^2\nu_\perp^2/4\xi^2$ and ν_⊥ is given in equation (45) for H || c and in equation (46) for H ⊥ c. Substituting this in system (68) we obtain after straightforward algebra:

$$\begin{eqnarray} &&a_{11}\Delta_1+a_{12}\Delta_2=0,\\ &&a_{21}\Delta_1+a_{22}\Delta_2=0, \end{eqnarray} \tag{ 72 }$$

with

$$\begin{equation*} a_{11} =\xi^2\ln t-\lambda_{11} {\cal J}_1,\tqs a_{12}=- \lambda_{12} {\cal J}_2, \end{equation*}$$

$$\begin{equation} a_{21}=- \lambda_{21} {\cal J}_1 ,\tqs a_{22}=\xi^2\ln t- \lambda_{22} {\cal J}_2, \end{equation} \tag{ 73 }$$

$$\begin{equation} {\cal J}_\alpha=\frac{n_\alpha}{2} \int_0^\infty \rmd\rho\,\rho\ln\tanh\frac{\pi T\rho}{\hbar} \langle \nu_\perp^2 \rme^{-\eta} \rangle_\alpha. \end{equation} \tag{ 74 }$$

Zero determinant of the linear system (72) gives ξ(t):

$$\begin{eqnarray} &&\fl \xi^4(\ln t)^2- \xi^2 \ln t (\lambda_{11} {\cal J}_1+\lambda_{22} {\cal J}_2) + \delta{\cal J}_1 {\cal J}_2=0.\qquad \end{eqnarray} \tag{ 75 }$$

For a single band n₂ = 0 and λ₁₁ = 1 and one obtains equation (44).

The order parameters Δ_1,2 at H_c2, as solutions of system (72), are determined only within an arbitrary factor, whereas their ratio is fixed: Δ₁/Δ₂ = −a₁₂/a₁₁. This means that Δ₁ and Δ₂ at H_c2 have the same coordinate dependences and, in particular, that they have coinciding zeros (this, of course, follows already from equation (70)). The gaps ratio is, in general, temperature dependent (for brevity, we use the term 'gap' instead of 'order parameter' although the latter is more accurate).

Introducing a dimensionless field h_c according to equation (49) one rewrites equation (75) as an equation for h_c(t):

$$\begin{eqnarray} &&\fl (\ln t)^2- 2h_c (n_1\lambda_{11} {\cal I}_1+n_2\lambda_{22} {\cal I}_2)\ln t \nonumber\\ &&+ 4h_c^2(n_1\lambda_{11}+n_2\lambda_{22} -1) {\cal I}_1 {\cal I}_2=0, \end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray} &&\fl {\cal I}_\alpha = \int_0^\infty \rmd s\,s\ln \tanh(st)\langle \mu_c \rme^{-\mu_c s^2h_c} \rangle_\alpha, \end{eqnarray} \tag{ 77 }$$

where μ_c is given by equation (49) for the corresponding band, and we took account of equations (66) and (65). As in the one-band case, this equation does not contain the anisotropy parameter γ and can be solved numerically for h_c(t). Equations of a structure similar to (76) have been employed in studies of H_c2 in two-band superconductors [8, 9].

Given h_c(t), one finds the upper critical field:

$$\begin{equation} H_{{\rm c2},c}= \frac{2\pi T_{\rm c}^2\phi_0}{\hbar^2v_0^2}\,h_c(t), \end{equation} \tag{ 78 }$$

where v₀ is expressed in terms of the Fermi energy and the total DOS in equation (47).

Writing the self-consistency condition for H ⊥ c, we obtain equation (76) in which, however, h_c(t) is now known and μ_{c, α} in integrals (77) is replaced with μ_{b, α}(γ) given by equation (51) for each band. Solving this numerically, one obtains γ(t).

The case of T → T_c is treated as was done for a one-band situation:

$$\begin{equation} {\cal I}_\alpha |_{T_{\rm c}} =-\frac{7\zeta(3)}{16} \langle \mu \rangle_\alpha \qquad \end{equation} \tag{ 79 }$$

(take equation (53), set Ω = 1 for the s-wave and add the band index). For H || c one takes μ = μ_c whereas μ = μ_b for H ⊥ c. The same argument which led to equation (56) gives for two bands at low temperatures:

$$\begin{equation} {\cal I}_\alpha(\mu ) |_{t\to 0} = \frac{\ln t}{2h_c} - \frac{{{\bit C}}+ \langle \ln (h_c\mu )\rangle_\alpha}{4h_c}. \end{equation} \tag{ 80 }$$

One can make progress analytically in looking for H_c2 near T_c and for T → 0. This calculation is useful for checking the numerical routine; for the sake of brevity we do not provide these somewhat cumbersome results.

5.1. Ratio Δ₂/Δ₁ at H_c2

It follows from system (72):

$$\begin{equation} \frac{\Delta_1}{\Delta_2} =-\frac{a_{22}}{a_{21}}=\frac{\ln t- 2n_2\lambda_{22}h_c{\cal I}_2}{2n_1\lambda_{12}h_c{\cal I}_1}, \end{equation} \tag{ 81 }$$

where

$$\begin{equation} \lambda_{12} =\sqrt{\frac{n_1n_2\lambda_{11}\lambda_{22}-n_1 \lambda_{11}-n_2 \lambda_{22}+1}{n_1n_2}}. \end{equation} \tag{ 82 }$$

We stress again that the coordinate independent ratio Δ₂(r)/Δ₁(r) makes sense only for order parameters having the same coordinate dependence (in particular, the same zeros and the same phases).

As T → 0, one can keep only the leading logarithmic term in ${\cal I}_\alpha$ of equation (80) to obtain¹

$$\begin{eqnarray} &&\frac{\Delta_1}{\Delta_2}\bigg|_{T=0} = \frac{1- n_2\lambda_{22}}{n_2\lambda_{12}}. \end{eqnarray} \tag{ 83 }$$

We point out that the zero-T gap ratio does not depend either on the Fermi surfaces involved or on the value of H_c2(0). Also, this ratio at t = 0 is the same for both field orientations; this is not the case for t ≠ 0.

Other than s-wave order parameters emerge if the coupling V(k, k') responsible for superconductivity is not a constant (or a 2 × 2 matrix of k independent constants). The formally simplest way to consider different from s-wave order parameters without going to details of microscopic interactions is to use a 'separable' potential:

$$\begin{equation} V_{\alpha\beta}({\bit k},{\bit k}^\prime) = V^{(0)}_{\alpha\beta}\Omega_\alpha({\bit k})\Omega_\beta({\bit k}^\prime), \end{equation} \tag{ 84 }$$

where $V^{(0)}_{\alpha\beta}$ is a k independent matrix, and look for the order parameters in the form

$$\begin{equation} \Delta_{\alpha} ({\bit r},T,{\bit k}) = \Psi_{\alpha}({\bit r},T)\Omega_\alpha({\bit k}) \end{equation} \tag{ 85 }$$

with the normalization 〈Ω²〉_α = 1 for each band. One can see that this leads to the self-consistency equation

$$\begin{eqnarray} &&\fl -\Psi_\alpha\ln t= 2\pi T \sum_{\omega}\left(\frac{\Psi_\alpha}{\hbar\omega} -\sum_\beta n_\beta \lambda_{\alpha\beta}\langle \Omega_\beta f_\beta \rangle_\beta\right), \end{eqnarray} \tag{ 86 }$$

where

$$\begin{equation} \lambda_{\alpha\beta} = V^{(0)}_{\alpha\beta}/V_0 , \end{equation} \tag{ 87 }$$

see also appendix B. The same algebra as in section 5 results in equation (76) for h_c(t) with

$$\begin{equation} {\cal I}_\alpha = \int_0^\infty \rmd s\,s\ln \tanh(st)\langle \Omega^2\mu_c \rme^{-\mu_c h_c s^2}\rangle_\alpha. \end{equation} \tag{ 88 }$$

As before, one calculates the anisotropy γ(t) replacing μ_c with μ_b(γ).

The Fermi surface as an ellipsoid of rotation is an interesting example in its own right and as a model system for calculating H_c2 in uniaxial materials with closed Fermi surfaces. Since H_c2 is weakly sensitive to fine details of Fermi surfaces, calculations performed for ellipsoids might be relevant for realistic shapes as well.

Similarly, open Fermi surfaces (extending to boundaries of the Brillouin zone) in uniaxial materials can be studied qualitatively by considering rotational hyperboloids. The formal treatment of these shapes is similar to that of ellipsoids. This work is still in progress, we show some of it in appendix D.

Consider an uniaxial superconductor with the electronic spectrum

$$\begin{equation} E({\bit k})=\hbar^2\left(\frac{k_x^2+k_y^2}{2m_{ab}}+\frac{k_z^2}{2m_c}\right), \end{equation} \tag{ 89 }$$

so that the Fermi surface is an ellipsoid of rotation with z being the symmetry axis.

In spherical coordinates (k, θ, φ) we have

$$\begin{eqnarray} &&\fl E({\bit k})=\frac{\hbar^2k^2}{2m_{ab}}\left( \sin^2\theta+\frac{m_{ab}}{m_c}\cos^2\theta \right)=\frac{\hbar^2k^2}{2m_{ab}}\Gamma(\theta), \end{eqnarray} \tag{ 90 }$$

so that

$$\begin{equation} k_{\rm F}^2(\theta)=\frac{2m_{ab} E_{\rm F}}{\hbar^2 \Gamma(\theta)}. \end{equation} \tag{ 91 }$$

The Fermi velocity is v(k) = ∇_kE(k), with the derivatives taken at k = k_F:

$$\begin{eqnarray} &&\fl v_x=\frac{v_{ab} \sin\theta\cos\phi}{\sqrt{\Gamma(\theta)}}, \tqs v_y=\frac{v_{ab}\sin\theta\sin\phi}{\sqrt{\Gamma(\theta)}},\nonumber\\ &&\fl v_z=\varepsilon \frac{v_{ab}\cos\theta}{\sqrt{\Gamma(\theta)}},\tqs \varepsilon=\frac{m_{ab}}{m_c}, \tqs v_{ab}= \sqrt{\frac{2E_{\rm F}}{m_{ab}}}.\nonumber\\ \end{eqnarray} \tag{ 92 }$$

The value of the Fermi velocity, $v =(v_x^2+v_y^2+v_z^2)^{1/2}$, is given by

$$\begin{equation} v= v_{ab}\sqrt{\frac{\sin^2\theta+\varepsilon^2 \cos^2\theta}{\sin^2\theta+\varepsilon\cos^2\theta}}=v_{ab}\sqrt{\frac{\Gamma_1(\theta)}{\Gamma(\theta)}}. \end{equation} \tag{ 93 }$$

The DOS N(0) is defined as an integral over the Fermi surface:

$$\begin{equation} N(0)=\int\frac{\hbar^2\rmd^2{\bit k}_{\rm F}}{(2\pi\hbar)^3v}. \end{equation} \tag{ 94 }$$

The integral over the Fermi surface can be done by integrating over the solid angle dΩ = sin θ dθ dφ:

$$\begin{equation} N(0)= \frac{m^2_{ab} v_{ab}}{2\pi^2\hbar^3} \int\frac{\rmd\Omega}{4\pi \sqrt{\Gamma(\theta)\Gamma_1(\theta)}}. \end{equation} \tag{ 95 }$$

The Fermi surface average of a function A(θ, φ) is

$$\begin{eqnarray} &&\fl \langle A( \theta,\phi) \rangle =\frac{1}{D} \displaystyle\int\frac{\rmd\Omega\,A ( \theta,\phi)}{4\pi \sqrt{\Gamma(\theta)\Gamma_1(\theta)}}, \end{eqnarray} \tag{ 96 }$$

$$\begin{eqnarray} &&\fl D = \int\frac{\rmd\Omega} {4\pi \sqrt{\Gamma(\theta,\varepsilon)\Gamma_1(\theta,\varepsilon)}} =\frac{F(\cos^{-1}\sqrt{\varepsilon},1-\varepsilon)}{\sqrt{1-\varepsilon}},\nonumber\\ \end{eqnarray} \tag{ 97 }$$

where F is an incomplete elliptic integral of the first kind. If the function A depends only on the polar angle θ, one can employ the variable u = cos θ:

$$\begin{eqnarray} &&\fl \langle A( \theta ) \rangle =\frac{1}{D(\varepsilon)} \int_0^1\frac{\rmd u\,A(u)}{\sqrt{\Gamma(u,\varepsilon)\Gamma_1(u,\varepsilon)}},\qquad\qquad \end{eqnarray} \tag{ 98 }$$

$$\begin{eqnarray} &&\fl \Gamma = 1+(\varepsilon-1)u^2,\tqs \Gamma_1 = 1+(\varepsilon^2-1)u^2 . \end{eqnarray} \tag{ 99 }$$

It is useful for the following to have a relation between v_ab and v₀ of equation (47) for a one-band situation:

$$\begin{equation} v_{ab}^3= D (\varepsilon)\, v_0^3. \end{equation} \tag{ 100 }$$

7.1. H || c

One obtains μ_c using equations (47), (49) and (92):

$$\begin{equation} \mu_ c =D^{2/3}(\varepsilon)\frac{\sin^2\theta}{\Gamma (\theta,\varepsilon)}. \end{equation} \tag{ 101 }$$

Hence, we can solve equation (48) for h_c(t) for a spheroid with a fixed m_ab/m_c = ε.

Examples of this calculation for the s-wave order parameter, Ω = 1, are shown in the upper panel of figure 1 for a prolate ellipsoid with ε = 0.2 and for an oblate one with ε = 5 (the latter corresponds to the ratio of spheroid semi-axes $\sqrt{5})$; equations (76) and (77) for h_c are solved numerically by employing any of the available packages, such as Mathematica, able to find roots of nonlinear equations and do multiple integrals. The plotted h_c(t) is H_c2,c normalized on $\phi_0/ (\hbar^2v_0^2/2\pi T_{\rm c}^2)$ with v₀ given in equation (47) in terms of the Fermi energy and the total DOS. For comparison, the same results are shown in the lower panel of figure 1 in the traditional HW normalization

$$\begin{equation} h_c^*= \frac{H_{{\rm c2},c}(T)}{T_{\rm c}H_{{\rm c2},c} ^\prime(T_{\rm c})}= \frac{h_{c}(t)}{h_{c}^\prime(1)}. \end{equation} \tag{ 102 }$$

It is seen, therefore, that for one-band s-wave materials, although the actual values of h_c(0) vary, the curves of H_c2,c(T) have qualitatively similar shapes for different Fermi surfaces.

Zoom In Zoom Out Reset image size

7.2. γ(t)

To solve equation (50) for H || a, we need

$$\begin{eqnarray} \fl \mu_b &=\frac{v_x^2+\gamma^2v_z^2}{v^2_0}\nonumber\\ \fl\quad & =D^{2/3}(\varepsilon)\frac{\sin^2\theta\cos^2\varphi+ \gamma^2 \varepsilon^2\cos^2\theta}{\Gamma (\theta,\varepsilon)}. \end{eqnarray} \tag{ 103 }$$

The anisotropy parameter γ is calculated with the help of equations (50) and (51) and shown in figure 2. It is worth observing that for the s-wave case, γ depends on the Fermi surface shape but is temperature independent.

Zoom In Zoom Out Reset image size

One can see that γ = 1 for a Fermi sphere with ε = 1, as it should. One can show that γ(ε) behave approximately as $1/\sqrt{\varepsilon}$. In particular, we observe that for oblate Fermi spheroids, γ < 1, i.e. H_c2,ab < H_c2,c.

7.3. d-wave on a one-band ellipsoid

For this case Ω = Ω₀ cos 2φ and one can verify that the condition 〈Ω²〉 = 1 yields $\Omega_0^2=2$, the same value for any Fermi spheroid. Equation (48) then can be solved numerically with the results shown in figure 3 for values of ε given in the caption. The anisotropy parameter for d-wave is shown in figure 2; unlike s-wave, it decreases on warming.

Zoom In Zoom Out Reset image size

It is worth observing that, for a fixed ε, h_c(t) for s- and d-wave order parameters are the same. This is because the Fermi surface average in equation (48) involves $(1/2)\int_0^{\pi}\sin\theta \rmd\theta$ for the s-wave whereas for the d-wave we have $(1/4\pi)\int_0^{\pi}\sin\theta \rmd\theta\int_0^{2\pi}2\cos^22\phi\,\rmd\phi$, which has the same value as for the s-case.

7.4. Order parameter modulated along k_z

The gap function suggested by the ARPES data for Ba_0.6K_0.4Fe₂As₂ has a general form [4]

$$\begin{equation} \Delta = \Delta_0(1+\eta\cos k_za). \end{equation} \tag{ 104 }$$

This order parameter varies along the Fermi surface with changing k_z; it does not have zeros if |η|<1. Depending on the sign of η, it has maximum or minimum at the 'equator' k_z = 0.

To apply this dependence for Fermi spheroids, we write

$$\begin{equation} k_z^2 =k_{\rm F}^2 \cos^2\theta = \frac{2m_{ab} E_{\rm F} \cos^2\theta}{\hbar^2\Gamma(\varepsilon,\theta)}, \end{equation} \tag{ 105 }$$

where k_F is taken from equation (91). Since $m_{ab}=2E_{\rm F}/v_{ab}^2$ and v_ab = v₀D^1/3, we obtain

$$\begin{equation} k_za = \frac{2 E_{\rm F} a\cos \theta}{\hbar v_0D^{1/3}(\varepsilon)\sqrt{\Gamma(\varepsilon,\theta)}}. \end{equation} \tag{ 106 }$$

We now choose the length scale

$$\begin{equation} a = \frac{\hbar v_0}{2 E_{\rm F}}, \end{equation} \tag{ 107 }$$

so that

$$\begin{equation} k_za = \frac{\cos \theta}{D^{1/3}(\varepsilon)\sqrt{\Gamma(\varepsilon,\theta)}}. \end{equation} \tag{ 108 }$$

Note that in the isotropic case, the length a = a₀/(3π²)^1/3 with a₀ being the interparticle spacing (the unit cell size).

To adopt the order parameter (104) for our formalism we define

$$\begin{equation} \Omega^2 = \frac{(1+\eta\cos (k_za))^2}{\langle(1+\eta\cos (k_za))^2\rangle}, \end{equation} \tag{ 109 }$$

as to satisfy the normalization 〈Ω²〉 = 1 (the average in the denominator is calculated according to equation (D10)).

Numerical evaluation of h_c(t) for parameters given in the caption of figure 4 results in a curve qualitatively similar to that of HW (the upper panel); however, the anisotropy parameter decreases on warming for η > 0 as shown in the middle panel. It is quite remarkable that changing the sign of η, i.e. going from order parameters with a maximum at the Fermi surface 'equator' k_z = 0 (η > 0) to ones with a minimum (η < 0), not only changes the temperature dependence of γ to the opposite but also affects its absolute values as well (for positive η, the anisotropy parameter is noticeably larger than for η < 0 for other parameters kept the same).

Zoom In Zoom Out Reset image size

One can readily evaluate how the anisotropy of the London penetration depth Λ, γ_Λ = Λ_c/Λ_ab, changes with temperature for η = −0.5 (for which the H_c2 anisotropy is shown in the lowest panel of figure 4). To this end we note that the GL theory requires the same values of these two anisotropies at T_c, so that γ_Λ(T_c) = γ(T_c) ≈ 2.17 according to figure 4. As T → 0, $\gamma_\Lambda^2(0)\to\langle v_x^2\rangle / \langle v_z^2\rangle$ (in the clean limit, the order parameter does not enter the anisotropy of the London depth) [20, 38, 39]. The calculation of these averages is straightforward for a spheroid with ε = 0.1: γ_Λ(0) ≈ 3.38. Thus, the Λ-anisotropy decreases on warming, unlike H_c2-anisotropy. This qualitative behavior of the Λ-anisotropy is, in fact, seen in experiments on Ba(Fe_1−xCo_x)₂As₂, (Ba_1−xK_x)Fe₂As₂ and NdFeAs(O_1−xF_x) [40].

In figure 5 we demonstrate that the effect of γ increasing on warming remains if the Fermi surface changes from a prolate spheroid with ε = 0.1 to a sphere ε = 1 and to oblate spheroid with ε = 5. In particular, these features challenge the common belief that temperature dependence of the anisotropy parameter is always related to a multi-band situations.

Zoom In Zoom Out Reset image size

We note in concluding this section that other possible anisotropic order parameters can be treated within our scheme in a similar manner.

To apply the theory developed for two-band materials one first should map actual band structure upon two spheroids, the procedure we demonstrate in some detail on the well-studied MgB₂. When calculating parameters μ_c (and μ_b) needed in this mapping for each band, one should bear in mind that in the two-band situation we have

$$\begin{equation} \mu_ {c,\alpha} = \left[\frac{D(\varepsilon_\alpha)}{n_\alpha}\right]^{2/3} \,\frac{\sin^2\theta}{\Gamma (\theta,\varepsilon_\alpha)}, \end{equation} \tag{ 110 }$$

see appendix C.

8.1. MgB₂

We take this example to demonstrate that our procedure yields H_c2,c(T) and the anisotropy γ(T) in agreement with existing data (see, e.g., [12, 41]) and with calculations of [2, 3, 9]. We stress that our calculations of H_c2 are done with the same set of coupling parameters as those used for the zero-field properties of this material as described in [38]. The four Fermi sheets of MgB₂ can be grouped into two effective bands with nearly constant zero-field gaps for each group [42]. The two effective bands are mapped here upon two ellipsoids [2]. We describe this procedure in some detail.

We take the following data from the band structure calculations [44]: the relative DOS of our model are n₁ ≈ 0.42 and n₂ ≈ 0.58 for σ and π bands, respectively [42, 43]. The band calculations [43] also provide the averages over separate Fermi sheets: $\langle v_a^2\rangle_1=23$, $\langle v_c^2\rangle_1=0.5$ and $\langle v_a^2\rangle_2=33.2$, $\langle v_c^2 \rangle_2=42.2 \times 10^{14} \, {\rm cm}^{2} \, {\rm s}^{-2}$.

To map this system onto two Fermi ellipsoids, we note that the averages over spheroids are given by

$$\begin{equation} \langle v_x^2\rangle = \frac{v_{ab}^2}{2D(\varepsilon)} \int_0^1\frac{\rmd u(1-u^2)}{\Gamma ^{3/2}(\varepsilon)\Gamma_ 1 ^{1/2} (\varepsilon)}, \end{equation} \tag{ 111 }$$

$$\begin{equation} \langle v_z^2\rangle = \frac{v_{ab}^2\varepsilon^2}{D(\varepsilon)} \int_0^1\frac{\rmd u \,u^2}{\Gamma ^{3/2}(\varepsilon)\Gamma_{1}^{1/2}(\varepsilon)} ,\end{equation} \tag{ 112 }$$

where u = cos θ. The integrals here can be expressed in terms of elliptic integrals, alternatively they can be evaluated numerically. Forming the ratio $\langle v_x^2\rangle/ \langle v_z^2\rangle$ we obtain an equation which can be solved for ε. This gives ε₁ = 0.028 67 and ε₂ = 1.273. For a given ε and, e.g., $\langle v_z^2\rangle\equiv\langle v_c^2\rangle$, we obtain $v_{ab}^2$ for two ellipsoids: $v_{ab,1}^2=0.6019\times 10^{14}$ and $v_{ab,2}^2=1.436\times 10^{16}\,({\rm cm}\,{\rm s}^{-1})^{2}$. Next, we write

$$\begin{eqnarray} \fl \frac{1}{v_0^3}&= \frac{\pi^2\hbar^3(N_1+N_2)}{2E_F^2}= \frac{1}{v_{01}^3}+ \frac{1}{v_{02}^3}\nonumber\\ \fl\quad &= \frac{D(\varepsilon_1)}{v_{ab,1}^3} + \frac{D(\varepsilon_2)}{v_{ab,2}^3}, \end{eqnarray} \tag{ 113 }$$

where equation (100) has been used; this gives $v_0^2=3.867\times 10^{14}\,({\rm cm}\,{\rm s}^{-1})^{2}$, the constant used in the field normalization.

To obtain normalized coupling constants λ_αβ = V_αβ/V₀ we use the effective values (calculated including Coulomb repulsion), which in our notation read: N₁V₁₁ = 0.807, N₂V₂₂ = 0.276, N₁V₂₁ = 0.118, N₂V₁₂ = 0.086 [44]. Using equation (64) one evaluates 1/V₀N(0) = 1.211, and obtains the normalized coupling

$$\begin{equation} \lambda_{11}= \frac{V_{11}}{V_0}=\frac{N_1V_{11}}{N(0)V_0}\frac{1}{n_1}\approx 2.328. \end{equation} \tag{ 114 }$$

Similarly, we find λ₂₂ = 0.5765, λ₁₂ = λ₂₁ = 0.340. With these input parameters we have solved equations (76) and (77) for h_c(t). The result is shown in the upper panel of figure 6.

Zoom In Zoom Out Reset image size

Given h_c(t), we rewrite the same equations where μ_c is replaced with μ_b(γ) and solve them for γ(t). The latter is shown in the lower panel of figure 6.

Hence, the calculation gives h_c(0) ≈ 0.292. We now use relation (78) between the dimensionless h_c and physical H_c2,c, where T_c = 39.5 K and v₀ has been estimated above, to obtain H_c2,c(0) ≈ 2.8 T, the value close to the observed [12, 41]. The general behavior of H_c2,c(T) is close to that of HW, the fact confirmed by experiment. It should be mentioned that our result is close to the calculations of Miranovic et al [2], Dahm and Schopohl [3] and Palistrant [9]. We, however, do not reproduce the calculated H_c2,c(T) with changing curvature and a substantial upturn at low Ts of [8].

The general formula for the gaps ratio at H_c2 is given in equation (81). Figure 7 shows this ratio as a function of temperature at H_c2,c. For comparison we show the gap ratio in zero field calculated with the help of the same coupling constants [38]. It is instructive to note that the gap ratio at H_c2 exceeds the zero-field value at all temperatures, which can be interpreted as a stronger suppression of the small gap by the magnetic field than that of the large one at the leading σ band. At first sight, one would expect the two ratios to coincide as T → T_c, which our results clearly do not show. Such an expectation, however, would not be justified: even when H_c2 → 0, the superconductor in the mixed state differs from the uniform state by an extra magnetic field suppression of the order parameter.

Zoom In Zoom Out Reset image size

One often finds in the literature the statement that a small gap in MgB₂ is substantially or even completely suppressed by a large enough field. This would correspond to a substantial increase in Δ₁/Δ₂ or even divergence of this ratio at some field under H_c2. Our result, however, shows that even at H_c2 the gap ratio is finite and of the same order at all Ts. We conclude that full suppression of the small gap does not happen at any field H ⩽ H_c2. On the other hand, assuming (as was done in [2]) that the gap ratio at H_c2 is the same as that calculated in zero field, is also incorrect.

8.2. λ₁₁ ∼ λ₂₂ ≪ |λ₁₂|

This case is close to theoretical models of pnictides in which the inter-band coupling is assumed dominant. We consider here a limiting situation λ₁₁ = λ₂₂ = 0 to simplify the algebra. Indeed, for the two field directions we have

$$\begin{equation} (\ln t)^2- 4h_c^2 \, {\cal I}_1(\mu_c )\, {\cal I}_2(\mu_c )=0, \end{equation} \tag{ 115 }$$

$$\begin{equation} (\ln t)^2- 4h_c^2 \, {\cal I}_1(\mu_b)\, {\cal I}_2(\mu_b)=0. \end{equation} \tag{ 116 }$$

The first equation here can be solved for h_c(t). Since μ_b depends on γ, the second gives an equation for γ(t). The latter can also be written in a different form if one subtracts equation (115) from (116):

$$\begin{equation} {\cal I}_1(\mu_c ) \,{\cal I}_2(\mu_c )= {\cal I}_1(\mu_b)\, {\cal I}_2(\mu_b). \end{equation} \tag{ 117 }$$

Figure 8 shows the numerical solution for γ(t) for a few representative parameter sets. It is worth noting the particularly informative feature of figure 8: it shows that the anisotropy γ(t) is not necessarily a monotonic function of temperature. We note that numerical calculations of γ(t) showing an extremum at 0 < t < 1 are quite robust.

Zoom In Zoom Out Reset image size

8.2.1. Δ₁/Δ₂.

The general result for the gap ratio (81) gives for λ₁₁ = λ₂₂ = 0:

$$\begin{equation} \frac{\Delta_1}{\Delta_2} =\frac{\ln t}{2n_1\lambda_{12}h_c{\cal I}_1}. \end{equation} \tag{ 118 }$$

Near T_c we use $\cal I$ s of equation (53) and h_c from equation (115) to obtain

$$\begin{equation} \frac{\Delta_1}{\Delta_2}\bigg|_{T_c} = \frac{1}{n_1\lambda_{12}}\sqrt{\frac{\langle \Omega^2\mu_c \rangle_2}{\langle \Omega^2\mu_c \rangle_1}} = \sqrt{\frac{n_2 \langle \Omega^2\mu_c \rangle_2} {n_1 \langle \Omega^2\mu_c \rangle_1}}, \end{equation} \tag{ 119 }$$

since for λ₁₁ = λ₂₂ = 0, the normalized $\lambda_{12}^2=1/n_1n_2$. As T → 0, we can keep only the leading term ∼ ln t in $\cal I$ of equation (56):

$$\begin{equation} \frac{\Delta_1}{\Delta_2}\bigg|_{T=0} =\sqrt{\frac{n_2}{n_1}}. \end{equation} \tag{ 120 }$$

It is instructive to note that the last relation in the form $n_1\Delta_1^2= n_2\Delta_2^2$ at T = 0 suggests the equipartition of condensation energy between the bands (provided the only non-zero coupling is λ₁₂).

8.2.2. Two bands with Δ(k_z).

This example is of interest because it may have implications for understanding the behavior of H_c2(T) and, in particular, its anisotropy in Fe-based materials. In figure 9 we show h_c(t) and γ(t) for two nearly cylindrical bands with order parameters modulated along k_z according to equation (104). Modulations are characterized by η₁ = η₂ = −0.2; other parameters are given in the caption. The main feature worth paying attention to is the anisotropy γ(t), which is increasing on warming for negative ηs.

Zoom In Zoom Out Reset image size

It is worth noting that experimental anisotropy of Fe-based materials behaves qualitatively similar to that shown by the lowest curve in the lower panel of figure 9 [13]. We, however, do not have enough information to fix the necessary parameters for realistic calculations (one needs partial DOS, Fermi surfaces characterized separately for relevant bands for evaluating the geometric parameters ε, and the order parameters Δ(k_F)). Hence, we take our results as having a generic qualitative value.

The upper critical field H_c2 and its anisotropy are among the easiest properties to examine when a new superconductor is discovered. This work provides a relatively straightforward scheme for evaluating the orbital H_c2(T) and its anisotropy γ(T) for single and two-band uniaxial materials. We reproduce here the main points of our approach.

The input parameters for two-band materials are (i) the coupling matrix V_αβ (or the normalized couplings λ_αβ and T_c), (ii) the symmetries of the order parameter on two bands given as Ω_β(θ, φ), β is the band index, with the normalization (3) and (iii) the characteristics of electron systems (Fermi surfaces, averages of squared Fermi velocities, DOS).

The equation to solve for the reduced upper critical field h_c(t) parallel to the c crystal axis of a two-band clean uniaxial material reads

$$\begin{eqnarray} &&\fl (\ln t)^2- 2h_c (n_1\lambda_{11} {\cal I}_1+n_2\lambda_{22} {\cal I}_2)\ln t \nonumber\\ &&+ 4h_c^2(n_1\lambda_{11}+n_2\lambda_{22} -1) {\cal I}_1 {\cal I}_2=0,\nonumber\\ &&\fl {\cal I}_\beta = \int_0^\infty \rmd s\,s\ln \tanh(st)\langle \Omega^2\mu_c \rme^{-\mu_c s^2h_c} \rangle_\beta, \end{eqnarray} \tag{ 121 }$$

where μ_c is given by equation (49) for the corresponding band, and Ωs describe the order parameter symmetries. After h_c(t) is found, one solves equation (121), in which μ_c is replaced by μ_b(γ) of equation (51), for γ(t). The one-band case is obtained by setting n₂ = 0, n₁ = 1 and λ₁₁ = 1. The case of two s-wave bands corresponds to Ω₁ = Ω₂ = 1.

Because H_c2 is determined by equations containing only integrals over the Fermi surface, it is insensitive to fine details of the Fermi surface shape. Therefore, one can replace actual Fermi surfaces with ellipsoids (or with spheroids, for uniaxial materials). Given the averages of the squared Fermi velocities over each band, one establishes the geometry of corresponding rotational ellipsoids (the squared ratio of semi-axes, ε). This procedure is described in section 8.1 on the well-studied example of MgB₂.

One numerically solves equation (121) by employing any of the available packages (such as Mathematica) able to find roots of nonlinear equations and to do multiple integrals. The scheme can also be applied to the case of two bands with order parameters of different symmetries.

By design, our method is applicable for clean materials with a moderate H_c2(0); paramagnetic limiting effects are out of the scope of this work [8]. The method differs from those previously employed by not involving explicit coordinate dependent Δ_1,2 and minimization relative to the vortex lattice structures in calculating H_c2(T) [2, 3, 18, 45]. The main feature of the two-band derivation is that the linearized GL equation (70) is assumed to hold at H_c2(T) at any T, the ansatz proven correct by satisfying the self-consistency equation of the theory. The method is tested on the well-studied example of MgB₂ where it shows satisfactory agreement with data and with other calculations.

Our main results are as follows:

• 1.  
  We find that in clean one-band s-wave materials, the T dependence of the anisotropy γ cannot be caused by the Fermi surface anisotropy (however, the paramagnetic limit, which is out of the scope of this paper, may suppress γ at low T and cause γ(T) to increase on warming).
• 2.  
  For other than s-wave symmetry, γ depends on temperature even for one-band materials. This dependence is pronounced for open Fermi surfaces as well as for order parameters depending on k_z. Thus, the common belief that the temperature dependence of the anisotropy parameter is always related to multi-band situations is incorrect.
• 3.  
  Our scheme of calculating H_c2 for two-band materials does not utilize any assumptions about the field and temperature dependences of the order parameters Δ_α in two bands [2]. In fact, the gap ratio is calculated self-consistently and in general turns out temperature dependent. Although both Δ₁ and Δ₂ turn zero at H_c2(T), their ratio is finite and in the examples examined is larger than at zero field.
• 4.  
  The case of exclusively inter-band coupling is discussed, which might be relevant while interpreting data on H_c2 and its anisotropy in Fe-based compounds.
• 5.  
  For order parameters of the form Δ = Δ₀(1 + η cos k_za) (one of the candidates suggested for pnictides), the anisotropy parameter γ(t) depends on the sign of η (or ηs for two bands). In particular, γ(T) increases on warming in a nearly linear fashion (as for pnictides) for both ηs negative.

Some ideas described in this text were conceived in discussions with Predrag Miranovich while working on H_c2(T) of MgB₂ in 2002; VK is grateful to Predrag for this experience. We thank Andrey Chubukov for turning our attention to order parameters of the form (104) and our Ames Lab colleagues John Clem, Andreas Kreyssig, Sergey Bud'ko, Makariy Tanatar and Paul Canfield for interest and critique. We are grateful to Erick Blomberg for reading the manuscript and useful remarks. Work at the Ames Laboratory is supported by the Department of Energy - Basic Energy Sciences under Contract No DE-AC02-07CH11358.

Both sides of equation (48) diverge logarithmically when t → 0, so that these divergences, in fact, cancel out. However, in numerical work this cancellation may not always be exact, which may cause the numerical solutions to become unreliable in this limit. Here, we provide an alternative form of this equation free of this shortcoming.

To this end, consider an identity:

$$\begin{eqnarray} &&\fl \ln t=\left \langle \Omega^2 2 h_c \mu \int_0^\infty \rmd s\, s \ln (st) \, \rme^{-\mu h_c s^2}\right\rangle \nonumber\\ &&\quad +\frac{{\bit C}+ \langle\ln(\Omega^2h_c\mu) \rangle}{2}, \end{eqnarray} \tag{ A1 }$$

which is verified by direct integration. We now combine it with equation (48):

$$\begin{eqnarray} &&\fl \frac{{\bit C}+ \langle\Omega^2\ln(h_c\mu) \rangle}{4 h_c} = \int_0^\infty \rmd s\, s \ln \frac{\tanh(st)}{st}\langle \Omega^2\mu \rme^{-\mu h_c s^2}\rangle.\nonumber\\ \end{eqnarray} \tag{ A2 }$$

As t → 0, the integral on the rhs goes to zero, and we immediately obtain the result (57) for h_c(0).

This question had been discussed, e.g., in [38]. Since our notation of normalized λ_αβ differs from that in the literature, we provide here corresponding relations. The s-wave self-consistency equation for H = 0,

$$\begin{equation} \Delta_ \nu =2\pi T\sum_{\mu,\omega} N_\mu V_{\nu\mu} f_\mu(\omega), \end{equation} \tag{ B1 }$$

gives near T_c where f_μ = Δ_μ/ℏω:

$$\begin{equation} \Delta_ \nu =\sum_{\mu} n_\mu N(0)V_0\lambda_{\nu\mu} \Delta_ \mu\sum_\omega^{\omega_D}\frac{2\pi T}{\hbar\omega},\qquad \end{equation} \tag{ B2 }$$

where V₀ is to be defined. We choose it so that

$$\begin{equation} \frac{1}{N(0)V_0}=\sum_\omega^{\omega_{\rm D}} \frac{2\pi T}{\hbar\omega}=\ln\frac{2\hbar\omega_{\rm D}}{1.76 T_{\rm c}}, \end{equation} \tag{ B3 }$$

or, which is the same, $1.76 T_{\rm c}=2\hbar\omega_{\rm D}\rme^{-1/N(0)V_0}$. We then obtain a linear and homogeneous system of equations Δ_ν = ∑_μn_μλ_νμΔ_μ, the zero determinant of which gives equation (66).

For other than s-wave order parameters on two bands, we take the coupling potential in the form (84) and the order parameters as in equation (85). We then obtain the self-consistency equation $\Psi_ \nu =2\pi T\sum_{\mu,\omega} N_\mu V^{(0)}_{\nu\mu} \langle \Omega_\mu f_\mu(\omega)\rangle$. We now denote $\lambda_{\nu\mu} = V^{(0)}_{\nu\mu}/V_0$, equation (87), and recall that f_μ = Ω_μ Ψ_μ/ℏω near T_c. This gives Ψ_ν = ∑_μn_μλ_νμΨ_μ, i.e. the same system of equations as above for the s-wave case and the same zero-determinant condition (66) albeit with renormalized couplings (87).

By definition, $\mu_{c,\beta}\propto v^2_{ab,\beta}/v_0^2$ with v₀ given in equation (47). For v_{ab, β} we have two relations: one with the Fermi energy, $m_{ab,\beta}v_{ab,\beta}^2=2E_{\rm F}$, and the other in terms of the band DOS, equation (95):

$$\begin{equation} N_\beta=\frac{m^2_{ab,\beta} v_{ab,\beta}}{2\pi^2\hbar^3} \,D(\varepsilon_\beta),\tqs \beta=1,2. \end{equation} \tag{ C1 }$$

One excludes m_{ab, β} from these two relations to obtain

$$\begin{equation} v_{ab,\beta}^3= \frac{2E_{\rm F}^2}{\pi^2\hbar^3N_\beta}. \end{equation} \tag{ C2 }$$

Hence, we have

$$\begin{equation} \mu_{c,\beta}=\frac{v_{ab,\beta}^2}{v_0^2}= \left[\frac{D(\varepsilon_\beta)}{n_\beta}\right]^{2/3} , \end{equation} \tag{ C3 }$$

a clear generalization of equation (100) for the one-band spheroid. This gives equation (110) used in the text.

The theory employed above is designed to model Fermi surfaces closed within the first Brillouin zone. Here we consider an example of the Fermi surface which crosses the zone boundary, i.e. it is open. Perhaps, the simplest shape to consider is a rotational hyperboloid which is a property of the carriers energy of the form

$$\begin{equation} E({\bit k})=\hbar^2\left(\frac{k_x^2+k_y^2}{2m_{ab}}-\frac{k_z^2}{2m_c}\right). \end{equation} \tag{ D1 }$$

A schematic picture is shown in figure 10. In spherical coordinates (k, θ, φ) we have

$$\begin{eqnarray} &&\fl E({\bit k})=\frac{\hbar^2k^2}{2m_{ab}}\left( \sin^2\theta-\frac{m_{ab}}{m_c}\cos^2\theta \right)=\frac{\hbar^2k^2}{2m_{ab}}\Gamma_2(\theta),\nonumber\\ &&\fl \Gamma_2 = \sin^2\theta-\epsilon\cos^2\theta,\tqs \epsilon = m_{ab}/m_c. \end{eqnarray} \tag{ D2 }$$

It is seen from the figure that in the first quadrant of the plane k_x, k_z, the Fermi surface is situated at $\theta>\theta_m>\tan^{-1}\sqrt{\epsilon}$, i.e. everywhere at the Fermi surface Γ₂ > 0.

Zoom In Zoom Out Reset image size

The angle θ_m corresponds to the crossing of the Fermi hyperboloid with the zone boundary k_z,BZ = 2π/c where c is the unit cell size along the $\hat c$ direction:

$$\begin{equation} \fl\tan^2\theta_m =\epsilon+\frac{m_{ab} c^2E_{\rm F}}{2\pi^2\hbar^2} =\epsilon +\alpha. \end{equation} \tag{ D3 }$$

The parameter $\alpha = k ^2_{\rm F}(\pi/2)/k_{z,{\rm BZ}}^2\,$ in most situations of interest is less than unity (k_F(π/2) is the radius of the hyperboloid neck).

The Fermi momentum is given by

$$\begin{equation} k_{\rm F}^2(\theta)=\frac{2m_{ab} E_{\rm F}}{\hbar^2 \Gamma_2(\theta)}. \end{equation} \tag{ D4 }$$

The Fermi velocity is ${\bit v}({\bit k})={\boldsymbol \nabla}_{\bit k}E({\bit k})_{{\bit k}_{\rm F}}$:

$$\begin{eqnarray} &&\fl v_x=\frac{v_{ab} \sin\theta\cos\phi}{\sqrt{\Gamma_2(\theta)}}, \tqs v_y=\frac{v_{ab}\sin\theta\sin\phi}{\sqrt{\Gamma_2(\theta)}},\nonumber\\ &&\fl v_z=\epsilon \frac{v_{ab}\cos\theta}{\sqrt{\Gamma_2(\theta)}},\tqs v_{ab}= \sqrt{2E_{\rm F}/m_\parallel}. \end{eqnarray} \tag{ D5 }$$

Further, we have for $v =(v_x^2+v_y^2+v_z^2)^{1/2}$:

$$\begin{equation}v= v_{ab}\sqrt{\frac{\sin^2\theta+\epsilon^2 \cos^2\theta}{\sin^2\theta-\epsilon\cos^2\theta}}=v_{ab}\sqrt{\frac{\Gamma_1(\theta)}{\Gamma_2(\theta)}}. \end{equation} \tag{ D6 }$$

The DOS N(0) is defined by the integral of equation (94), which can be written as an integral over the solid angle dΩ = sin θ dθ dφ:

$$\begin{equation} N(0)= \frac{m^2_{ab} v_{ab}}{2\pi^2\hbar^3} \int\frac{\rmd\Omega}{4\pi \sqrt{\Gamma_2(\theta)\Gamma_1(\theta)}}, \end{equation} \tag{ D7 }$$

where the integration over θ is extended from θ_m to π − θ_m. The Fermi surface average of a function A(θ, φ) is

$$\begin{equation} \langle A \rangle =\frac{1}{D} \displaystyle\int \frac{\rmd\Omega\,A( \theta,\phi)} {4\pi \sqrt{\Gamma_2(\theta)\Gamma_1(\theta)}}, \end{equation} \tag{ D8 }$$

$$\begin{equation} D = \int\frac{\rmd\Omega} {4\pi \sqrt{\Gamma_2(\theta,\epsilon)\Gamma_1(\theta,\epsilon)}}. \end{equation} \tag{ D9 }$$

For A depending only on θ, one can employ u = cos θ:

$$\begin{eqnarray} &&\fl \langle A \rangle =\frac{1}{D(\epsilon)} \int_0^{u_m}\frac{\rmd u\,A(u)} {\sqrt{\Gamma_2(u,\epsilon)\Gamma_1(u,\epsilon)}}, \end{eqnarray} \tag{ D10 }$$

$$\begin{eqnarray} &&\fl \Gamma_2 = 1-(\epsilon+1)u^2,\tqs \Gamma_1 = 1+(\epsilon^2-1)u^2, \end{eqnarray} \tag{ D11 }$$

where the upper limit is $u_m=\cos \theta_m=1/\sqrt{1+\epsilon+\alpha}$. In particular, one obtains

$$\begin{equation} D(\epsilon,\alpha) = \frac{F(\tan^{-1}\sqrt{(1+\epsilon)/\alpha},1-\epsilon)}{\sqrt{1+\epsilon}} \tqs \end{equation} \tag{ D12 }$$

where F is an incomplete elliptic integral of the first kind.

As for ellipsoids, the relation between v_ab and v₀ defined in equation (47) for a one-band situation holds

$$\begin{equation} v_{ab}^3= D \, v_0^3, \end{equation} \tag{ D13 }$$

however, with a different D.

For ${\bit H}\parallel \hat c$, the relevant electron orbits are circular, as for spheres and rotational ellipsoids, and we do not expect qualitative deviations from the latter. Figure 11 shows h_c(t) calculated with the help of equations (48) and (49) for both s- and d-wave order parameters.

Zoom In Zoom Out Reset image size

Numerical evaluation of the anisotropy shows that γ(t) decreases on warming, figure 11. Reduction of the anisotropy on warming is substantial for a single band open Fermi surface, an interesting observation given the common belief that the temperature dependence of anisotropy is a multi-band property. We also note that this reduction is stronger than that of open Fermi surfaces.