Dark Andreev States in Superconductors

Bardeen-Cooper-Schrieffer theory of superconductivity (Bardeen et al., 1957) and Bogoliubov-de Gennes equations have proven to be relatively simple and reliable tools to describe a variety of conventional superconductors. A key simplifying assumption of this approach is that the superconducting order parameter is energy-independent (Altland and Simons, 2010). While it is clearly not the case in any real superconductor, many features such as the quasiparticle spectrum, thermodynamic and electromagnetic properties (Altland and Simons, 2010; Bardeen et al., 1957) appear insensitive to this approximation. It is reasonable, because the omitted energy dependence normally does not affect low-energy physics. In contrast, accurate determination of the transition temperature is sensitive to the details of the phonon dispersion, the dynamical screening of Coulomb interaction, and the structure of the superconducting gap. The latter can be obtained from the Migdal-Eliashberg (ME) equations (Eliashberg, 1960; Marsiglio, 2020) which are integral equations in both energy and momentum.

This Migdal-Eliashberg theory Marsiglio et al. (1988) predicts a non-trivial structure of the superconducting gap as a function of real frequency. In particular, it was shown that for superconductivity mediated by the Einstein phonon modes (Marsiglio, 2018; Mirabi et al., 2020), the gap function has sharp resonance-like features. In the weak coupling regime, these features can be well approximated by the Lorentz function centered at the Debye frequency. Another interesting real-frequency behavior of the gap function is discussed in (Christensen and Chubukov, 2021) where the authors predict the possibility of formation of frequency-domain vortices in the presence of phonon-induced attraction and Coulomb repulsion.

In this work, we study how a frequency-dependent order parameter affects the quasiparticle spectrum and tunneling properties of a superconductor. We assume that the frequency dependence has a Lorentz shape, which corresponds to optical-phonon-mediated pairing (Marsiglio, 2020). We solve the corresponding generalized BdG equations and show that an additional gap emerges in the qusiparticle spectrum at higher energies. In the case of an SNS junction, we also find that additional Andreev in-gap high-energy states can form. We use analogy with quantum optics to interpret these high-energy peaks in terms of “dark” resonance features (Fleischhauer and Lukin, 2000; Lukin et al., 1999; Fleischhauer and Lukin, 2002) of the BdG Hamiltonian with a frequency-dependent order parameter. This suggests a speculation that these finite-energy states could potentially be used to reliably store quantum information.

In this work we consider the conventional s-wave superconductor but keeping a complete realistic frequency dependence of the order parameter. We restrict our discussion to the mean-field level and rely on the Bogoliubov-de Gennes approach, which requires a straightforward generalization. Consider an electron gas with the creation (annihilation) operators $\psi^{\dagger}_{{\bf k,\sigma=\shortuparrow,\downarrow}}$ $(\psi_{{\bf k},\sigma})$, where ${\bf k}$ denotes the electron momentum. For the two-component spinor $\Psi_{{\bf k},n}=\{\psi_{{\bf k},\downarrow}\left(i\varepsilon_{n}\right),\psi_{-{\bf k},\uparrow}^{\dagger}\left(-i\varepsilon_{n}\right)\}^{T}$, the Green function is defined as: $\hat{{\cal G}}_{{\bf k}}(i\varepsilon_{n})\equiv-\langle\Psi_{{\bf k},n}\otimes\Psi_{{\bf k},n}^{\dagger}\rangle$, where $\varepsilon_{n}=\left(2n+1\right)\pi/\beta$ is the standard the fermionic Matsubara frequency, $\beta$ is the inverse temperature, and $n\in\mathbb{Z}$. Neglecting the momentum dependence of the self-energy, the Green function of an interacting Fermi gas can be written as (Marsiglio, 2020; Schrieffer, 2018):

$$\hat{{\cal G}}_{{\bf k}}^{-1}\left(i\varepsilon_{n}\right)=i\varepsilon_{n}\mathcal{Z}_{n}\hat{\tau}_{0}-\phi_{n}\hat{\tau}_{1}-\xi_{k}\hat{\tau}_{3},$$

(1)

where $\xi_{k}=k^{2}/2m-\mu$, $m$ is the electron mass, $\mu$ is the chemical potential, $\tau_{i}$ are Pauli matrices. $\mathcal{Z}_{n}$ denotes the inverse quasiparticle residue obtained from the odd part of the normal-state self-energy. The off-diagonal matrix element $\phi_{n}$ stands for the anomalous self energy. $\mathcal{Z}_{n}$ and $\phi_{n}$ naturally appear in both intrinsic and proximity-induced superconductors (Chubukov et al., 2020; Liu et al., 2019). The BdG equations correspond to $\hat{{\cal G}}_{{\bf k}}^{-1}\left(i\varepsilon\right)\chi=0$, which parameterically determines the sought-after quasiparticle dispersion (here, $\chi$ is a Nambu spinor).

To calculate the quasiparticle spectrum, we need an explicit form of the frequency dependence of the anomalous self-energy. In what follows, we consider the aforementioned Lorentz-shaped form (parameterized by its amplitude $\phi_{0}$ and the characteristic frequency $\omega_{0}$)

$$\phi_{n}=\phi_{0}\,\frac{\omega_{0}^{2}}{\varepsilon_{n}^{2}+\omega_{0}^{2}}.$$

(2)

As demonstrated in Refs. (Marsiglio, 2018; Mirabi et al., 2020) this solution naturally appears in superconductors, where pairing is induced by an optical phonon mode with the Einstein spectrum at weak coupling. We now consider the spectrum of quasiparticles in Eq. (1).

To develop intuition, it is instructive to consider an auxiliary setup, which involves a fictitious flat-band superconductor with a frequency-independent gap proximity coupled to a normal metal as shown in Fig. 1. As we show, a proper choice of parameters in this setup gives rise to a Green function, which replicates Eq. 1 with the order parameter (2). The advantage of this construction is that its BdG Hamiltonian below involves only standard, frequency-independent parameters:

	$\displaystyle\hat{{\cal G}}_{k}^{\prime-1}\left(i\varepsilon_{n}\right)$	$\displaystyle=i\varepsilon_{n}\hat{\tau}_{0}\hat{\sigma}_{0}-\hat{H}_{\text{BdG}}$		(3)
	$\displaystyle H_{\text{BdG}}=$	$\displaystyle\xi_{k}\hat{\tau}_{3}\frac{\left(\hat{\sigma}_{3}+\hat{1}\right)}{2}-\omega_{0}\hat{\tau}_{1}\frac{\left(\hat{1}-\hat{\sigma}_{3}\right)}{2}+t\hat{\tau}_{3}\hat{\sigma}_{1},$		(4)

Here $\hat{\tau}_{i}$ represents Pauli matrices in the Nambu space and $\hat{\sigma}$ parametrizes an effective two-band model: $\left(\hat{\sigma}_{3}+\hat{1}\right)/2$ projects on the normal metal fermion modes and $\left(\hat{1}-\hat{\sigma}_{3}\right)/2$ projector on the fictitious flat-band superconductor. The order parameter in the latter is set the characteristic phonon frequency, $\omega_{0}$. We note that we assumed no intrinsic order parameter for the original fermions. The coefficient $t$ denotes the tunneling amplitude between the superconductor and the metal. The auxiliary superconductor can be integrated-out generating both the anomalous and normal self energies. By choosing $t=\sqrt{\phi_{0}\omega_{0}}$ we can exactly match frequency dependence of the order parameter to reproduce Eq. 2. The corresponding $\mathcal{Z}$-factor in Eq. (1) is equal to $\mathcal{Z}_{n}=\left[1+\phi_{0}\omega_{0}/(\omega_{0}^{2}+\varepsilon_{n}^{2})\right]\approx 1$, for $\phi_{0}\ll\omega_{0}$. This procedure effectively replaces the integrating out the bosonic Einstein-phonon degree of freedom (which generates the frequency dependence of the gap in the physical setup) with the integrating out the degrees of freedom of the auxiliary flat-band superconductor. Note that this picture is introduced for the purposes of illustration only. All results can reproduced in the original model directly. As we discuss in the supplementary material, the inverse quasiparicle residue, $\mathcal{Z}_{n}$ cannot be set identically to one because it would result in an unstable spectrum.

The quasiparticle spectrum can now be found by diagonalizing the static BdG Hamiltonian $H_{\text{BdG}}$ given in Eq. (4). The result is shown on Fig.  1 (b) and it has two band gaps. First, we observe the conventional BCS-like band gap at low energies $[-\phi_{0},\phi_{0}]$. It can be obtained by e.g. neglecting the frequency-dependence of the self-energy in Eq. (1), which reduces it to the textbook case. The second band gap at high energies is a specific feature of the two-band system as defined in Eq. (3). As a result of the flatness of one of the dispersion relations, the avoided crossing forms a band gap close to the frequency, $\omega_{0}$. The value of this second band gap can be readily obtained analytically from Eq. (3):

$$\phi_{2}=\sqrt{\frac{\omega_{0}}{2}\left(\sqrt{\omega_{0}(4\phi_{0}+\omega_{0})}+\phi_{0}+\frac{\omega_{0}}{2}\right)}-\omega_{0}\approx\phi_{0},$$

(5)

where the approximate sign corresponds to the limit $\phi_{0}\ll\omega_{0}$, which must hold for weak coupling.

We now define the local density of states (LDOS) of the electron gas as $\frac{-1}{\pi}\int d\xi_{k}\text{Im}\left(G_{{\bf k}}\left(\omega+i0^{+}\right)\right)_{1,1}$, where $G_{{\bf k}}$ is the Green function 3 analytically continued to real frequencies. As the direct consequence of the band gap the LDOS is strongly depleted at energies $\approx[\omega_{0},\omega_{0}+\phi_{0}]$ as shown on inset in Fig. 2. We note that the exactly zero density of states is a feature of the flat-band dispersion of the auxiliary superconductor/Einstein phonons. However, as we discuss in the SM, introduction of a finite curvature to the phonon dispersion would still lead to a significant depletion of the density of states. As we discuss below this secondary gap can host additional Andreev (Sauls, 2018) reflection peaks, observable in metal-superconductor heterostructures.

We now explore how the additional sharp Lorentz-like features of the gap function affect the superconducting proximity effect. In order to describe the transmission and reflection of quasiparticles, we employ the Blonder-Tinkham-Klapwijk (BTK) formalism (Blonder et al., 1982). We consider a heterostructure consisting of normal and superconducting metals (NS). Following (Blonder et al., 1982), we consider the scattering of an incident electron off of the barrier. The strength of the proximity effect can be characterized by the probability for an electron to scatter into a hole-type excitation.

Within the BTK theory the boundary condition is given by:

	$\displaystyle\vec{\psi}_{N}\left(0\right)$	$\displaystyle=\vec{\psi}_{S}\left(0\right),$		(6)
	$\displaystyle\frac{\partial_{z}}{2m_{N}}\vec{\psi}_{N}\left(0\right)$	$\displaystyle=\frac{\partial_{z}}{2m_{S}}\vec{\psi}_{S}\left(0\right)+H\vec{\psi}_{S}(0).$		(7)

where $m_{S}$ and $m_{N}$ are the effective electron masses and $H$ is the $\delta$-barrier height. Performing the analytic continuation and replacing $i\epsilon_{n}\rightarrow\omega+i0^{+}$, the wavefunctions on superconducting side satisfy the equation $\hat{G}_{{\bf k}}^{-1}(\omega+i0^{+})\vec{\psi}_{S}=0$. On the normal side the equation is the same with the substitution $\phi(\omega+i0^{+})\rightarrow 0$ and $\mathcal{Z}(\omega+i0^{+})\rightarrow 1$. In the following we do not explicitly write $0^{+}$ for shortness. The normal-state solution representing an incident electron and reflected electron and hole components is:

$$\vec{\psi}_{N}=\begin{bmatrix}1\\ 0\end{bmatrix}e^{ik_{\text{e}}z}+\begin{bmatrix}r_{N}\\ 0\end{bmatrix}e^{-ik_{\text{e}}z}+\begin{bmatrix}0\\ r_{A}\end{bmatrix}e^{ik_{\text{h}}z},$$

(8)

where $r_{N}$ and $r_{A}$ denote the reflection amplitude in the electron and hole channels respectively and the electron/hole momenta are given by $k_{\text{e}/\text{h}}=\sqrt{2m\left(\mu\pm\omega\right)}$. Analogously we find the solution for the quasi-electrons and quasi-holes propagating in the superconductor:

$\displaystyle\vec{\psi}_{S}$

$\displaystyle\approx C_{\text{qe}}\begin{bmatrix}1\\ \eta_{+}\end{bmatrix}e^{ik_{\text{qe}}z}+C_{\text{qh}}\begin{bmatrix}1\\ \eta_{-}\end{bmatrix}e^{-ik_{\text{qh}}z},$

(9)

where $C_{\text{qe}}$ and $C_{\text{qh}}$ are the corresponding amplitudes of the quasi-electron and quasi-holes and we denoted the corresponding coherence factors as $\eta_{\pm}=\phi(\omega)/\left(\mathcal{Z}(\omega)\omega\pm\sqrt{\mathcal{Z}^{2}(\omega)\omega^{2}-\phi^{2}(\omega)}\right)$. In the quasiclassical limit the quasielectron and quasihole momenta $k_{\text{qe}},k_{\text{qh}}$ can be taken to be equal to the corresponding Fermi momenta.

Upon solving the set of equations Eqs. (6-9) in the quasi-classical Belzig et al. (1999) limit and assuming $m_{S}=m_{N}$ we find the Andreev reflection coefficient to be:

	$\displaystyle r_{\text{A}}$	$\displaystyle=\frac{\eta_{-}\eta_{+}}{\eta_{-}+Z^{2}\left(\eta_{-}-\eta_{+}\right)},$		(10)
	$\displaystyle r_{N}$	$\displaystyle=\frac{-Z\left(i+Z\right)}{\eta_{-}+Z^{2}\left(\eta_{-}-\eta_{+}\right)}\left(\eta_{-}-\eta_{+}\right)$		(11)

where the normalized barrier height is defined as $Z\equiv mH/\sqrt{2m\mu}$. We note that the conventional Andreev reflection can be obtained from Eq. (10) by simply assuming frequency-independent $\mathcal{Z}$ and $\phi$. In the limit of $\mathcal{Z}=0$ the Andreev reflection coefficient is given by $r_{A}=\eta_{+}$. We therefore find that in order to have a strong Andreev reflection the condition $\phi\gg\mathcal{Z}\omega$ should be satisfied. The latter condition is always satisfied at very low frequencies leading to the conventional Andreev reflection (Sauls, 2018) result $\left|r_{\text{A}}\right|\approx 1$. However it can also be satisfied at large frequencies if the frequency-dependent order parameter $\Delta\left(\omega\right)\equiv\phi\left(\omega\right)/\mathcal{Z}\left(\omega\right)$ is larger than the frequency $\Delta\left(\omega\right)\gg\omega$. In this case the reflection is up to a possible phase factor identical to the low-energy case. As shown in Fig. 2, this scenario is realized for the Lorentz-like order parameter introduced above in the frequency range close to the resonance $\omega\sim\omega_{0}$.