Bound states around impurities in a superconducting bilayer

Due to the presence of the layer degree of freedom, the normal-state Hamiltonian Eq. 2 describes a two-band electronic system. Expressed in momentum space, the band energies are given by

where $a$ is the lattice constant. The multiband nature has an important consequence for the superconductivity: Pairing can occur between electrons in both the same and different bands. Of the pairing potentials listed in Tab. 1, only the uniform singlet state pairs electrons in time-reversed states, which guarantees purely intraband pairing; the other superconducting states will in general involve both intra- and interband pairing.

The concept of the superconducting fitness has recently been developed as a way to quantify the degree of interband pairing Ramires and Sigrist (2016); Ramires et al. (2018). For the pairing potential $\hat{\Delta}=\Delta_{0}\hat{\Delta}_{\nu}$, where $\Delta_{0}$ is the amplitude and $\hat{\Delta}_{\nu}$ is a dimensionless matrix for pairing in channel $\nu$, we define the fitness function

which is only nonzero if the interband pairing is present. For the two-band system considered here the fitness can be related to the gap opened at the Fermi energy

where

Eq. 13 takes values between $0$ and $1$, with $\tilde{F}_{\nu,{\bm{k}}}=0$ indicating a perfectly “fit” pairing state where there is only intraband pairing and thus the gap is maximal, whereas $\tilde{F}_{\nu,{\bm{k}}}=1$ implies that only interband pairing is present and no gap is opened at the Fermi energy at weak coupling.

The concept of superconducting fitness has been extended to the impurity potential Cavanagh and Brydon (2021). In a two-band system, whether or not the impurity is pair-breaking in the channel $\nu$ is given by the solution of

where $\lambda=-1$ ($1$) implies that the impurity is (is not) pair-breaking. In the case of on-site singlet pairing this is equivalent to whether or not the impurity breaks time-reversal symmetry, and is thus a restatement of Anderson’s theorem Anderson (1959); for the nontrivial $s$-wave pairing potentials Eq. 14 therefore establishes a generalized Anderson’s theorem. An important difference compared to the usual Anderson’s theorem is that the critical temperature of the odd-parity $s$-wave states is generally suppressed even when the impurity is not pair-breaking if the fitness $0<\tilde{F}_{\nu,\bm{k}}<1$, albeit more slowly compared to the usual Gor’kov theory Mineev and Samokhin (1999).

We solve the self-consistency equations Eq. 9 for the pairing potential in real-space by diagonalizing the BdG Hamiltonian Eq. 8 implemented on a finite lattice with periodic boundary conditions Zhu (2016). For the calculations of the impurity bound state spectrum we choose a $31\times 31$ lattice with the impurity located at the center.

The four parameters $(t,\mu,t_{\perp},\alpha)$ of our tight-binding model describe a large variety of different band structures. Although our primary aim is to understand the influence of the intertwined spin and layer degrees of freedom on the impurity bound state spectrum, the bound states can also be influenced by the shape of the Fermi surface and the Fermi velocity Ortuzar et al. (2022); Uldemolins et al. (2022). To control for this variation we observe that in the absence of intralayer hopping $t=0$, the iso-energy lines of the bands Eq. 10 in the Brillouin zone are the same for all parameter choices so long as $\alpha\neq 0$. Thus, for fixed filling, the shape of the Fermi surface is the same for all values of $t_{\perp}$ and $\alpha$. This is equivalent to setting the chemical potential to be

for some fixed $k_{\text{F}}$; in the following we choose $k_{\text{F}}a=0.9$ which gives the Fermi surface shown in Fig. 1. Since the intralayer hopping appears with the identity matrix in Eq. 3, it does not influence the gap at the Fermi surface, and we thus do not expect that setting $t=0$ will qualitatively affect our results. This assumption will later be verified using the analytic $T$-matrix theory in Sec. II.3.

The Fermi velocity in our model is given by

where $k_{x}$ and $k_{y}$ lie on the Fermi surface $\sin^{2}(k_{x}a)+\sin^{2}(k_{y}a)=\sin^{2}(k_{\text{F}}a)$. The magnitude of the Fermi velocity decreases with increasing $t_{\perp}$, which should reduce the coherence length even for fixed pairing potential, and may affect the impurity bound state spectrum. We keep the Fermi velocity constant by scaling $\alpha$ as

where $\alpha_{0}$ is the value of $\alpha$ at $t_{\perp}=0$. In the following we set $\alpha_{0}$ as the unit of energy. For fixed $k_{\text{F}}a$, the normal-state properties of the bilayer model now only depend on the ratio $t_{\perp}/\alpha_{0}$.

The effective gap of the six $s$-wave pairing states at the Fermi energy is given in Table. 1. All but the uniform singlet depend on the parameters of the model, reflecting the pair-breaking by the spin-layer texture of the normal-state band wavefunctions. Importantly, we see that the $A_{\text{1g}}$, $A_{\text{1u}}$, and $A_{\text{2u}}$ produce full isotropic gaps at the Fermi energy. In contrast, the two $E_{\text{u}}$ states have point nodes. Since impurity bound states are only well-defined in the presence of a full gap, we will not examine the $E_{\text{u}}$ states within the mean-field theory; the impurity virtual bound states which appear in these cases will be analyzed using the $T$-matrix theory developed below. Moreover, we will also not present results for the nontrivial $A_{\text{1g}}$ state, as this is essentially equivalent to the trivial $A_{\text{1g}}$ state when only a single band crosses the Fermi energy Cavanagh and Brydon (2021). In solving the mean-field equations, the effective coupling constant $g$ in Eq. 7 is chosen such that the pairing amplitude $\Delta_{0}=\frac{g}{2}\langle\mathbf{c_{r_{j}}}\hat{\Delta}_{\nu}\mathbf{c_{r_{j}}}\rangle=0.08\alpha_{0}$ in the absence of the impurity.

To complement the numerical mean-field calculation we also study the appearance of bound states at the impurity site within the $T$-matrix approximation. This has the benefit that it can be performed analytically if we relax the self-consistency requirement Eq. 9, i.e. we treat $\Delta_{0}=\frac{g}{2}\langle\mathbf{c_{r_{j}}}\hat{\Delta}_{\nu}\mathbf{c_{r_{j}}}\rangle$ as independent of $\mathbf{r_{j}}$.

The non-self-consistent $T$ matrix is defined as

where $\hat{V}_{\text{imp}}$ is the impurity potential in Nambu grading and

is the momentum-average of the Green’s function $G_{0}({\bf k},i\omega_{n})$ of the superconductor in the absence of the impurity, with $N$ the number of lattice points.

Analytic expressions for the Green’s functions are complicated and the exact evaluation of the momentum sum in Eq. 18 is difficult. To make progress, we replace the sum with an integral over energy near the Fermi surface and a Fermi-surface average

where ${N}(\xi)$ is the density of states (DOS) and $\Lambda$ is an energy cutoff. Although analytic expressions for the DOS of the tight-binding model exist, here we pursue a more general approach by adopting a linear approximation to the DOS

where ${N}_{0}$ is the DOS at the Fermi energy and ${N}_{0}^{\prime}$ is the first-order derivative with respect to energy. The parameter ${N}_{0}^{\prime}$ reflects the particle-hole asymmetry of the DOS. Although the presence of particle-hole asymmetry is not essential for the formation of bound states, it is necessary for achieving quantitative agreement with the tight-binding model.

Using these approximations, the momentum-averaged Green’s functions for the three pairing states are

where $\tau_{j}$ are the Pauli matrices in Nambu space and its product with $\eta$ and $\sigma$ matrices implies a Kronecker product, $\tilde{t}_{\perp}=t_{\perp}/\mu$, and

The parameter $y\approx{N}^{\prime}_{0}\Lambda$ encodes the particle-hole asymmetry, and we treat it as a fitting parameter. A sketch of the derivation of Eq. 21 is provided in appendix A. In particular, we show that Eq. 21 is valid for any band structure where only a single band crosses the Fermi energy. This gives us confidence that our results hold beyond the restricted parameter space chosen for the lattice model where we set $t=0$.

The relatively simple form of the integrated Green’s function allows us to evaluate the $T$ matrix Eq. 18 analytically. The impurity bound states are then determined by making the analytic continuation $i\omega_{n}\rightarrow\omega+i0^{+}$ and solving for the poles at subgap energies $|\omega|<\Delta_{\text{eff}}$. In the following, we express our results in turns of dimensionless quantities $\tilde{\omega}=\omega/|\Delta_{\text{eff}}|$, $\tilde{V}=\frac{\pi}{2}N_{0}V$ and $\tilde{y}=y/(\frac{\pi}{2}N_{0})$.

The intralayer singlet $A_{\text{1g}}$ state pairs electrons in time-reversed states, and is thus subject to Anderson’s theorem. Since the normal state Hamiltonian has time-reversal symmetry, this state is perfectly fit in the clean limit and there is only intraband pairing. The only pair-breaking impurities are those that break time-reversal symmetry, and we accordingly only find bound states around the impurities in these cases, as shown as a function of impurity strength and interlayer coupling in Fig. 2 and 3, respectively. These figures also show the excellent quantitative agreement between the mean-field and the $T$-matrix theories. Analytic expressions for the bound states obtained via the $T$-matrix method are complicated and are given in the appendix B.1.

Our results reveal two distinct classes of bound states: class (a) impurities where the bound states are two-fold degenerate, and class (b) impurities where the bound states are not degenerate. The degeneracy of the states in case (a) is ensured by a generalized Kramers’ theorem based on the preserved $\mathcal{PT}$ symmetry; in contrast, in case (b) both $\mathcal{T}$ and $\mathcal{PT}$ symmetries are broken, and there is no symmetry enforcing degeneracy. We note that even in case (b) the spectra is two-fold degenerate in the case of vanishing interlayer coupling.

We now turn to the bound state spectrum in the case of the fully-gapped odd-parity states, which is considerably richer than the $A_{\text{1g}}$ state. The classification of the spectra is the same for the two odd-parity states, reflecting the fact that in the clean limit the BdG Hamiltonian is the same up to unitary transformation, specifically we have

where ${\cal H}_{\text{BdG},A_{\text{1u}}(A_{\text{2u}})}$ is the BdG Hamiltonian matrix for the $A_{\text{1u}}$ ($A_{\text{2u}}$) state in the absence of the impurity, and

The impurity potential is generally not invariant under this unitary transform: this allows us to map each impurity potential in the $A_{\text{1u}}$ pairing state to another impurity potential in the $A_{\text{2u}}$ pairing state. It also follows from Eq. 23 that the fitness of the $A_{\text{1u}}$ and $A_{\text{2u}}$ states are the same in the clean limit, with both pairing states unfit with respect to the interlayer coupling $t_{\perp}$. Both states are perfectly fit when $t_{\perp}=0$, and in this limit there also exists a unitary transformation that maps the odd-parity state to the $A_{\text{1g}}$ intralayer singlet state Fu and Berg (2010). We therefore expect that the bound states for the odd-parity $s$-wave pairing will strongly depend on the interlayer coupling.

We find six different classes of bound state spectra, which are shown in Figs. 4 and 5 as a function of $\tilde{V}$ and $\tilde{t}_{\perp}$, respectively. There is again excellent agreement between the self-consistent real-space lattice theory and the continuum $T$-matrix approximation, using the same particle-hole asymmetry parameter $\tilde{y}=0.50$ as for the $A_{\text{1g}}$ results. The dependence on the impurity strength shows a clear difference between the fit impurities (cases (a), (c), and (e)) and the unfit impurities (cases (b), (d), and (f)): although bound states are generally present, only for the unfit impurities do the bound states cross zero energy as a function of impurity strength. The closing of the gap by the impurity evidence the pair-breaking nature of the unfit impurities. The role of the fitness with respect to the normal-state Hamiltonian is clearly revealed by the dependence on the interlayer coupling $t_{\perp}$ shown in Fig. 5: bound states for the fit impurities are always absent for vanishing intralayer coupling, as anticipated by the mapping between the odd-parity and the $A_{\text{1g}}$ states.

Greater insight is provided by the analytic expressions for the bound states found in the $T$-matrix theory. These expressions are generally rather complicated, so to more clearly reveal the relationship we set $\tilde{y}=0$; general expressions with $\tilde{y}\neq 0$ are provided in Appendix B.2. In the particle-hole-symmetric limit, the bound states are

where $\lambda=\pm 1$ is the solution to Eq. 14 (i.e. the fitness of the impurity), and the parameter $S_{t_{\perp}}=+1$ ($-1$) according as the impurity commutes (anticommutes) with the interlayer hopping term. The combinations of $(S_{t_{\perp}},\lambda)$ give rise to four distinct bound state spectra: $(1,1)$ corresponding to cases (a) and (e), $(1,-1)$ corresponding to (b) and (f), $(-1,1)$ corresponding to (c), and $(-1,-1)$ corresponding to (d). The difference between the (a) and (e) (or (b) and (f)) cases is only apparent for nonzero asymmetry parameter $\tilde{y}$. In particular, a nonzero value of the asymmetry parameter is necessary to lift the two-fold degeneracy of the bound states, and is also responsible for the asymmetric dependence on the impurity potential seen in cases (a) and (b), and the interlayer coupling strength in case (b). These effects of the asymmetry parameter are consistent with previous work Wang et al. (2012); Mashkoori et al. (2017).

Despite these limitations, Eq. 25 nevertheless reveals the role of fitness. We first consider the different impurity fitness cases, $\lambda=\pm 1$: in the case of the fit impurity $\lambda=1$, the numerator of Eq. 25 is always positive and so the impurity cannot close the gap; on the other hand, for unfit impurities $\lambda=-1$ the numerator vanishes at a critical value $\tilde{V}=1/\sqrt{1-\tilde{t}_{\perp}^{2}}$, thus closing the gap. We also observe that the bound-state energies depend nontrivially on the interlayer coupling, reflecting the interplay of the impurity with the spin and layer degrees of freedom of the clean Hamiltonian, and thus the normal-state fitness of the pairing potential. The dependence on $S_{t_{\perp}}$ is generally more subtle, but in the case $S_{t_{\perp}}=-1$ it ensures that the term under the square root is a complete square. For $\lambda=1$ this cancels the numerator, and thus there are no subgap states due to the impurity.

The parameters $S_{t_{\perp}}$ and $\lambda$ also enter the Born approximation result for the suppression of the critical temperature by the presence of a finite concentration $n_{\text{imp}}$ of impurities Cavanagh and Brydon (2020, 2021). Here the critical temperature is given by the solution of the equation

where $T_{c0}$ is the critical temperature in the clean limit and the effective scattering rate is

The first two terms inside the brackets give the normal-state scattering rate, while the bulk fitness parameter Eq. 13 controls the effective scattering rate in the superconducting state. For the unfit impurities with $\lambda=-1$ the effective scattering rate is higher than the normal-state scattering rate, and the critical temperature is thus suppressed by increasing impurity concentration faster than the Abrikosov-Gor’kov predictions for an odd-parity state Mineev and Samokhin (1999); Conversely, the fit impurities with $\lambda=1$ have a lower effective scattering rate than in the normal-state, and so the critical temperature displays a weaker suppression than predicted by the usual theory. Intriguingly, in case (c) with $(S_{t_{\perp}},\lambda)=(-1,1)$, the scattering rate vanishes, implying that there is no suppression of the critical temperature. Our results are consistent with the observation that the presence of impurity bound states, even when produced by the fit impurities, indicates some degree of pair-breaking.

A general impurity potential may involve a linear combination of all the matrices listed in Tab. 2, and hence includes both fit and unfit components. Since the formation of bound states is a nonlinear effect we cannot directly infer the bound state spectrum of a general impurity potential from the spectra produced by each component. A physically-interesting case of a multi-component impurity has the impurity potential localized on a single sublattice, i.e. arising from the replacement of an atom at a particular sublattice site by another species. The corresponding impurity potential has the matrix form $(\eta_{0}\pm\eta_{z})\sigma_{\mu}$, which involves two matrices which generally belong to different classes in Tab. 2.

In Fig. 6 we plot the subgap spectrum for sublattice-localized potential ($\mu=0$) and magnetic ($\mu=x$, $y$, $z$) impurities in the two odd-parity states. For both pairing states the subgap spectrum is independent of the polarization of the magnetic impurity. For the $A_{2u}$ state, the two components of the potential (magnetic) impurity are both fit (unfit). Consistent with our argument above, only the the subgap states of the unfit magnetic impurity cross zero energy.

The situation with the $A_{1u}$ state is more subtle, as for both potential and magnetic impurities one component is fit and the other is unfit. In the particle-hole-symmetric limit (i.e. $\tilde{y}=0$), the subgap states asymptotically go to zero at infinite impurity strength; switching on particle-hole asymmetry shifts the zero-crossing to a finite value of $\tilde{V}$. Thus, the presence of the unfit impurity potential is still consistent with the appearance of a zero-crossing of the subgap states.

Close to the $\Gamma$ point, and keeping terms up to linear order in momentum, a minimal three-dimensional generalization of our normal-state Hamiltonian is

where $m$ is the mass term, $v$ is the in-plane velocity, and $v_{z}$ is the out-of-plane velocity. The presence of the out-of-plane velocity $v_{z}$ extends the model to three dimensions; scaling the $z$-component of the momentum allows us to eliminate the velocity anisotropy, i.e. in the following we set $v_{z}=v$. In terms of the parameters of the tight-binding model we have $m=t_{\perp}$ and $v=\alpha a$. We note that a term proportional to $\eta_{z}\sigma_{z}$ is allowed by symmetry, but we neglect it as it appears with much higher power of momenta: for the tetragonal system considered here it requires taking the momentum expansion to the fifth order, whereas for trigonal or hexagonal systems it appears at the third order.

The momentum-averaged Green’s function for both the $A_{\text{1g}}$ and $A_{\text{1u}}$ pairing states in the three-dimensional model is the same as in Eq. 21, and thus the conclusions of Sec. III remain valid. We note that including the $\eta_{z}\sigma_{z}$ term in the Hamiltonian does introduce some gap anisotropy for the $A_{\text{1u}}$ state and may therefore modify our results. So long as this term is small compared to the other terms in Eq. 28, however, we expect that the effects on the bound state spectrum will be small.

The gap opened by the $A_{\text{2u}}$ pairing state on a three-dimensional Fermi surface is required by symmetry to have a node along the $k_{z}$ axis Mineev and Samokhin (1999). In our model, the gap at the Fermi energy is given by

where $\tilde{v}=v/\mu$. Restricted to the $k_{z}=0$ plane we recover the two-dimensional theory, with a uniform gap $\Delta_{0}\sqrt{1-\tilde{m}^{2}}$.

The derivation of the momentum-averaged Green’s function $G_{0}(\omega)$ is quite similar to that in Appendix A, except that accounting for the variation of the gap over the Fermi surface gives a more complicated frequency-dependence:

The virtual bound states are most clearly resolved by examining the local density of states (LDOS) at the impurity site; the deviation of the LDOS from the background contribution from the bulk electronic structure is given by

where $T(\omega)$ is the $T$-matrix, and the factor of $\frac{1}{2}(\tau_{0}+\tau_{z})$ in the trace selects the electron-like components of the Green’s function.

In Fig. 7 we plot the deviation of the LDOS Eq. 31 for impurity potentials corresponding to the same classification as in Fig. 4. We immediately observe that the spectrum is very similar to the fully-gapped two-dimensional model, but with the impurity bound states broadened into virtual bound states due to the hybridization with bulk quasiparticles. This indicates that the impurity physics is still largely controlled by the fitness parameters as found above. Case (c) appears to be an exception to this rule, as we observe subgap features in the three-dimensional system whereas there are no bound states in two dimensions. This discrepancy is expected since the cancellation in Eq. 27 no longer holds: specifically the Fermi-surface average of the fitness is now less than $1-\tilde{t}_{\perp}^{2}$ due to the nodal gap, and so there is a finite effective scattering rate, and thus the impurity is pair-breaking in three dimensions.

A number of other features of the LDOS are notable. The subgap resonances become sharper as one approaches the middle of the gap, reflecting the reduced density of bulk quasiparticle states, which vanish at zero energy. The LDOS also has a pronounced particle-hole asymmetry which arises from the particle-hole asymmetry of the normal-state Kariyado and Ogata (2010); Beaird et al. (2012). We also observe that the LDOS features fade out with increasing impurity potential strength. This is likely due to the shift of wavefunction weight away from the impurity site as the potential increases.

The gaps opened by the $E_{\text{u}}$ pairing states in the three-dimensional model have a similar form compared to the $A_{\text{2u}}$ state, albeit with point nodes along the $x$- and $y$-axes. Indeed, the momentum-space BdG Hamiltonian for the clean $A_{\text{2u}}$ state can be mapped to that for each of the $E_{\text{u}}$ states; specifically, for the two $E_{\text{u}}$ states we have

where

Note that the swapping of the momentum components in Eqs. 32 and 33 does not change the momentum-averaged Green’s function, and thus does affect our results in the $T$-matrix theory. Thus, as for the $A_{\text{1u}}$ and $A_{\text{2u}}$ states in the two-dimensional model, our results for the $A_{\text{2u}}$ state in the three-dimensional model can be mapped to each $E_{\text{u}}$ state individually, with due alteration of the impurity potential by the unitary transform. The impurity potentials corresponding to each of the subgap spectra in Fig. 7 are listed in Tab. 2.