Higgs modes in proximized superconducting systems

TL;DR

This paper investigates how the proximity effect in superconductor-normal metal hybrids influences Higgs modes, revealing damping of the standard mode and the emergence of two new modes affecting the system's dynamical behavior.

Contribution

It demonstrates the existence of two new Higgs modes in hybrid structures and analyzes their impact on the long-time dynamics of the superconducting order parameter.

Findings

Standard Higgs mode at 2Δ is damped by quasiparticle leakage.

Two new Higgs modes depend on both primary and induced gaps.

These modes influence the dynamical properties of superconducting devices.

Abstract

The proximity effect in hybrid superconducting - normal metal structures is shown to affect strongly the coherent oscillations of the superconducting order parameter $\Delta$ known as the Higgs modes. The standard Higgs mode at the frequency $2\Delta$ is damped exponentially by the quasiparticle leakage from the primary superconductor. Two new Higgs modes with the frequencies depending on both the primary and induced gaps in the hybrid structure are shown to appear due to the coherent electron transfer between the superconductor and the normal metal. Altogether these three modes determine the long-time asymptotic behavior of the superconducting order parameter disturbed either by the electromagnetic pulse or the quench of the system parameters and, thus, are of crucial importance for the dynamical properties and restrictions on the operating frequencies for superconducting devices based on the proximity effect used, e.g., in quantum computing, in particular, with topological low-energy excitations.

Figures

Equations

$$\begin{gathered} \delta \Delta \sim \frac{\cos (2{\Delta }_{0}t)}{\sqrt{{\Delta }_{0}t}}{e}^{-2{\Gamma }_{s}t}+ \\ D\frac{\cos [({\Delta }_0+{\Delta }_i)t]}{({\Delta }_{0}t{)}^p}{e}^{-{\Gamma }_{s}t}+D^2\frac{\cos (2{\Delta }_{i}t)}{({\Delta }_{0}t{)}^q}, \end{gathered}$$

$$\begin{gathered} \hat{H}=\sum _{k\sigma }{\xi }_k^s{\hat{a}}_{k\sigma }^{†}{\hat{a}}_{k\sigma }+\sum _k\left(\Delta {\hat{a}}_{k\uparrow }^{†}{\hat{a}}_{\bar{k}\downarrow }^{†}+{\Delta }^{\ast }{\hat{a}}_{\bar{k}\downarrow }{\hat{a}}_{k\uparrow }\right)+ \\ +\sum _{l\sigma }{\xi }_l^n{\hat{b}}_{l\sigma }^{†}{\hat{b}}_{l\sigma }+\sum _{kl\sigma }\left({\gamma }_{kl}{\hat{a}}_{k\sigma }^{†}{\hat{b}}_{l\sigma }+{\gamma }_{kl}^{\ast }{\hat{b}}_{l\sigma }^{†}{\hat{a}}_{k\sigma }\right) \end{gathered}$$

$$\Delta =\frac{\lambda }{V_s}\sum _k⟨{\hat{a}}_{k\uparrow }{\hat{a}}_{\bar{k}\downarrow }⟩,$$

$$\begin{array}{c}i\frac{\partial {\breve{G}}_k^s}{\partial t}-{\breve{H}}_k^s{\breve{G}}_k^s-{\breve{\Sigma }}^s\ast {\breve{G}}_k^s=\delta (t-t^{\prime }),\\ i\frac{\partial {\breve{G}}_l^n}{\partial t}-{\breve{H}}_l^n{\breve{G}}_l^n-{\breve{\Sigma }}^n\ast {\breve{G}}_l^n=\delta (t-t^{\prime }),\end{array}$$

$${\breve{G}}_{k/l}^{s/n}=\left(\begin{array}{cc}{\check{G}}_{k/l}^{s/n(R)}& {\check{G}}_{k/l}^{s/n(K)}\\ 0& {\check{G}}_{k/l}^{s/n(A)}\end{array}\right),$$

$${\breve{H}}_{k/l}^{s/n}=\left(\begin{array}{cc}{\check{H}}_{k/l}^{s/n}& 0\\ 0& {\check{H}}_{k/l}^{s/n}\end{array}\right),$$

$${\check{H}}_k^s=\left(\begin{array}{cc}{\xi }_k^s& \Delta \\ {\Delta }^{\ast }& -{\xi }_k^s\end{array}\right),\text{leavevmode}{\check{H}}_l^n=\left(\begin{array}{cc}{\xi }_l^n& 0\\ 0& -{\xi }_l^n\end{array}\right).$$

$$\Delta (t)={\Delta }_0+\delta \Delta (t)+i{\Delta }_0\delta \theta (t).$$

$$\begin{array}{c}\delta \Delta (\omega )=K_{\Delta }(\omega )\delta \Delta (\omega )+K^{\prime }(\omega ){\Delta }_0\delta \theta (\omega )\\ \delta \theta (\omega )=K^{\prime \prime }(\omega )\frac{\delta \Delta (\omega )}{{\Delta }_0}+K_{\theta }(\omega )\delta \theta (\omega )\end{array}$$

$$\delta \Delta (t)=\frac{1}{2\pi }\int \frac{f(\omega )e^{-i\omega t}d\omega }{1-K_{\Delta }(\omega )}.$$

$$\begin{gathered} \frac{1-K_{\Delta }(\omega )}{\lambda {\nu }_s}=i\frac{\sqrt{4{\Delta }_0^2-{\omega }_2^2}}{2{\omega }_2}\cdot \ln \left[\frac{F({\omega }_1,{\omega }_2)}{F(i{\Gamma }_s,{\omega }_2)}\right]+ \\ i\frac{\sqrt{4{\Delta }_0^2-{\omega }^2}}{2\omega }\ln \left[\frac{F(-i{\Gamma }_s,\omega )}{F({\omega }_1,\omega )}\cdot \frac{\omega +i\sqrt{4{\Delta }_0^2-{\omega }^2}}{-\omega +i\sqrt{4{\Delta }_0^2-{\omega }^2}}\right], \end{gathered}$$

$$\begin{gathered} \delta \Delta (t)\approx -i\sum _j\frac{e^{-i{\omega }_{j}t}{\mathrm{res}}_{{\omega }_j}f(\omega )}{1-K_{\Delta }({\omega }_j)}+ \\ +\frac{a\cos (2{\Delta }_{0}t-\alpha )}{\sqrt{{\Delta }_{0}t}}{e}^{-2{\Gamma }_{s}t}+\frac{{\Gamma }_s}{{\Delta }_0}\frac{b\cos ({\Delta }_{0}t-\beta )}{({\Delta }_{0}t{)}^{3/2}}{e}^{-{\Gamma }_{s}t}. \end{gathered}$$

$${\hat{A}}_k=\left(\begin{array}{c}{\hat{a}}_{k\uparrow }\\ {\hat{a}}_{\bar{k}\downarrow }^{†}\end{array}\right),{\hat{B}}_l=\left(\begin{array}{c}{\hat{b}}_{l\uparrow }\\ {\hat{b}}_{\bar{l}\downarrow }^{†}\end{array}\right).$$

$$\begin{array}{c}{\check{G}}_{k{k}^{\prime }}^{ss(R)}(t,t^{\prime })=-i\Theta (t-t^{\prime })⟨{\hat{A}}_k(t){\hat{A}}_{k^{\prime }}^{†}(t^{\prime })+{\hat{A}}_{k^{\prime }}^{†T}(t^{\prime }){\hat{A}}_k^T(t)⟩\\ {\check{G}}_{k{k}^{\prime }}^{ss(A)}(t,t^{\prime })=i\Theta (t^{\prime }-t)⟨{\hat{A}}_k(t){\hat{A}}_{k^{\prime }}^{†}(t^{\prime })+{\hat{A}}_{k^{\prime }}^{†T}(t^{\prime }){\hat{A}}_k^T(t)⟩\\ {\check{G}}_{k{k}^{\prime }}^{ss(K)}(t,t^{\prime })=-i⟨{\hat{A}}_k(t){\hat{A}}_{k^{\prime }}^{†}(t^{\prime })-{\hat{A}}_{k^{\prime }}^{†T}(t^{\prime }){\hat{A}}_k^T(t)⟩\end{array}$$

$${\breve{G}}_{k{k}^{\prime }}^{\alpha \beta }(t,t^{\prime })=\left(\begin{array}{cc}{\check{G}}_{k{k}^{\prime }}^{\alpha \beta (R)}(t,t^{\prime })& {\check{G}}_{k{k}^{\prime }}^{\alpha \beta (K)}(t,t^{\prime })\\ 0& {\check{G}}_{k{k}^{\prime }}^{\alpha \beta (A)}(t,t^{\prime })\end{array}\right),$$

$$\begin{array}{c}i\frac{\partial }{\partial t}{\breve{G}}_{k{k}^{\prime }}^{ss}-{\breve{H}}_k^s{\breve{G}}_{k{k}^{\prime }}^{ss}-\sum _l{\gamma }_{kl}{\breve{\tau }}_3{\breve{G}}_{l{k}^{\prime }}^{ns}=\delta (t-t^{\prime }){\delta }_{k{k}^{\prime }}\\ i\frac{\partial }{\partial t}{\breve{G}}_{lk}^{ns}-{\breve{H}}_l^n{\breve{G}}_{lk}^{ns}-\sum _{k^{\prime }}{\gamma }_{k^{\prime }l}^{\ast }{\breve{\tau }}_3{\breve{G}}_{k^{\prime }k}^{ss}=0\end{array},$$

$$\begin{array}{c}i\frac{\partial }{\partial t}{\breve{G}}_{l{l}^{\prime }}^{nn}-{\breve{H}}_n^s{\breve{G}}_{l{l}^{\prime }}^{nn}-\sum _k{\gamma }_{kl}^{\ast }{\breve{\tau }}_3{\breve{G}}_{k{l}^{\prime }}^{sn}=\delta (t-t^{\prime }){\delta }_{l{l}^{\prime }}\\ i\frac{\partial }{\partial t}{\breve{G}}_{kl}^{sn}-{\breve{H}}_k^s{\breve{G}}_{kl}^{sn}-\sum _{l^{\prime }}{\gamma }_{k{l}^{\prime }}{\breve{\tau }}_3{\breve{G}}_{l^{\prime }l}^{nn}=0\end{array},$$

$${\breve{H}}_k^s=\left(\begin{array}{cc}{\check{H}}_k^s& 0\\ 0& {\check{H}}_k^s\end{array}\right),\text{leavevmode}{\check{H}}_k^s=\left(\begin{array}{cc}{\xi }_k^s& \Delta \\ {\Delta }^{\ast }& -{\xi }_k^s\end{array}\right),$$

$${\breve{H}}_l^n=\left(\begin{array}{cc}{\check{H}}_l^n& 0\\ 0& {\check{H}}_l^n\end{array}\right),\text{leavevmode}{\check{H}}_l^n=\left(\begin{array}{cc}{\xi }_l^n& 0\\ 0& -{\xi }_l^n\end{array}\right),$$

$${\breve{\tau }}_3=\left(\begin{array}{cc}{\check{\tau }}_3& 0\\ 0& {\check{\tau }}_3\end{array}\right),\text{leavevmode}{\check{\tau }}_3=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).$$

$$\Delta (t)=\frac{i\lambda }{4V}\sum _k\mathrm{Tr}\left[\left({\check{\tau }}_1+i{\check{\tau }}_2\right){\check{G}}_{kk}^{ss(K)}(t,t)\right],$$

$${\check{\tau }}_1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\check{\tau }}_2=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right).$$

$${\breve{G}}_{lk}^{ns}=\sum _{k^{\prime }}{\gamma }_{k^{\prime }l}^{\ast }{\breve{\mathcal{G}}}_l^n{\breve{\tau }}_3\ast {\breve{G}}_{k^{\prime }k}^{ss},$$

$${\breve{G}}_{kl}^{sn}=\sum _{l^{\prime }}{\gamma }_{k{l}^{\prime }}{\breve{\mathcal{G}}}_k^s{\breve{\tau }}_3\ast {\breve{G}}_{l^{\prime }l}^{nn},$$

$$(X\ast Y)(t,t^{\prime })=\int X(t,t^{\prime \prime })Y(t^{\prime \prime },t^{\prime })d{t}^{\prime \prime }$$

$$i\frac{\partial }{\partial t}{\breve{\mathcal{G}}}_k^s-{\breve{H}}_k^s{\breve{\mathcal{G}}}_k^s=\delta (t-t^{\prime }),$$

$$i\frac{\partial }{\partial t}{\breve{\mathcal{G}}}_l^n-{\breve{H}}_l^n{\breve{\mathcal{G}}}_l^n=\delta (t-t^{\prime }).$$

$$i\frac{\partial }{\partial t}{\breve{G}}_{k{k}^{\prime }}^{ss}-{\breve{H}}_k^s{\breve{G}}_{k{k}^{\prime }}^{ss}-\sum _{k^{\prime \prime }}{\breve{\Sigma }}_{k{k}^{\prime \prime }}^s\ast {\breve{G}}_{k^{\prime \prime }{k}^{\prime }}^{ss}=\delta (t-t^{\prime }){\delta }_{k{k}^{\prime }},$$

$$i\frac{\partial }{\partial t}{\breve{G}}_{l{l}^{\prime }}^{nn}-{\breve{H}}_l^n{\breve{G}}_{l{l}^{\prime }}^{nn}-\sum _{l^{\prime \prime }}{\breve{\Sigma }}_{l{l}^{\prime \prime }}^n\ast {\breve{G}}_{l^{\prime \prime }{l}^{\prime }}^{nn}=\delta (t-t^{\prime }){\delta }_{l{l}^{\prime }},$$

$${\breve{\Sigma }}_{k{k}^{\prime \prime }}^s=\sum _l{\gamma }_{kl}{\gamma }_{k^{\prime \prime }l}^{\ast }{\breve{\tau }}_3{\breve{\mathcal{G}}}_l^n{\breve{\tau }}_3,$$

$${\breve{\Sigma }}_{l{l}^{\prime \prime }}^n=\sum _k{\gamma }_{kl}^{\ast }{\gamma }_{k{l}^{\prime \prime }}{\breve{\tau }}_3{\breve{\mathcal{G}}}_k^s{\breve{\tau }}_3.$$

$$\begin{array}{c}i\frac{\partial {\breve{G}}_k^s}{\partial t}-{\breve{H}}_k^s{\breve{G}}_k^s-{\breve{\Sigma }}^s\ast {\breve{G}}_k^s=\delta (t-t^{\prime }),\\ i\frac{\partial {\breve{G}}_l^n}{\partial t}-{\breve{H}}_l^n{\breve{G}}_l^n-{\breve{\Sigma }}^n\ast {\breve{G}}_l^n=\delta (t-t^{\prime }),\end{array}$$

Full text

The progress of modern nanotechnology opens new horizons for engineering superconducting correlations in various hybrid structures and creating, in fact, novel types of artificial superconducting materials with controllable properties Oreg et al. (2010); Qi and Zhang (2011); Sato and Ando (2017); Heersche et al. (2007); Sato et al. (2008); Kopnin et al. (2013a, b); Lee and Lee (2018); Lutchyn et al. (2018); Aasen et al. (2016); Alicea (2012). The proximity phenomenon arising in a non-superconducting material from the electron exchange with a primary superconductor can generate the induced superconducting ordering in a wide class of materials, including unconventional ones Oreg et al. (2010); Qi and Zhang (2011); Sato and Ando (2017); Heersche et al. (2007); Sato et al. (2008); Lee and Lee (2018); Lutchyn et al. (2018); Aasen et al. (2016); Alicea (2012). The resulting superconducting state in these materials can controllably reveal the exotic properties very rarely found in natural metals or alloys and strongly different from the ones of the primary superconductor. The induced Cooper pairs can change, e.g., their spin structure from the singlet to a triplet one in the presence of strong spin-orbit coupling and Zeeman (or exchange) field Oreg et al. (2010); Lutchyn et al. (2018); Aasen et al. (2016) This spin transformation affects, of course, the momentum space structure of pairs: the routine s-wave condensate can turn into an exotic p-wave one. The resulting Cooper pair structure leads to the formation of topological low-energy excitations such as Majorana fermions Oreg et al. (2010); Lutchyn et al. (2018); Aasen et al. (2016); Alicea (2012) and possesses a high potential for the development of new types of nanoelectronic devices perspective for applications in quantum computing, quantum information processing, quantum annealing, quantum memory and others Lutchyn et al. (2018); Aasen et al. (2016).

No wonder that the study of both equilibrium and nonequilibrium spectral and transport properties of these systems with engineered superconducting state has become recently one of the central research directions in condensed matter physics. While dc properties of these structures have been investigated in numerous theoretical and experimental works, the dynamic effects and, in particular, high frequency response remains an appealing problem which definitely deserves deeper understanding. Indeed, the limitations on the operating frequencies for above mentioned proximized devices Khaymovich et al. (2017); Semenoff and Sodano (2007); Tewari et al. (2008); Fu (2010); Bolech and Demler (2007); Nilsson et al. (2008) can be solely given by the dynamic characteristics of their induced superconducting ordering. The nonlinear dynamic effects are also known to provide a new route to the fascinating physics of coherent modification of the density of states etc. induced by the microwave irradiation Semenov et al. (2016); Devyatov and Semenov (2019)

Clearly, as typical quantum computing devices operate at the temperatures well below the gap of the primary superconductor, the study of the relaxation dynamics of the order parameter close to the critical temperature of the superconducting transition is irrelevant for their description. More adequate theory can be obtained by considering the so-called coherent quantum-mechanical dynamics of the system, which neglects inelastic scattering processes. In addition, in superconducting systems with the unconventional pairing the dynamic response is known to provide important information about the order parameter structure Bittner et al. (2015); Krull et al. (2016) working as a spectroscopic tool Fauseweh et al. (2017). The analogous method can provide an insight about the internal structure of the primary and induced Cooper pairing in superconducting hybrids.

Indeed, even the linear dynamic response of the superconductor near equilibrium provides a detection method of the superconducting gap structure via the coherent order parameter oscillations known also as Higgs modes Volkov and Kogan (1974); Yuzbashyan et al. (2005, 2006a); Yuzbashyan and Dzero (2006); Yuzbashyan et al. (2006b); Higgs (1964). The name is given due to the analogy to the Higgs boson in particle physics Higgs (1964). In the low-temperature limit these near-equilibrium Higgs oscillations of the the order parameter magnitude are described by the asymptotic long-time expression $\Delta-\Delta_{0}\sim\cos(2\Delta_{0}t)/\sqrt{\Delta_{0}t}$, where $\Delta_{0}$ is the superconducting gap in equilibrium Volkov and Kogan (1974); Yuzbashyan et al. (2006a, b); Barankov et al. (2004); Yuzbashyan et al. (2005); Yuzbashyan and Dzero (2006). The Higgs mode has been first detected using the Raman spectroscopy in the superconductors with charge density wave ordering Sooryakumar and Klein (1980); Méasson et al. (2014) as a peak at the frequency $2\Delta_{0}$ in the Raman spectrum of $2H$-$\mathrm{NbSe}_{2}$ below the superconducting transition temperature. Recent progress of the THz experimental techniques allowed the direct observation of the order parameter oscillations using the pump–probe method Matsunaga et al. (2013). The broad band pump excites a power–law relaxing Higgs mode at the frequency $2\Delta_{0}$ while the narrow band pump with a well defined frequency $\omega$ induces the oscillations of the order parameter at the frequency $2\omega$. The resonant third harmonic generation Matsunaga et al. (2014); Silaev (2019) provides an evidence of the Higgs mode excitation in the superconductors. An alternative method of the Higgs mode detection through the second harmonic measurements has been recently for the current carrying states Moor et al. (2017).

In this Letter we address an important effect of the quasiparticle spectrum of the superconducting system on long-time dynamic properties of Higgs modes and apply it to a system with proximity induced superconducting gap. Key physical phenomena related to the presence of the induced gap are demonstrated on the example of a junction of the superconductor (S) and normal metal (N) coupled via an insulating barrier (I) with a finite transparency. We analyze the distinctive features of the Higgs modes in this structure and make the predictions for the experimentally accessible relevant quantities.

Let’s first consider the qualitative picture of the Higgs dynamics in SIN system, Fig. 1, before proceeding with the further microscopic calculations. The Higgs modes in superconductors can be interpreted as a coherent splitting-recovery process of Cooper pairs. The energy difference between the ground state without quasiparticles and the excited state with two quasiparticles governs the frequency $2\Delta_{0}$ of this coherent superposition of above states as each unpaired quasiparticle brings an additional energy $\Delta_{0}$. Due to this simple qualitative reasoning we may expect the frequencies of the Higgs modes to be determined by the quasiparticle spectrum of the system. In SIN-structure the superconducting correlations penetrate to the normal metal and induce the hard gap $\Delta_{i}$ in the whole system McMillan (1968). $\Delta_{i}$ depends on the transparency of the barrier and the size of the normal subsystem.

We claim that there are three Higgs modes in the SIN-structure corresponding to three possible processes shown in Fig. 1. First, as in the isolated superconductor the Cooper pair may coherently split into two electrons both being located in the superconductor, Fig. 1(a). The energy $2\Delta_{0}$ of the these unpaired electrons determines the frequency of this process. However, in the presence of the normal metal the Cooper pair splitting can be accompanied by the coherent tunneling process of either or both electrons, see Fig. 1(b) and (c), respectively. The minimal energy of each electron which tunnels to the normal metal should be $\Delta_{i}$, so the frequency of the corresponding Higgs mode is given by $\Delta_{0}+\Delta_{i}$ and $2\Delta_{i}$ for $n_{e}=1$ and $2$ electrons tunneling to the normal metal, respectively. The amplitudes $A$ of these modes are expected to be reduced by the factor of $D^{n_{e}}$, where $D$ is transparency of the barrier (provided $D\ll 1$). On top of that in the first two processes, Fig. 1(a) and (b), the coherent superposition can be destroyed by the incoherent decay of each electron located in the superconductor to the normal metal (see green wavy lines in Fig. 1) because the hard gap in the spectrum of the whole system $\Delta_{i}$ is below the energy of quasiparticle in the superconductor $\Delta_{0}$. This effect results in an exponential damping of the Higgs modes with rate $\Gamma_{s}$ to each $\Delta_{0}$ frequency contribution. The value $\Gamma_{s}$ is an inverse lifetime of the electron in the superconductor determined by the tunneling rate from the superconductor to the normal metal. To sum up, the main result of our work can be written from the qualitative perspective as the following structure of the gap oscillations in SIN system:

[TABLE]

with certain power-law decay rates $p$ and $q$.

This qualitative picture can be confirmed by the direct microscopic calculations. The considered SIN system can be described by the following Hamiltonian:

[TABLE]

where $\hat{a}_{k\sigma}$ and $\hat{a}_{k\sigma}^{\dagger}$ are the electron annihilation and creation operators in the superconducting layer, $k$ is the index of the single–electron state and $\sigma$ is the projection of the electron spin, $\bar{k}$ denotes the index of the state obtained from the state $k$ by the time inversion operation. The operators $\hat{b}_{l\sigma}$ and $\hat{b}_{l\sigma}^{\dagger}$ are the electron annihilation and creation operators in the normal layer. The last term describes the tunneling between the superconductor and the normal metal. Assuming the insulating interlayer to be dirty so that the electron momentum is not conserved we consider the tunneling matrix elements $\gamma_{kl}$ to be the Gaussian uncorrelated random values $\langle\gamma_{kl}\gamma_{k^{\prime}l^{\prime}}^{\ast}\rangle=\gamma^{2}\delta_{kk^{\prime}}\delta_{ll^{\prime}}$ McMillan (1968); S-i . This model approximately describes tunneling junction if the magnitude of the tunneling matrix element $\gamma^{2}$ is proportional to $S/(V_{s}V_{n})$, where $S$ is the junction area, $V_{s}$ and $V_{n}$ are the volumes of the superconducting and normal subsystems, respectively. The superconducting order parameter $\Delta$ is given by the self-consistency equation:

[TABLE]

where $\lambda$ is the pairing constant and $V_{s}$ is the volume of the superconductor, the brackets $\langle\ldots\rangle$ denote a quantum-mechanical averaging. Here we assume the thickness of the superconducting subsystem to be small compared to the superconducting coherence length so we can consider $\Delta$ to be homogeneous within the sample.

We study the dynamics of the system using non-equilibrium Keldysh technique. Following approach developed in the Refs. McMillan (1968); Volkov et al. (1995); Kopnin and Melnikov (2011); Kopnin et al. (2013a, b) we treat the tunneling operator within the self-consistent Born approximation and write the following set of the equations for the Green’s functions in the superconductor and normal metal:

[TABLE]

where the Green’s function of the superconductor $\breve{G}_{k}^{s}$ and the Green’s function of the normal metal $\breve{G}_{l}^{n}$ are the $4\times 4$ matrices in the Keldysh–Nambu space

[TABLE]

and $\breve{H}_{k}^{s}$ and $\breve{H}_{l}^{n}$ are the single mode Hamiltonians of the superconductor and the normal metal:

[TABLE]

The self-energies of the superconductor $\breve{\Sigma}^{s}$ and the normal metal $\breve{\Sigma}^{n}$ describe the tunneling between two subsystems and self-consistently expressed through the Green’s functions $\breve{G}_{l}^{n}$ and $\breve{G}_{k}^{s}$, respectively. These self-energies are proportional to the tunneling rates of electrons from the superconductor $\breve{\Sigma}_{s}\propto\Gamma_{s}=\pi\gamma^{2}\nu_{n}V_{n}$ and the normal metal $\breve{\Sigma}_{n}\propto\Gamma_{n}=\pi\gamma^{2}\nu_{s}V_{s}$, where $\nu_{s}$ and $\nu_{n}$ are the densities of normal states per the unit volume in the superconductor and the normal metal at the Fermi level, respectively. The derivation of the Dyson equations (32) and the explicit form of the self-energies are given in Appendix.

In order to study the near–equilibrium dynamics of the system we expand the order parameter near $\Delta_{0}$ assuming its phase to be zero in equilibrium without loss of the generality:

[TABLE]

Here $\delta\Delta$ and $\delta\theta$ are the perturbations of the magnitude and the phase of $\Delta$, respectively. One can introduce the corresponding perturbations to the Green’s functions and the self-energies and linearize the system of the equations for the Green’s functions with respect to the perturbation of the superconducting order parameter. Along with the linearized self-consistency equation (21) this system can be reduced to the equations for the eigenmodes of the order parameter in the Fourier form (see Appendix for details):

[TABLE]

The off-diagonal kernels $K^{\prime}$ and $K^{\prime\prime}$ are equal exactly to zero in the systems which have electron–hole symmetry and can be neglected if the Fermi level both in the superconductor and the normal metal is far from the van Hove singularities in the density of states von . The singular points of $[1-K_{\Delta}(\omega)]^{-1}$ and $[1-K_{\theta}(\omega)]^{-1}$ correspond to the frequencies of the Higgs and Anderson-Bogoliubov modes of the superconductor. The perturbations of the magnitude and the phase of the order parameter are completely independent so hereafter we focus only on the study of the Higgs modes.

The linear response of the superconducting order parameter to an external force $f(\omega)$ takes the following form:

[TABLE]

Previous studies Tsuji and Aoki (2015); Matsunaga et al. (2013, 2014); Silaev (2019) show that this force can originate, e.g., from the pulses of the external electromagnetic field. In Fig. 4 the typical frequency dependencies of the real and imaginary parts of $[1-K_{\Delta}(\omega)]^{-1}$ are shown for some particular values of the tunneling rates. The spectrum of the Higgs modes appears to be consistent with the picture shown in Fig. 1. The broadened features at the frequencies $\omega\approx\pm 2\Delta_{0}$ and $\omega\approx\pm(\Delta_{0}+\Delta_{i})$ are seen clearly, while the singularity at $\omega=\pm 2\Delta_{i}$ can be seen only in the zoomed inset. The latter singularity corresponds to the low frequency Higgs mode, Fig. 1(c). This mode has no exponential damping as $\Delta_{i}$ is a hard gap of the whole system (the singular point of the kernel is exactly at the real axis). It means that this low-frequency mode gives the major contribution to the oscillations of the order parameter in the long time limit, $t\gg\Gamma_{s}^{-1},\Gamma_{n}^{-1}$, however the amplitude of this mode is few orders of magnitude lower than the amplitude of the usual Higgs mode due to the low transparency of the insulating barrier.

The kernel $K_{\Delta}$ can be evaluated analytically in the zero temperature limit and $\Gamma_{n}=0$ which corresponds to the bulk normal metal $V_{n}\to\infty$ with the suppressed induced superconducting ordering and vanishing induced gap $\Delta_{i}=0$:

[TABLE]

where $\omega_{1}=\omega+i\Gamma_{s}$, $\omega_{2}=\omega+2i\Gamma_{s}$, and $F(x,y)=2\Delta_{0}^{2}-xy+\sqrt{\Delta_{0}^{2}-x^{2}}\sqrt{4\Delta_{0}^{2}-y^{2}}$. In the isolated superconductor, $\Gamma_{s}=0$, the kernel as a function of the complex frequency $\omega$ has two branch points at $\omega=\pm 2\Delta_{0}$ which give the usual polynomially damped Higgs mode in the superconductor. In the presence of tunneling, $\Gamma_{s}>0$, these branch points shift to the points $\omega=\pm 2\Delta_{0}-2i\Gamma_{s}$ corresponding to exponential damping of this mode. Moreover, two additional branch points at $\omega=\pm\Delta_{0}-i\Gamma_{s}$ appear triggering a new Higgs mode, which has already revealed itself in Fig. 4. The absence of the low frequency mode at $2\Delta_{i}$ is rather expected for the considered case $\Gamma_{n}=0$ because the electrons of the normal metal are not affected by the proximity with the superconductor and, thus, cannot form a Cooper pair as in Fig. 1(c). As a result they do not contribute to the order parameter oscillations. The asymptotic behavior of the order parameter perturbation in the intermediate-time limit $\Delta_{0}^{-1}\ll t\ll\Gamma_{s}^{-1}$ reads as

[TABLE]

The details of calculations, restrictions on the analytic properties of the exciting force $f(\omega)$, and the expressions of parameters $a$, $b$, $\alpha$ and $\beta$ depending on the exciting force are given in Appendix. The first term in (12) describes the forced oscillations of the order parameter which occur at the frequencies corresponding to the poles of the Fourier spectrum of the external force $f(\omega)$. The finite tunneling rate $\Gamma_{s}$ leads to the exponential damping of the oscillations of the order parameter and appearance of a new Higgs mode at the frequency $\Delta_{0}$ with the magnitude suppressed by factor $\Gamma_{s}/\Delta_{0}$. This mode corresponds to the middle-frequency mode $\Delta_{0}+\Delta_{i}$ shown in Fig. 1(b) in the limit $\Delta_{i}=0$. Let’s emphasize once again that all the above analysis was related to the case of the small tunneling rates $\Gamma_{s}\ll\Delta_{0}$. In the opposite limit $\Gamma_{s}\gg\Delta_{0}$ we get the gapless superconductivity with a pure imaginary frequency of the Higgs mode: $\omega=-i\Delta_{0}^{2}/\Gamma_{s}$. The relaxation of the order parameter in this limit is described by the decaying exponent $\delta\Delta(t)\propto\exp(-\Delta_{0}^{2}t/\Gamma_{s})$.

The Higgs modes of the SIN system can be studied using a pump–probe experiment similar to the one developed in the Ref. Matsunaga et al. (2013). The measured $\delta\Delta(t)$ can be analyzed using the Fourier transform. The Fourier spectrum is expected to have two peaks at $\omega\approx 2\Delta_{0}$ and $\omega\approx\Delta_{0}+\Delta_{i}$ corresponding to the above Higgs modes of the SIN system. In the limit of low transparency $D\ll 1$ the mode with the frequency $2\Delta_{i}$ has too low magnitude so it’s experimental observation can be hampered, however it may be visible at intermediate transparencies $D\sim 1$, $\Delta_{i}\lesssim\Delta_{0}$. Another way to detect the Higgs modes is the experimental studying of the frequency dependence of the third harmonic generation Silaev (2019); Matsunaga et al. (2014). The electromagnetic wave with the frequency $\Omega$ excites the Higgs mode with the frequency $2\Omega$. The magnitude of the generated nonlinear signal with the frequency $3\Omega$ should depend on the magnitude of the oscillations of the order parameter and therefore one can expect the appearance of the broadened resonances in the third harmonic response if the frequency $2\Omega$ is close to the frequency of any of the Higgs modes. However, the effect of generation of non-equilibrium quasiparticles by the electromagnetic radiation with the frequency $\Omega>2\Delta_{i}$ may complicate the observation of the resonant effect in the third harmonic generation. The overheating can be significantly reduced at intermediate transparencies $D\sim 1$ when the induced gap is high enough $3\Delta_{i}>\Delta_{0}$. This improves the observability of the resonance at $2\Omega=\Delta_{0}+\Delta_{i}$ predicted above.

To sum up, the effect of the quasiparticle spectrum of the superconducting system on the low-temperature dynamics of the order parameter has been studied on an example of a hybrid superconductor-insulator-normal metal system. Two new Higgs modes in such system have been discovered. The frequencies of these modes are formed by sums of the two quasiparticle energies, which are in a good correspondence with the qualitative interpretation of the Higgs modes as of coherent processes of splitting and recovery of the Cooper pairs, Fig. 1. The proposals of experimental observation of these new Higgs modes in hybrid SIN system have been developed basing on the existing THz techniques used for study of the Higgs modes in the pure superconductors.

This work was supported, in part, by the Russian Foundation for Basic Research under Grant No. 17-52-12044, the Foundation for the Advancement to Theoretical Physics and Mathematics “BASIS” Grant No. 17-11-109, and the German-Russian Interdisciplinary Science Center (G-RISC) funded by the German Federal Foreign Office via the German Academic Exchange Service (DAAD). In the part of numerical calculations, the work was supported by Russian Science Foundation (Grant No. 17-12-01383). I. M. K. acknowledges the support of the German Research Foundation (DFG) Grant No. KH 425/1-1. V. L. V. and A. S. M. appreciate warm hospitality of the Max-Planck Institute for the Physics of Complex Systems, Dresden, Germany, extended to them during their visits when this work was done.

Appendix A Derivation of the Dyson equations

In this section we derive the equations for the Green’s functions of the superconductor and the normal metal. We define the Nambu pseudospinor operators as follows:

[TABLE]

Using these operators we define the retarded, advanced and Keldysh Green’s functions of the superconductor in the following way:

[TABLE]

The Green’s functions of the normal metal $G_{ll^{\prime}}^{ll(RAK)}$ and the tunneling Green’s functions $G_{kl}^{sn(RAK)}$ and $G_{lk}^{ns(RAK)}$ are defined in the same way replacing the operators $\hat{A}$ with $\hat{B}$. One can construct $4\times 4$ matrix Green’s functions in the Keldysh–Nambu space in the usual way:

[TABLE]

where the indices $\alpha$ and $\beta$ denote $s$ and $n$. The equations for the Green’s functions read as:

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

The selfconsistency condition reads as

[TABLE]

where

[TABLE]

Using the Green’s functions of the isolated superconductor and the normal metal $\breve{\mathcal{G}}_{k}^{s}$ and $\breve{\mathcal{G}}_{l}^{n}$ we may eliminate the tunneling Green’s functions:

[TABLE]

where $\ast$ denotes the convolution operation

[TABLE]

and the Green’s functions $\breve{\mathcal{G}}_{k}^{s}$ and $\breve{\mathcal{G}}_{l}^{n}$ satisfy the following equations:

[TABLE]

Thus, we can write two independent equations for the Green’s functions in the superconductor and the normal metal:

[TABLE]

[TABLE]

where the self-energies of the superconductor and the normal metal are:

[TABLE]

[TABLE]

These equations written in the diagram form are shown in the Fig. 3(a). These equations are not practical to use as they contain the tunneling matrix elements $\gamma_{kl}$ which are the random numbers. One can average the equations over the random matrix elements, the average of the product of the matrix elements can be expanded as sums of correlators due to the Wick theorem, thus, the averaged Green’s function can be written as a sum of diagrams. We omit the diagrams with the intersecting correlators (thus, we neglect the vertex corrections) and take account only of the diagrams with the consequent and nested correlators. Such approach is known as selfconsistent Born approximation. The diagrammatic form of the Dyson equation for this approximation is shown in the Fig. 3(b). After the averaging the self-energies and the Green’s function appear to be diagonal in the normal mode picture and obey the following equations:

[TABLE]

[TABLE]

In the wide band approximation the sum over the normal modes can be replaced with the integral over the normal energy $\sum_{k}\to\nu_{s}V_{s}\int\;d\xi_{k}^{s}$ ($\sum_{l}\to\nu_{n}V_{n}\int\;d\xi_{l}^{n}$):

[TABLE]

Here we have introduced the tunneling rates $\Gamma_{s}=\pi\gamma^{2}\nu_{n}V_{n}$ and $\Gamma_{n}=\pi\gamma^{2}\nu_{s}V_{s}$.

Appendix B Derivation of the equations for the eigenmodes

Let us introduce the perturbations of the Green’s functions and the self-energies with respect to the perturbation of the magnitude $\delta\Delta$ and the phase $\delta\theta$ of the superconducting order parameter:

[TABLE]

where $\breve{G}_{0k}^{s}$, $\breve{G}_{0l}^{n}$ and $\breve{\Sigma}^{s/n}$ are the equilibrium Green’s functions and self-energies of the superconductor and the normal metal with the account of tunneling. The closed set of equations for the linear perturbations of the Green’s functions and the self-energies reads as:

[TABLE]

The perturbation to the single mode Hamiltonian of the superconductor is $\delta\breve{H}^{s}=\delta\Delta\breve{\tau}_{1}-i\Delta_{0}\delta\theta\breve{\tau}_{2}$. One can easily solve the equations for $\delta G$ in the Fourier form:

[TABLE]

where the Fourier transform of the Green’s functions is defined as follows:

[TABLE]

We introduce the quasiclassic Green’s function $\delta\breve{g}^{s}=\int\delta\breve{G}_{k}^{s}\;d\xi_{k}^{s}$ and write an algebraic equation for it:

[TABLE]

The above equation can be considered as a system of 12 linear equations for the 12 components of matrix Green’s function $\breve{g}(\omega,\omega^{\prime})$ (retarded, advanced and Keldysh, each of them is $2\times 2$ matrix). Each integral can be evaluated analytically because the equilibrium Green’s functions are rational functions of the normal energies $\xi_{k}^{s}$ and $\xi_{l}^{n}$. The solution can be written in the following form:

[TABLE]

One should use the selfconsistency equation (21)

[TABLE]

and obtain the expression for all the kernels in the equation (7) of the main text of the paper:

[TABLE]

Appendix C Bulk normal metal $\Gamma_{n}=0$

In the case $\Gamma_{n}=0$ the self-energies of the superconductor are determined by the equilibrium Green’s functions of the normal metal, and, assuming zero temperature limit we have:

[TABLE]

The self-energies are constant $\delta\Sigma^{s}=0$ so the solution of the equations (36)

[TABLE]

using the selfconsistency equation (21) one can obtain an expression for the kernel of the Higgs mode:

[TABLE]

This integral diverges logarithmically, however it can be regularized using the equilibrium selfconsistency equation for $\Delta_{0}$. Finally, one can obtain the Eq. (11) of the main text of the paper:

[TABLE]

The plots of $[1-K_{\Delta}(\omega)]^{-1}$ are shown in Fig. 4 for the cases of the (a) low and (b) high tunneling rates.

C.1 Derivation of the long time asymptotics of the order parameter

We suppose that the external force $f(\omega)$ is a meromorphic function of the frequency $\omega$ one can close up the integral over the real freqeuncies into a contour integral as it is shown in the Fig. 5:

[TABLE]

where $\mathcal{C}$ is the integration contour, $\omega_{j}$ are the poles of the external force $f(\omega)$ within the contour $\mathcal{C}$ and $\mathop{\mathrm{res}}_{\omega_{j}}f(\omega)$ is the residue of $f(\omega)$ at the pole $\omega_{j}$, thus, the integral in the equation (8) of the main paper can be expressed as a sum of the integrals along the branch cuts and the terms with the residues of the external force:

[TABLE]

The integrals along the branch cuts can be evaluated approximately assuming $f(\omega)$ are regular near the branch points of $K_{\Delta}(\omega)$ and that the main contribution to these integrals comes from the from the vicinity of the singularities. The expansion of the kernel $K_{\Delta}(\omega)$ near its branch points at $\omega=\pm 2\Delta_{0}-2i\Gamma_{s}$ and $\omega=\pm\Delta_{0}-i\Gamma_{s}$ in the limit $\Gamma_{s}\ll\Delta_{0}$ reads as follows:

[TABLE]

Using these expansions one can finally obtain an expression for the near equilibrium oscillations of the superconducting gap:

[TABLE]