text:
# Mixing of 𝑑_{𝑥²-𝑦²} and 𝑑_{𝑥𝑦} superconducting states for different filling and temperature PACS # 74.20.Fg, 74.62.-c, 74.25.Bt. KEYWORDS mixed-symmetry state, high-𝑇_𝑐 superconductor, 𝑑_{𝑥²-𝑦²} and 𝑑_{𝑥𝑦} order parameters

## I Introduction

Many years after the discovery of the unconventional high-$`T_c`$ cuprates  with a high critical temperature $`T_c`$ the symmetry of its order parameter continues the subject of active investigation. From both experimental observations  and related theoretical analyses  there is a consensus that immediately below $`T_c`$, the symmetry of the order parameter is of the $`d_{x^2y^2}`$ type. However, many experiments  and theoretical analyses  suggest that at lower temperatures the order parameter has a mixed symmetry of the $`d_{x^2y^2}+\mathrm{exp}(i\theta )\chi `$ type, where $`\chi `$ represents a minor-component state with a distinct symmetry superposed on the major component $`d_{x^2y^2}`$. From theoretical analysis the mixing angle $`\theta `$ can have the values $`0`$, $`\pi /2,\pi ,`$ and $`3\pi /2`$. For mixing angles $`\pi /2`$ and $`3\pi /2`$ the time-reversal symmetry is broken for these states . Four possible candidates for the minor symmetry state $`\chi `$ are the $`d_{xy}`$, $`s`$, $`s_{x^2+y^2}`$ and $`s_{xy}`$ states. From a group theoretical point of view these states belong to the same irreducible representation of the orthorhombic point group. However, there is still controversy about the nature of the minor component and the value of the mixing angle in different high-$`T_c`$ materials. There have been extensive discussions in recent literature about the nature of the minor component and the value of the mixing angle. Without repeating those discussions in detail here we present a brief summary of these analyses which suggest a mixed-symmetry state and refer the interested reader to these references .

There are experimental evidences based on the Josephson supercurrent for tunneling between a conventional $`s`$-wave superconductor (Pb) and twinned or untwinned single crystals of YBa2Cu3O7 (YBCO) that YBCO has a mixed-symmetry state with a minor $`s`$-wave component at lower temperatures . Kouznetsov et al.  performed $`c`$-axis Josephson tunneling experiments by depositing a conventional superconductor (Pb) across the single twin boundary of a YBCO crystal and from a study of the critical current under an applied magnetic field they also found evidence of a mixed-symmetry state with a minor $`s`$-wave component. Similar conclusion was reached by Sridhar et al.  from the measurement of the microwave complex conductivity of high quality YBa2Cu3O7-δ single crystals at 10 GHz using a high-$`Q`$ Nb cavity.

From an analysis of angle-resolved photoemission spectroscopy experiment, Ma et al.  detected a temperature-dependent gap anisotropy in Bi2Sr2CaCu2O8+x which can be associated with a mixed-symmetry state with a minor $`d_{xy}`$ component. More recently, Krishana et al.  reported a phase transition to a mixed-symmetry state with a minor $`d_{xy}`$ component in Bi2Sr2CaCu2O8 induced by a magnetic field from a study of thermal conductivity as a function of temperature and applied field. Later, the experimental study of thermal conductivity in Bi2Sr2CaCu2O8 by Movshovich et al.  also leads to the same conclusion.

Different theoretical investigations  provided explanations of the observation by Krishana et al.  that for weak magnetic field possibly a minor time-reversal symmetry-breaking state is induced on a dominating $`d_{x^2y^2}`$ state in Bi2Sr2CaCu2O8. This mixed-symmetry state is likely to be a minor $`d_{xy}`$ component on a $`d_{x^2y^2}`$ state at lower temperatures. The above investigations in Bi2Sr2CaCu2O8 suggest a minor $`d_{xy}`$ component in the order parameter. In view of this we perform a theoretical study of this problem using a tight-binding model including orthorhombic distortion and second-nearest-neighbor hopping. Both orthorhombicity and filling are expected to play a crucial role in the evolution of these mixed-symmetry states and we study in this paper the properties of the mixed-symmetry states of the type $`d_{x^2y^2}+\chi d_{xy}`$ at different temperature, filling, orthorhombicity, and second-neighbor hopping. A similar study of the mixture of the $`d_{x^2y^2}`$ and different types of $`s`$ states will be reported elsewhere.

There is no suitable microscopic theory for high-$`T_c`$ cuprates and there is controversy about a proper description of the normal state and the pairing mechanism for such materials. Also, for underdoped systems the existence of a pseudo gap above the superconducting critical temperature can not yet be understood within a microscopic theory. However, with the increase of doping the superconducting critical temperature increases and the pseudo gap is reduced and eventually it disappears at optimal doping. In the absence of a microscopic theory, we use a phenomenological two-dimensional tight-binding model with appropriate lattice symmetry for studying some of the general features of the mixed-symmetry states of overdoped high-$`T_c`$ superconductors where there is no pseudo gap. This model has been used successfully in describing many properties of high-$`T_c`$ materials .

In high-$`T_c`$ cuprates, it is generally believed that superconductivity resides mainly in the CuO2 planes with nearly tetragonal symmetry. Hence, quite expectedly, the values of critical temperature and superconducting gap of high-$`T_c`$ materials are very sensitive to the level of doping which determines the number of electrons available for conduction in the two-dimensional CuO2 plane of the cuprates. However, a large anisotropy in the penetration depth between the $`a`$ and $`b`$ directions has been observed experimentally  in YBCO and this cannot be understood in a purely tetragonal system. Therefore, in YBCO there is the evidence of a weak orthorhombic distortion. Hence we shall consider orthorhombic distortion in the present study of mixed-symmetry states.

One of the objectives of the present study is to investigate how the superconductivity in the uncoupled and mixed-symmetry states changes with doping or filling. The present single-band model can acomodate a maximum filling of two electrons (spin up and down) per unit cell. In addition to studying the evolution of the mixed-symmetry states with respect to filling or doping in presence of second-nearest-neighbor hopping for tetragonal and orthorhombic lattice, we also study the temperature dependence of the order parameters for a fixed filling.

In Sec. II we describe the formalism. In Sec. III we present our numerical results and finally in Sec. IV we give a brief summary.

## II Theoretical Formulation

We base the present study on the two-dimensional tight binding model which we describe below. This model is sufficiently general for considering mixed angular momentum states on tetragonal and orthorhombic lattice, employing nearest and second-nearest-neighbor hopping integrals. In the present tight-binding model the effective interaction $`V_{\mathrm{𝐤𝐪}}`$ for transition from a momentum $`𝐪`$ to $`𝐤`$ is taken to be separable, and is expanded in terms of some general basis functions $`\eta _{i𝐤}`$, labeled by index $`i`$, so that

$$V_{\mathrm{𝐤𝐪}}=\underset{i}{}V_i\eta _{i𝐤}\eta _{i𝐪}$$
(1)
The separable nature of the interaction facilitates the solution of the gap equation. The orthogonal functions $`\eta _{i𝐤}`$ are associated with a one-dimensional irreducible representation of the point group of square lattice $`C_{4v}`$ and are appropriate generalizations of the circular harmonics incorporating the proper lattice symmetry. Here $`V_i`$ is the coupling of effective interaction in the specific angular momentum state specified by the index $`i`$ above. In general we have the following $`\eta `$’s

$`\eta _{1𝐪}`$ $`=`$ $`\mathrm{cos}q_x\beta \mathrm{cos}q_y,d_{x^2y^2}\text{-wave},`$ (2)
$`\eta _{2𝐪}`$ $`=`$ $`2\mathrm{sin}q_x\mathrm{sin}q_y,d_{xy}\text{-wave},`$ (3)
$`\eta _{3𝐪}`$ $`=`$ $`1,s\text{-wave},`$ (4)
$`\eta _{4𝐪}`$ $`=`$ $`\mathrm{cos}q_x+\beta \mathrm{cos}q_y,s_{x^2+y^2}\text{-wave},`$ (5)
$`\eta _{5𝐪}`$ $`=`$ $`2\mathrm{cos}q_x\mathrm{cos}q_y,s_{xy}\text{-wave},`$ (6)
etc. Here $`\beta =1`$ corresponds to square lattice and $`\beta 1`$ represents orthorhombic distortion. The orthogonality property of functions $`\eta `$’s is given by

$`{\displaystyle \underset{𝐪}{}}\eta _{i𝐪}\eta _{j𝐪}=0,\mathrm{},ij.`$ (7)
In the present investigation of mixed-symmetry states involving $`d_{x^2y^2}`$ and $`d_{xy}`$, it is appropriate to consider only two terms in the sum (1) with $`\eta _1`$ and $`\eta _2`$ above. We shall consider singlet Cooper pairing and subsequent condensation in $`d_{x^2y^2}`$ and $`d_{xy}`$ waves, denoted by indices 1 and 2, respectively, and the mixed-symmetry state formed by these two. Consequently, the potential couplings $`V_1`$ and $`V_2`$ will be taken to be positive corresponding to attractive interactions.

We consider a single tight-binding two-dimensional band with electron dispersion relation including second-nearest-neighbor hopping. In this case the quasiparticle dispersion relation relating the electronic energy $`ϵ_𝐤`$ and momentum $`𝐤`$ is taken as

$$ϵ_𝐤=2t(\mathrm{cos}k_x+\beta \mathrm{cos}k_y2\gamma \mathrm{cos}k_x\mathrm{cos}k_y)\mu ,$$
(8)
where $`t`$ and $`\beta t`$ are the nearest-neighbor hopping integrals along the in-plane $`a`$ and $`b`$ axes, respectively, and $`\gamma t`$ is the second-nearest-neighbor hopping integral. In Eq. (8) $`\mu `$ is the chemical potential measured with respect to the Fermi energy and is determined once the filling $`n`$ is specified. The nearest-neighbor hopping parameter $`t`$ is typically taken to be $`0.1`$ eV. The parameter $`\beta `$ destroys the symmetry between the $`a`$ and $`b`$ directions in the CuO2 planes in this simple model. The potential $`V_{\mathrm{𝐤𝐪}}`$ above also possesses such a symmetry-breaking term. The energy $`ϵ_𝐤`$ is measured with respect to the Fermi surface. Such a one-band model with different first-neighbor-hopping parameters in the $`a`$ and $`b`$ directions is the simplest approximate way of including in the theoretical description the effect of orthorhombicity.

At a finite temperature $`T`$, one has the following gap equation

$`\mathrm{\Delta }_𝐤`$ $`=`$ $`{\displaystyle \underset{𝐪}{}}V_{\mathrm{𝐤𝐪}}{\displaystyle \frac{\mathrm{\Delta }_𝐪}{2E_𝐪}}\mathrm{tanh}{\displaystyle \frac{E_𝐪}{2k_BT}}`$ (9)
with $`E_𝐪=[(ϵ_𝐪\mu )^2+|\mathrm{\Delta }_𝐪|^2]^{1/2},`$ and $`k_B`$ the Boltzmann constant.

It has been observed that the critical temperature $`T_c`$ is sensitive to the level of doping which determines the number of available conduction electrons in the CuO2 plane. In this model the chemical potential $`\mu `$ and the the filling $`n`$ are determined by the number equation

$$n=1\underset{𝐪}{}\frac{ϵ_𝐪\mu }{E_𝐪}\mathrm{tanh}\frac{E_𝐪}{2k_BT}.$$
(10)
The filling $`n`$ can be related to the the experimental doping $`\delta `$ in the three dimensional Brillouin zone by $`n=1\delta `$. In this work we study the variation of $`n`$ from 0 to 1 (half filling).

The order parameter $`\mathrm{\Delta }_𝐪`$ has the following anisotropic form:

$$\mathrm{\Delta }_𝐪=\mathrm{\Delta }_1\eta _{1𝐪}+C\mathrm{\Delta }_2\eta _{2𝐪},$$
(11)
where $`C\mathrm{exp}(i\theta )`$ is a complex number of unit modulus $`|C|^2=\mathrm{cos}^2\theta +\mathrm{sin}^2\theta =1`$. We substitute Eqs. (1) and (11) into the gap equation (9) and using the orthogonality property (7) obtain the two following coupled equations for $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$:

$`\mathrm{\Delta }_1`$ $`=`$ $`{\displaystyle \underset{𝐪}{}}V_1\eta _{1𝐪}{\displaystyle \frac{\mathrm{\Delta }_1\eta _{1𝐪}+C\mathrm{\Delta }_2\eta _{2𝐪}}{2E_𝐪}}\mathrm{tanh}{\displaystyle \frac{E_𝐪}{2k_BT}},`$ (12)
$`C\mathrm{\Delta }_2`$ $`=`$ $`{\displaystyle \underset{𝐪}{}}V_2\eta _{2𝐪}{\displaystyle \frac{\mathrm{\Delta }_1\eta _{1𝐪}+C\mathrm{\Delta }_2\eta _{2𝐪}}{2E_𝐪}}\mathrm{tanh}{\displaystyle \frac{E_𝐪}{2k_BT}}.`$ (13)
In order to proceed further we have to know the nature of $`C`$, e. g., whether $`C`$ is real, imaginary or complex. As $`|C|^2=1`$, a purely real $`C`$ implies $`C=\pm 1`$ and in this case no further simplification of the coupled Eqs. (12) and (13) is possible and one has

$`\mathrm{\Delta }_1`$ $`=`$ $`{\displaystyle \underset{𝐪}{}}V_1\eta _{1𝐪}{\displaystyle \frac{\mathrm{\Delta }_1\eta _{1𝐪}\pm \mathrm{\Delta }_2\eta _{2𝐪}}{2E_𝐪}}\mathrm{tanh}{\displaystyle \frac{E_𝐪}{2k_BT}},`$ (14)
$`\pm \mathrm{\Delta }_2`$ $`=`$ $`{\displaystyle \underset{𝐪}{}}V_2\eta _{2𝐪}{\displaystyle \frac{\mathrm{\Delta }_1\eta _{1𝐪}\pm \mathrm{\Delta }_2\eta _{2𝐪}}{2E_𝐪}}\mathrm{tanh}{\displaystyle \frac{E_𝐪}{2k_BT}},`$ (15)
with

$$E_𝐪=[(ϵ_𝐪\mu )^2+(\mathrm{\Delta }_1\eta _{1𝐪}\pm \mathrm{\Delta }_2\eta _{2𝐪})^2]^{1/2},$$
(16)
The upper sign in Eqs. (14), (15), and (16) corresponds to a mixing angle $`\theta =0`$ ($`C=1`$) and the lower sign to $`\theta =\pi `$ ($`C=1`$).

However, Eqs. (12) and (13) get simplified when $`C`$ is imaginary, e. g., $`C=\pm i`$. Consequently,

$$E_𝐪=[(ϵ_𝐪\mu )^2+\mathrm{\Delta }_1^2\eta _{1𝐪}^2+\mathrm{\Delta }_2^2\eta _{2𝐪}^2]^{1/2},$$
(17)
and Eqs. (12) and (13) reduce to

$`\mathrm{\Delta }_1`$ $`=`$ $`{\displaystyle \underset{𝐪}{}}V_1{\displaystyle \frac{\mathrm{\Delta }_1\eta _{1𝐪}^2}{2E_𝐪}}\mathrm{tanh}{\displaystyle \frac{E_𝐪}{2k_BT}},`$ (18)
$`\mathrm{\Delta }_2`$ $`=`$ $`{\displaystyle \underset{𝐪}{}}V_2{\displaystyle \frac{\mathrm{\Delta }_2\eta _{2𝐪}^2}{2E_𝐪}}\mathrm{tanh}{\displaystyle \frac{E_𝐪}{2k_BT}},`$ (19)
since in this case

$`{\displaystyle \underset{𝐪}{}}{\displaystyle \frac{\eta _{1𝐪}\eta _{2𝐪}}{2E_𝐪}}\mathrm{tanh}{\displaystyle \frac{E_𝐪}{2k_BT}}=0,`$ (20)
which follows from the definitions of $`\eta _{1𝐪}`$ and $`\eta _{2𝐪}`$ and Eq. (17). However, Eq. (20), which is responsible for the simplification for $`C=\pm i`$ ($`\theta =\pi /2`$ and $`3\pi /2`$), does not hold for $`C=\pm 1`$ or for a general complex $`C`$.

Finally, for a general complex $`C`$ one can separate in a straightforward fashion Eqs. (12) and (13) into their real and imaginary parts. In that case Eq. (20) is not valid, and this results in four equations for the two unknowns $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$. These four equations are consistent only if $`\mathrm{\Delta }_1=0`$ or $`\mathrm{\Delta }_2=0`$, which means that there could not be any mixing between the two components. So mixed-symmetry states are allowed only for mixing angles $`\theta =0,\pi /2,\pi ,`$ and $`3\pi /2`$, or for $`C=\pm 1,`$ and $`\pm i`$ and we shall consider only these cases in the following. Essentially the same conclusion was reached in a different context in Ref.  in a discussion of mixed-symmetry s-d states on the continuum. In that case mixing was also possible only for mixing angles $`0,\pi /2,\pi `$ and $`3\pi /2`$. However, in Ref.  the discrete sum over $`𝐪`$ was replaced by a continuum integral over angle. We note that the $`C=\pm i`$ cases are physically equivalent and lead to the same gap equations (18) and (19). On the other hand, $`C=\pm 1`$ cases corresponding to $`\theta =0`$ and $`\pi `$, respectively, lead to two different sets of gap equations (14) and (15).

The ultraviolet momentum-space divergence of the Bardeen-Cooper-Schrieffer equation was originally neutralized by a physically-motivated Debye cut off . This procedure had the advantage of reproducing the experimentally observed isotope effect. It can also be handled by using the technique of renormalization . Here we introduce a cut off in the momentum sums of the gap equation. As there is no pronounced isotope effect in the high-$`T_c`$ cuprates, the present cut off is merely a mathematical one without any reference to the phonon-induced Debye cut off. In Eqs. (12) and (13) both the interactions $`V_1`$ and $`V_2`$ are assumed to be energy-independent constants for $`|ϵ_𝐪\mu |<k_BT_D`$ and zero for $`|ϵ_𝐪\mu |>k_BT_D`$, where $`k_BT_D`$ is the present cut off.

## III Numerical Result

We solve the coupled set of equations (14) and (15) or (18) and (19) in conjunction with the number equation (10). This model gives us the opportunity to study the mixed-symmetry state composed of $`d_{x^2y^2}`$ and $`d_{xy}`$ for different filling and temperature for both tetragonal and orthorhombic lattice. This study will also provide the interesting variation of the chemical potential and critical temperature with filling. We solve the above sets of equations numerically by the method of iteration and calculate the gaps $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$ at various fillings and temperatures.

We perform calculations on a perfect square lattice with second-nearest-neighbor hopping contribution. In addition to the superconductivity in coupled $`d_{x^2y^2}+d_{xy}`$ and $`d_{x^2y^2}+id_{xy}`$ waves we also consider the uncoupled $`d_{x^2y^2}`$ and $`d_{xy}`$ waves. The effect of orthorhombic distortion on our investigation is also studied. The results of our study have interesting variation as the second nearest hopping parameter is varied and this is studied in detail in the following for mixing angles 0 and $`\pi /2`$. Throughout the calculation we consider the cut off $`k_BT_D=0.1t`$ eV with the parameter $`t=0.2586`$ eV. This corresponds to a cut-off of $`T_D=300`$ K. The results for the chemical potential $`\mu `$ and order parameters $`\mathrm{\Delta }_{x^2y^2}`$ and $`\mathrm{\Delta }_{xy}`$ presented in this work are all in units of $`t`$.

First, we study the superconductivity in uncoupled $`d_{x^2y^2}`$ and $`d_{xy}`$ waves by setting $`\mathrm{\Delta }_2=0`$ and $`\mathrm{\Delta }_1=0`$, respectively, in the equations above at $`T=0`$ on square lattice ($`\beta =1`$). For $`d_{x^2y^2}`$ wave we use potential parameters $`V_1=0.73t`$ and $`V_2=0`$, and for $`d_{xy}`$ wave we use $`V_1=0`$ and $`V_2=2.1t`$. The resulting $`\mathrm{\Delta }`$’s with the variation of filling $`n`$ are plotted in Figs. 1(a) and 1(b), respectively, for three different values of the second-nearest-neighbor hopping parameter $`\gamma `$: $`\gamma =`$ 0, 0.1 and 0.2. In the present study we shall always use $`V_1=0.73t`$ on square lattice; for $`\gamma =0`$ this leads to a critical temperature of 71 K . The value of $`V_2`$ would be chosen to have a reasonable mixture between the two components.

From Fig. 1(a) we find that for $`\gamma =0`$, $`\mathrm{\Delta }_1\mathrm{\Delta }_{x^2y^2}`$ has a pronounced maximum at half filling and it reduces quickly to zero as $`n`$ decreases. As second nearest hopping $`\gamma `$ increases, the position of this maximum shifts to smaller $`n`$. For $`\gamma =0.2`$, this maximum is at $`n0.8`$. In Fig. 1(b), we see that for $`\gamma =0`$, $`\mathrm{\Delta }_2\mathrm{\Delta }_{xy}`$ also has a maximum at $`n=1`$; however, this maximum is much wider than that in the case of $`\mathrm{\Delta }_1`$. The position of this maximum also shifts towards smaller $`n`$ as the parameter $`\gamma `$ increases. However, the change in the nature of the curves as $`\gamma `$ increases is more pronounced for $`d_{x^2y^2}`$ in Fig. 1(a) than for $`d_{xy}`$ in Fig. 1(b). The variation of $`\mathrm{\Delta }`$ with $`n`$ in Fig. 1(a) is also quite distinct from that in Fig. 1(b).

Now we study superconductivity in coupled $`d_{x^2y^2}+d_{xy}`$ wave at $`T=0`$ on a square lattice, which corresponds to mixing angle 0. In this case Eqs. (14) and (15) are applicable. It is found that with $`\gamma =0,\beta =1,`$ and $`V_1=0.73t`$, mixing between the two components occur for $`V_2>1.7t`$. Here we employ a $`V_2`$ slightly larger than this critical value to have adequate mixing. The parameters for this model are the following: $`\beta =1`$, $`V_1=0.73t`$, and $`V_2=2.1t`$. In Fig. 2 we plot $`\mathrm{\Delta }_1\mathrm{\Delta }_{x^2y^2}`$ and $`\mathrm{\Delta }_2\mathrm{\Delta }_{xy}`$ in this case for different filling $`n`$ and for $`\gamma =0`$, 0.05, 0.1, and 0.2. Although the nature of $`\mathrm{\Delta }_1`$ of Fig. 2 is similar to that of $`\mathrm{\Delta }_1`$ of Fig. 1(a), the present maxima are narrower. The order parameter $`\mathrm{\Delta }_2`$ in the present mixed-symmetry case has a different nature for larger $`n`$ compared to the uncoupled case. Here, in the presence of second-neighbor hopping, as $`n`$ decreases, $`\mathrm{\Delta }_2`$ first decreases rapidly to a minimum and then increases and finally decreases to zero. Moreover, the order parameter $`\mathrm{\Delta }_2`$ for this mixed wave has a minimum at the same position of the maximum of $`\mathrm{\Delta }_1`$.

Next we study superconductivity in coupled $`d_{x^2y^2}+id_{xy}`$ wave at $`T=0`$ on a square lattice, governed by Eqs. (18) and (19), which corresponds to the mixing angle $`\pi /2`$ or $`3\pi /2`$. The parameters for this model are the following: $`\beta =1`$, $`V_1=0.73t`$, and $`V_2=1.7t`$. In this case to have mixing between the two components for $`\gamma =0`$ the critical value of $`V_2`$ is $`1.0t`$. In Fig. 3 we plot $`\mathrm{\Delta }_1\mathrm{\Delta }_{x^2y^2}`$ and $`\mathrm{\Delta }_2\mathrm{\Delta }_{xy}`$ in this case for different filling $`n`$ and $`\gamma =0`$, 0.05, 0.1, and 0.2. The qualitative nature of $`\mathrm{\Delta }_1`$ of Fig. 3 is similar to that of Figs. 1(a) and 2. The order parameter $`\mathrm{\Delta }_2`$ for this coupled wave also has a minimum at the same position of the maximum of $`\mathrm{\Delta }_1`$. The minimum in $`\mathrm{\Delta }_{xy}`$ is less pronounced in this case than in the case of $`d_{x^2y^2}+d_{xy}`$ considered in Fig. 2.

Next we study the behavior of the mixed-symmetry states in the presence of orthorhombic distortion ($`\beta =0.95`$). First we consider the $`d_{x^2y^2}+d_{xy}`$ wave. In this case we needed to increase $`V_1`$ to $`0.97t`$ in order to have a critical temperature of 71 K for $`\gamma =0`$ . The parameters for this model are the following: $`\beta =0.95`$, $`V_1=0.97t`$ and $`V_2=2.1t`$. To have mixing in this case the critical value of $`V_2`$ is $`1.5t`$ for $`\gamma =0`$. In Fig. 4 we plot $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$ in this case for different filling $`n`$ and second neighbor parameter $`\gamma =0`$, 0.05, 0.1, and 0.2. We find significant change in the behavior of the order parameters in this case compared to the case of square lattice discussed above. In this case $`\mathrm{\Delta }_{xy}`$ could become negative for non-zero values of the second-neighbor hopping parameter $`\gamma `$ for values of $`n`$ close to half filling and the phase of the mixing changes with a variation of $`n`$ in the presence of second-nearest hopping $`\gamma `$. For a non-zero $`\gamma `$ the mixing angle $`\theta `$ changes from 0 to $`\pi `$ as $`n`$ increases. Also, with the increase of $`\gamma `$ the zero of $`\mathrm{\Delta }_{xy}`$ shifts to a lower value of $`n`$,

Finally, we consider the $`d_{x^2y^2}+id_{xy}`$ wave on orthorhombic lattice with the following parameters: $`\beta =0.95`$, $`V_1=0.97t`$ and $`V_2=1.7t`$. For $`\gamma =0`$, the critical $`V_2`$ for mixing in this case is $`1.2t`$. In Fig. 5 we plot $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$ in this case for different filling $`n`$ for $`\gamma =0`$, 0.05, 0.1, and 0.2. Here, the maximum of $`\mathrm{\Delta }_{x^2y^2}`$ is also accompanied by the minimum of $`\mathrm{\Delta }_{xy}`$. The order parameter $`\mathrm{\Delta }_1`$ behaves in a similar way as in Fig. 4, but $`\mathrm{\Delta }_2`$ does not become negative. Both in Figs. 4 and 5 we observe a dip in $`\mathrm{\Delta }_1`$ near the maximum. This behavior was not found in Figs. 2 and 3 for square lattice. In Figs. 2 $``$ 5 we do not find any superconductivity for very low values of filling.

Now we investigate the temperature dependence of the order parameters in different cases. We studied several cases for different values of $`n`$, $`\gamma `$ and $`\beta (=1,0.95)`$. The qualitative nature of the temperature dependence in different cases is quite similar, e. g., the order parameters start from a large non-zero value at $`T=0`$. With the increase of temperature they first remain essentially constant and then decrease very slowly. With further increase of temperature they decrease more rapidly and eventually become zero at a large enough temperature. Figure 6 illustrates a typical example of this behavior for $`\beta =1`$, $`n=0.8`$, and $`\gamma =0.2`$ for the $`d_{x^2y^2}+d_{xy}`$ and $`d_{x^2y^2}+id_{xy}`$ cases. The potential parameters in these cases are the same as in Figs. 2 and 3. The present temperature dependence of the order parameters is quite similar to our earlier studies for $`n=1`$ for $`d_{x^2y^2}+d_{xy}`$ and $`d_{x^2y^2}+id_{xy}`$ waves .

In Fig. 7 we plot the chemical potential $`\mu `$ versus $`n`$ for the four different models considered in Figs. 2 $``$ 5 for the $`d_{x^2y^2}+d_{xy}`$ and $`d_{x^2y^2}+id_{xy}`$ cases on tetragonal and orthorhombic lattice. The parameters for the model are the same as in Figs. 2 $``$ 5. Surprisingly enough we find that $`\mu `$ is essentially the same for the two types of mixing $`d_{x^2y^2}+d_{xy}`$ and $`d_{x^2y^2}+id_{xy}`$, independent of the lattice symmetry. However, $`\mu `$ depends on the parameter $`\gamma `$ denoting second-neighbor hopping. So in Fig. 7 we have four different curves for four different values of second-neighbor hopping. The nature of variation of $`\mu `$ with $`n`$ changes with $`\gamma `$. The fall in $`\mu `$ with $`n`$ is most rapid in the case of $`\gamma =0`$ and slows down as we increase $`\gamma `$.

In Figs. 8 (a) and 8 (b) we plot $`T_c`$ versus $`n`$ for the models of Figs. 2 and 3 for three values of the second-neighbor hopping parameter $`\gamma `$: 0, 0.1, and 0.2. We find that $`T_c`$ has a pronounced maxima as a function of $`n`$. The maximum value of $`T_c`$ increases with the increase of $`\gamma `$ for both $`d_{x^2y^2}+d_{xy}`$ and $`d_{x^2y^2}+id_{xy}`$ waves. A similar increase in the superconducting critical temperature with second-nearest-neighbor hopping term was also observed in a recent study of temperature and filling dependence of the superconducting phase in the Penson-Kolb-Hubbard model.

Finally, we study the interesting ratio $`2\mathrm{\Delta }_{\text{max}}/T_c`$ as function of filling $`n`$ in different cases. The maximum value $`\mathrm{\Delta }_{\text{max}}`$ of $`\mathrm{\Delta }_𝐤`$ was found numerically in the $`𝐤`$ plane and the ratio $`2\mathrm{\Delta }_{\text{max}}/T_c`$ was calculated in different cases as a function of filling $`n`$. As a function of $`n`$, the quantity $`\mathrm{\Delta }_{\text{max}}`$ has a rapid variation (not shown explicitly here) similar to $`T_c`$ of Fig. 8. Although $`T_c`$ as well as $`\mathrm{\Delta }_{\text{max}}`$ are found to exhibit rapid variation as function of $`n`$, the ratio $`2\mathrm{\Delta }_{\text{max}}/T_c`$ is found to be a reasonably smooth function of $`n`$. This ratio varies from 4.0 for $`n=1`$, to a maximum of about 4.7 as $`n`$ decreases; the maximum is obtained for $`n0.4`$. The average value of $`2\mathrm{\Delta }_{\text{max}}/T_c`$ is found to be about 4.3.

## IV Summary

In this work we have studied different properties of the mixed $`d_{x^2y^2}+d_{xy}`$ and $`d_{x^2y^2}+id_{xy}`$ superconducting states in two-dimensional cuprates using a tight-binding model with appropriate lattice symmetry. Recent experiments and other theoretical investigations suggested the formation of such mixed-symmetry states in cuprates on tetragonal and orthorhombic lattice. We have included the effect of orthorhombic distortion of the lattice as well as of second nearest hopping on the order parameter of these states in the tight-binding model. The allowed mixing angles are 0, $`\pi /2,\pi ,`$ and $`3\pi /2`$. We studied the variation of order parameters $`\mathrm{\Delta }_{x^2y^2}`$ and $`\mathrm{\Delta }_{xy}`$ with filling $`n`$ for tetragonal and orthorhombic lattices for different second-neighbor hopping for pure and mixed-symmetry states. These results are reported in Figs. 1 $``$ 5. In the case of mixed-symmetry $`d_{x^2y^2}+d_{xy}`$ wave on orthorhombic lattice in the presence of second-neighbor hopping, the mixing angle changes from 0 to $`\pi `$ as filling $`n`$ is varied. We have also studied the temperature dependence of the order parameter under different situations. This temperature dependence is quite similar to the $`n=1`$ case studied by us previously . We also studied the dependence of the critical temperature $`T_c`$ for the mixed-symmetry state on filling $`n`$ under different situations and find that this dependence exhibits a maximum and rapid variation. The maximum of $`T_c`$ increases with the increase of second-nearest hopping. The ratio $`2\mathrm{\Delta }_{\text{max}}/T_c`$ is a slowly varying function of $`n`$ and varies from 4 for $`n=1`$ to a maximum of about 4.7 as $`n`$ decreases.

We thank Conselho Nacional de Desenvolvimento Científico e Tecnológico and Fundação de Amparo à Pesquisa do Estado de São Paulo for financial support.

Figure Captions:

1. The order parameters (a) $`\mathrm{\Delta }_{x^2y^2}`$ and (b) $`\mathrm{\Delta }_{xy}`$ at different filling $`n`$ for second nearest hopping $`\gamma =0`$ (full line), 0.1 (dashed-dotted line), and 0.2 (dashed line) for square lattice. The potential parameters for (a) are $`V_1=0.73t`$, $`V_2=0`$, and for (b) are $`V_1=0`$, $`V_2=2.1t`$.

2. The order parameters $`\mathrm{\Delta }_{x^2y^2}`$ (dashed line), $`\mathrm{\Delta }_{xy}`$ (full line) for different filling $`n`$ for mixed $`d_{x^2y^2}+d_{xy}`$ wave for second nearest hopping $`\gamma =0`$, 0.05, 0.1, and 0.2 on square lattice. The parameters for the model are $`\beta =1`$, $`V_1=0.73t`$, $`V_2=2.1t`$.

3. Same as in Fig. 2 for mixed $`d_{x^2y^2}+id_{xy}`$ wave. The parameters for the model are $`\beta =1`$, $`V_1=0.73t`$, $`V_2=1.7t`$.

4. Same as in Fig. 2 on orthorhombic lattice. The parameters for the model are $`\beta =0.95`$, $`V_1=0.97t`$, $`V_2=2.1t`$.

5. Same as in Fig. 3 on orthorhombic lattice. The parameters for the model are $`\beta =0.95`$, $`V_1=0.97t`$, $`V_2=1.7t`$.

6. Temperature dependence of the order parameters $`\mathrm{\Delta }_{x^2y^2}`$ and $`\mathrm{\Delta }_{xy}`$ for mixed-symmetry $`d_{x^2y^2}+d_{xy}`$ (dashed line) and $`d_{x^2y^2}+id_{xy}`$ (full line) waves for $`\beta =1`$, $`n=0.8`$, and $`\gamma =0.2`$. The larger order parameters represent the $`d_{x^2y^2}`$ wave. The parameters of the models are the same as in Figs. 2 and 3, respectively.

7. Chemical potential $`\mu `$ versus filling $`n`$ for the four models of Figs. 2 $``$ 5 for $`\gamma =0`$ (full line), 0.05 (short dashed line), 0.1 (long dashed line), and 0.2 (dashed-dotted line). Four different models lead to identical result.

8. Critical temperature $`T_c`$ versus $`n`$ for the models of (a) Fig. 2 and (b) Fig. 3 for $`\gamma =0`$ (full line), 0.1 (short dashed line), and 0.2 (long dashed line).


format: text

text:
# 1 Introduction

## 1 Introduction

In this paper we consider Slater-Jastrow trial wave functions, $`\mathrm{\Psi }=De^J`$, where $`D`$ is a Slater determinant and $`J`$ the Jastrow factor that takes account of the electronic correlations. Such wave functions are employed in variational quantum Monte Carlo (VQMC) simulations of solids ,, where an approximate many-electron ground state is obtained by numerically optimizing an explicit parametrized trial wave function. Remarkably better results are obtained from VQMC calculations using Slater-Jastrow functions than from Hartree-Fock (HF) or density-functional theory within the local density approximation (LDA). For example, cohesive energies obtained using VQMC are typically an order of magnitude more accurate than those obtained employing HF or LDA -.

We use a Slater determinant of LDA orbitals and a Jastrow factor that includes pairwise correlation terms $`u(𝐫_𝐢,𝐫_𝐣)`$. In simulation of solids $`u`$ is usually taken to be both homogeneous and isotropic, i.e., $`u(𝐫_𝐢,𝐫_𝐣)`$ is assumed to depend only on the interelectronic distance $`r_{ij}=|𝐫_𝐢𝐫_𝐣|`$.

The long-range behaviour of $`u(𝐫_𝐢,𝐫_𝐣)`$ is determined by the zero-point motion of the plasmons, both in homogeneous  and inhomogeneous systems  while the short-range behaviour is governed by the so-called cusp condition . Gaskell minimized the variational energy of the uniform electron gas in the random-phase approximation and found a Jastrow function which combined both the short-range and the long-range behaviour. Ceperley used this Jastrow function in his variational Monte Carlo calculations of the spin-unpolarized homogeneous electron gas .

Here we generalize Gaskell’s approach to inhomogeneous systems and derive an inhomogeneous and anisotropic Jastrow factor that depends on the non-interacting structure factor and density of the system. The Jastrow factor is represented as a double Fourier integral. The full interacting many-electron Hamiltonian $`\widehat{H}`$ is expressed as a function of the plane-wave expansion coefficients of the Jastrow function and, within the random phase approximation (RPA) $`\widehat{H}`$ is minimized with respect to these coefficients. This way the approximated ground state energy is obtained, and the trial wavefunction of the system can be constructed from the resulting inhomogeneous and anisotropic Jastrow function. In the homogeneous limit, our approach reduces to that of Gaskell.

The rest of the paper is organized as follows. In Sec.2.1. we describe the formalism used in Gaskell’s original derivation of the Jastrow function for the homogeneous electron gas (HEG). Sec. 2.2 presents our generalization to systems with density inhomogeneities where, for the sake of simplicity, we consider systems with density variations in only one direction. Sec. 3 discusses the numerical results obtained for the IHEG, and Sec.4 concludes.

## 2 Formalism

A Slater-Jastrow trial wave function is the product of a totally antisymmetric Slater determinant $`D`$ and a totally symmetric Jastrow factor $`e^J`$: $`\mathrm{\Psi }=De^J`$. The orbitals used in D are obtained from LDA calculations. Here we consider spin unpolarized systems. The Slater determinant builds in exchange effects but neglect the electronic correlations caused by the Coulomb interactions. The most important correlation effects occur when pairs of electrons approach each other, and these may be included by choosing pairwise Jastrow factors of the form:

$$J=\underset{i,j}{}u(𝐫_i,𝐫_j),$$
(1)
where $`u(𝐫_i,𝐫_j)=u(𝐫_j,𝐫_i)`$. So the trial wave function can be written as:

$$\mathrm{\Psi }(𝐑)=D\mathrm{exp}\left[\underset{i,j}{}u(𝐫_i,𝐫_j)\right].$$
(2)
Note that we include the $`i=j`$ terms in the sum over $`i`$ and $`j`$, which are omitted by many authors as they give a constant contribution for the HEG. (Unless otherwise is stated, we use atomic units throughout, i.e., $`e^2=\mathrm{}=m_e=1`$.)

The full interacting many-electron Hamiltonian is given by

$$\widehat{H}=\widehat{T}+\widehat{V}_{ee}+\widehat{V}_{ext},$$
(3)
where $`\widehat{T}`$ and $`\widehat{V}_{ee}`$ are the kinetic energy and the electron-electron interaction operators respectively and $`\widehat{V}_{ext}`$ is an external potential. The expectation value of the energy with respect to the above trial wavefunction can be written as (the normalized $`\mathrm{\Psi }`$ is taken):

$$E=<\mathrm{\Psi }|H|\mathrm{\Psi }>=\underset{i=1}{\overset{N}{}}ϵ_i+𝑑𝐫n(𝐫)V_{ext}(𝐫)𝑑𝐫n(𝐫)V_s(𝐫)+<T_2>+<V_{ee}>,$$
(4)
where $`V_s`$ is the effective single-particle potential associated with the orbitals in the Slater determinant, $`ϵ_i`$ are the energy eigenvalues which correspond to such one electron orbitals and $`n(𝐫)`$ is the electron density. $`<T_2>`$ and $`<V_{ee}>`$ are defined as:

$$<T_2>=𝑑R\mathrm{\Psi }^{}(𝐑)\mathrm{\Psi }(𝐑)\left[\underset{m=1}{\overset{N}{}}_m(\underset{i,j}{}U(𝐫_i,𝐫_j))_m(\underset{i,j}{}U(𝐫_i,𝐫_j))\right]$$
(5)
and

$$<V_{ee}>=𝑑𝐑\mathrm{\Psi }^{}(𝐑)\mathrm{\Psi }(𝐑)(\frac{1}{2}\underset{ij}{}\frac{1}{|𝐫_i𝐫_j|}).$$
(6)
Here $`𝐑=(𝐫_\mathrm{𝟏},𝐫_\mathrm{𝟐},\mathrm{}𝐫_𝐍)`$ represents all electron coordinates. Several different types of Jastrow factors are commonly used, the best known are based on the RPA in the form due to Gaskell and the form based on the RPA of Bohm and Pines . Here we will briefly describe the Jastrow factor due to Gaskell for the homogeneous electron gas before presenting its generalization to inhomogeneous systems.

### 2.1 The homogeneous RPA Jastrow factor

For homogeneous systems the correlation term $`u`$ depends only on the interelectronic distance $`r_{ij}`$ and may be expressed as follows:

$$u(𝐫_𝐢,𝐫_𝐣)=u(r_{ij})=\underset{k}{}u(k)e^{\mathrm{i}kr_{ij}},$$
(7)
where $`k`$ is the modulus of the total wavevector. Also, the single-particle orbitals are a set of planewaves. Inserting the above expression in Eq.(4), the total energy in the random-phase approximation of $`T_2`$ can be written as :

$$E=\underset{kk_F}{}\frac{1}{2}k^2+2N^2\underset{𝐤}{}k^2u^2(k)S(k)+2\pi n_0\underset{𝐤}{}\frac{1}{k^2}[S(k)1]$$
(8)
where $`NS(k)=<\mathrm{\Psi }|\widehat{n}_𝐤\widehat{n}_𝐤^{}|\mathrm{\Psi }>`$ is the static structure factor of the uniform interacting electron gas and $`\widehat{n}_𝐤=_ie^{\mathrm{i}𝐤𝐫_i}`$ are the number density operators. $`k_F`$ and $`n_0=N/V`$ are the Fermi energy and the average electron densityi, respectively, $`N`$ being the total number of electrons in the system and $`V`$ the volume. Gaskell assumed that the structure factor for the interacting electron gas is close to that of the non-interacting electron gas and, for a small Jastrow function $`u(𝐫)`$ 111This corresponds to the high density limit. can be approximated by the perturbation formula :

$$S(k)=\frac{S_0(k)}{1+4u(k)S_0(k)},$$
(9)
where $`S_0(k)`$ is the structure factor of the non-interacting electron gas . The variational energy in the random-phase approximation (Eq. 3) can then be minimized with respect to $`u(k)`$ to yield:

$$u(k)=\frac{1}{4S_0(k)}+\frac{1}{4}\left[\frac{1}{4(S_0(k))^2}+\frac{16\pi n_0}{k^4}\right]^{1/2}.$$
(10)
### 2.2 The inhomogeneous RPA Jastrow factor

Here we consider the case of an inhomogenous system translationally invariant in the plane perpendicular to the $`z`$ axis. In this case $`u`$ can be written as:

$$u(𝐫_\mathrm{𝟏},𝐫_\mathrm{𝟐})=\underset{𝐪}{}\underset{k_{z_1},k_{z_2}}{}u_{k_{z_1},k_{z_2}}e^{\mathrm{i}𝐪(𝐫_{}\mathrm{𝟏}𝐫_\mathrm{𝟐})}e^{\mathrm{i}k_{z_1}z_1}e^{\mathrm{i}k_{z_2}z_2},$$
(11)
where $`𝐫=(𝐫_{},z)`$ and $`𝐤=(𝐪,k_z)`$. $`u_{k_{z_1},k_{z_2}}`$ are the two-dimensional Fourier coefficients of the correlation function $`u`$:

$$u_{k_{z_1},k_{z_2}}(q)=\frac{1}{V}𝑑𝐫_{}e^{\mathrm{i}𝐪𝐫_{}}𝑑z_1𝑑z_2e^{\mathrm{i}k_{z_1}z_1}e^{\mathrm{i}k_{z_2}z_2}u(𝐫_{},z_1,z_2).$$
(12)
We use the following definition of the number density operator: $`\widehat{n}_{𝐪,k_z}=_ie^{i𝐪𝐫_{}^{}{}_{i}{}^{}}e^{ik_zz_i}`$, so that the Slater-Jastrow trial wave function is given by:

$$\mathrm{\Psi }=D\mathrm{exp}\left[\underset{𝐪}{}\underset{k_{z_1},k_{z_2}}{}u_{k_{z_1},k_{z_2}}(𝐪)\widehat{n}_{𝐪,k_{z_1}}\widehat{n}_{𝐪k_{z_2}}^{}\right].$$
(13)
#### 2.2.1 Deriving an approximation for the inhomogeneous structure factor

In this section we develop an approximation based on the RPA to obtain the interacting structure factor in terms of Fourier coefficients of the above Jastrow factor. Let us define:

$$\mathrm{\Psi }(\lambda )=D\mathrm{exp}\left[\lambda \underset{𝐪}{}\underset{k_{z_1},k_{z_2}}{}u_{k_{z_1},k_{z_2}}(𝐪)\widehat{n}_{𝐪,k_{z_1}}\widehat{n}_{𝐪k_{z_2}}^{}\right].$$
(14)
The expectation value we need to calculate is of the form:

$$I(\lambda )\frac{<\mathrm{\Psi }(\lambda )|\widehat{n}_{𝐪,k_{z_1}}\widehat{n}_{𝐪k_{z_2}}^{}|\mathrm{\Psi }(\lambda )>}{<\mathrm{\Psi }(\lambda )|\mathrm{\Psi }(\lambda )>}.$$
(15)
with $`\lambda =1`$. With the above trial function this becomes

$$I(\lambda )=\frac{DD^{}\widehat{n}_{𝐪_\mathrm{𝟏},k_{z_1}}\widehat{n}_{𝐪_\mathrm{𝟐}k_{z_2}}^{}\mathrm{exp}\left[2\lambda _{𝐪_\mathrm{𝟐}}_{k_{z_3},k_{z_4}}u_{k_{z_3},k_{z_4}}(𝐪_\mathrm{𝟐})\widehat{n}_{𝐪_\mathrm{𝟐},k_{z_3}}\widehat{n}_{𝐪_\mathrm{𝟐}k_{z_4}}^{}\right]}{DD^{}\mathrm{exp}\left[2\lambda _{𝐪_\mathrm{𝟐}}_{k_{z_3},k_{z_4}}u_{k_{z_3},k_{z_4}}(𝐪_\mathrm{𝟐})\widehat{n}_{𝐪_\mathrm{𝟐},k_{z_3}}\widehat{n}_{𝐪_\mathrm{𝟐}k_{z_4}}^{}\right]}.$$
(16)
Following Gaskell we use the identity:

$$\frac{1}{I(\lambda )}\frac{1}{I(0)}=_0^\lambda \frac{}{\lambda ^{}}\left[\frac{1}{I(\lambda ^{})}\right]𝑑\lambda ^{}$$
(17)
to obtain a perturbative approximation for $`I(1)`$. Taking the derivative inside the integral with respect to $`\lambda `$ and expanding in a Taylor series the leading term is

$$\begin{array}{cc}\frac{}{\lambda ^{}}(\frac{1}{I(\lambda ^{})})|_{\lambda =0}=\hfill & 4u_{k_{z_1},k_{z_2}}(𝐪_\mathrm{𝟏})\left[\frac{<0|(\widehat{n}_{𝐪_\mathrm{𝟏},k_{z_1}}\widehat{n}_{𝐪_\mathrm{𝟏}k_{z_2}}^{})^2|0><0|\widehat{n}_{𝐪_\mathrm{𝟏},k_{z_1}}\widehat{n}_{𝐪_\mathrm{𝟏}k_{z_2}}^{}|0>^2}{<0|\widehat{n}_{𝐪_\mathrm{𝟏},k_{z_1}}\widehat{n}_{𝐪_\mathrm{𝟏}k_{z_2}}^{}|0>^2}\right]\hfill \\ & \\ & +2_{𝐪_\mathrm{𝟏},k_{z_3},k_{z_4}(𝐪_\mathrm{𝟐}\pm 𝐪_\mathrm{𝟏},k_{z_3}k_{z_1},k_{z_4}k_{z_2})}u_{k_{z_3},k_{z_4}}(𝐪_\mathrm{𝟐})\hfill \\ & \\ & \times \left[\frac{<0|\widehat{n}_{𝐪_\mathrm{𝟐},k_{z_3}}\widehat{n}_{𝐪_\mathrm{𝟐}k_{z_4}}^{}\widehat{n}_{𝐪_\mathrm{𝟏},k_{z_1}}\widehat{n}_{𝐪_\mathrm{𝟏}k_{z_2}}^{}|0>}{<0|\widehat{n}_{𝐪_\mathrm{𝟏},k_{z_1}}\widehat{n}_{𝐪_\mathrm{𝟏}k_{z_2}}^{}|0>^2}\frac{<0|\widehat{n}_{𝐪_\mathrm{𝟐},k_{z_3}}\widehat{n}_{𝐪_\mathrm{𝟐}k_{z_4}}^{}|0><0|(\widehat{n}_{𝐪_\mathrm{𝟏},k_{z_1}}\widehat{n}_{𝐪_\mathrm{𝟏}k_{z_2}}^{}|0>)}{<0|\widehat{n}_{𝐪_\mathrm{𝟏},k_{z_1}}\widehat{n}_{𝐪_\mathrm{𝟏}k_{z_2}}^{}|0>^2}\right]\hfill \end{array}$$
(18)
where $`<0|\widehat{A}|0>`$ denotes the expectation value of an operator $`\widehat{A}`$ with respect to $`\mathrm{\Psi }(0)=D`$. Following Gaskell, we use the following approximation:

$$<0|(\widehat{n}_{𝐪,k_{z_1}}\widehat{n}_{𝐪k_{z_2}}^{})^2|0>=2[<0|(\widehat{n}_{𝐪,k_{z_1}}\widehat{n}_{𝐪k_{z_2}}^{})|0>]^2+(1/N),$$
(19)
and ignore correlations between different $`\widehat{n}^{}s`$ in the second term. The resulting approximation for $`I(\lambda )`$ becomes

$$I(\lambda )=\frac{I(0)}{1+\lambda I(0)4u_{k_{z_1},k_{z_2}}(𝐪)}.$$
(20)
The corresponding approximation for the Fourier coeffecients of the structure factor is obtained form the above equation for $`\lambda =1`$, yielding

$$S_{k_{z_1},k_{z_2}}(𝐪)=\frac{S_{k_{z_1},k_{z_2}}^0(𝐪)}{1+4u_{k_{z_1},k_{z_2}}(𝐪)S_{k_{z_1},k_{z_2}}^0(𝐪)}.$$
(21)
#### 2.2.2 Expectation value of the energy

In order to obtain the expression for the total energy of an inhomogeneous system we use the RPA approximation

$$<\widehat{n}_{𝐪_1,k_{z_1}}\widehat{n}_{𝐪_2,k_{z_2}}\widehat{n}_{𝐪_1+𝐪_2,k_{z_3}}>=<\widehat{n}_{𝐪_1,k_{z_1}}\widehat{n}_{𝐪_2,k_{z_2}}><\widehat{n}_{𝐪_1+𝐪_2,k_{z_3}}>.$$
(22)
In the case of a system with density variation along $`z`$ we have:

$$<\widehat{n}_{𝐪_1+𝐪_2,k_z}>=\{\begin{array}{cc}V\rho _{k_z}\hfill & \text{ if }𝐪_1+𝐪_2=0\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(23)
where $`\rho _{k_z}`$ represent the Fourier coefficients of the electronic density.

Making these assumptions and employing the expression for the structure factor obtained in the previous section, the terms $`<T_2>`$ and $`<V_{ee}>`$ which appear in the expression for the total energy Eq.(4) can be written as:

$$\begin{array}{cc}\hfill <T_2>=& 2VN_{𝐪_1,k_{z_1},k_{z_2},k_{z_3},k_{z4}}\widehat{u}_{k_{z_1},k_{z_2}}(q_1)\widehat{u}_{k_{z_3},k_{z4}}(q_1)\left[(q_1)^2k_{z_1}k_{z_3}\right]\hfill \\ & \\ & \times \rho _{k_{z_1}+k_{z_3}}\left(S_{k_{z_2},k_{z_4}}(q_1)\right)\hfill \\ & \\ \hfill <V_{ee}>=& n_0_{𝐪,k_z}\frac{2\pi }{(q^2+k_z^2)^2}\left[\left(S_{k_z,k_z}(q)\right)1\right].\hfill \end{array}$$
(24)
The above equation combined with our approximation (Eq.21) for the structure factor leads to an expression for the total energy in terms of the Fourier coefficients of the $`u`$ function, the electron density and the non-interacting structure factor. We take both the density and the non-interacting structure factor from LDA calculations and minimize $`E`$ with respect to the $`\widehat{u}_{k_z,k_z^{}}(q)`$ coefficients, which then can be used to construct the trial wave functions for VQMC calculations.

## 3 Results

This section presents the numerical results obtained for uniform and non-uniform spin-unpolarized $`N`$electron systems ($`N=38`$) contained in a simple cubic (SC) simulation cell and satisfying periodic boundary conditions within the cell. Density modulations in $`z`$ are induced by subjecting the system to and external potential of the from

$$V_{ext}(𝐫)=V_0cos(qz)$$
(25)
with $`V_0=0.1`$ a.u.. $`q=0.6k_F`$ is a reciprocal lattice vector of the SC simulation cell. The average density of the system corresponds to $`r_s=2`$. The electron density is uniform in the $`xy`$ plane and varies in th $`z`$ direction.

The inhomogeneous $`u(𝐫_\mathrm{𝟏},𝐫_\mathrm{𝟐})`$ term is represented by fixing the position $`𝐫_\mathrm{𝟏}`$ of the first electron, while moving the second electron, $`𝐫_\mathrm{𝟐}`$, along the $`z`$ direction on a line passing through $`𝐫_\mathrm{𝟏}`$.

### 3.1 Homogeneous system

In Eq.(7) we saw that $`u(z_1,z_2)`$ may be expanded as a Fourier series. In Fig. 1 the convergence of the energy with respect to the cutoff wavevector, $`k_{cutoff}`$, used in the Fourier expansion of the $`u`$ factor is shown. The typical cutoff needed is small enough to be computationally feasible, so that such a representation for $`u(z_1,z_2)`$ is found to be very useful.

Looking at the expression for the energy we see that for $`u(k)=0`$, $`<T_2>=0`$ and $`<V_{ee}>`$ is just the exchange energy of a HEG. For non-zero $`u(k)`$ factors $`<T_2>+<V_{ee}>`$ includes both exchange and correlation energies. The correlation energy is defined as:

$$E_c=[<T_2>+<V_{ee}>][<T_2>+<V_{ee}>]_{u(k)=0}$$
(26)
Table 1 presents the values for the correlation energy obtained in this work in comparison with those obtained by Gaskell . Differences due to finite size effects have to be taken into account, as we have considered a very small system with only $`N=38`$ electrons.

The results for the energy obtained with the uniform system previously described and using the homogeneous RPA u factor given in Section. 2.1. are presented in Table 2.

### 3.2 Inhomogeneous systems

The inhomogeneous RPA jastrow factor considered here is the one derived in Section. 2.2. The energies derived are presented in Table 2. Exchange-correlation energy is defined as:

$$E_{xc}=<T_2>+<V_{ee}>E_H,$$
(27)
where $`E_H`$ represents Hartree energy.

In Fig.2a we plot the inhomogeneous RPA $`u`$ term, for the three different positions of $`𝐫_\mathrm{𝟏}`$ indicated in Fig. 2b which presents the LDA density and the average electron density. The homogeneous Jastrow factor is shown as well, which depends only on the relative distance between the electrons and not their individual coordinates. Anisotropy of the Jastrow factor can be observed. This is such that the $`u`$ function is stronger on the side where the density is lower. This can be explained by the fact that the RPA screening is more effective where the electron density is high, so that the $`u`$ function is weaker on the high density side .

## 4 Conclusions

We generalize the RPA treatment of the homogeneous electron gas by Gaskell to the inhomogeneous case. As a result we obtain a Slater-Jastrow trial wave functions containing an anisotropic and inhomogeneous Jastrow factor. The latter may be obtained from the non-interacting structure factor and density of the system.

Results for the energy and the inhomogeneous correlation $`u`$ factor using an electronic system subject to a sinusoidal external potential are presented. We compare these results with those obtained using the homogeneous $`u`$ factor by Gaskell, and we clearly observe the anisotropy introduced by the inhomogeneous $`u`$ factor. Future work will focus on testing the performance of our inhomogeneous Jastrow factor using variational Monte Carlo and in further refining our approximation of the interacting structure factor.

## 5 Acknowledgements

We acknowledge partial support by the Basque Unibertsitate eta Ikerketa Saila.


text:
# Andreev levels in a single-channel conductor

## I Introduction

Several recent works have identified and studied deviations from mean-field theory in the sub-gap density of states of a normal metal in contact with a superconductor. The excitation spectrum below the gap of the bulk superconductor consists of a coherent superposition of electron and hole excitations, coupled by Andreev reflection at the normal-metal–superconductor (NS) interface. The energy of these Andreev levels fluctuates from sample to sample, but such mesoscopic fluctuations are ignored in mean-field theory. Because of these fluctuations, the ensemble averaged density of states $`\nu (\epsilon )`$ acquires a tail that extends below the mean-field gap, vanishing only at the Fermi level (zero excitation energy $`\epsilon `$). The fluctuations become particularly large if the size of the normal metal is greater than the localization length.

The purpose of this paper is to analyze an extreme case of complete breakdown of mean-field theory, which is still sufficiently simple that it can be solved exactly. This is the case of single-mode conduction through a disordered normal-metal wire attached to a superconductor. The localization length in this geometry is equal to the elastic mean free path $`\mathrm{}`$, so that the wire crosses over with increasing length $`L`$ from the ballistic regime directly into the localized regime — without an intermediate diffusive regime. Perturbation theory is possible in the quasiballistic regime $`\mathrm{}L`$, but for $`\mathrm{}<L`$ an essentially non-perturbative approach is required. We will use an approach based on a scaling equation (also known as invariant embedding), that has proved its use before in different contexts.

We will contrast the quasiballistic and localized regimes, as well as the two geometries with a single superconducting contact (NS junction) or with two superconducting contacts at both ends of the normal metal wire (SNS junction). If we assume that the two superconductors have the same phase, so that there is no supercurrent flowing through the normal metal, then the Andreev levels of the SNS junction are doubly degenerate in the absence of disorder. This degeneracy is broken by disorder. We find that for weak disorder the probability distribution of the splitting is given precisely by Wigner’s surmise from random-matrix theory. (The spectra of chaotic systems have spacings described by Gaudin’s distribution, which is close to, but not identical with Wigner’s surmise.)

In the localized regime the fluctuations of the Andreev levels become greater than their spacing, and they can no longer be distinguished in the mean density of states, which decreases smoothly to zero on approaching the Fermi level. The energy scale for this soft gap is exponentially small because of localization, given by $`\epsilon _g=(\mathrm{}v_F/\mathrm{})e^{L/\mathrm{}}`$. The decay of $`\nu (\epsilon )`$ for $`\epsilon \epsilon _g`$ has a log-normal form $`\mathrm{exp}\left[\frac{\mathrm{}}{4L}\mathrm{ln}^2\left(\epsilon _g/\epsilon \right)\right]`$. Such log-normal tails are characteristic of rare fluctuations in the localized regime and have appeared recently in the context of the superconductor proximity effect.

## II Quasiballistic regime

### A NS junction

The NS junction consists of a piece of normal metal of length $`L`$ connected at one end to a superconductor (see Fig. 1a). The width of the normal metal is of the order of the Fermi wave length $`\lambda _F`$, such that there is a single propagating mode at the Fermi energy $`E_F`$. We assume an ideal junction, without any tunnel barrier and with $`E_F`$ much greater than the superconducting gap $`\mathrm{\Delta }_0`$. An electron incident on the superconductor with energy $`\epsilon <\mathrm{\Delta }_0`$ above the Fermi level is then Andreev reflected as a hole at energy $`\epsilon `$ below the Fermi level, with the phase shift

$$\varphi _A=\mathrm{arccos}(\epsilon /\mathrm{\Delta }_0),\pi /2<\varphi _A<0.$$
(1)
We wish to know at which $`\epsilon `$ a bound state (Andreev level) will form in the normal metal.

The electron and hole components of the wave function $`\psi (x)=\mathbf{(}u(x),v(x)\mathbf{)}`$ satisfy the Bogoliubov-de Gennes (BdG) equation

$$\left(\begin{array}{cc}_0& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& _0^{}\end{array}\right)\psi =\epsilon \psi ,$$
(2)
where $`_0=(\mathrm{}^2/2m)^2/x^2+V(x)`$ is the Hamiltonian of the normal metal (with disorder potential $`V`$) and $`\mathrm{\Delta }(x)=\mathrm{\Delta }_0\theta (x)`$ is the superconducting gap (which vanishes in the normal-metal region $`x>0`$). For narrow junctions (width much less than the superconducting coherence length $`\xi _0=\mathrm{}v_F/\mathrm{\Delta }_0`$) the depletion of $`\mathrm{\Delta }(x)`$ on the superconducting side may be neglected, hence the step function $`\theta (x)`$. At the closed end $`x=L`$ of the normal metal we impose the boundary condition $`\psi (L)=0`$.

In this section we address the quasiballistic regime of mean free path $`\mathrm{}L`$. We can then treat $`V`$ as a small perturbation on the ballistic bound states

$`\psi (x)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{Z}}}\left(\begin{array}{c}\mathrm{sin}[(k_F+k)(xL)]\\ \mathrm{sin}[(k_Fk)(xL)\pi n]\end{array}\right),0<x<L,`$ (6)
$`\psi (x)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{Z}}}\left(\begin{array}{c}\mathrm{sin}[k_Fx(k_F+k)L]\\ \mathrm{sin}[k_Fx(k_Fk)L\pi n]\end{array}\right)\mathrm{exp}\left({\displaystyle \frac{x}{\xi _0}}\mathrm{sin}\varphi _A\right),x<0.`$ (9)
The normalisation constant is $`Z=L\frac{1}{2}\xi _0/\mathrm{sin}\varphi _A`$ for $`k_FL1`$. (We denote $`k_F=mv_F/\mathrm{}=2\pi /\lambda _F`$.) The wave number $`k=\epsilon /\mathrm{}v_F`$ should satisfy the quantization condition

$$2kL+\varphi _A=\pi n,n=0,1,2,\mathrm{}.$$
(10)
The total number of Andreev levels within the gap is $`2L/\pi \xi _0`$ for $`L\xi _0`$. (There remains one level if $`L\xi _0`$.)

To first order in $`V`$ the energy level is shifted by the matrix element

$$\delta \epsilon =_0^L𝑑xV(x)\left[u(x)^2v(x)^2\right].$$
(11)
We assume a potential with a short-range correlation, expressed by

$$V(x)=0,V(x)V(x^{})=\frac{\mathrm{}^2v_F^2}{\mathrm{}}\delta (xx^{}),$$
(12)
where $`\mathrm{}`$ stands for the disorder average. It follows that the distribution of an Andreev level around its ballistic value is a Gaussian with zero mean, $`\delta \epsilon =0`$, and variance

$$\delta \epsilon ^2=\frac{\mathrm{}^2v_F^2\left(2L+\xi _0\mathrm{sin}\varphi _A\right)}{2\mathrm{}\left(2L\xi _0/\mathrm{sin}\varphi _A\right)^2}.$$
(13)
By way of illustration, we show in Fig. 2 the mean density of states of an NS junction containing three Andreev levels ($`\xi _0/L=0.24`$) with mean free path $`\mathrm{}=12L`$. The Gaussian given by Eq. (13) agrees very well with the numerical solution of the BdG equation (data points).

We briefly explain the numerical method. The BdG equation is solved numerically on a one-dimensional grid (lattice constant $`a`$) by replacing the Laplacian by finite differences and truncating the Hamiltonian matrix in the superconducting region, where the wave function is evanescent for energies in the superconducting gap. The resulting tight-binding model has nearest-neighbor coupling $`\gamma =\mathrm{}^2/2ma^2`$ (band width $`4\gamma `$). We set $`E_F=\gamma `$ and $`\mathrm{\Delta }_0=0.1\gamma `$, corresponding to $`\lambda _\mathrm{F}=6a`$ and $`\xi _0=10\sqrt{3}a`$. The disorder is modelled by a random on-site potential which is uniformly distributed in the interval $`(W,W)`$. The mean free path from the Born approximation, $`l=3E_F(4\gamma E_F)a/W^2`$, was found to fit well to the prediction of one-dimensional scaling theory for the mean inverse transmission probability, $`T^1=\frac{1}{2}[1+\mathrm{exp}(2L/l)]`$, in the complete range from the quasiballistic to the localized regime. (The localization length $`\xi `$ is related to the mean free path by $`\xi =2\mathrm{}`$, cf. Ref. .) This allows for a parameter-free comparison of the analytical and numerical results for the ensemble-averaged density of states.

### B SNS junction

The quasiballistic regime in an SNS junction (Fig. 1b) is qualitatively different from the NS case of the previous subsection. The reason is the double degeneracy of the unperturbed Andreev levels. This degeneracy exists if the phase of the order parameter in the two superconductors is the same, which is what we assume in this paper. Let us examine the splitting of the Andreev levels by the disorder potential.

The SNS junction has energy gap

$$\mathrm{\Delta }(x)=\mathrm{\Delta }_0\theta (x)+\mathrm{\Delta }_0\theta (xL).$$
(14)
The quantization condition reads

$$kL+\varphi _A=n\pi ,n=0,1,2,\mathrm{}.$$
(15)
There are $`L/\pi \xi _0`$ Andreev levels (for $`L\xi _0`$), each level being doubly degenerate. We choose the two independent eigenfunctions $`\psi _\pm (x)`$ such that they carry zero current. They are given by

$`\psi _+(x)={\displaystyle \frac{1}{\sqrt{Z^{}}}}\left(\begin{array}{c}\mathrm{cos}(k_Fx)\\ \mathrm{cos}(k_Fx\varphi _A)\end{array}\right)\mathrm{exp}\left({\displaystyle \frac{x}{\xi _0}}\mathrm{sin}\varphi _A\right),x<0,`$ (19)
$`\psi _+(x)={\displaystyle \frac{1}{\sqrt{Z^{}}}}\left(\begin{array}{c}\mathrm{cos}[(k_F+k)x]\\ \mathrm{cos}[(k_Fk)x\varphi _A]\end{array}\right),0<x<L,`$ (22)
$`\psi _+(x)={\displaystyle \frac{1}{\sqrt{Z^{}}}}\left(\begin{array}{c}\mathrm{cos}(k_Fx+kL)\\ \mathrm{cos}(k_Fx+\pi n)\end{array}\right)\mathrm{exp}\left({\displaystyle \frac{xL}{\xi _0}}\mathrm{sin}\varphi _A\right),x>L,`$ (25)
and $`\psi _{}(x)`$ is obtained by replacing cosine by sine. The normalization constant is now $`Z^{}=L\xi _0/\mathrm{sin}\varphi _A`$.

To first order in $`V`$ the levels are splitted symmetrically around the ballistic value, by an amount $`\pm \frac{1}{2}s`$. The basis (II B) is chosen in such a way that the off-diagonal elements of the perturbation vanish. The shift of each level can then be calculated from Eq. (11) using the corresponding eigenfunction. We again calculate the probability distribution $`P(s)`$ of the level splitting using Eq. (12). The result is

$$P(s)=\frac{\pi s}{2s^2}\mathrm{exp}\left(\frac{\pi s^2}{4s^2}\right),$$
(26)
with average splitting

$$s=\mathrm{\Delta }\sqrt{\frac{\pi }{2\mathrm{}}}\frac{\xi _0\sqrt{L+\xi _0\mathrm{sin}\varphi _A}}{L\xi _0/\mathrm{sin}\varphi _A}.$$
(27)
We recognize Eq. (26) as Wigner’s surmise of random matrix theory.

In Fig. 3 we compare Eq. (26) with numerical data. The agreement is excellent for a range of mean free paths in the quasiballistic regime. The mean position of the splitted levels fluctuates only to higher orders in $`L/\mathrm{}`$. This makes it possible to resolve the splitting in the mean density of states (see inset in Fig. 3).

## III Localized regime

### A NS junction

In order to go beyond the quasiballistic regime into the localized regime $`L\mathrm{}`$ we write the quantization condition for the Andreev levels in an NS junction in the form

$$r(\epsilon )r(\epsilon )^{}e^{2i\varphi _A}=1,$$
(28)
where $`r(\epsilon )=e^{i\varphi (\epsilon )}`$ is the reflection amplitude of the disordered normal metal. \[The hole has reflection amplitude $`r^{}(\epsilon )`$.\] In terms of the phase shifts we have

$$\mathrm{\Phi }(\epsilon )\frac{\varphi (\epsilon )\varphi (\epsilon )}{2}+\varphi _A(\epsilon )=\pi n,n=0,1,\mathrm{}.$$
(29)
The density of states $`\nu (\epsilon )=_n\delta (\epsilon \epsilon _n)`$ is related to the scattering phase shifts by

$$\nu (\epsilon )=\frac{1}{\pi }\frac{d}{d\epsilon }\mathrm{Im}\mathrm{ln}\mathrm{sin}\mathrm{\Phi }(\epsilon +i0^+),$$
(30)
where $`0^+`$ denotes a positive infinitesimal. The imaginary part of the logarithm jumps by $`\pi `$ whenever $`\mathrm{sin}\mathrm{\Phi }(\epsilon )`$ changes sign, hence it counts the number of levels below $`\epsilon `$. The derivative with respect to $`\epsilon `$ then gives the density of states. It is convenient to write Eq. (30) as a Taylor series,

$$\nu (\epsilon )=\frac{1}{\pi }\frac{d}{d\epsilon }\left(\mathrm{\Phi }+\mathrm{Im}\underset{m=1}{\overset{\mathrm{}}{}}\frac{1}{m}e^{2im\mathrm{\Phi }}\right),$$
(31)
which converges because $`\mathrm{\Phi }(\epsilon +i0^+)`$ is equivalent to $`\mathrm{\Phi }(\epsilon )+i0^+`$.

We seek the disorder averaged density of states $`\nu (\epsilon )`$. One way to proceed is by means of the Berezinskii technique. An alternative way, that we will follow here, is to start from the scaling equation for the probability distribution $`P(\varphi _N)`$ of the phase shift $`\varphi _N=\frac{1}{2}[\varphi (\epsilon )\varphi (\epsilon )]`$. This equation has the form

$$\frac{P}{L}=\frac{}{\varphi _N}\left(\frac{2\epsilon }{\mathrm{}v_F}+\frac{1}{\mathrm{}}\frac{}{\varphi _N}\mathrm{sin}^2\varphi _N\right)P.$$
(32)
The initial condition is $`lim_{L0}P(\varphi _N)=\delta (\varphi _N)`$.

The first moment satisfies $`\varphi _N/L=2\epsilon /\mathrm{}v_F`$, hence

$$\varphi _N=\frac{2\epsilon L}{\mathrm{}v_F}.$$
(33)
Multiplication of Eq. (32) by $`\mathrm{exp}(2im\varphi _N)`$ and integration over $`\varphi _N`$ from $`0`$ to $`\pi `$ yields a set of recursive differential equations for the moments $`R_m=e^{2im\varphi _N}`$,

$$\frac{R_m}{L}=\frac{m^2}{\mathrm{}}(R_{m+1}+R_{m1}2R_m)+\frac{4i\epsilon }{\mathrm{}v_F}mR_m,$$
(34)
with the initial condition $`R_m(0)=1`$. We solve this set of equations by truncating the vector $`(R_1,R_2,\mathrm{}R_M)`$ at a sufficiently large value of $`M400`$ and diagonalizing the corresponding tri-diagonal matrix. From Eq. (31) we then find the mean density of states.

The result is shown in Fig. 4 for $`\xi _0/L=0.24`$ and ratios $`\mathrm{}/L`$ ranging from the quasiballistic regime to the localized regime. Agreement with the numerical solution of the BdG equation is excellent over the whole range.

In the localized regime $`L\mathrm{}`$ the individual Andreev levels can no longer be distinguished in the mean density of states, because the broadening of the levels becomes greater than the spacing. In this regime we distinguish two energy ranges, $`\epsilon \epsilon _g`$ and $`\epsilon \epsilon _g`$, where $`\epsilon _g=(\mathrm{}v_F/\mathrm{})e^{L/\mathrm{}}`$.

For energies higher than $`\epsilon _g`$ we may use the $`L\mathrm{}`$ limit of the distribution $`P(\varphi _N)`$, obtained by setting the left-hand-side of Eq. (32) equal to zero. The resulting moments are

$$\underset{L\mathrm{}}{lim}R_m=\underset{0}{\overset{\mathrm{}}{}}𝑑\sigma e^\sigma \left(\frac{\sigma }{\sigma i\omega }\right)^m,\omega =\frac{4\epsilon \mathrm{}}{\mathrm{}v_F}.$$
(35)
We then calculate the mean density of states from Eq. (31), with the result

$`\nu (\epsilon )={\displaystyle \frac{2L}{\pi \mathrm{}v_F}}+{\displaystyle \frac{1}{\pi \sqrt{\mathrm{\Delta }_0^2\epsilon ^2}}}+f(\epsilon ),\epsilon \epsilon _g,`$ (36)
$`f(\epsilon )={\displaystyle \frac{}{\epsilon }}\mathrm{Im}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d\sigma }{\pi }}\left[{\displaystyle \frac{e^\sigma }{\sigma i\omega }}{\displaystyle \frac{e^\sigma (1e^{2i\varphi _A})}{\sigma (1e^{2i\varphi _A})i\omega }}\right].`$ (37)
The first term on the right-hand-side of Eq. (36) is the energy independent density of states $`\nu _0`$ in an isolated normal metal. The main effect of the superconductor for $`\epsilon \epsilon _g`$ is an enhancement of the density of states close to the gap $`\mathrm{\Delta }_0`$ of the bulk superconductor (second term). The third term is negative for sufficiently small $`\epsilon `$ and is a precursor of the soft gap near the Fermi level. For $`\xi _0\mathrm{}`$ and $`\epsilon \mathrm{}v_F/\mathrm{}`$ the reduction term $`f(\epsilon )`$ can be simplified as

$$f(\epsilon )=\frac{2\mathrm{}}{\pi \mathrm{}v_F}\left(\mathrm{ln}\frac{\mathrm{}v_F}{8\epsilon \mathrm{}}\gamma \right),\epsilon _g\epsilon \frac{\mathrm{}v_F}{\mathrm{}},$$
(38)
where $`\gamma 0.58`$ is Euler’s constant.

Near the Fermi level, for $`\epsilon \epsilon _g`$, the mean density of states vanishes as a result of the proximity to the superconductor. This “soft gap” appears no matter how strongly localized the normal metal is. The coefficients $`R_m`$ may now be treated as analytical functions of the parameter

$$z=\frac{4i\epsilon \mathrm{}m}{\mathrm{}v_F},R_m=R(z).$$
(39)
Taking the limit $`\epsilon 0`$ we deduce from Eq. (32) the partial differential equation

$$\mathrm{}\frac{R}{L}=z^2\frac{^2R}{z^2}zR,$$
(40)
with initial condition $`lim_{L0}R(z)=1`$. This differential equation has been studied before in the theory of one-dimensional localization, but not in connection with the proximity effect. The result for the mean density of states, derived in the Appendix, is given by

$$\nu (\epsilon )=\frac{2\mathrm{}}{\pi ^{3/2}\mathrm{}v_F}\mathrm{exp}\left[\frac{\mathrm{}}{4L}\mathrm{ln}^2\frac{\pi \epsilon _g}{\epsilon }\frac{u\mathrm{}}{2L}\left(\mathrm{ln}\frac{u\mathrm{}}{2L}1\right)\right],$$
(41)
where $`u=\mathrm{ln}\pi \mathrm{}v_F/\epsilon \mathrm{}=\mathrm{ln}\pi \epsilon _g/\epsilon +L/\mathrm{}`$. The leading logarithmic asymptotic of this expression in the limit $`\epsilon \epsilon _g`$ has the log-normal tail

$$\nu (\epsilon )\mathrm{exp}\left[\frac{\mathrm{}}{4L}\mathrm{ln}^2\frac{\pi \epsilon _g}{\epsilon }\right],\epsilon \epsilon _g.$$
(42)
The same log-normal tail was found in Ref.  for a many-channel diffusive conductor. In that case the factor $`\mathrm{}/L`$ is replaced by the Drude conductance of the normal metal and the energy scale $`\epsilon _g`$ is replaced by the Thouless energy $`\mathrm{}D/L^2`$ (with $`D`$ the diffusion constant). In our single-channel localized conductor nether the Drude conductance nor the Thouless energy play a role.

### B SNS junction

In contrast to the quasiballistic regime, the NS and SNS junctions are similar in the localized regime. (At least for the case of zero current through the SNS junction considered here.) Unfortunately, there exists no simple scaling equation as Eq. (32) that can describe the density of states of the SNS junction. We therefore rely on the numerical solution of the BdG equation. In Fig. 5 we show that the mean density of states of an NS junction of length $`L`$ is close to that of an SNS junction of length $`2L`$. This factor of $`2`$ has an obvious explanation in the ballistic regime \[compare Eqs. (10) and (15)\], but it is remarkable that it still applies to the localized regime.

## IV Conclusion

In summary, we have calculated the effect of disorder on the spectrum of Andreev levels in single-channel NS and SNS junctions. The non-perturbative effects of localization in the one-dimensional case can be studied exactly, at least in the NS geometry. Our research is of theoretical interest in view of recent studies of the sub-gap density of states beyond mean field theory, but may also be of experimental interest in view of recent progress made in superconductor–carbon-nanotube devices.

###### Acknowledgements.

We thank Piet Brouwer for a crucial discussion at the initial stage of this project. This research was supported by the “Nederlandse organisatie voor Wetenschappelijk Onderzoek” (NWO) and by the “Stichting voor Fundamenteel Onderzoek der Materie” (FOM). M.T. and N.A.M. thank the visitors program at the Max-Planck-Institut für Physik komplexer Systeme, Dresden. N.A.M. also acknowledges support by the “Ingeniørvidenskabelig Fond og G.A. Hagemanns Mindefond”.

## A Derivation of the log-normal tail

The differential operator on the right-hand-side of Eq. (40) has eigenfunctions

$$f_p(z)=2\sqrt{z}K_p(2\sqrt{z}),$$
(A1)
where $`K_p(z)`$ is the modified Bessel function, such that

$$\left(z^2\frac{^2}{z^2}z\right)f_p(z)=\frac{p^21}{4}f_p(z).$$
(A2)
The solution to Eq. (40) with the initial condition $`lim_{L0}R(z)=1`$ is

$$R(z)=f_1(z)+\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\nu \frac{\nu \mathrm{sinh}(\pi \nu /2)}{\pi (\nu ^2+1)}f_{i\nu }(z)e^{(\nu ^2+1)L/4\mathrm{}}.$$
(A3)
To obtain the density of states of the NS junction it is convenient to define the inverse Laplace transform

$$F(\lambda )=\frac{1}{2\pi i}\underset{i\mathrm{}+0^+}{\overset{i\mathrm{}+0^+}{}}\frac{dz}{(4\lambda )^2}R(z)\mathrm{exp}\left(\frac{z}{4\lambda }\right).$$
(A4)
From Eq. (31) we find for $`\epsilon \epsilon _g`$ the mean density of states in terms of the function $`F`$,

$$\nu (\epsilon )=\frac{4\mathrm{}}{\pi \mathrm{}v_F}F\left(\frac{\epsilon \mathrm{}}{\pi \mathrm{}v_F}\right).$$
(A5)
Our aim is to find the asymptotic form of $`F(\lambda )`$ in the limit $`\lambda 0`$. The inverse Laplace transform of the modified Bessel functions in Eq. (A3) can be found in Ref. . We obtain

$`F(\lambda )=F_0(\lambda ){\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑\nu \lambda ^{(i\nu +1)/2}e^{(\nu ^2+1)L/4\mathrm{}}`$ (A6)
$`\times {\displaystyle \frac{{}_{1}{}^{}F_{1}^{}(\frac{3}{2}+\frac{i\nu }{2},1+i\nu ,4\lambda )}{2\sqrt{\pi }(1i\nu )\mathrm{\Gamma }\left(i\nu /2\right)}},`$ (A7)
where $`F_0(\lambda )=\mathrm{exp}(4\lambda )`$. The integrand has a single pole $`\nu =i`$ in the lower half of the complex plane and the residue from this pole cancels the term $`F_0`$. Let us shift the contour by the transformation $`\nu \nu (i\mathrm{}/L)\mathrm{ln}(1/\lambda )`$ and consider the limit $`\lambda e^{L/\mathrm{}}`$. In this limit the contour is shifted through the pole so that the term $`F_0`$ is cancelled. Moreover, the hypergeometric function $`{}_{1}{}^{}F_{1}^{}`$ can be replaced by unit in this limit. Thus, we end up with the integral

$`F(\lambda )={\displaystyle \frac{1}{2\sqrt{\pi }}}e^{\frac{\mathrm{}}{4L}\left(\mathrm{ln}\frac{1}{\lambda }\frac{L}{\mathrm{}}\right)^2}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑\nu e^{\nu ^2L/4\mathrm{}}`$ (A8)
$`\times \left[\left(i\nu 1{\displaystyle \frac{\mathrm{}}{L}}\mathrm{ln}\lambda \right)\mathrm{\Gamma }\left({\displaystyle \frac{i\nu }{2}}{\displaystyle \frac{\mathrm{}}{2L}}\mathrm{ln}\lambda \right)\right]^1.`$ (A9)
The asymptotic form of this integral in the limit $`\lambda e^{L/\mathrm{}}`$ can be found by evaluation of the expression in square brackets in the point $`\nu =0`$ and calculation of the Gaussian integral. Using the asymptotic formula for the Euler gamma function one obtains the mean density of states given in Eq. (41).