text:
# Edge states and determination of pairing symmetry in superconducting 𝐒𝐫_𝟐𝐑𝐮𝐎_𝟒

## I Introduction

Numerous experiments suggest that the superconducting state of $`\mathrm{Sr}_2\mathrm{RuO}_4`$  is unconventional (see review ). Strong suppression of $`T_c`$ by non-magnetic impurities  and absence of the Hebel-Slichter peak in NMR experiments  indicate that it is not $`s`$-wave. The absence of spin susceptibility reduction below $`T_c`$ in the Knight shift measurements indicates the spin-triplet pairing . The $`\mu `$SR experiment suggests a superconducting state with broken time-reversal symmetry . Small-angle neutron scattering reveals a square vortex lattice, which is interpreted as the indication for a two-component order parameter . These experiments initially led to suggestion of the 2D chiral (time-reversal non-invariant), isotropic, nodeless $`p`$-wave pairing potential . However, the power-law temperature dependences found in specific heat , nuclear relaxation rate , penetration depth , thermal conductivity , and ultrasonic attenuation  indicate nodes in the energy gap. In response to these experiments, the following alternative order parameters were proposed : anisotropic $`p`$-wave , $`p`$-wave with horizontal lines of nodes , and $`f`$-waves with vertical  or horizontal lines of nodes . It was also proposed that the $`\alpha `$ and $`\beta `$ bands of $`\mathrm{Sr}_2\mathrm{RuO}_4`$ are either not superconducting  or have horizontal lines of nodes . Big in-plane anisotropy of ultrasound attenuation  may support the vertical lines of nodes. However, thermal conductivity depends very little on the orientation of an in-plane magnetic field , which is against the vertical lines .

Electron tunneling between a normal metal and a superconductor proved to be an important tool in determining superconducting symmetry. Observation of the zero-bias conductance peak (ZBCP)  due to the formation of the midgap Andreev bound states  confirmed the $`d`$-wave symmetry of the high-$`T_c`$ cuprates. Electron tunneling into a three-dimensional $`p`$-wave superconductor with the pairing potential corresponding to the B-phase of 3He was considered in Ref. . Electron tunneling in $`\mathrm{Sr}_2\mathrm{RuO}_4`$ was studied theoretically in Refs.  and  for the 2D isotropic chiral unitary and non-unitary $`p`$-waves. In this paper, we calculate the tunneling conductance curves for the alternative order parameters listed above and compare them with experiments .

## II General formalism

We model the tunneling contact by two semi-infinite regions, normal (N) and superconducting (S), with a flat interface (I) perpendicular to the $`a`$ axis of $`\mathrm{Sr}_2\mathrm{RuO}_4`$ (see Fig. 1). The $`x`$ and $`y`$ axes are selected along the crystal axes $`a`$ and $`b`$, respectively. The tunneling conductance is calculated at zero temperature by solving the Bogoliubov-de Gennes equations in the ballistic regime following Refs. . The tunneling barrier is modeled by a delta-function potential of strength $``$, so the boundary conditions at the interface are $`\psi _n|_I=\psi _s|_I`$ and $`\widehat{v}_n\psi _n|_I=\widehat{v}_s\psi _s|_I2i\psi _n|_I`$. Here $`\psi _n`$ and $`\psi _s`$ are the electron wave functions in the normal and superconducting regions, and $`\widehat{v}_{n,s}`$ are the velocity components perpendicular to the interface. The pairing potential of a triplet superconductor can be expressed as $`\widehat{\mathrm{\Delta }}=i\widehat{\sigma }_y(\widehat{𝝈}𝐝)\mathrm{\Delta }(𝐤,x)`$, where $`\widehat{𝝈}`$ are the Pauli matrices operating on the electron spin indices, the vector $`𝐝`$ indicates the direction of spin polarization, and $`𝐤`$ is the relative momentum of electrons in a Cooper pair. In this paper, we consider the cases where $`𝐝`$ does not depend on $`𝐤`$ and is pinned to the $`c`$ axis of $`\mathrm{Sr}_2\mathrm{RuO}_4`$. The spatial dependence of the pairing potential is taken to be step-like: $`\mathrm{\Delta }(𝐤,x)=\mathrm{\Delta }(𝐤)\mathrm{\Theta }(x)`$. Assuming that the electron momentum $`k_y`$ parallel to the interface is conserved, the tunneling conductance per one layer, $`\overline{G}(V)`$, can be written as an integral over $`k_y`$ :

$`\overline{G}(V)`$ $`=`$ $`{\displaystyle \frac{2e^2}{h}}l{\displaystyle _{k_y^{\mathrm{max}}}^{k_y^{\mathrm{max}}}}{\displaystyle \frac{dk_y}{h}}G(V,k_y),`$ (1)
$`G(V,k_y)`$ $`=`$ $`D{\displaystyle \frac{1+|\mathrm{\Gamma }|^2D(1D)|\mathrm{\Gamma }|^4}{\left|1(1D)\mathrm{\Gamma }^2\mathrm{exp}(i\mathrm{\Phi })\right|^2}}.`$ (2)
In Eq. (1), $`e`$ is the electron charge, $`h`$ is the Planck constant, $`l`$ is the length of the interface, and $`V`$ is the bias voltage. We consider the case where the Fermi surfaces of both the normal metal and superconductor are circular, and their radii are $`k_F^{(n)}>k_F^{(s)}`$, as shown in Fig. 1. The limits of integration are set by the smaller Fermi momentum: $`k_y^{\mathrm{max}}=k_F^{(s)}`$. In Eq. (2), $`D(k_y)=4v_n(k_y)v_s(k_y)/\{[v_n(k_y)+v_s(k_y)]^2+4^2\}`$ is the normal-state transmission coefficient, where $`v_{n,s}(k_y)`$ are the Fermi velocities components perpendicular to the interface. $`\mathrm{\Phi }(k_y)=\varphi _A(k_y)\varphi _B(k_y)`$ is the phase difference of the superconducting pairing potentials at the points on the Fermi surface connected by electron reflection from the interface (points A and B in Fig. 1) . $`\mathrm{\Gamma }(V,k_y)`$ is

$`\mathrm{\Gamma }`$ $`=`$ $`{\displaystyle \frac{eV\mathrm{sgn}(eV)\sqrt{(eV)^2|\mathrm{\Delta }(k_y)|^2}}{|\mathrm{\Delta }(k_y)|}},|eV||\mathrm{\Delta }|,`$ (3)
$`\mathrm{\Gamma }`$ $`=`$ $`{\displaystyle \frac{eVi\sqrt{|\mathrm{\Delta }(k_y)|^2(eV)^2}}{|\mathrm{\Delta }(k_y)|}},|eV||\mathrm{\Delta }|.`$ (4)
Denoting $`G(V,k_y)`$ as $`𝒢(V,k_y)`$ for the subgap voltage $`|eV||\mathrm{\Delta }|`$, where $`|\mathrm{\Gamma }(V)|=1`$, we rewrite Eq. (2) as

$`𝒢(V,k_y)={\displaystyle \frac{D^2(k_y)/2[1D(k_y)]}{D^2(k_y)/4[1D(k_y)]+F^2(V,k_y)}},`$ (5)
$`F=\sqrt{1|eV/\mathrm{\Delta }|^2}\mathrm{cos}(\mathrm{\Phi }/2)(eV/|\mathrm{\Delta }|)\mathrm{sin}(\mathrm{\Phi }/2).`$ (6)
For a high barrier, $`Dv_nv_s/^21`$. Then Eq. (5) becomes $`𝒢(V,k_y)\pi D(k_y)\delta [F(V,k_y)]`$, and Eq. (1) gives

$$\overline{𝒢}(V)C_{k_y^{\mathrm{max}}}^{k_y^{\mathrm{max}}}𝑑k_yv_n(k_y)v_s(k_y)\delta [F(V,k_y)],$$
(7)
where $`C=2\pi e^2l/h^2^2`$. The delta-function contributes to the integral (7) when $`F(V,k_y)=0`$. Using Eq. (6), this condition can be written as $`eV=E(k_y)`$, where

$$E(k_y)=|\mathrm{\Delta }(k_y)|\mathrm{cos}[\mathrm{\Phi }(k_y)/2]$$
(8)
for $`\mathrm{sin}(\mathrm{\Phi }/2)0`$, i.e. $`0\mathrm{\Phi }2\pi `$, which is the appropriate interval for $`\mathrm{\Phi }`$. $`E(k_y)`$ (8) is nothing but the energy of a surface Andreev state with the momentum $`k_y`$ obtained for the impenetrable barrier ($`\mathrm{}`$) . The energy density of these surface states (DOS) is

$$\rho (ϵ)=\underset{k_F^{(s)}}{\overset{k_F^{(s)}}{}}\frac{dk_y}{h}\delta [ϵE(k_y)]=\underset{j}{}\frac{1}{h|_{k_y^{(j)}}E(k_y^{(j)})|},$$
(9)
where $`k_y^{(j)}`$ is the $`j`$-th root of the equation $`E(k_y)=ϵ`$.

The delta-function in Eq. (7) indicates that the subgap tunneling takes place when the bias voltage matches the energy of a surface state: $`eV=E(k_y)`$. Denoting the $`j`$-th root of this equation as $`k_y^{(j)}`$ and resolving the delta-function in Eq. (7), we find

$$\overline{𝒢}(V)C\underset{j}{}\frac{v_n(k_y^{(j)})v_s(k_y^{(j)})\sqrt{|\mathrm{\Delta }(k_y^{(j)})|^2(eV)^2}}{|_{k_y^{(j)}}E(k_y^{(j)})|}.$$
(10)
Because the denominators in Eqs. (10) and (9) are the same, both DOS and tunneling conductance exhibit peaks (the inverse-square-root van Hove singularities) at the turning points of the energy dispersion, where $`_{k_y}E(k_y)=0`$. The peak positions depend solely on the pairing potential $`\mathrm{\Delta }(k_y)`$ through Eq. (8) and not on the band structure details. The number of such turning points does not change upon small continuous deformation of $`\mathrm{\Delta }(k_y)`$, so it is a topological feature. However, because of the additional factors in Eq. (10), tunneling conductance is not simply proportional to DOS, as often assumed. Eq. (10) vanishes at $`k_yk_F^{(s)}`$, where the normal component $`v_s`$ of the Fermi velocity goes to zero. It also vanishes when $`\mathrm{\Delta }(k_y)0`$ or $`E(k_y)\pm \mathrm{\Delta }(k_y)`$ (i.e. $`\mathrm{\Phi }(k_y)0,\mathrm{\hspace{0.25em}2}\pi `$). In these cases, the quasiparticle localization length perpendicular to the interface, $`\lambda =\mathrm{}v_s/\sqrt{|\mathrm{\Delta }|^2E^2}`$, diverges ; thus, the probability to find a quasiparticle on the surface vanishes.

## III Application to $`\mathrm{𝐒𝐫}_\mathrm{𝟐}\mathrm{𝐑𝐮𝐎}_\mathrm{𝟒}`$

In Figs. 2, 3, 4, and 5, we show how this general formalism applies to the alternative superconducting order parameters proposed for $`\mathrm{Sr}_2\mathrm{RuO}_4`$ in the literature. In order to focus on the features of superconductor, rather than normal metal, we took $`k_F^{(n)}=5k_F^{(s)}`$. Then, the perpendicular velocity of the normal metal is approximately constant for all momenta $`|k_y|k_F^{(s)}`$: $`v_nv_F^{(n)}`$, while in the superconductor it is $`v_s=v_F^{(s)}\sqrt{1(k_y/k_F^{(s)})^2}`$, where $`v_F^{(n)}`$ and $`v_F^{(s)}`$ are the corresponding Fermi velocities. The numerical calculations are done for the high barrier $`Z=/\sqrt{v_F^{(n)}v_F^{(s)}}=3`$. As the figures illustrate, in this limit, the approximate curves given by Eq. (10), which includes only the contribution of the localized surface states, agree well with the exact curves obtained from Eqs. (1) and (2), which also include the contribution of the extended bulk states.

### A Nodeless pairing potential

For the 2D isotropic chiral $`p`$-wave, the pairing potential and the energy dispersion of the edge states are

$$\mathrm{\Delta }_1(𝐤)=\mathrm{\Delta }_0(k_x+ik_y)/k_F^{(s)},E_1(k_y)=\mathrm{\Delta }_0k_y/k_F^{(s)}.$$
(11)
Substituting Eq. (11) into Eqs. (9) and (10), we find

$`\rho (ϵ)=k_F^{(s)}/h\mathrm{\Delta }_0\rho _0,`$ (12)
$`\overline{𝒢}(V)2G_0[1(eV/\mathrm{\Delta }_0)^2],`$ (13)
where $`G_0=\pi e^2lk_F^{(s)}v_F^{(n)}v_F^{(s)}/h^2^2`$ is the normal-state conductance. The energy dispersion $`E_1(k_y)`$ (11) has no turning points, and the corresponding DOS (12) is flat: $`\rho (ϵ)=\mathrm{const}`$. The tunneling conductance curve $`\overline{G}(V)`$ (13) is parabolic for $`|eV|<\mathrm{\Delta }_0`$, attaining the maximal value $`2G_0`$ at $`V=0`$ and vanishing at $`|eV|=\mathrm{\Delta }_0`$, as shown in Fig. 2. For $`|eV|>\mathrm{\Delta }_0`$, $`\overline{G}(V)`$ increases with increasing $`|V|`$ and saturates at the normal-state value $`G_0`$ at $`|eV|\mathrm{\Delta }_0`$. This behavior is in agreement with previous numerical calculations . One can recognize that the curves shown in Fig. 4 of Ref.  and the corresponding figures of Ref.  for $`|eV|<\mathrm{\Delta }_0`$ are actually distorted parabolas. The deviations are due to the finite values of $`Z`$ and to the finite size of the acceptance cone utilized in the numerical calculations . Our analytical expression (13) applies in the case $`Z1`$ and for the full acceptance cone.

### B Pairing potentials with vertical lines of nodes

In this subsection, we consider pairing potentials with vertical lines of nodes, which originate from modulation of the gap function $`\mathrm{\Delta }(𝐤)`$ in the $`(k_x,k_y)`$ plane.

For the anisotropic chiral $`p`$-wave proposed in Ref. , the pairing potential and the energy dispersion are

$`\mathrm{\Delta }_2(𝐤)`$ $`=`$ $`\mathrm{\Delta }_0[\mathrm{sin}(k_xa/\mathrm{})+i\mathrm{sin}(k_ya/\mathrm{})],`$ (14)
$`E_2(k_y)`$ $`=`$ $`\mathrm{\Delta }_0\mathrm{sin}(k_ya/\mathrm{}),`$ (15)
where $`a`$ is the lattice constant, and $`k_F^{(s)}a=0.9\pi \mathrm{}`$. Substituting Eqs. (14) and (15) into Eqs. (1), (2), (4), (10) and (9), we calculate the tunneling conductance and DOS curves shown in Fig. 3. They exhibit peaks at $`eV=\pm \mathrm{\Delta }_0`$ originating from the turning points in the edge states dispersion. Equation $`eV=E(k_y)`$ has one solution $`k_y^{(j)}`$ for $`|eV|<0.31\mathrm{\Delta }_0`$ and two solutions for $`0.31\mathrm{\Delta }_0<|eV|<\mathrm{\Delta }_0`$. The switch causes the discontinuity of DOS and the slope change of tunneling conductance at $`eV=\pm 0.31\mathrm{\Delta }_0`$.

For the $`f`$-wave proposed in Ref. , we have

$`\mathrm{\Delta }_3(𝐤)`$ $`=`$ $`\mathrm{\Delta }_0(k_x^2k_y^2)(k_x+ik_y)/(k_F^{(s)})^3,`$ (16)
$`E_3(k_y)`$ $`=`$ $`\mathrm{\Delta }_0|k_x^2k_y^2|k_y/(k_F^{(s)})^3.`$ (17)
As shown in Fig. 4, the edge states dispersion has cusps at $`k_y=\pm k_F^{(s)}/\sqrt{2}`$, where $`\mathrm{\Delta }_3(𝐤)`$ vanishes, and two turning points in between. Consequently, the tunneling conductance and DOS curves have peaks at $`eV=\pm 0.26\mathrm{\Delta }_0`$.

For another $`f`$-wave proposed in Refs. , we have

$`\mathrm{\Delta }_4(𝐤)`$ $`=`$ $`2\mathrm{\Delta }_0k_xk_y(k_x+ik_y)/(k_F^{(s)})^3,`$ (18)
$`E_4(k_y)`$ $`=`$ $`2\mathrm{\Delta }_0k_x^2k_y/(k_F^{(s)})^3.`$ (19)
Eqs. (18) and (19) transform into Eqs. (16) and (17) upon a $`\pi /4`$ rotation in the $`x`$-$`y`$ plane. As shown in Fig. 5, the turning points in the energy dispersion and the peaks in the tunneling conductance and DOS curves occur at $`eV=\pm 0.76\mathrm{\Delta }_0`$. Their position is approximately the same as in Fig. 3 for the anisotropic $`p`$-wave. However, in Fig. 5, unlike in all other figures, tunneling conductance vanishes at zero bias, even though DOS remains finite. This characteristic feature of the $`f`$-wave potential (18) is a consequence of the square-root factor in Eq. (10).

Thus, tunneling conductance $`\overline{G}(V)`$ exhibits two peaks located symmetrically around $`V=0`$ for all pairing potentials with vertical lines of nodes described in this subsection. This conclusion is in agreement with recent numerical calculations .

### C Pairing potentials with horizontal lines of nodes

Thus far we considered pairing potentials that do not depend on $`k_z`$, the momentum component perpendicular to the $`\mathrm{Sr}_2\mathrm{RuO}_4`$ planes. Ref.  also proposed the $`p`$-wave with horizontal lines of nodes at certain $`k_z`$. Its pairing potential can be obtained from Eq. (11) by the following substitution

$$\mathrm{\Delta }_0\stackrel{~}{\mathrm{\Delta }}_0(p_z)=\mathrm{\Delta }_0\mathrm{cos}(p_z),$$
(20)
where we switched to the dimensionless variable $`p_z=k_zc/\mathrm{}`$ ($`c`$ is the interplane distance).

Now DOS is obtained by averaging $`\rho (ϵ)`$ from Eq. (12) over $`p_z`$ using Eq. (20):

$`\rho (ϵ)`$ $`=`$ $`{\displaystyle \underset{|\stackrel{~}{\mathrm{\Delta }}_0(p_z)|>|ϵ|}{}}{\displaystyle \frac{dp_z}{2\pi }}{\displaystyle \frac{k_F}{h|\stackrel{~}{\mathrm{\Delta }}_0(p_z)|}}`$ (21)
$`=`$ $`{\displaystyle \frac{\rho _0}{\pi }}\mathrm{ln}{\displaystyle \frac{1+\sqrt{1(ϵ/\mathrm{\Delta }_0)^2}}{1\sqrt{1(ϵ/\mathrm{\Delta }_0)^2}}}.`$ (22)
In Eq. (21), $`\stackrel{~}{\mathrm{\Delta }}_0(p_z)`$ vanishes linearly with $`p_z`$ at the nodes, which results in a logarithmic divergence of $`\rho (ϵ)`$ at $`ϵ0`$:

$$\rho (ϵ)\mathrm{ln}|\mathrm{\Delta }_0/ϵ|,\mathrm{for}|ϵ|\mathrm{\Delta }_0.$$
(23)
The logarithmic divergence is cut off at $`|\stackrel{~}{\mathrm{\Delta }}_0(p_z)|=|ϵ|`$, because Eq. (12) applies only to the subgap states with $`|ϵ|<|\stackrel{~}{\mathrm{\Delta }}_0(p_z)|`$. The plot of $`\rho (ϵ)`$ is shown in Fig. 2.

Tunneling conductance $`\overline{𝒢}(V)`$ is obtained by substituting Eq. (20) into Eq. (13) and averaging over $`p_z`$:

$`\overline{𝒢}(V)`$ $`=`$ $`2G_0{\displaystyle \underset{|\stackrel{~}{\mathrm{\Delta }}_0(p_z)|>|eV|}{}}{\displaystyle \frac{dp_z}{2\pi }}\left(1\left|{\displaystyle \frac{eV}{\stackrel{~}{\mathrm{\Delta }}_0(p_z)}}\right|^2\right)`$ (24)
$`=`$ $`{\displaystyle \frac{4}{\pi }}G_0\left(\mathrm{arccos}\left|{\displaystyle \frac{eV}{\mathrm{\Delta }_0}}\right|\left|{\displaystyle \frac{eV}{\mathrm{\Delta }_0}}\right|\sqrt{1\left|{\displaystyle \frac{eV}{\mathrm{\Delta }_0}}\right|^2}\right).`$ (25)
In Eq. (24), the divergence of $`1/|\stackrel{~}{\mathrm{\Delta }}_0(p_z)|^2`$ at the node is cut off at $`|eV|`$, because Eq. (13) applies only for $`|eV|<|\stackrel{~}{\mathrm{\Delta }}_0(p_z)|`$. To treat the divergence, it is convenient to change the variable of integration from $`p_z`$ to $`\mathrm{\Delta }`$:

$$\overline{𝒢}(V)=2G_0\left(1\frac{4}{2\pi }_0^{|eV|}\frac{d\mathrm{\Delta }}{\mathrm{\Delta }^{}}\frac{4|eV|^2}{2\pi }_{|eV|}^{\mathrm{\Delta }_0}\frac{d\mathrm{\Delta }}{\mathrm{\Delta }^{}\mathrm{\Delta }^2}\right).$$
(26)
Here the factors of 4 come from the two nodal lines and the integration over the two sides of each line. The slope $`\mathrm{\Delta }^{}=|d\stackrel{~}{\mathrm{\Delta }}_0(p_z)/dp_z|`$ is approximately constant at the nodes: $`\mathrm{\Delta }^{}\mathrm{\Delta }_0`$. Then the integrals in Eq. (26) can be taken, and we find

$$\overline{𝒢}(V)2G_0\left(1\frac{4}{\pi }\left|\frac{eV}{\mathrm{\Delta }_0}\right|\right)\mathrm{for}|eV|\mathrm{\Delta }_0.$$
(27)
Thus, tunneling conductance has a cusp at zero bias. The curve for $`\overline{G}(V)`$ calculated exactly starting from Eqs. (1) and (2) is shown in Fig. 2 and indeed demonstrates the cusp in agreement with Eq. (27).

Notice the difference in Fig. 2 between the singular logarithmical divergence of DOS and the triangular-shaped cusp in tunneling conductance at zero bias. Both features result from the pile-up of surface Andreev states at zero energy caused by vanishing gap at the nodes. These are robust features of the chiral pairing potentials with horizontal lines of nodes, independent of the band structure details. While Eqs. (22) and (25) are derived specifically for the gap (20) being a cosine function, the asymptotic Eqs. (23) and (27) are valid for any generic chiral pairing potential with horizontal lines of nodes. For example, the $`f`$-wave with horizontal lines of nodes ,

$`\mathrm{\Delta }_6(𝐤)`$ $`=`$ $`i\mathrm{\Delta }_0\mathrm{sin}(p_z)(k_x+ik_y)^2/(k_F^{(s)})^2,`$ (28)
$`E_6(k_y)`$ $`=`$ $`\mathrm{sgn}(k_y)\mathrm{\Delta }_0[2(k_y/k_F^{(s)})^21]|\mathrm{sin}(p_z)|,`$ (29)
produces similar curves for the in-plane tunneling. However, being an odd function of $`p_z`$, it also exhibits a ZBCP in the $`c`$-axis tunneling.

## IV Comparison with experiment

Point-contact tunneling spectroscopy  found a ZBCP in $`\mathrm{Sr}_2\mathrm{RuO}_4`$. The data were fitted using the 2D isotropic $`p`$-wave (11) and assuming a narrow acceptance cone, i.e. a small value of $`k_y^{\mathrm{max}}`$. In this situation, only the states with $`k_y0`$ and hence $`E(k_y)0`$ are probed; thus all pairing potentials, except the $`f`$-wave (18), would show a ZBCP. However, such a ZBCP, in contrast to nonchiral cuprates  and organic superconductors , does not identify the pairing symmetry uniquely.

The more recent tunneling experiments  were performed on cleaved junctions of $`\mathrm{Sr}_2\mathrm{RuO}_4`$ with inclusions of Ru. Although nominally their setup corresponds to tunneling along the $`c`$ axis, the experimentalists believe that the actual tunneling takes place in the $`(a,b)`$ plane via the inclusions of Ru. They found a cusp-like peak at zero bias in the $`1.4`$-K phase and a sharper ZBCP in the $`3`$-K phase. The latter was attributed to a nonchiral superconductivity developing at the grain boundaries . An experimental curve from Ref.  representing the $`1.4`$-K phase is shown in Fig. 6. It has the characteristic triangular shape with a cusp at zero bias and agrees very well with our theoretical curve calculated for the pairing potential (20). To making the comparison, we treated the horizontal and vertical scales and the vertical offset (the background) of the theoretical curve as fitting parameters. The best fit corresponds to the gap $`\mathrm{\Delta }_0=0.8`$ meV.

The good agreement between the theoretical and experimental curves is a strong indication in favor of the $`p`$\- or $`f`$-wave pairing potential with horizontal lines of nodes in the 1.4-K phase. Further measurements showing the absence or presence of a ZBCP in the $`c`$-axis tunneling could discriminate between the $`p`$\- and $`f`$-wave cases. Experiment  found that an additional sharp ZBCP with the width of the order of 0.1 meV develops on top of the cusp at lower temperatures 0.5 and 0.32 K. It may be attributed to a contribution from the $`c`$-axis tunneling, if the pairing potential is the $`f`$-wave (28).

The described above scenario corresponds to the horizontal nodes in the main electron band $`\gamma `$ of $`\mathrm{Sr}_2\mathrm{RuO}_4`$. In the alternative scenario proposed in Ref. , the main band $`\gamma `$ is nodeless, but the two other bands $`\alpha `$ and $`\beta `$ have horizontal lines of nodes. In this case, one would expect a superposition of a parabolic curve for the nodeless band and triangular curves for the bands with nodes. Perhaps the sharp ZBCP developing at low temperatures 0.5 and 0.32 K could be interpreted as the triangular curve corresponding to the gap of the order of 0.1 meV in the $`\alpha `$ and $`\beta `$ bands.

## V Conclusions

We derived the analytical formula (10) for subgap tunneling at low transparency of the barrier, which takes into account only the contribution of the surface Andreev states. This formula produces tunneling curves in a simple and physically transparent way. For the chiral pairing potentials with vertical lines of nodes, tunneling curves show double peaks, which originate from the turning points in the energy dispersion of the surface Andreev states. On the other hand, for the chiral pairing potentials with horizontal lines of nodes, we find a single triangular-shaped peak with a cusp at zero bias, which results from the pile-up of the surface Andreev state at zero energy caused by vanishing gap at the nodes.

Double peaks in tunneling conductance were observed experimentally in some point contacts with $`\mathrm{Sr}_2\mathrm{RuO}_4`$ in Ref. . However, such double-peaked spectra were associated with contacts with high transparency of the barrier , a regime that is not addressed in our paper. On the other hand, no double peaks were observed in experiment , which only found a single triangular-shaped peak at zero bias. As shown in Fig. 6, this peak is very well fitted by our calculations for horizontal lines of nodes. It is also visually similar to the single peak found for contacts with low transparency in Fig. 3 of experiment  (which shows resistance, rather than conductance). Thus, we conclude that the superconducting pairing potential in $`\mathrm{Sr}_2\mathrm{RuO}_4`$ most likely has horizontal lines of nodes.

*Note added in proof.* Quantitative comparison shows that our calculation for horizontal lines of nodes also fits the experimental data of Ref. 29 very well. The curves are shown at
http://www2.physics.umd.edu/~yakovenk/talks/Sr2RuO4/

We are grateful to Ying Liu and Gernot Goll for sending us the experimental data of Refs. 30 and 29, and to Yoshiteru Maeno, Ying Liu, and Igor Z̆utić for useful discussions. This work was supported by the Packard Foundation and by the NSF Grant DMR-9815094.


format: text

text:
# Ordering of the three-dimensional Heisenberg spin glass in magnetic fields

## Abstract

Spin and chirality orderings of the three-dimensional Heisenberg spin glass are studied under magnetic fields in light of the recently developed spin-chirality decoupling-recoupling scenario. It is found by Monte Carlo simulations that the chiral-glass transition and the chiral-glass ordered state, which are essentially of the same character as their zero-field counterparts, occur under magnetic fields. Implication to experimental phase diagram is discussed.

In the studies of spin glasses, much effort has been devoted to the properties under magnetic fields. Unfortunately, our understanding of them still has remained unsatisfactory. Most of the numerical studies have focused on the properties of the simple Ising model, especially the three-dimensional (3D) Edwards-Anderson (EA) model. While the existence of a true thermodynamic spin-glass (SG) transition has been established for this model in zero field, the question of its existence or nonexistence in magnetic fields has remained unsettled.

If one tries to understand real experimental SG ordering, one has to remember that many of real SG materials are more or less Heisenberg-like rather than Ising, in the sense that the random magnetic anisotropy is considerably weaker than the isotropic exchange interaction. Numerical simulations have indicated that the isotropic 3D Heisenberg SG with finite-range interaction does not exhibit the conventional SG order at finite temperature in zero field. Since applied fields generally tend to suppress the SG ordering, a true thermodynamic SG transition is even more unlikely under magnetic fields in case of the Heisenberg model.

Experimentally, however, a rather sharp transition-like behavior has been observed under magnetic fields in typical Heisenberg-like SG magnets, e.g., canonical SG like AuFe and CuMn, although it is not completely clear whether the observed anomaly corresponds to a true thermodynamic transition. The situation is in contrast to the zero-field case where the existence of a true thermodynamic SG transition has been established experimentally. Set aside the question of the strict nature of the SG “transition”, it is experimentally observed that a weak applied field lowers the zero-field SG transition temperature rather quickly . For higher fields, the SG “transition” becomes much more robust to fields, where the “transition temperature” shows much less field dependence . Such behaviors of the SG transition temperature under magnetic fields $`T_g(H)`$ were often interpreted in terms of the mean-field model . Indeed, the mean-field Sherrington-Kirkpatrick (SK) model with an infinite-range Heisenberg exchange interaction with weak random magnetic anisotropy exhibits a transition line similar to the experimental one : i.e., the so-called de Almeida-Thouless (AT) line $`H(T_g(0)T_g(H))^{3/2}`$ in weak-field regime where the anisotropy is important, and the Gabay-Toulouse (GT) line $`H(T_g(0)T_g(H))^{1/2}`$ in strong-field regime where the anisotropy is unimportant. Nevertheless, if one notes that the true finite-temperature transition under magnetic fields, though possible in the infinite-range SK model, is unlikely to occur in a more realistic finite-range Heisenberg model, an apparent success of the mean-field model in explaining the experimental phase diagram should be taken with strong reservation.

Recently, one of the present authors has proposed a scenario, the spin-chirality decoupling-recoupling scenario, aimed at explaining some of the puzzles concerning the experimentally observed SG transition in zero field . In this scenario, chirality, which is a multispin variable representing the sense or the handedness of local noncoplanar spin structures induced by spin frustration, plays an essential role. In a fully isotropic Heisenberg SG, in particular, this scenario claims the occurrence of a novel chiral-glass ordered state in which only the chirality exhibits a glassy long-range order (LRO) while the spin remains paramagnetic. At the chiral-glass transition, among the global symmetries of the Hamiltonian, only the $`Z_2`$ spin reflection (inversion) symmetry is broken spontaneously with keeping the $`SO(3)`$ spin rotation symmetry preserved. Note that this picture entails the spin-chirality (or $`SO(3)Z_2`$) decoupling on long length and time scales: Namely, although the chirality is not independent of the spin on microscopic length scale, it eventually exhibits a long-distance behavior entirely different from the spin. Such a chiral-glass transition was indeed observed in zero field in a recent Monte Carlo simulation by Hukushima and Kawamura . It was also found there that the critical properties associated with the chiral-glass transition were different from those of the Ising SG, and that the chiral-glass ordered state exhibited a one-step-like novel RSB.

In the chirality scenario of Ref., experimental SG transition in real Heisneberg-like SG magnets is regarded essentially as a chiral-glass transition “revealed” via the magnetic anisotropy. Weak but finite random magnetic anisotropy inherent to real magnets “recouples” the spin to the chirality, and the chiral-glass transition shows up as an experimentally observable spin-glass transition.

The purpose of the present Letter is to reexamine the SG ordering of the 3D isotropic Heisenberg SG under applied fields in light of the above chirality scenario. We first argue some of the possible consequences of the chirality picture on the finite-field SG properties, and then perform extensive Monte Carlo (MC) simulations to check how the spin and chirality really order in fields. We show that the chiral-glass phase, essentially of the same character as the zero-field one, remains stable as a true thermodynamic phase in applied fields. It is also found that the associated chiral-glass transition line possesses some of the character of the GT line of the mean-field model, yet its physical origin entirely different.

Though we expect that our argument holds quite generally, we fix here our model Hamiltonian, which is just the one used in our MC simulation below. We consider the isotropic classical Heisenberg model on a 3D simple cubic lattice,

$$=\underset{<ij>}{}J_{ij}\stackrel{}{S}_i\stackrel{}{S}_jH\underset{i}{}S_i^z,$$
(1)
where $`H`$ is the intensity of the magnetic field applied along the $`z`$ direction. The nearest-neighbor coupling $`J_{ij}`$ is assumed to take the value $`J`$ or $`J`$ with equal probability ($`\pm J`$ distribution). Local scalar chirality is defined for three neighboring spins: Here we define it at the $`i`$-th site and in the $`\mu `$-th direction ($`\mu =x,y,z`$) by $`\chi _{i\mu }=\stackrel{}{S}_{i\widehat{e}_\mu }\stackrel{}{S}_i\times \stackrel{}{S}_{i+\widehat{e}_\mu }`$, $`\widehat{e}_\mu `$ being a unit lattice vector in the $`\mu `$-th direction.

Applied fields reduce the global symmetry of the Hamiltonian (1) from the zero-field one $`Z_2\times SO(3)`$ to $`Z_2\times SO(2)`$, where $`Z_2`$ refers to the spin reflection with respect to an arbitrary plane including the field axis, and $`SO(2)`$ to the spin rotation around the field axis. Note that, even in fields, the chiral $`Z_2`$ reflection symmetry, characterized by the mutually opposite signs of the chirality, is kept intact. The chirality is a pseudoscalar invariant under $`SO(2)`$ (or $`SO(3)`$ in case of $`H=0`$) rotations but changes sign under $`Z_2`$ reflections. The same chiral $`Z_2\times SO(2)`$ symmetry also appears in the $`XY`$ SG with two-component spins.

Since the chiral $`Z_2`$ is supposed to be decoupled from the spin $`SO(3)`$ already in zero field, and the applied field serves only to reduce the decoupled $`SO(3)`$ to $`SO(2)`$, one naturally expects that the $`Z_2`$ chiral-glass transition, essentially of the same type as the zero-field one, would persist even under magnetic fields . More specifically, the chiral-glass transition in a field should lie in the same universality class as the zero-field one, characterized by the same set of exponents, and the chiral-glass ordered state in a field exhibits the same kind of one-step-like RSB as in zero field.

Since the chiral-glass order in a field is expected to be essentially of the same character as in zero field, the chiral-glass transition temperature under magnetic fields $`T_{\mathrm{CG}}(H)`$ should be a regular function of $`H`$. If one takes account of the obvious symmetry $`HH`$, the chiral-glass transition line at low enough fields should behave as

$$T_{\mathrm{CG}}(0)T_{\mathrm{CG}}(H)=cH^2+c^{}H^4+\mathrm{}.$$
(2)
Generally, the coefficient $`c`$ could be either positive or negative. Interestingly, the above form of the chiral transition line is similar to the so-called GT line of the mean-field Heisenberg SK model. We emphasize, however, that their physical origin is entirely different. The quadratic dependence of the chiral-glass transition line is simply of regular origin, whereas that of the GT-line in the SK model cannot be regarded so.

Concerning the spin order, applied fields trivially induce a nonzero longitudinal (parallel to the field) spin order at any temperature. The behavior of the transverse (perpendicular to the field) component could be more nontrivial. If one recalls the recent numerical results on the 3D XY SG which indicate the absence of the standard SG order at least just below the chiral-glass transition , the transverse spin order is also unlikely to arise in the present case, at least just below the $`Z_2`$ chiral-glass transition.

In order to examine whether the above expectation really holds or not, we next perform extensive MC simulations on the isotropic $`\pm J`$ Heisenberg SG model (1). Simulations are performed for a variety of fields $`H/J=0.05,0.1,0.5,2.0,3.0,5.0`$. The lattices studied are simple-cubic lattices with $`L^3`$ sites with $`L=616`$ with periodic boundary conditions. The system is fully equilibrated with use of the temperature-exchange method. We show here the data explicitly for the particular field value of $`H/J=0.5`$, where the sample average is taken over 128-800 bond realizations.

By running two independent sequences of systems (replica 1 and 2) in parallel, we compute a scalar chiral overlap $`q_\chi `$ between the chiralities of the two replicas by $`q_\chi =\frac{1}{3N}_{i\mu }\chi _{i\mu }^{(1)}\chi _{i\mu }^{(2)}`$, as well as a spin-overlap tensor $`q_{\mu \nu }`$ between the transverse $`\mu `$ and $`\nu `$ components ($`\mu `$, $`\nu `$=$`x,y`$) of the spin by $`q_{\mu \nu }=\frac{1}{N}_iS_{i\mu }^{(1)}S_{i\nu }^{(2)}`$. Then, in terms of these overlaps, we calculate the Binder ratios of the chirality $`g_\chi `$, and of the transverse ($`XY`$) components of the spin $`g_T`$ defined in the standard manner: See Ref. for detailed definition. The results are shown in Fig.1. The Binder ratio of the chirality $`g_\chi `$ exhibits a negative dip which, with increasing $`L`$, tends to deepen and shift toward lower temperature. Furthermore, $`g_\chi `$ of various $`L`$ cross at a temperature slightly above the dip temperature $`T_{\mathrm{dip}}`$ on negative side of $`g_\chi `$, eventually merging at temperatures lower than $`T_{\mathrm{dip}}`$. The observed behavior of $`g_\chi `$ is similar to the one observed in zero field. As argued in Ref., the persistence of a negative dip and the crossing occurring at $`g_\chi <0`$, are strongly suggestive of the occurrence of a finite-temperature transition where $`g_\chi (T_{\mathrm{CG}}^{})`$ and $`g_\chi (T_{\mathrm{CG}})`$ take negative values in the $`L\mathrm{}`$ limit. In the inset of Fig.1(a), we plot the negative-dip temperature $`T_{\mathrm{dip}}(L)`$ versus $`1/L`$. The data lie on a straight line fairly well, and its extrapolation to $`1/L=0`$ gives an estimate of the bulk chiral-glass transition temperature, $`T_{\mathrm{CG}}/J0.25`$. (More precisely, $`T_{\mathrm{dip}}(L)`$ should scale with $`L^{1/\nu }`$ where $`\nu `$ is the chiral-glass correlation-length exponent. As shown below, our estimate of $`\nu 1.3`$ comes close to unity, more or less justifying the linear extrapolation employed here. Extrapolation with respect to $`L^{1/1.3}`$ yields $`T_{\mathrm{CG}}/J0.23`$.)

In sharp contrast to $`g_\chi `$, Binder ratio of the transverse component of the spin $`g_T`$ decreases monotonically toward zero with increasing $`L`$, without a negative dip nor a crossing, suggesting that the transverse component of spin remains disordered even below $`T_{\mathrm{CG}}`$.

Validity of our estiamte of $`T_{\mathrm{CG}}`$ has been checked also from the behavior of other quantities, e.g., the temporal-decay of the equilibrium chirality autocorrelation function (data not shown here), which yields $`T_{\mathrm{CG}}/J=0.23(2)`$ in accord with the above estimate. In order to probe the possible RSB in the chiral-glass ordered state, we display in Fig.2 the distribution function of the chiral-overlap defined by $`P(q_\chi ^{})=[\delta (q_\chi q_\chi ^{})]`$ where $`<\mathrm{}>`$ and \[$`\mathrm{}`$\] represent the thermal average and the sample average, respectively. The existence of a growing “central peak” at $`q_\chi =0`$ for larger $`L`$, in addition to the standard “side-peaks” corresponding to $`\pm q_{\mathrm{CG}}^{\mathrm{EA}}`$, suggests the occurrence of a one-step-like peculiar RSB in the chiral-glass ordered state. Similar behavior was observed in the chiral-glass state in zero field. The existence of a persistent negative dip in the Binder ratio $`g_\chi `$ is also consistent with the occurrence of such a one-step-like RSB.

With setting $`T_{\mathrm{CG}}/J=0.23`$ as determined above, we perform the standard finite-size scaling of the chiral-glass order parameter $`[<q_\chi ^2>]`$ and of the chiral autocorrelation function, to estimate various chiral-glass exponents. We then get $`\nu =1.3(2)`$, $`\eta =0.6(3)`$, $`z=5.3(5)`$, which turns out to agree within errors with the corresonding zero-field exponents of Ref. . The results seem consistent with a common universality class occurring both in zero-field and finite-field chiral-glass transitions.

Similar calculations and analysis are repeated for other field values as well. For $`H/J=0.1`$, in particular, we have performed the same scale of intensive calculation as was done for $`H/J=0.5`$, to find that all qualitative features are similar. The chiral-glass transition with a one-step-like RSB occurs at $`T_{\mathrm{CG}}/J=0.21(2)`$, with the exponents $`\nu =1.3(2)`$, $`\eta =0.6(3)`$, $`z=4.9(5)`$, which agrees within errors with our estimates for $`H/J=0.5`$.

By collecting the $`T_{\mathrm{CG}}`$ values determined for other field values, we construct a phase diagram in the temperature-magnetic field plane as shown in Fig.3. The chiral-glass state remains quite robust against magnetic fields. Indeed, $`T_{\mathrm{CG}}(H)`$ is not much reduced from the zero-field value even at field as large as ten times of $`T_{\mathrm{CG}}(0)`$. This somewhat surprising property probably arises from the fact that the magnetic field couples in the Hamiltonian directly to the spin, not to the chirality, and the effective coupling between the field and the chirality is rather weak. At lower fields, the chiral-glass transition line is almost orthogonal to the $`H=0`$ axis, consistent with the expected behavior Eq.(2). Our data are even not inconsistent with the coefficient $`c`$ in Eq.(2) being slightly negative so that $`T_{\mathrm{CG}}(H)`$ initially increases slightly with $`H`$, though it is difficult to draw a definite conclusion due to the scatter of our estimate of $`T_{\mathrm{CG}}(H)`$.

Finally, we wish to discuss implications of our results to real experimental SG. In real Heisenberg-like SG magnets, weak but finite random magnetic anisotropy neglected here plays a role. Among other things, the anisotropy recouples the spin to the chirality, transforming the chiral-glass state into the spin-glass state. Furthermore, in the presence of both random magnetic anisotropy and magnetic field, all global symmetries of the Hamiltonian will be lost. Nevertheless, we expect that the chiral-glass transition should still persist in fields as a pure RSB transition, not accompanying the global $`Z_2`$-symmetry breaking. Difference in the broken symmetries in zero- and in finite fields, however, causes a singular crossover behavior in fields, in sharp contrast to the fully isotropic case. We believe that this singular crossover line expected in the weak-field regime of anisotropic system is nothing but the AT-like transition line ubitiously observed experimentally. Further detailed nature of this crossover expected in the low-field regime will be discussed in a separate paper.

Meanwhile, in the strong-field regime where the applied field dominates the anisotropy, main features of the experimental phase diagram would be described by our present results on isotropic system. Thus, the experimental observation of the GT-like field-insensitive transition line in the high-field regime is fully consistent with our observation of the field-insensitive chiral-glass transition line. Other interesting possibility revealed by our analysis is that the SG transition line might extend to higher fields than hitherto suspected. For example, for canonical SG AuFe, Campbell et al recently determined the SG phase boundary in fields by torque measurements: At lower fields, applied fields rapidly suppress the SG transition giving rise to the standard AT-like behavior, whereas around the maximal fields of the measurements ($`7`$T) the SG transition line becomes almost-field independent. If one roughly estimates the effective $`J`$ of this system from its zero-field SG transition temperature, it is of order 50K. If one assumes that the energy scale of our present model calculation could roughly be applicable to AuFe, the SG ordered state should extend to fields much higher than $`10`$T without much reduction, or even with a slight increase, in $`T_g`$, although, considering the difference in microscopic details between the present model and real AuFe, one cannot expect a truely quantitative correspondence. Anyway, further high-field experiments on AuFe and on other Heisenberg-like SG magnets might be worthwhile to determine the SG phase boundary in the high-field regime.

The numerical calculation was performed on the Hitachi SR8000 at the supercomputer center, ISSP, University of Tokyo. The authors are thankful to Dr.K. Hukushima for useful discussion.


text:
# Suppression of the “quasiclassical” proximity gap in correlated-metal–superconductor structures

## Abstract

We study the energy and spatial dependence of the local density of states in a superconductor–correlated-metal–superconductor Josephson junction, where the correlated metal is a non-Fermi liquid (described by the Falicov-Kimball model). Many-body correlations are treated with dynamical mean-field theory, extended to inhomogeneous systems. While quasiclassical theories predict a minigap in the spectrum of a disordered Fermi liquid which is proximity-coupled within a mesoscopic junction, we find that increasing electron correlations destroy any minigap that might be opened in the absence of many-body correlations.

Fermi-liquid metals have excitation spectra that typically vary on energy scales of electron volts. Metals that become superconducting, have all low-energy electrons form pairs. Since it takes an energy equal to the superconducting gap $`\mathrm{\Delta }`$ to break a pair of electrons and form an excitation, there is a low-energy gap (on the order of meV) in the single-particle density of states (DOS). The original states at the Fermi level are “pushed” to excitations near $`\pm \mathrm{\Delta }`$, which yields a singularity in the DOS at zero temperature (and a large peak at finite $`T`$). When a superconductor ($`S`$) is connected to a normal metal ($`N`$) to form a $`SNS`$ Josephson junction, the superconductivity leaks into the normal metal via the proximity effect, and a weak link is established between the two $`S`$ through the $`N`$. What happens to the low-energy electrons in the $`N`$ is quite interesting. An electron near the Fermi level of the $`N`$ is localized within the $`N`$ because there are no single-particle states at low energy for it to scatter into within the $`S`$. Instead, the electron is retroreflected into a hole in the $`N`$, and creates a superconducting pair in the $`S`$ via a process called Andreev reflection. This reflection occurs at the $`SN`$ interface and the $`NS`$ interface, creating Andreev bound states with well-defined energy levels. These states are doubly degenerate, one carrying supercurrent to the right and one to the left. The states are broadened into peaks in the DOS when one averages over all different perpendicular momenta. Thus, Andreev reflection mixes the electron and hole states in the same proportion that they are mixed to form Bogoliubov quasiparticles in the $`S`$, with weights determined by the self-consistency condition. Such partially superconducting properties of the proximity-coupled normal metal are responsible for Josephson effect, as well as other peculiar phenomena in inhomogeneous systems which have been drawing increased attention over the past decade due to advances in mesoscopic superconductivity .

Another aspect of the proximity effect is the modification of the DOS in both the $`S`$ (“inverse proximity effect” ) and the $`N`$ side of a $`SN`$ boundary, which becomes most conspicuous in mesoscopic  confined geometries . For example, in a sufficiently long $`SNS`$ junction at low enough temperature, the proximity of the superconductor induces a minigap in the local density of states (LDOS) inside $`N`$ interlayer  that has chaotic classical dynamics . The minigap is of the order of the Thouless energy $`E_{\mathrm{Th}}=\mathrm{}/t_{\mathrm{dwell}}`$, where $`t_{\mathrm{dwell}}`$ is the typical time spent by an electron during its diffusive motion  ($`t_{\mathrm{dwell}}L^2/𝒟`$ for a $`N`$ strongly coupled to a $`S`$, where $`𝒟`$ is the diffusion constant) or during its chaotic ballistic motion  through the $`N`$ region of size $`L`$, before escaping into the superconductor (for integrable classical dynamics in the $`N`$, the LDOS is nonzero, but vanishes nearly linearly at the Fermi level ). These results were obtained using either quasiclassical approaches  or mean-field treatment by random matrix theory . However, recent calculations , which include additional quantum effects through a supersymmetric non-linear $`\sigma `$-model (NLSM) , show that mesoscopic fluctuations cause an exponentially small smearing of the quasiclassical gap in a diffusive $`SNS`$ junction (the DOS tails appearing below the quasiclassical gap edge are due to prelocalized states  which couple weakly to the $`S`$ leads). Using tunneling spectroscopy, the proximity affected LDOS can be measured as a function of the distance from the $`SN`$ interface in both the normal metal  and superconductor .

While both quasiclassical and NLSM calculations rely on the picture of well-defined (noninteracting) quasiparticles, little attention has been paid to proximity effects in systems where the $`N`$ layers are dominated by strong many-body correlations (except in one dimension where anomalously enhanced DOS have been found in a Luttinger liquid coupled to a superconductor ). Here we explore the LDOS, in both the superconducting and correlated metal sides of a superconductor–correlated-metal–superconductor ($`SCmS`$) Josephson junction. The

$`Cm`$ region is a non-Fermi liquid modeled by a Falicov-Kimball (FK) Hamiltonian . We find that increasing electron correlations completely destroy any initially open “quasiclassical” minigap in the DOS of the correlated metal, as shown in Fig. 1; this occurs due to the extensive broadening of the “Andreev bound states” by the scattering in the $`Cm`$. The appearance of a nonzero LDOS inside the $`S`$ region within a distance on the order of the superconducting coherence length $`\xi _S`$ is plotted in Fig. 2. Our analysis is fully self-consistent (i.e., we take into account the suppression of the superconducting order parameter inside the $`S`$ leads) at zero Josephson current across the junction. The calculation includes all (many-body) quantum effects encompassed in the

dynamical mean-field theory , which has only recently been generalized to treat inhomogeneous normal systems  and Josephson junctions .

Thus, our principal result is substantially different from the standard lore of a proximity-induced “hard minigap” (i.e., no states inside an energy interval $`E_gE_{\mathrm{Th}}`$), which is supported by both quasiclassical calculations  and the picture of bound states induced by Andreev reflection . Such qualitative considerations  give an estimate for $`E_g`$, which comes close to the values of the quasiclassical minigap obtained from the solution of the Usadel equation: $`E_g=0.78E_{\mathrm{Th}}`$ in the diffusive $`N`$ layer of an $`INS`$ structure  ($`I`$ is an insulator that specularly reflects), or $`E_g=3.12E_{\mathrm{Th}}`$ in a diffusive $`SNS`$ junction . The perturbative analysis of the electron-electron interaction, within quasiclassics, only generates a slightly smaller $`E_g`$ . The NLSM calculations  for a $`SNS`$ junction, with $`N`$ being a disordered noninteracting electron system, find nonzero DOS at all $`E<E_g`$, but the tail of subgap states is small for good metals (except near $`E_g`$). Our results are complementary to these theories, showing how the minigap can disappear at a critical correlation strength, and are of special interest in understanding the limitations of a phenomenological application of standard proximity-effect theory to experiments dealing with unconventional inhomogeneous structures, such as high-$`T_c`$ Josephson junctions with underdoped cuprates (a strongly correlated electron system) playing the role of the “normal region” .

The $`SCmS`$ Josephson junction is modeled by a Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{ij\sigma }{}}t_{ij}c_{i\sigma }^{}c_{j\sigma }+{\displaystyle \underset{i}{}}U_i\left(c_i^{}c_i{\displaystyle \frac{1}{2}}\right)\left(c_i^{}c_i{\displaystyle \frac{1}{2}}\right)`$ (1)
$`+`$ $`{\displaystyle \underset{i\sigma }{}}U_i^{FK}c_{i\sigma }^{}c_{i\sigma }\left(w_i{\displaystyle \frac{1}{2}}\right),`$ (2)
on an infinite set of stacked square lattice planes, whose connectivity is the same as a simple cubic lattice (with lattice constant $`a`$). Here $`c_{i\sigma }^{}`$ ($`c_{i\sigma }`$) creates (destroys) an electron of spin $`\sigma `$ at site $`i`$, $`t_{ij}=t`$ (the energy unit) is the hopping integral between nearest neighbor sites $`i`$ and $`j`$ (both within the planes and between planes), $`U_i<0`$ is the attractive Hubbard interaction for sites within the superconducting planes, $`U_i^{FK}`$ is the FK interaction for planes within the $`Cm`$ region, and $`w_i`$ is a classical variable that equals 1 if an $`A`$ ion occupies site $`i`$ and is zero if a $`B`$ ion occupies site $`i`$. The chemical potential $`\mu `$ is set equal to zero to yield half filling in the $`S`$ and $`Cm`$. The negative-$`U`$ Hubbard term describes the real-space pairing of electrons due to a local instantaneous attractive interaction. This generates a superconducting order in the $`S`$ leads which, when treated in the Hartree-Fock approximation, is equivalent to conventional BCS theory, except that here the DOS is non-constant and provides the energy cutoff.

The superconducting layers have $`U_i=2`$ and $`w_i=0`$ for all sites. Such a homogeneous bulk superconductor is characterized by the usual BCS parameters: the transition temperature $`T_c=0.11t`$, the zero-temperature order parameter $`\mathrm{\Delta }=0.198t`$, and the coherence length $`\xi _S=\mathrm{}v_F/(\pi \mathrm{\Delta })4a`$. The Cm interlayer is described by a half-filled FK model in the symmetric limit of half filling for the “ions” $`w_i=0.5`$. One can view the FK metal as a binary alloy of $`A`$ and $`B`$ ions at 50% concentration with $`U_{FK}`$ being the difference in site energy between the $`A`$ and $`B`$ ionic sites. The many-body problem is solved by taking an annealed average that yields the coherent potential approximation. This is the simplest many-body problem which, nevertheless, mimics a metal-insulator transition (MIT) of the type seen in the repulsive Hubbard model (except the metallic phase is not a Fermi liquid here). In the bulk, the FK correlated metal undergoes a MIT at $`U_{FK}4.9t`$ (which is close to half of the bandwidth $`6t`$). This is illustrated in Fig. 3 which shows the DOS in a bulk $`Cm`$ as a function of $`U_{FK}`$. The DOS is independent of temperature . Since the system is not a Fermi liquid for nonzero $`U_{FK}`$, the DOS first develops a pseudogap, and then is suppressed entirely to zero as the correlations increase, eventually driving the system into a correlated insulator. The opening of the gap is continuous. In order to focus only on the modification of the DOS induced by the proximity effect, we choose $`U_{\mathrm{FK}}2`$ for the strength of Coulomb interaction in the FK correlated metal, which ensures that the bulk DOS around the band center is essentially constant.

The problem of inhomogeneous superconductivity

is solved by employing the Nambu-Gorkov matrix formulation for Green functions with a local self-energy . We treat the problem self-consistently in the complex order parameter $`\mathrm{\Delta }_i`$ for the part of the junction comprised of the normal region (containing 5 or 10 planes) and 30 superconducting planes on each side of the $`Cm`$ interlayer. Inside this “self-consistent-part” of the infinite $`SCmS`$ junction, the superconducting gap $`\mathrm{\Delta }_i`$ heals to its bulk value since all signatures of the inverse proximity effect are gone on the length scale of few $`\xi _S`$ away from the $`SCm`$ interface. The calculation is performed at the temperature $`T=0.09T_c`$ where the BCS gap is fully developed. Details of our computational algorithm have been given elsewhere. The final result is the self-consistent Green function which allows us to compute the (many-body) LDOS as a function of plane position $`z_i`$ and frequency $`\omega `$ from the real-axis analytic continuation, $`\rho (\omega ,z_i)=\mathrm{Im}G(\omega +i\delta ;z_i,z_i)/\pi `$.

The strength of the superconductivity in the $`N`$ interlayer is quantified by the nonzero pair amplitude $`F(z_i)=\mathrm{\Delta }_i/|U_i|`$ (a two-particle property) which decays exponentially, due to the absence of an attractive interaction, on the length scale $`\xi _N=\sqrt{\mathrm{}𝒟/2\pi k_BT}`$  (or as a power law at zero temperature in a clean normal metal ). This is the length scale over which two thermal electrons with energy $`\omega \pi k_BT`$ in the $`N`$, correlated by Andreev reflection over the length scale $`L_\omega =\sqrt{\mathrm{}𝒟/2\omega }`$, lose their relative phase coherence (which then determines the thermodynamic critical Josephson current $`I_c`$). However, single-particle properties (like the DOS) , or kinetic quantities , can be influenced on a much longer length scale (where Josephson coupling vanishes) $`L_\omega \xi _N`$ for low-energy electrons $`\omega kT`$, which is ultimately limited by the mesoscopic phase-breaking length $`L_\varphi `$ . Therefore, the anomalies in the DOS for low $`\omega `$ extend up to the energy dependent distances $`L_\omega `$  from the $`SN`$ interface (which is smeared upon approaching $`L_\varphi `$ ). In a “closed” geometry, where a finite-size $`N`$ is disconnected from any electron reservoirs , this leads to a position independent minigap edge at $`E_gE_{\mathrm{Th}}`$ for sufficiently long $`E_{\mathrm{Th}}\mathrm{\Delta }`$ diffusive interlayer ($`\mathrm{}L`$, where $`\mathrm{}`$ is the mean free path). Thus, both the quasiclassical and NLSM descriptions of the proximity effect rely on the essential concepts in disordered Fermi-liquid physics, such as $`E_{\mathrm{Th}}`$ governing thermodynamic and quantum transport phenomena in mesoscopic systems. Since our $`Cm`$ layer is a non-Fermi liquid, these concepts are not directly transferable. Therefore, to compare our findings with standard notions, we proceed along a phenomenological route frequently (but unwarrantedly) employed in experiments on $`SNS`$ Josephson junctions with an unconventional $`N`$ interlayer . Namely, we extract an effective diffusion constant from the Kubo conductivity of the FK model , using the Einstein relation $`\sigma _{\mathrm{FK}}=2e^2N(0)𝒟`$ and the DOS at half-filling $`N(0)`$ from Fig. 3, and then compute $`E_{\mathrm{Th}}=\mathrm{}𝒟/L^2`$. It is interesting to check if such an energy scale provides any heuristic guidance in interpreting our results at small $`U_{\mathrm{FK}}`$ where the hard minigap is present in the $`Cm`$ spectrum. For example, in the largest $`U_{\mathrm{FK}}=2.0`$ sample, the resistivity of the FK correlated metal is $`\rho _{\mathrm{FK}}240`$ $`\mu \mathrm{\Omega }`$cm (assuming $`a=3`$ Å). From here we get $`𝒟2ta^2/\mathrm{}`$ and $`\xi _N5.6a`$. This is surprisingly close to the true $`\xi _N6.7a`$ extracted from the decay of $`I_c`$ in the $`SCmS`$ junction as a function of the $`Cm`$ layer thickness  (the agreement improves for smaller $`U_{FK}`$). The quasiclassical analysis for a mesoscopic diffusive junction of the same resistivity and with thickness $`L=10a`$ would give $`E_{\mathrm{Th}}0.1\mathrm{\Delta }`$ and $`E_g0.32\mathrm{\Delta }`$. However, at $`U_{\mathrm{FK}}=2.0`$ no gap is found in the $`Cm`$ spectrum, while only a small dip (Fig. 1) in the LDOS persists as a remnant of the minigap opened for $`U_{\mathrm{FK}}1.0`$. Moreover, the ratios of the minigap sizes $`E_g^a/E_g^b`$ in the $`SCmS`$ junctions with two different thicknesses $`L_a=5a`$ and $`L_b=10a`$ are: (i) $`U_{\mathrm{FK}}=0.1E_g^a/E_g^b=2.9`$, (ii) $`U_{\mathrm{FK}}=0.25E_g^a/E_g^b=3.1`$, (iii) $`U_{\mathrm{FK}}=0.5E_g^a/E_g^b=4.4`$, is a function of $`U_{\mathrm{FK}}`$. These are different from the expected results of $`E_g^a/E_g^b=(L_b/L_a)^2=4`$ or $`E_g^a/E_g^b=L_b/L_a=2`$, which would follow from $`E_gE_{\mathrm{Th}}`$ analogy with the quasiclassical description of the proximity effect in a diffusive Fermi-liquid metal of the same resistivity , or a clean but chaotic interlayer , respectively.

The last two pieces of information needed to characterize proximity-induced effects in correlated metals are: (1) the LDOS is position dependent, while $`E_g`$ is spatially constant (which is the same as the quasiclassical phenomenology ), and does not change upon lowering the temperature below our reference $`T/T_c=0.091`$ for Fig. 1; (2) the minigap is open for small enough $`U_{\mathrm{FK}}`$, where Andreev bound states also clearly coexist with it, as shown by the peaks below $`\mathrm{\Delta }`$ in the $`U_{\mathrm{FK}}=0.1`$ case in Fig. 1. In fact, we find the largest minigap in the limit $`U_{\mathrm{FK}}0`$ (corresponding to a clean $`SNS`$ junction), which is generated by the normal reflection at the $`SN`$ interface due to a non-negligible $`\mathrm{\Delta }/\mu 0.03`$ determining the amplitude of the scattering  ($`\mu `$ is the Fermi energy measured from the bottom of the band). It appears that increasing of electron correlations by increasing $`U_{\mathrm{FK}}`$ then just leads to a monotonic vanishing of any initially open minigap in the noninteracting case \[$`E_g(U_{\mathrm{FK}}=0)`$ depends on $`U_i`$, as illustrated in the inset of Fig. 1, and is outside of the quasiclassical approximation , but belongs to the realm of “noninteracting quasiparticle” physics\]. This should be contrasted with the quasiclassical minigap as a function of quenched disorder : the minigap forms for arbitrarily small concentration of impurities, increases with $`1/\mathrm{}`$ to a maximum value when $`\mathrm{}L`$, and then decreases in the diffusive limit as $`E_{\mathrm{Th}}v_F\mathrm{}`$. On the superconducting side of our $`SCmS`$ junction we find that position dependent LDOS is nonzero in the energy range $`2\mathrm{\Delta }`$ (Fig. 2), and decays to zero on the length scale of a few $`\xi _S`$ from the $`SCm`$ boundary .

Acknowledgments: Support from ONR grant number N00014-99-1-0328 is gratefully acknowledged. Computer calculations were partially supported by HPC time from the Arctic Region Supercomputer Center.


text:
# Determination of superconducting anisotropy from magnetization data on random powders as applied to LuNi2B2C, YNi2B2C and MgB2.

## Abstract

The recently discovered intermetallic superconductor MgB2 appears to have a highly anisotopic upper critical field with $`H_{c2}^{max}/H_{c2}^{min}\gamma >5`$. In order to determine the temperature dependence of both $`H_{c2}^{max}`$ and $`H_{c2}^{min}`$ we propose a method of extracting the superconducting anisotropy from the magnetization $`M(H,T)`$ of randomly oriented powder samples. The method is based on two features in $`(M/T)_H`$: the onset of diamagnetism at $`T_c^{max}`$, that is commonly associated with $`H_{c2}`$, and a kink in $`M/T`$ at a lower temperature $`T_c^{min}`$. Results for LuNi2B2C and YNi2B2C powders are in agreement with anisotropic $`H_{c2}`$ obtained from magneto-transport measurements on single crystals. Using this method on four different types of MgB2 powder samples we are able to determine $`H_{c2}^{max}(T)`$ and $`H_{c2}^{min}(T)`$ with $`\gamma 6`$.

Anisotropic type-II superconductors have been studied for a number of years with a revisited interest after the discovery of strongly anisotropic high-$`T_c`$ cuprates. The possibility has recently been raised that MgB2 ($`T_c40`$ K) may also be highly anisotropic. The quantitative characterisation of the anisotropy is usually one of the first questions to be addressed after the material is synthesized. One of the most accurate methods of extracting the anisotropy parameter $`\gamma =H_{c2}^{max}/H_{c2}^{min}`$ is by measuring the torque upon the single crystal sample in intermediate fields $`H_{c1}HH_{c2}`$ inclined relative to the crystal axes (we are interested here in uniaxial materials with large Ginzburg-Landau (GL) parameter $`\kappa =\lambda /\xi `$; $`\lambda `$ and $`\xi `$ are the penetration depth and the coherence length, $`H_{c1}`$ and $`H_{c2}`$ are the lower and the upper critical fields). The method is based on the existence of the transverse magnetization in anisotropic superconductors which causes the torque, the angular dependence of which depends only on $`\gamma `$.

Often, materials are first synthesized in a polycrystalline form. Anisotropy estimates of such samples can be inferred by a number of measurements. In one of the common schemes,  a powder of the material is mixed with a low magnetic background epoxy and aligned in a high magnetic field. The alignment is made permanent by curing the epoxy. The necessary requirements include a powder that consists of single crystalline grains with a considerable normal state magnetic anisotropy. The samples obtained in this way are suitable for anisotropic magnetization measurements. Usually, uncertainties in a degree of alignment lead to underestimates of $`\gamma `$.

We describe here a simple method to assess the anisotropy and to delineate $`H_{c2}^{max}(T)`$ and $`H_{c2}^{min}(T)`$ by analysis of the magnetization $`M(T)`$ taken for several values of the applied field $`H`$. We establish the validity of this method by comparison of anisotropic $`H_{c2}(T)`$ measured directly (by magneto-transport) on single crystals of quaternary borocarbide superconductors LuNi2B2C and YNi2B2C with the values extracted from the $`M(H,T)`$ data on randomly oriented powders. We then use the method to determine the curves of $`H_{c2}^{max}(T)`$ and $`H_{c2}^{min}(T)`$ for MgB2 and to obtain $`\gamma 6`$, close to the values given in Ref. .

Single crystals of LuNi2B2C and YNi2B2C were grown using the Ni2B flux growth technique. Four different samples of MgB2 were used in this work: (i) polycrystalline sintered Mg11B2, (ii) MgB2 wire segments, and commercial 99.5% purity MgB2 powders from two sources: (iii) Accumet Materials Co. and (iv) Alfa Aesar. Samples (i) and (ii) were synthesized by reacting elemental boron with Mg at $`950^{}`$ C.

The anisotropic $`H_{c2}(T)`$ of LuNi2B2C and YNi2B2C single crystals were measured directly by monitoring resistance $`R(H,T)`$ near the superconducting transition. For the magnetization measurements, single crystals of borocarbides and bulk samples of MgB2 were ground to powder, mixed with Epocast 121 epoxy (the epoxy has weak featureless diamagnetism in the $`HT`$ range of interest), stirred to assure random orientation of particles and cured at $`120^{}`$ C for 1 hour. DC magnetization measurements were performed in Quantum Design MPMS-5 or MPMS-7 SQUID magnetometers. Resistance of borocarbides were measured using LR-700 AC resistance bridge ($`f`$ = 16 Hz, $`I`$ = 1-3 mA) and the $`HT`$ environment of the MPMS instruments. Contacts to the samples were made with Epo-tek H20E silver epoxy. Currents were flowing in the $`ab`$ plane, whereas the field was applied perpendicular to the current, along $`c`$ axis or in the $`ab`$ plane.

The data on anisotropic $`H_{c2}`$ for non-magnetic borocarbides were reported in several publications. Since the $`H_{c2}`$ anisotropy appears to be somewhat sample-dependent, we performed direct magneto-transport measurements on samples from the same batches that were used for random powder magnetization measurements. Figure 1 shows representative resistance, $`R(H)`$, data for LuNi2B2C measured at constant temperatures for two field orientations. The resistively measured superconducting transition broadens slightly in field; this results in criterion-dependent $`H_{c2}(T)`$ curves. Still, the anisotropy of $`H_{c2}(T)`$ turns out practically criterion-independent. In the following we use the criterion of the maximum slope intersection with the $`R=0`$ line, as shown in Fig.1.

An example of magnetization data, $`M(T)`$, in a constant field, $`H=30`$ kG, for a powder sample of LuNi2B2C is shown in Fig. 2 along with the temperature derivative $`M/T`$. We interpret the behavior of $`M/T`$ as follows. Upon cooling a powdered sample in a fixed $`H`$, there is a deviation from a roughly $`T`$ independent normal state magnetization (for non-magnetic materials as is the case for the three compounds discussed here) to an increasingly diamagnetic signal at $`T=T_c^{max}`$. A second sharp feature will occur in $`M/T`$ when the sample temperature passes through $`T_c^{min}`$. This can be most clearly understood by considering what happens to the sample upon warming. For $`T<T_c^{min}`$ all grains are superconducting whereas for $`T>T_c^{min}`$ part of them will become normal, depending upon their orientation with respect to the applied magnetic field. Therefore, upon warming through $`T_c^{min}`$ there will be a kink in $`M/T`$ associated with the onset of normal state properties in an increasing number of appropriately oriented grains. In Fig. 2 this can be seen in the $`M/T`$ plot with $`T_c^{max}=10.6`$ K and $`T_c^{min}=8.7`$ K. The temperatures $`T_c^{min}`$ and $`T_c^{max}`$ are marked with vertical arrows in Fig. 2. It is worth noting that when a single $`H_{c2}(T)`$ value is obtained from measurements on polycrystalline samples it is actually the maximum upper critical field $`H_{c2}^{max}(T)`$.

The $`M(T)`$ data taken for a number of fixed fields were analyzed in a similar manner; the resulting curves $`H_{c2}^{min}(T)`$ and $`H_{c2}^{max}(T)`$ are presented in the upper panel of Fig. 4 (open symbols). These are plotted together with the directly measured anisotropic $`H_{c2}(T)`$ from magneto-transport measurements on LuNi2B2C single crystals from the same batch. It is seen that both sets of data are in agreement. We obtain $`\gamma 1.21.3`$. Albeit slightly higher, this value of $`\gamma `$ is consistent with values reported in literature: e.g. $`\gamma 1.2`$.  Similar measurements and analysis were performed on YNi2B2C. The results are shown in lower panel of Fig. 4. The slightly higher mismatch between the magneto-transport $`H_{c2}`$’s of single crystals and those extracted from powder magnetization measurements than for the case of LuNi2B2C might be due to arbitrary in-plane field orientation in magneto-transport measurements, choice of criteria and/or error bars in the measurements.

It is worth noting that the method of extracting the anisotropic $`H_{c2}(T)`$ and the anisotropy parameter $`\gamma `$ from the powder data presented here constitutes a robust procedure independent of a particular model for describing the anisotropy. Moreover, the method is just based on the existence of a kink in $`M/T`$ located at $`T=T_c^{min}`$, a feature which should be present for any angular distribution of the grains (as long as it is continuous, but not necessarily random). The analysis of the $`M(H,T)`$ data can be pushed further to to relate $`\gamma `$ to other material characteristics ($`\kappa `$, $`\lambda `$, etc).  This, however, would have involved more information (e.g., randomness, assumptions on linear dependences of $`M`$ on $`H_{c2}H`$ and $`T_cT`$) and we will not pursue this here.

With the validity of this method established for borocarbides, we can use it for MgB2. An example of magnetization curves for two values of $`H`$ and their temperature derivatives are shown in Fig. 5. As in borocarbides, a clear feature at $`T_c^{min}`$ is seen in $`M/T`$. With increasing $`H`$, this feature moves down in temperature faster than $`T_c^{max}`$ and disappears for $`H>25`$ kG.

$`H_{c2}^{min}(T)`$ and $`H_{c2}^{max}(T)`$ curves were deduced from magnetization data collected in fixed fields for four different samples of MgB2. These data are presented in Fig. 6. Data for the sintered powder and the powdered wire segments are very similar. The $`H_{c2}^{max}(T)`$ curves are consistent with the reported  polycrystalline magneto-transport $`H_{c2}(T)`$. The extracted anisotropy is $`\gamma 56`$. Both commercial MgB2 powders have $`T_c`$ values that are 1-1.5 K lower than sintered powder or wire segments (see Fig. 6), presumably due to a higher levels of impurities. The features in $`M/T`$ are less pronounced in commercial powders, however the resulting anisotropies are similar to that observed in clean MgB2 (Fig. 6). For these samples $`\gamma 67`$.

There have been some initial reports of the $`H_{c2}`$ anisotropy for MgB2. The anisotropy $`\gamma `$ for separate particles settled on a flat surface is 1.73.  For a hot-pressed bulk sample, $`\gamma 1.1`$.  Measurements on $`c`$-axis oriented thin films  gave $`\gamma =1.82`$, with the higher $`\gamma `$ for films with higher resistivity and lower $`T_c`$. Recently, there were several announcements on growth and anisotropic properties of sub-mm sized single crystals of MgB2.  In all three cases the residual resistance ratio was $`57`$ and the transition temperature was $`0.51`$ K lower than for polycrystalline materials.  The $`H_{c2}`$ anisotropy was found to be in the range $`2.63`$. A significantly higher $`H_{c2}`$ anisotropy ($`\gamma 69`$) was inferred from conduction electron spin resonance measurements and from fits to the reversible part of $`M(H)`$ on high purity and high residual resistance ratio samples . The results for MgB2 obtained in this work are consistent with the latter and are at the higher end of the wide-spread of values in the literature.

The method of extracting the superconducting anisotropy suggested in this paper works well for random powders of LuNi2B2C and YNi2B2C, for which the results were compared with direct measurements of $`H_{c2}`$. Moreover, the $`\gamma `$ values so obtained are in agreement with the microscopic theory according to which

$$\gamma ^2=\frac{\mathrm{\Delta }(𝐤_F)v_{ab}^2}{\mathrm{\Delta }(𝐤_F)v_c^2},$$
(1)
where $`v_i`$ are the Fermi velocities and $`\mathrm{}`$ stand for Fermi surface averages. For the isotropic gap function $`\mathrm{\Delta }(𝐤_F)=`$ constant, the band structure estimates give:

$$\gamma =\sqrt{v_{ab}^2/v_c^2}1.2.$$
(2)
This, however, is not the case for MgB2. The ratio $`v_{ab}^2/v_c^2`$ averaged over the whole Fermi surface for this material is close to unity, whereas our values of $`\gamma ^2`$ are in the range $`2550`$. This is indicative of a strong anisotropy of $`\mathrm{\Delta }(𝐤_F)`$. There are arguments  that the electron-phonon interaction is particularly strong on the Fermi surface sheets shaped as slightly distorted cylinders along the $`c`$ crystal direction. The ratio $`v_{ab}^2/v_c^2`$ averaged only over these cylinders is estimated as $`40`$. If the gap $`\mathrm{\Delta }`$ on the remaining Fermi surface sheets is negligible, we expect the superconducting anisotropy as $`\gamma \sqrt{40}6`$, the value close to that extracted from our analysis of the magnetization.

In conclusion, we suggest a simple method of evaluation of the anisotropy of the upper critical field from the analysis of the temperature dependent magnetization of randomly oriented powders. In the case of non-magnetic borocarbides the results are in good agreement with the direct measurements on single crystals and with the band-structure calculations provided the gap is isotropic. For MgB2 the estimated anisotropy is $`\gamma 57`$, which can be reconciled with the band calculations only if the superconducting gap on the cylindrically shaped sheets of the Fermi surface is dominant.

We thank K. D. Belashchenko and V. P. Antropov for discussions and for providing the $`v_{ab}^2/v_c^2`$ value for cylindrical parts of MgB2 Fermi surface. Ames Laboratory is operated for the U. S. Department of Energy by Iowa State University under contract No. W-7405-ENG-82. This work was supported by the Director of Energy Research, Office of Basic Energy Sciences.


text:
# Vortex nucleation in Bose-Einstein condensates in an oblate, purely magnetic potential

## Abstract

We have investigated the formation of vortices by rotating the purely magnetic potential confining a Bose-Einstein condensate. We modified the bias field of an axially symmetric TOP trap to create an elliptical potential that rotates in the radial plane. This enabled us to study the conditions for vortex nucleation over a wide range of eccentricities and rotation rates.

The existence of quantised vortices is one of the most striking and fascinating signatures of superfluidity. There has been theoretical interest in vortices in a dilute gas BEC for many years (for a review see for example ) and the recent experimental observations  have stimulated much further work. Vortices have been nucleated using phase imprinting  and by using a laser beam ’stirrer’. We have used a new, purely magnetic excitation scheme, which like the optical stirrer has a close analogy to the ’rotating bucket’ experiment with liquid He . Vortices provide dramatic visual evidence for the single condensate wavefunction, the detailed mechanisms for vortex formation, stabilisation and decay reveal new information about the dynamics of the excited states of the condensate and its interaction with the thermal cloud. The excited states and hence the conditions for nucleation depend on the geometry of the confining potential . Our work, in an oblate geometry compliments the work done to date in prolate and spherical traps. The range of eccentricity available in our purely magnetic rotating trap  is much greater than in previous experiments, where the radial symmetry was broken by ‘stirring’ the condensate with the dipole force from off-resonant laser beams.

We have measured the relationship between the eccentricity of the trap in the plane of rotation and the angular frequency $`\mathrm{\Omega }`$ that must be satisfied to nucleate vortices in our trap. The trap has a harmonic potential with oscillator frequencies $`\omega _x<\omega _y<\omega _z`$. The shape of the trap is characterised by the deformation parameter $`ϵ`$, given by

$$ϵ=\frac{\omega _y^2\omega _x^2}{\omega _x^2+\omega _y^2}$$
(1)
We show that when the trap is adiabatically ramped from circular, $`ϵ=0`$, to a given eccentricity, there is both an upper and lower critical rotation rate for nucleation. We present the threshold conditions for vortex nucleation in our geometry and finally discuss the role of the thermal cloud in both the nucleation and stabilization processes.

A detailed description of the rotating elliptical trap is given in  and only a summary of the technique is presented here. Our TOP trap consists of a spherical quadrupole field and a rapidly rotating bias field . The zero of the magnetic field describes a circle in the radial (XY) plane. The atoms experience a time-averaged, axially symmetric harmonic potential. If the X and Y components of the bias field ($`B_x`$ and $`B_y`$) are not of equal amplitude, then the axial symmetry is broken and the potential becomes elliptical in the XY plane. The ratio of the radial trap frequencies $`\omega _x/\omega _y`$, can be calculated numerically from $`B_y/B_x`$. For small eccentricities $`\omega _x/\omega _y1(B_y/B_x1)/4`$. To rotate the elliptical potential we modulated $`B_x`$ and $`B_y`$ at $`\mathrm{\Omega }`$, a frequency much lower than the TOP frequency $`\omega _0`$. The final bias field has the form of a rotation matrix through angle $`\mathrm{\Omega }t`$ multiplying the X and Y components of an ellipse:

$$\left(\begin{array}{c}B_x\\ B_y\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\mathrm{\Omega }t& \mathrm{sin}\mathrm{\Omega }t\\ \mathrm{sin}\mathrm{\Omega }t& \mathrm{cos}\mathrm{\Omega }t\end{array}\right)\left(\begin{array}{c}EB_t\mathrm{cos}\omega _0t\\ B_t\mathrm{sin}\omega _0t\end{array}\right)$$
(2)
During any particular evaporative cooling run the value of $`\mathrm{\Omega }`$, the angular frequency at which the elliptical trap rotates, was fixed (for technical reasons) and so the trap was always spinning. To create a ’static’ trap for efficient evaporative cooling, the bias field ratio, E, is set to 1 so that the trap is symmetric around the rotation axis. This symmetry was checked to be accurate to $`ϵ=0\pm 0.005`$ using measurements of the frequency of dipole oscillations in the X and Y directions. Evaporative cooling followed by an adiabatic expansion gave a condensate of $`2\times 10^4`$ atoms, at a temperature of $`0.5T/T_c`$. At this stage the trap is axially symmetric, $`\omega _x=\omega _y`$, with trap frequencies typically $`\omega _x/2\pi =62`$ Hz and $`\omega _z/2\pi =175`$ Hz.

To create vortices, the value of E was ramped linearly over 200ms from 1 to its final value, to give a trap that was elliptical in the rotating frame. If one creates an eccentric trap by increasing (decreasing) the TOP bias field then all three trap frequencies are reduced (increased). Since vortex lifetimes depend on the mean trap frequency , we adjusted the quadrupole field in the adiabatic expansion stage, to ensure that all traps had the same average radial trap frequency $`\omega _{}`$, defined as

$$\omega _{}=\sqrt{\frac{\omega _x^2+\omega _y^2}{2}}.$$
(3)
The condensate was then held in the rotating anisotropic trap  for a further 800 ms before being released. After 12ms of free expansion the cloud was imaged along the axis of rotation using an absorption imaging system with 3 $`\mu `$m resolution. Figure 1 shows images of the expanded condensate at different stages during the nucleation process. Initially the cloud elongates, confirming that nucleation is being mediated by excitation of a quadrupole mode (fig.1(a)) . Then finger-like structures appear on the outside edge of the condensate which eventually close round and produce vortices, $`800`$ms after rotation began (fig. 1(b)). Approximately 200 ms later, these have moved to equilibrium positions within the bulk of the condensate. Figures 1(c) and (d) show typical, single-shot images of stable vortex arrangements. The depth of each vortex (in the integrated absorption profile) is approximately $`95\%`$ of the surrounding condensate. The core diameters of $`5\mu `$m after 12 ms of free expansion are consistent with the predictions in . The maximum number of vortices we observed was seven, limited by the number of atoms in our condensate.

Our first study of the nucleation process involved counting the number of vortices as a function of the normalised trap rotation rate, $`\overline{\mathrm{\Omega }}=\mathrm{\Omega }/\omega _{}`$, for a fixed eccentricity. Results for trap deformations $`ϵ`$ = 0.084 and 0.041 are given in fig. 2. These graphs show a maximum and minimum value of $`\overline{\mathrm{\Omega }}`$ for nucleation at a given eccentricity. Increasing $`ϵ`$ increases the range of $`\overline{\mathrm{\Omega }}`$ over which vortices may be nucleated, both by lowering $`\overline{\mathrm{\Omega }}_{min}`$ and raising $`\overline{\mathrm{\Omega }}_{max}`$. In the limit of small eccentricities, $`\sqrt{2}\omega _{}`$ is the frequency of the m=2 quadrupole mode, which has been shown elsewhere to play a critical role in the nucleation process . Rotation of the ellipsoidal trap at half this frequency, i.e. $`\overline{\mathrm{\Omega }}_c=1/\sqrt{2}0.71`$, excites this mode. The plots in fig. 2 confirm that the nucleation depends on resonant excitation of the quadrupole mode as seen previously . The resonance is broader at higher eccentricity, as intuitively expected for stronger driving.

Our second study involved holding $`\overline{\mathrm{\Omega }}`$ constant and counting the number of vortices as a function of the trap deformation $`ϵ`$. Figure 3 shows our results for two cases: (a) $`\overline{\mathrm{\Omega }}>\overline{\mathrm{\Omega }}_c`$ and (b) $`\overline{\mathrm{\Omega }}<\overline{\mathrm{\Omega }}_c`$. Interestingly we were able to nucleate vortices under adiabatic conditions when $`\overline{\mathrm{\Omega }}<\overline{\mathrm{\Omega }}_c`$. We employ the adiabaticity criterion $`\dot{ϵ}/ϵ\omega _{}`$. As a further check of the adiabaticity we varied the time for the eccentricity ramp between $`200`$ ms and $`1`$ s, and detected no difference in the number of vortices formed.

The critical values of $`\overline{\mathrm{\Omega }}`$ and $`ϵ`$ for nucleation were extracted from plots such as fig. 2 and fig. 3 for many values of $`\overline{\mathrm{\Omega }}`$ and compiled on fig. 4. The data points show the minimum eccentricity required for a given rotation rate and map out region 2, within which vortices nucleate. $`\overline{\mathrm{\Omega }}_c`$ appears to be a critical rotation frequency at which vortices can be nucleated with minimum eccentricity, as predicted in . Changing $`\overline{\mathrm{\Omega }}`$ in either direction requires a more elliptical trap for nucleation, although different physical processes control the upper and lower limits as explained below.

At $`\overline{\mathrm{\Omega }}>\overline{\mathrm{\Omega }}_c`$, the condensate follows a particular quadrupole mode at small eccentricity, region 3 in fig. 4. This mode is referred to as the ’overcritical branch’ in  and has an elliptical density distribution which is orthogonal to the trap potential. It then nucleates vortices when the eccentricity is too large for the quadrupole mode to be a solution of the hydrodynamic equations. The boundary of the region in the $`ϵ`$ versus $`\overline{\mathrm{\Omega }}`$ plot where this quadrupole mode exists is given by

$$ϵ=\frac{2}{\overline{\mathrm{\Omega }}}\left(\frac{2\overline{\mathrm{\Omega }}^21}{3}\right)^{3/2}.$$
(4)
We determined this relation from the solutions of the hydrodynamic equations for superfluids and it is plotted as a solid line in fig. 4. This line agrees well with the experimental data for the critical conditions for nucleation for $`\overline{\mathrm{\Omega }}>\overline{\mathrm{\Omega }}_c`$ and a wide range of $`ϵ`$.

Below $`\overline{\mathrm{\Omega }}_c`$, the deformation needed to nucleate vortices appears to increase linearly with $`\overline{\mathrm{\Omega }}`$. This boundary cannot be explained in terms of the stability limit of a quadrupole mode - the ’normal branch’ is stable on both regions 1 and 2 and the ’overcritical branch’ is stable in neither. Our data appears to be at variance with the results in , where no vortices were seen when the eccentricity was increased adiabatically and $`\overline{\mathrm{\Omega }}<\overline{\mathrm{\Omega }}_c`$.

A mechanism for the creation of vortices at a frequency below $`\overline{\mathrm{\Omega }}_c`$ has been proposed in . They have shown that there are regions in the plot of $`ϵ`$ versus $`\overline{\mathrm{\Omega }}`$, both above and below $`\overline{\mathrm{\Omega }}_c`$, where the quadrupole solutions of the hydrodynamic equations become dynamically unstable. For frequencies above $`\overline{\mathrm{\Omega }}_c`$ the predicted instability domains coincide with the experimentally observed vortex domains in , and this work, thus indicating a link between their instability analysis and vortices. To make a quantitative prediction for the boundary between regions 1 and 2 shown in fig. 4, will require further detailed work for our specific case. We also note that it is possible that there exist non-negligible terms of higher order than quadratic in the rotating magnetic potential. These have been demonstrated to excite high order rotating surface modes, e.g. hexapole, which lead to vortex nucleation at frequencies less than $`\overline{\mathrm{\Omega }}_c`$ .

Another possible mechanism for observation of vortices below $`\overline{\mathrm{\Omega }}_c`$ is that the thermal cloud plays an important role. Transfer of angular momentum to the condensate from the spinning thermal cloud may provide a mechanism for vortices to form at $`\overline{\mathrm{\Omega }}<\overline{\mathrm{\Omega }}_c`$. However the transfer rate of angular momentum must be greater than any loss rate due to residual trap anisotropy . In , gravity produces a small static eccentricity in the trap in the plane of rotation. The ‘spin down’ time for a rotating thermal cloud in the presence of a deformation parameter $`ϵ`$ of only $`0.01`$ is very short, 0.5 s, compared to the spin up time of 15 s and hence the cloud may never gain significant angular momentum. However, in our experiment, gravity acts along the rotation axis and hence the trap is symmetric in the plane of rotation, giving a more favourable ratio of spin-up to spin-down times.

With this hypothesis in mind, we tested our nucleation curve at lower temperature to see if there was any change when the amount of thermal cloud was reduced. When acquiring the data of fig. 4 the rf was turned off after condensation and some heating was observed during the nucleation procedure, resulting in a temperature around $`0.8T_c`$. To achieve a lower temperature we left on the so-called ‘rf shield’ so as to give an effective trap depth of $`800`$ nK during the nucleation process. This resulted in a temperature of $`0.5T_c`$. No significant change was observed in the nucleation curve at this lower temperature. However, this does not totally rule out a role for the thermal cloud in the nucleation process since even at $`0.5T_c`$, there was still a significant $`20\%`$ of atoms in the thermal cloud.

Although the exact amount of thermal cloud seemed to have little effect on the nucleation conditions, it had a striking effect on the behaviour of vortices after nucleation. Without the rf shield during the nucleation process, vortices were only occasionally found in an equilibrium configuration (i.e. 1 vortex in the centre, 3 vortices in an equilateral triangle as in fig. 1) and $`400`$ms after forming they had already moved to the edge of the condensate before disappearing. However with the rf shield present, the vortices were normally found in equilibrium positions $`200`$ ms after formation and had a lifetime of $`4s`$, only limited by the decay of the condensate itself.

In summary, we have used a purely magnetic rotating trap to investigate conditions for vortex nucleation (after the formation of the condensate) over a wide range of trap eccentricities. For a given eccentricity, we observe both an upper and lower limit to the rotation rate for nucleation. The upper limit confirms the predictions in , but over a much wider range of parameters. However the lower limit to the rotation rate is different to that reported elsewhere. Further theoretical work is required to explain the linear dependence of $`ϵ`$ on $`\overline{\mathrm{\Omega }}`$ shown in fig. 4. The thermal cloud is shown to destabilise vortex arrays, thus measurements of the vortex lifetime in an oblate geometry were made with an rf shield to prevent heating.

The rotating anisotropic magnetic trap is well suited to experiments that require a large trap eccentricity. It has recently been used to investigate the irrotational behaviour of vortex-free condensates at low rotation rates . A large eccentricity is also needed to spin up a thermal cloud , hence this trap may be useful in attempts to condense directly into a vortex state from a spinning thermal cloud in thermodynamic equilibrium.

We would like to thank D. Feder for useful discussions. This work was supported by the EPSRC and the TMR program (No. ERB FMRX-CT96-0002). O.M. Maragò acknowledges the support of a Marie Curie Fellowship, TMR program (No. ERB FMBI-CT98-3077).


text:
# Building dwarf galaxies out of merged young star clusters

## 1 Introduction

Systems of interacting galaxies (like The Antennae or Stephan’s Quintet) show intense star formation. In the tidal tails one finds especially young massive star clusters ($`R_{\mathrm{eff}}210`$ pc; $`M10^410^7`$ M; see e.g. Whitmore et al. (1999)). These star clusters are not distributed uniformly but are clustered into groups, some containing more than a hundred clusters. With this project we investigate the future fate of these clusters of star clusters (what we call “super-cluster”). Following the arguments of Kroupa (1998) these super-clusters may be bound entities, meaning that the star clusters in the super-cluster are not only bound to the tidal field of the parent galaxy but also bound to each other. Therefore, the star clusters are very likely to merge with each other on a few dynamical times of the super-cluster. This merging process and the resulting merger object are the major topic of this project.

## 2 Setup

We simulate super-clusters containing $`N_0=20`$ star clusters. Each of them is represented by a Plummer sphere, containing 100,000 particles, a total mass of $`10^6`$ M, a scale-length (Plummer radius) of $`R_{\mathrm{pl}}=6`$ pc and a crossing time of 1.4 Myr. All star clusters are distributed into an encompassing Plummer sphere with varied scale-length, $`R_{\mathrm{pl}}^{\mathrm{sc}}=50,75,150,300`$ pc, and this super-cluster is placed on a circular orbit at distance $`D=5,10,20,30,50,100`$ kpc around an analytic galactic potential with an asymptotically flat rotation curve with $`v_\mathrm{c}=220`$ km/s and $`R_{\mathrm{gal}}=4`$ kpc. For the discussion of the results we define two dimensionless parameters for our systems. $`\alpha =R_{\mathrm{pl}}/R_{\mathrm{pl}}^{\mathrm{sc}}=0.02,0.04,0.08,0.12`$ is the ratio between the scale-length of the single cluster and the scale-length of the super-cluster. $`\beta =R_{\mathrm{cut}}^{\mathrm{sc}}/R_{\mathrm{tidal}}=0.2`$$`2.7`$ is the ratio between the size and the tidal radius of the super-cluster. Several random realizations of each parameter combination are performed to reduce stochastic errors. The simulations are carried out with the particle mesh code Superbox (Fellhauer et al. 2000) which has the ability to trace many objects with different levels of high resolution sub-grids. While the outermost grid covers the whole orbit around the galactic centre with the lowest resolution, the medium sized grids of each star cluster cover the size of the super-cluster and all intermediate interactions between the clusters in medium resolution (10–50 pc), and the innermost high resolution (1 pc) grids cover each single cluster. Both medium and high resolution grids stay focused on the simulated objects and travel with them through the simulation area.

## 3 Compact Merger Objects

A detailed description of the merging time-scales and the efficiency of the merger process will appear in Fellhauer et al. (2001).

The resulting merger objects are compact, spherical and dense objects. They have half-mass radii from 40–150 pc and sizes of 120–1,000 pc. The mass of these objects resembles between 75 and 95 % of the total mass of the initial super-cluster. Central densities are derived by fitting exponential profiles to the inner parts. The range of central densities is 100–1,000 M/pc3 with exponential scale lengths of about 10 pc. In the outer parts a steep power law ($`r^{3.5}`$) fits the density best. The surface-density can be well fitted with a de-Vaucouleur profile. The velocity dispersion in the central region is about 30 km/s but measured along the line-of-sight through the centre one derives about 20 km/s. All this data suggests that our merger objects have similar properties to nucleated dwarf galaxies (dE,N). Because the merger process lasts for about a dozen crossing times of the super-cluster (for small $`\alpha `$’s this can sum up to more than 1 Gyr) the merger object is surrounded by several star clusters which are still in the process of merging or are on near-circular orbits around the super-cluster and slowly decaying through dynamical friction. Therefore, our model can account for high specific cluster frequencies ($`S_\mathrm{N}`$), which are found for dE,N-galaxies.

Furthermore, our simple models do not account for destructive processes like tidal heating on eccentric orbits, dynamical friction in the parent galaxy and decay of orbit, shocks due to disc or bulge passages. All these processes would lead to a mass-loss of the merger object and a decrease in its central density. Translated to observations, our objects would transform into dE,noN with time, while the number of remaining clusters is also decreasing. Observations show that dE,N have higher $`S_\mathrm{N}`$ than dE,noN.

Our objects also show extended tidal tails along the whole orbit around the galactic centre. Escaped star clusters are found in these tidal features.

## 4 Dissolved Merger Object

In two test calculations the super-cluster was set up on an eccentric orbit ($`D_{\mathrm{apo}}=60`$ kpc, $`D_{\mathrm{peri}}=30`$ kpc) so that tidal heating could act on the merger objects. The merger objects get dissolved with time (mass-loss $`5`$ % per perigalacticon passage). Due to a disruptive interaction with a second object the central densities drop down from several hundred to 0.3–0.002 M/pc3, but can still be fitted with an exponential in the innermost part. But in the outer parts the material is not bound to the object anymore. It has a “fluffy” structure extending a few kpc. The surface-density can not be fitted with a de-Vaucouleur profile anymore. The measured central velocity dispersions are in the order of 1–3 km/s, but measured along the line-of-sight through the centre and the outer unbound parts one derives about 10 km/s. Compared with observations these objects look like dwarf spheroidal galaxies (c.f. Klessen & Kroupa 1998). Star clusters remaining around these objects can account for the very high $`S_\mathrm{N}`$ of some dSph-galaxies.

## 5 Conclusions

With our simulations we are able to show, that the star clusters in a super-cluster merge efficiently. Comparing our parameters with observations we predict an exponential decrease with time of the number of clusters if the observed super-clusters are bound entities.

The resulting merger object has similar properties like dE,N-galaxies, and if destructive processes are taken into account the next step in the evolution of our objects would be dE,noN. This also gives a natural explanation why dE,N have higher $`S_\mathrm{N}`$ than dE,noN.

The final stage of the evolution would be a dSph-galaxy, suggesting that at least some of these objects may not have a predicted high dark matter content.

References:

Fellhauer, M., Kroupa, P., Baumgardt, H., Bien, R., Boily, C.M., Spurzem, R., & Wassmer, N. 2000, NewA, 5, 305
Fellhauer, M., & Kroupa P. 2000, ASP Conf. Ser., 211, 2411
Fellhauer, M., Baumgardt, H., Kroupa, P., & Spurzem, R. 2001, Cel. Mech. & Dyn. Astron. in press; astro-ph/0103052
Gallagher, S.C., Hunsberger, S.D., Charlton, J.C., & Zaritsky, D. 2000, ASP Conf. Ser., 211, 2471
Grebel, E.K. 2000, ASP Conf. Ser., 211, 2621
Hunsberger, S.D., Gallagher, S.D., Charlton, J.C., & Zaritsky, D. 2000, ASP Conf. Ser., 211, 2541
Klessen, R.S., & Kroupa, P. 1998, ApJ, 498, 143
Kroupa, P. 1998, MNRAS, 300, 200
Mateo, M. 1998, ARAA, 36, 435
Whitmore, B.C., Zhang, Q., Leitherer, C., & Fall, S.M. 1999, AJ, 118, 1551
Zhang, Q., & Fall, S.M. 1999, ApJL, 527, 81L

1: ASP Conf. Ser. 211: Proceedings of the workshop ‘Massive Stellar Clusters’ held in Strasbourg Nov. 1999; eds. A. Lancon, C.M. Boily