Interplay of Fermi velocities and healing lengths in two-band superconductors

Characteristic length scales associated with different contributing condensates constitute one of the cornerstone features of multiband superconductors. Multiple condensates in one system interfere, which results in unconventional coherent phenomena Milorad2015 . Effects of such interference are most pronounced when the spatial lengths of the contributing condensates are notably different. Various definitions of such lengths are in use, including those related to the gap function slope in the vortex core Caroli1964 ; Gygi1991 , the maximum density of the supercurrent Sonier2004 , the radius of a cylinder containing energy equal to the condensation energy Gennes1966 ; Tinkham1996 , or the healing length along which the condensate reaches $60$-$80\%$ (there are different choices) of its bulk value Komendova2011 ; Komendova2012 ; Saraiva2017 . All such definitions produce similar results (except of the slope definition that fails at nearly zero temperatures due to the Kramer-Pesch collapse Kramer1974 ; Chen2014 ) and either of them can be employed to characterize the condensate spatial scales.

Since 1970s it is well known Geilikman1967 that the spatial lengths of different band condensates in multiband materials are the same in the Ginzburg-Landau (GL) domain, see also Refs. Koshelev2004, ; Koshelev2005, ; Geyer2010, ; Kogan2011, ; Shanenko2011, ; Vagov2012, . However, using the perturbative expansion of the microscopic equations in the small deviation from the critical temperature $T_{c}$ to one order beyond the GL theory (extended GL), one finds that the band-dependent condensate lengths can be different Shanenko2011 ; Vagov2012 ; Saraiva2017 . This conclusion was confirmed by numerically solving the two-band Bogoliubov-de Gennes (BdG) equations Komendova2012 . Moreover, it was also demonstrated Komendova2012 that the condensate characteristic length associated with a weaker band notably increases when approaching the critical temperature of this band taken as a separate superconductor (the hidden critical point). As the length of the stronger band condensate remains unaffected in this case, one can get an increased difference between the two lengths governed by the hidden criticality. On the other hand, the condensate lengths for sufficiently strong interband couplings tend to be nearly the same, as demonstrated in Refs. Saraiva2017, and Ichioka2017, . This can possibly explain the recent scanning tunneling microscopy measurements Fente2016 and can be referred to as the “length-scales locking” Ichioka2017 .

Though the disparity between the condensate lengths is more pronounced for weaker interband couplings, the “length-scales locking” regime can be shifted toward larger values of the interband couplings. Indeed, it was recently shown within the extended GL formalism Saraiva2017 that the difference between the condensate healing lengths in a two-band superconductor is very sensitive to the ratio of the band Fermi velocities $v_{Fi}\;(i=1,2)$ which varies within a wide range, see, for example, Table 1 illustrating some experimental results.

Furthermore, this ratio can be altered by doping in superconductors Lee2006 , changing the topology of the Fermi surface Cappelluti2007 , engineering the interface of the system Peng2014 , applying the pressure Tissen2013 ; Suderow2005 , and changing the characteristic size of nanoscale superconductors via the quantum-size effects Blatt1963 ; Shanenko2006a ; Guan2010 ; Chen2010 ; Shanenko2008 ; Chen2012a . As the impact of the Fermi velocities on the condensate lengths was investigated by means of the extended GL formalism and in the absence of the magnetic effects, it is of importance to compliment the conclusions of Ref. Saraiva2017, by investigating temperatures far below $T_{c}$ and including the magnetic field.

In this work we explore the interplay between the band Fermi velocities and the condensate healing lengths by numerically solving the two-band BdG equations for a single vortex solution in the entire range of the temperatures below $T_{c}$. As the local magnetic field is not neglected, the BdG equations are supplemented by Ampere’s law introducing an additional magnetic coupling between the contributing condensates. The special attention is also given to the effect of the hidden criticality at which the disparity between the healing lengths is most pronounced.

The paper is organized as follows. In Sec. II, we outline the formalism of the BdG equations for a single-vortex state in a two-band condensate. The numerical results and related discussions are given in Sec. III including three subsections. The first subsection presents results for the zero temperature $T=0$ and zero external field $H=0$. Here one can find the healing lengths $\xi_{1}$ and $\xi_{2}$ as functions of the Fermi velocities ratio $v_{F2}/v_{F1}$ and the interband coupling $g_{12}$. For illustration, we also show how the healing lengths are extracted from the spatially dependent gap functions. The results for $T\neq 0$ and $H=0$ are discussed in the second subsection. Here we investigate the healing lengths as functions of $T$ for different parameters $v_{F2}/v_{F1}$ and $g_{12}$. In the third subsection we discuss the results for $T\neq 0$ and $H\neq 0$. Conclusions are given in Sec. IV.

To investigate how the spatial scales of the partial band condensates in a two-band superconductor are sensitive to the band Fermi velocities, a single vortex solution of the two-band BdG equations is considered in a cylinder with the vortex line parallel to the $z$ axis of this cylinder. We utilize the standard microscopic model of a two-band superconductor Suhl1959 ; Moskal1959 with the conventional $s$-wave pairing in both bands, controlled by the symmetric coupling matrix $g_{ii^{\prime}}$ ($i,i^{\prime}=1,2$). The intraband couplings $g_{11}$ and $g_{22}$ are chosen so that the critical temperature of band $1$, taken as a separate superconductor, is larger than the critical temperature in the decoupled band $2$, i.e., we have stronger band $1$ and weaker band $2$. The two condensates are coupled through the Josephson-like transfer of Cooper pairs controlled by $g_{12}$. The parabolic single-particle energy dispersion is assumed for charge carriers in both bands. For our calculation we choose quasi-2D bands, as multiband materials often exhibit quasi-2D Fermi surfaces, see e.g. Ref. Paglione2010, . An external magnetic field is applied along the $z$ axis of the cylinder while the dependence of the quasi-2D band dispersions on the $z$ projection of the single-particle momentum is minor and neglected in our calculations. The superconductor is in the clean limit.

The corresponding BdG equations read Komendova2012 ; Araujo2009

where $u_{i\nu}(\bf r)$ and $v_{i\nu}({\bf r})$ are the electron-like and hole-like wave functions associated with band $i$ ($\nu$ is the set of the relevant quantum numbers); $E_{i\nu}$ and $\Delta_{i}({\bf r})$ are the corresponding quasiparticle energy and the spatial pair potential (gap function); and the single-particle Hamiltonian for the charge carriers in band $i$ is given by

with $m_{i}$ the electron band mass, $\mu_{i}=m_{i}v_{Fi}^{2}/2$ the chemical potential measured from the lower edge of the corresponding band, ${\bf D}=\bm{\nabla}-\mathbbm{i}\frac{e}{\hbar\mathbbm{c}}{\bf A}$, and ${\bf A}({\bf r})$ the vector potential.

As the problem is solved in a self-consistent manner, the band gap functions and the vector potential depend on the solutions of Eqs. (1) as

and

where $f(E_{i^{\prime}\nu})$ is the Fermi-Dirac distribution and the supercurrent density is given by

The summation in Eqs. (3) and (5) goes over the quasiparticle states with positive energies. In addition, Eq. (3) includes only the states for which the averaged single-electron energy taken at zero field Shanenko2008 $\langle\hat{H}_{e,i}\rangle|_{{\bf A}=0}$ falls into the range $[-\hbar\omega_{D},~{}\hbar\omega_{D}]$, with $\omega_{D}$ the Debye frequency assumed the same for both contributing bands. Similar results (with deviations of about $1$-$2\%$) can be obtained when selecting $E_{i\nu}<\hbar\omega_{D}$ in Eq. (3), see e.g. Ref. Gygi1991, .

The Josephson-like coupling between the two contributing bands is not explicitly present in the Bogoliubov-de Gennes Eqs. (1), appearing in the self-consistency gap equation Eq. (3). The magnetic coupling between the condensates manifests itself through the presence in Eqs. (1) of the vector potential that is related to the both contributing condensates by means of Ampere’s law Eq. (4). We remark that to go beyond the adopted model, the pairing of electrons from different bands should be taken into consideration, i.e. in addition to the transfer of the Cooper pairs from one band to another, one accounts for an extra coupling through the interband Cooper pairs, including one electron from band $1$ and another from band $2$, see e.g. Ref. Shanenko2015, . In this case the coupling between bands appears in the Bogoliubov-de Gennes equations Shanenko2015 ; Vargas2020 . However, in most cases the interband pairing is suppressed due to incommensurability of the Fermi momenta in different bands and can be neglected.

Considering a single vortex oriented along the $z$ direction, we follow the previous studies of a single vortex solution within the single-band Gygi1991 ; Bardeen1969 ; Hayashi1998 and two-band BdG equations Komendova2012 ; Araujo2009 . Due to the cylindrical geometry, we can write

and

where $\rho,\theta$ and $z$ are the cylindrical coordinates, $L$ is the unit cell of the periodic boundary conditions in the $z$-direction, and $\nu=\{j,m,k_{z}\}$ with $j$ the radial quantum number, $m$ the azimuthal quantum number being half an odd integer, and $k_{z}$ the wavenumber in the $z$-direction. As mentioned above, the dependences of the quasi-2D band dispersions on $k_{z}$ are neglected and so, the wave functions $u_{i\nu}({\bf r})$ and $v_{i\nu}({\bf r})$ are not dependent on $k_{z}$ either.

For the chosen gauge and symmetry, ${\bf A}({\bf r})=A_{\theta}(\rho){\bf e}_{\theta}$, with ${\bf e}_{\theta}$ the azimuthal unit vector. The two boundary conditions for $A_{\theta}(\rho)$ are set as: (1) the magnetic field approaches the external one $H{\bf e}_{z}$ far away from the cylinder; (2) the magnetic field is finite at the origin of the coordinates (the vortex center). The latter assumes $A_{\theta}(0)=0$ to avoid the divergence of the field. In addition, the transverse quantum confinement requires the boundary conditions $u_{ijm}(\rho=R)=0$ and $v_{ijm}(\rho=R)=0$, where $R$ is the radius of the cylinder.

To represent the BdG equations in the matrix form, we expand the radial parts of the particle-like and hole-like wave functions $u_{ijm}(\rho)$ and $v_{ijm}(\rho)$ in terms of the normalized Bessel functions of the first kind

where superscripts “-” and “+” are for $u$ and $v$ functions, respectively, and $\alpha_{i,\eta}$ is the $i$th zero of the corresponding Bessel function, i.e., $J_{\eta}(\alpha_{i,\eta})=0$. The expansion writes

where $N$ should be chosen sufficiently large to capture the essential features of the vortex solution. As a result, the BdG equations are reduced to the matrix ($2N\times 2N$) equation with the eigenvectors given by $\{c_{ijj^{\prime}m}\}$ (upper half of the column) and $\{d_{ijj^{\prime}m}\}$ (lower half). Then, the numerical solution of the problem is calculated in the self consistent manner. First, we choose some initial gap functions $\Delta_{i}(\rho)$ and vector potential $A_{\theta}(\rho)$ and find the corresponding eigenvalues and eigenstates of the matrix BdG equations. Second, we use the obtained sets $\{c_{ijj^{\prime}m}\}$ and $\{d_{ijj^{\prime}m}\}$ and the related quasiparticle energies $E_{ijm}$ to calculate the new position dependent gaps and vector potential by means of the equations Eqs. (3)-(5), (7), and (9). Third, the BdG equations are solved again with the calculated gap functions and vector potential. The procedure is repeated until the convergence.

Below the effective band-dependent electron masses $m_{i}$ are set to the free electron mass $m_{e}$, for simplicity. The intraband couplings are chosen such that $g_{11}N_{1}=0.3$ and $g_{22}N_{2}=0.24$, where $N_{i}$ is the normal density of states (DOS) per spin projection of band $i$. In the case of interest $N_{1}=N_{2}=(m_{e}/2\pi\hbar^{2}L)\sum_{k_{z}}\theta(k_{\rm max}-|k_{z}|)$, with $\theta(k_{\rm max}-|k_{z}|)$ the step function and $k_{\rm max}$ the maximal wavenumber in the $z$ direction. One can estimate $k_{\rm max}=\pi/a_{z}$, where $a_{z}$ is the corresponding lattice constant. Then, using the periodic boundary conditions for the motion in the $z$ direction, one gets $(1/L)\sum_{k_{z}}\theta(k_{\rm max}-k_{z})\sim 1/a_{z}$. Keeping in mind typical values for the lattice constant, one concludes that $1/a_{z}$ is of the order of $1$-$3{\rm nm}^{-1}$. For our calculations we choose $N_{1}=N_{2}=\tilde{N}m_{e}/2\pi\hbar^{2}$, with $\tilde{N}=1{\rm nm}^{-1}$ (similarly to Ref. Gygi1991, ). Notice that this choice and also the use of $m_{i}=m_{e}$ do not influence our conclusions because any changes in $N_{i}$ result simply in the renormalization of the intraband couplings $g_{11}$ and $g_{22}$, as we keep the same dimensionless couplings $g_{11}N_{1}$ and $g_{22}N_{2}$. Notice that the chosen values for the intraband couplings are in the typical range for multiband materials, see Ref. Vagov2016, and references therein. The interband coupling $g_{12}$ is varied in our study, in order to investigate effects of the interaction between the two contributing condensates.

To have different Fermi velocities $v_{F1}$ and $v_{F2}$, we choose different $\mu_{1}$ and $\mu_{2}$. For the stronger band we adopt $\mu_{1}=30$ meV, based on conservative estimates of the Fermi energy in emergent multiband superconductors, see e.g. Ref. Lubashevsky2012, . The chemical potential relative to the lower edge of the weaker band $\mu_{2}$ is varied in our calculations so that the ratio of the band Fermi velocities $v_{F2}/v_{F1}$ is altered by this variation.

To avoid effects of quantum confinement, the radius $R$ of the cylinder should be chosen sufficiently large. When taking the Debye frequency as $\hbar\omega_{D}=15$ meV (in the range of conventional values, see e.g. Ref. Fetter2003, ), one finds that for the zero temperature the healing lengths $\xi_{1}$ and $\xi_{2}$ are not sensitive to the cylinder radius for $R\gtrsim 100\,{\rm nm}$. For example, the calculations yield $\xi_{1}=19.2$ nm and $\xi_{2}=29.3$ nm when employing $g_{12}=0.05g_{11}$ and $v_{F2}/v_{F1}=1$ for $T,H=0$. Then, choosing $R=300\,{\rm nm}$, we safely have $R>\xi_{1,2}$ up to the temperatures $T\approx 0.99T_{c}$. Since the healing lengths increase with the temperature approximately as Fetter2003 $\propto\tau^{-1/2}$, with $\tau=1-T/T_{c}$, they approach $R$ at $T\approx 0.99T_{c}$. Only in this case $\xi_{1}$ and $\xi_{2}$ are affected by the geometry of the sample.

We also note that the presence of the boundary conditions $u_{ijm}(\rho=R)=0$ and $v_{ijm}(\rho=R)=0$ introduces an additional condensate length near the boundary. Indeed, here $\Delta_{i}(\rho)$ exhibits a series of the Friedel-like oscillations Troy1995 ; Martin1997 with the period of a half of the band-dependent Fermi wavelength $\lambda_{Fi}/2$. For the chosen parameters we have $\lambda_{F1}/2=1.1\,{\rm nm}$ and $\lambda_{F2}\sim\lambda_{F1}$. One sees that $\lambda_{F1,2}/2$ is much smaller than $\xi_{1,2}$ and the presence of the Friedel-like oscillations can in no way distort our results.

In this section we discuss the results of numerically solving the BdG equations for a single vortex in the two-band superconducting condensate within the model outlined in the previous section.

Before the discussion, we need to stress that the case of the zero external field $H=0$ does not assume the absence of any magnetic effects. The local field $B$ is always nonzero in the vortex core and the magnetic coupling between the two band condensates is present even for $H=0$ due to Ampere’s law given by Eq. (4). [Obviously, its impact on the healing lengths can be neglected only in deep type II.] Then, the question arises which boundary conditions for the magnetic field far beyond the vortex core we should use, to obtain relevant information about the condensate healing lengths in the mixed state. We recall that for a single-vortex solution in bulk we have $B\to 0$ at infinity, see e.g. Ref. Vagov2016, . Furthermore, near the lower critical field $H_{c1}$, an Abrikosov lattice exhibits a significant distance separating neighbouring vortices so that the single-vortex state is a good approximation for such a dilute lattice. In this case the local field $B$ is indeed exponentially small between vortices, being far smaller than the external field. Clearly, to describe this case, the boundary condition $B\to H=0$ should be applied far beyond the vortex core in our calculations.

Near the upper critical field $H_{c2}$ the local field $B$ approaches the external magnetic field between vortices in the vortex matter. We can model this situation by invoking the boundary condition $B\to H\not=0$ far beyond the vortex core. However, it is necessary to keep in mind that the healing lengths for the single-vortex state can deviate from the corresponding lengths in a dense Abrikosov lattice appearing close to $H_{c2}$. Below we investigate both $H=0$ and $H\not=0$. We expect that the former case gives the healing lengths in the two-band superconductor near the lower critical field while the latter case is more suitable to consider the mixed state near the upper critical field.

We have studied the effect of the band Fermi velocities on the healing lengths in a two-band superconductor by numerically solving the Bogoliubov-de Gennes equation for a single-vortex solution. Our results demonstrate that near the lower critical field the healing lengths of the two contributing condensates can be significantly different for sufficiently large values of the ratio of the band Fermi velocities $v_{F2}/v_{F1}$. This occurs far beyond the regime of nearly decoupled bands, at the interband couplings up to $g_{12}\sim g_{11},g_{22}$. The most pronounced difference between the healing lengths is observed in the vicinity of or below the hidden critical temperature. The “length-scales locking” regime takes place near the upper critical field and/or near the critical temperature.

Our study is connected with the long discussion about the possibility to have two coupled condensates with significantly different spatial profiles in the presence of the magnetic effects. Our work clearly demonstrates this possibility for a wide range of the physical parameters. The presence of different healing lengths can significantly change the magnetic response of multiband superconductors as compared to that of single-band ones. For example, it is known that the switching between superconductivity types I and II occurs through the finite intertype domain in the $\kappa$-$T$ plane ($\kappa$ is the Ginzburg-Landau parameter), see e.g. Refs. Brandt2011, and Vagov2020, . It has been proved Vagov2016 ; Cavalcanti2020 that this domain significantly enlarges in two-band and multiband superconductors (with respect to the single-band case) if the healing lengths of different contributing condensates are significantly different. Our present study is a solid compliment to these previous investigations based on the perturbation theory in the vicinity of $T_{c}$. We confirm that multiband materials with significantly different band Fermi velocities are most promising in searching for unconventional superconducting magnetic properties because of the presence of multiple condensates governed by different spatial scales.