Post-quench gap dynamics of two-band superconductors

TL;DR

This paper develops a new analytical method to study the non-adiabatic gap dynamics in two-band superconductors after a quench, revealing faster relaxation and beating patterns due to interband dephasing.

Contribution

The authors generalize the Volkov-Kogan method to analyze two-band superconductor dynamics, extending beyond single-band models and capturing interband effects.

Findings

Dephasing accelerates gap relaxation to a $t^{-3/2}$ power law.

Gap oscillations exhibit beating patterns from multiple asymptotic values.

Method accurately reproduces single-band solutions and aligns with numerical two-band results.

Abstract

Recent experimental progress in the fields of cold quantum gases and ultrafast optical spectroscopy of quantum materials allows to controllably induce and probe non-adiabatic dynamics of superconductors and superfluids. The time-evolution of the gap function before relaxation with the lattice is determined by the superposition of coherently evolving individual Cooper pairs within the manifold of the Bardeen-Cooper-Schrieffer (BCS) wavefunction. While dynamics following an abrupt quench of the pairing interaction strength in the single-band BCS model has been exactly solved due to the integrability of the model, the dynamics of post-quench multi-band superconductors remain under scrutiny. Here, we develop a generalization of the Volkov-Kogan Laplace-space perturbative method that allows us to determine the non-adiabatic gap dynamics of two-band fully gapped superconductors for a wide range of quench amplitudes. Our approach expands the long-time dynamics around the steady-state asymptotic value of the gap, which is self-consistently determined, rather than around the equilibrium value of the gap. We explicitly demonstrate that this method recovers the exact solution of the long-time gap dynamics in the single-band case and perfectly agrees with a numerical solution of the two-band model. We discover that dephasing of Cooper pairs from different bands leads to faster collisionless relaxation of the gap oscillation with a power-law of $t^{-3/2}$ instead of the well-known $t^{-1/2}$ behavior found in the single-band case. Furthermore, the gap oscillations display beating patterns arising from the existence of two different asymptotic gap values. Our results have important implications to a variety of two-band superconductors driven out of equilibrium, such as iron-based superconductors, MgB$_{2}$, and SrTiO$_{3}$.

Equations

$$H_{\text{BCS}}$$

$${\Delta }_{\alpha }={\Delta }_{\alpha }^{\prime }+i{\Delta }_{\alpha }^{\prime \prime }=-\frac{1}{N}\sum _{\mathbf{p},\beta }{V}_{\alpha \beta }⟨c_{-\mathbf{p},\downarrow ,\beta }{c}_{\mathbf{p},\uparrow ,\beta }⟩$$

$$\left(\begin{array}[]{c}\Delta_{1}\\ \Delta_{2}\end{array}\right)=\hat{\gamma}v\left(\begin{array}[]{c}\int_{-\Lambda}^{\Lambda}d\varepsilon\frac{\Delta_{1}}{2E_{1}}\tanh\left(\frac{E_{1}}{2T}\right)\\ \int_{-\Lambda}^{\Lambda}d\varepsilon\frac{\Delta_{2}}{2E_{2}}\tanh\left(\frac{E_{2}}{2T}\right)\end{array}\right)$$

$$\hat{\gamma}=\left(\begin{array}[]{cc}r&-\eta\\ -1&r\eta\end{array}\right)$$

$$H_{\text{BCS}}=-\sum _{\mathbf{k},\alpha }{\mathbf{B}}_{\mathbf{k},\alpha }\cdot {\hat{\mathbf{S}}}_{\mathbf{k},\alpha }+\text{const.}$$

$${\hat{S}}_{\mathbf{k},\alpha }^{-}$$

$${\hat{S}}_{\mathbf{k},\alpha }^{+}$$

$${\hat{S}}_{\mathbf{k},\alpha }^z$$

$${\Delta }_{\alpha }=-\frac{1}{N}\sum _{\mathbf{k},\beta }{V}_{\alpha \beta }{S}_{\mathbf{k},\beta }^{-}$$

$$S_{\mathbf{k},\alpha }^x$$

$$S_{\mathbf{k},\alpha }^y$$

$$S_{\mathbf{k},\alpha }^z$$

$$\frac{d}{dt}{\mathbf{S}}_{\mathbf{k},\alpha }\left(t\right)={\mathbf{S}}_{\mathbf{k},\alpha }\left(t\right)\times {\mathbf{B}}_{\mathbf{k},\alpha }\left(t\right)$$

$${\Delta }_{\alpha }\left(t\right)=v\left(t\right)\sum _{\beta }{\gamma }_{\alpha \beta }\int d\epsilon S_{\beta }^{-}\left(\epsilon ,t\right)$$

$$S_{\alpha }^z\left(\epsilon ,t\right)$$

$$S_{\alpha }^{-}\left(\epsilon ,t\right)$$

$${\Delta }_{\alpha }\left(t\right)$$

$${\delta }_{\alpha }\left(t\right)=v_f\sum _{\beta }{\gamma }_{\alpha \beta }{\int }_{-\Lambda }^{\Lambda }d\epsilon f_{\beta }\left(\epsilon ,t\right)$$

$${\dot{f}}_{\alpha }^{\prime }$$

$${\dot{f}}_{\alpha }^{\prime \prime }$$

$${\dot{g}}_{\alpha }$$

$${\ddot{f}}_{\alpha }^{\prime \prime }$$

$${\stackrel{...}{f}}_{\alpha }^{\prime }$$

$$f_{\alpha }^{\prime \prime }\left(s\right)+\frac{2s{S}_{\alpha ,\infty }^z}{s^2+4E_{\alpha ,\infty }^2}{\delta }_{\alpha }\left(s\right)$$

$$f_{\alpha }^{\prime }\left(s\right)-\frac{-4\epsilon S_{\alpha ,\infty }^z}{s^2+4E_{\alpha ,\infty }^2}{\delta }_{\alpha }\left(s\right)$$

$$\left({\hat{\Phi }}^{\infty }\left(s\right)+\hat{\mathcal{M}}\right)\vec{\delta }\left(s\right)=\frac{\vec{I}\left(s\right)}{s},$$

$${\hat{\Phi }}_{\alpha \beta }^{\infty }\left(s\right)$$

$${\hat{\mathcal{M}}}_{\alpha \beta }$$

$$⟨\dots ⟩=v_f\int d\epsilon \left(\dots \right).$$

$${\Phi }_{\alpha }^{\infty }(s)=\left(s^2+4{\Delta }_{\alpha ,\infty }^2\right)⟨\frac{S_{\alpha ,\infty }^x/{\Delta }_{\alpha ,\infty }}{s^2+4E_{\alpha ,\infty }^2}⟩$$

$$I_{\alpha }\left(s\right)$$

$$\times \left[{\Phi }_{\alpha }^i\left(s\right)-\frac{v_f}{v_i}\sum _{\beta }{\left({\hat{\gamma }}^{-1}\right)}_{\alpha \beta }\frac{{\Delta }_{\beta ,i}}{{\Delta }_{\alpha ,i}}\right],$$

$${\Phi }_{\alpha }^i\left(s\right)=\left(s^2+4{\Delta }_{\alpha ,\infty }^2\right)⟨\frac{S_{\alpha ,i}^x/{\Delta }_{\alpha ,i}}{s^2+4E_{\alpha ,\infty }^2}⟩$$

$$\vec{\delta }\left(s\right)={\left({\hat{\Phi }}^{\infty }\left(s\right)+\hat{\mathcal{M}}\right)}^{-1}\frac{\vec{I}\left(s\right)}{s}$$

$${\Phi }_{\alpha }^i\left(s\right)={\rm Y}\left({\tilde{\Delta }}_{\alpha ,i},\frac{s}{2{\Delta }_{\alpha ,\infty }}\right)$$

$${\rm Y}\left(\Delta ,x\right)=v_f\frac{\sqrt{\frac{x^2+1}{{\Delta }^2}}\arccos \left(\sqrt{\frac{x^2+1}{{\Delta }^2}}\right)}{\sqrt{1-\frac{1+x^2}{{\Delta }^2}}}$$

Full text

I Introduction

Superconductors that are perturbed into a state away from equilibrium display an extremely rich and interesting dynamical behavior. This originates from the interplay between the dynamics of its fermionic quasi-particle excitations and that of the superconducting order parameter, as expressed, for example, in the superconducting gap equation. Close to equilibrium, various collective modes emerge such as the Anderson-Bogoliubov phase mode Anderson (1958a) and the longitudinal Schmid (or Higgs) amplitude mode Schmid (1966); Pekker and Varma (2015), which describe phase and amplitude fluctuations of the order parameter. The transverse Carlson-Goldman mode describes the coupled oscillations of normal currents and supercurrents Carlson and Goldman (1975); Pals et al. (1982), whereas in multi-gap superconductors additional Leggett phase modes appear Leggett (1966), corresponding to oscillations of the relative phases of the different gaps. Interesting dynamics also occurs farther away from equilibrium, where one observes, for example, intriguing non-linear behaviors such as dynamic instabilities towards slowly damped Volkov and Kogan (1974); M.Gal’perin et al. (1981) or even undamped order parameter oscillations Barankov et al. (2004); Yuzbashyan et al. (2005a, 2006); Barankov and Levitov (2006).

Generally, the dynamic response of a superconductor depends on the type of perturbation that is applied, for example, whether it is adiabatic or non-adiabatic, linear or non-linear, and whether it is charge neutral or charged. It also depends on the hierarchy of a number of important timescales such as the quasi-particle energy relaxation time $\tau_{\varepsilon}$, the dynamical scale of the superconducting order parameter $\tau_{\Delta}$, the timescale of the external perturbation $\tau_{\text{pert}}$ and the characteristic observation time $t$ Aronov et al. (1981); Langenberg and Larkin (1986); Kopnin (2001). Here, we are interested in the case of $\tau_{\Delta}\ll\tau_{\varepsilon}$ and in fast, non-adiabatic perturbation occurring on a timescale $\tau_{\text{pert}}\ll\tau_{\Delta}\approx t$. This non-adiabatic, collisionless regime has been explored in a linearized approach close to equilibrium in the seminal work by Volkov and Kogan Volkov and Kogan (1974), who studied the gap dynamics of a single-band superconductor following a small and instantaneous perturbation. They found coherent gap oscillations that are only algebraically damped $\propto t^{-1/2}$, analogous to Landau damping in a collisionless plasma Landau (1946); Kamenev (2011). More recently, experimental progress on two distinct fronts have brought renewed interest to this field: (1) Ultrafast optical studies in the Terahertz regime have unveiled non-adiabatic, coherent gap dynamics in thin superconducting films Shimano and Tsuji (2019), for example, in NbN Matsunaga and Shimano (2012); Matsunaga et al. (2013, 2014); Beck et al. (2011) and Nb3Sn Yang et al. (2018); Cui et al. (2019); Yang et al. (2019); (2) Cold-atom realizations of superfluids and Bose-Einstein condensates have provided a fruitful avenue to induce non-adiabatic dynamics by performing rapid parameter changes such as quenching the pairing interaction strength Langen et al. (2015); Bloch et al. (2008).

The situation of a rapid parameter quench is theoretically particularly interesting as it is amenable to analytical approaches. Going beyond the linear analysis of Volkov and Kogan and exploiting the integrability of the Bardeen-Cooper-Schrieffer (BCS) Hamiltonian Richardson and Sherman (1964); Gaudin and Caux (2014); Dukelsky et al. (2004), a number of works have explored post-quench non-adiabatic dynamics of single-band BCS superconductors far-away from equilibrium Barankov et al. (2004); Yuzbashyan et al. (2005a, b); Barankov and Levitov (2006); Yuzbashyan et al. (2015, 2005c). It was discovered that non-equilibrium dynamics at times $\tau_{\Delta}\ll t\ll\tau_{\varepsilon}$ fall into one of three distinct classes (or ”phases”) Barankov and Levitov (2006); Yuzbashyan and Dzero (2006), which can be topologically distinguished by the number of complex roots of the spectral polynomial Yuzbashyan et al. (2005a, 2006): Phase I, where the gap decays exponentially to zero; Phase II, where the gap oscillates with frequency $2\Delta_{\infty}$ and decays algebraically $\propto t^{-1/2}$ to a finite value $\Delta_{\infty}$; and phase III, where persistent undamped gap oscillations occur. Phase II in the non-equilibrium quench phase diagram Barankov and Levitov (2006); Yuzbashyan and Dzero (2006) contains the linear regime around equilibrium studied by Volkov and Kogan Volkov and Kogan (1974). Finally, we note that the topological classification explains why terahertz induced gap dynamics is qualitatively similar to the case of a parameter quench, as has been observed in various numerical studies Chou et al. (2017); Papenkort et al. (2007, 2009); Krull et al. (2014); Akbari et al. (2013); Cui et al. (2019).

In this paper, we extend these previous studies by addressing numerically and analytically the gap dynamics of two-band superconductors following an interaction quench. Our motivation is on the fact that multi-band superconductivity is realized in a variety of materials with conventional and unconventional pairing mechanisms. Primary examples are MgB2 Bud’ko and Canfield (2015), the iron-based superconductors Chubukov et al. (2008), Sr2RuO4 Mackenzie and Maeno (2003), heavy fermions Wirth and Steglich (2016), strontium titanate Fernandes et al. (2013); Trevisan et al. (2018), and oxide heterostructures such as LaAlO3/SrTiO3 Scheurer and Schmalian (2015). Unconventional multi-orbital superfluidity has also been reported in cold-atom setups on the honeycomb lattice Soltan-Panahi et al. (2012). While different superconducting gap symmetries are possible in the presence of multiple Fermi surfaces, we will focus on the simplest case of $s$-wave superconductivity. As the quench dynamics is identical for $s^{+-}$ and $s^{++}$ pairing, corresponding to gaps with opposite or same signs on the two Fermi surfaces, respectively, our results apply to both cases. Quenches in two-band $s$-wave superconductors have so far only been studied numerically Akbari et al. (2013); Krull et al. (2016); Dzero et al. (2015), focusing on the coupling between the Higgs and the Leggett mode Krull et al. (2016) or the competition between superconductivity and spin-density wave order Dzero et al. (2015). Generalizations to quenches between other pairing symmetries such as time-reversal symmetry breaking $s+is$ or $s+id$ pairing are interesting avenues for further work. Indeed, a recent numerical study of Terahertz induced gap dynamics for $s+is$ pairing has revealed an unusual coupling between the Higgs amplitude and the Leggett (relative) phase mode Müller et al. (2018).

Exact solutions of the time-dependent two-band BCS model only exist for special fine tuned values of the intra- and inter-band interaction parameters, where the problem effectively reduces to the single-band case (see below and Ref. Dzero et al., 2015). It is an open question whether the generic two-band BCS model is integrable. Here, we develop a generalization of the Volkov-Kogan Laplace-space analysis in order to investigate the non-adiabatic post-quench dynamics in generic two-band BCS models. Like Volkov and Kogan we solve linearized equations of motion in Laplace space, but an important distinction of our work is that we expand around the long-time steady-state of the system instead of the equilibrium state. This allows us to explore the gap dynamics away from the weak-quench limit in a larger region of the non-equilibrium phase diagram. We achieve this methodological advancement by self-consistently solving for the steady-state value of the superconducting gap $\Delta_{\infty}$. We show in detail that our method reproduces the exact solution in phase II of the single-band model. For the two-band model we carefully check our analytical results by comparing to the numerical solution of the dynamics. We find that the oscillatory gap dynamics exhibits pronounced beating behavior due to the presence of two asymptotic gap values $\Delta_{1,\infty}$ and $\Delta_{2,\infty}$, which has been previously reported in a numerical investigation of Terahertz driven gap oscillations in two-band superconductors Akbari et al. (2013). A central new result of our work is that the decay of the gap oscillations due to Landau damping in two-band superconductors is governed by a power-law $\propto t^{-3/2}$ that is different from the one found in the single band case, where it is $\propto t^{-1/2}$ (see Fig. 1). Earlier numerical studies of multi-gap superconductors have reported power-law decays of $t^{-1/2}$, although in that case the dynamics was driven not a by an interaction quench, but by laser pulses Akbari et al. (2013). Interestingly, faster than $t^{-1/2}$ decay was also seen in the case of superconducting nanowires, where electronic subbands arise due to confinement  Zachmann et al. (2013). Finally, a similar $t^{-3/2}$ decay of the pairing amplitude has been found in quenches into the strong pairing (Bose-Einstein condensation (BEC)) regime in three dimensions, but by a different microscopic mechanism Gurarie (2009); Yuzbashyan et al. (2015).

The remainder of the paper is organized as follows: in Sec. II, we define the two-band BCS model and formulate it in terms of Anderson pseudospins. We then derive equations of motion of the pseudospins that govern the non-adiabatic dynamics of individual Cooper pairs and the gap following an instantaneous quench of the BCS coupling strength. In Sec. III, we present numerical solutions of the gap dynamics in the regime of weak quenches, which show the main features of oscillatory beating and algebraic decay $\propto t^{-3/2}$. In Sec. IV, we present our main analytical calculation and flesh out the details of our method to find the long-time dynamics of the gap using a self-consistent Laplace analysis. In Sec. IV.1, we derive linearized equations of motion around the long-time steady-state. We present the solution of these equations in Laplace space in Sec. IV.2, which depends on the steady-state values of the gap $\Delta_{\alpha,\infty}$ and the pseudo-spins $S_{\alpha,\infty}^{i}$. These values are determined in Sec. IV.3 by solving self-consistent equations via an ansatz for the non-equilibrium distribution function in the steady-state. We first show that our method yields the exact solution in the single-band model, and then apply it to the two-band case, where only numerical solutions are available. Finally, in Sec. IV.4, we discuss the long-time gap dynamics in real-time by performing an inverse Laplace transformation. We explicitly show how the new power-law decay exponent emerges from a distinct analytical structure of the gap in Laplace space and demonstrate how one re-obtains the single-band result. We conclude in Sec. V, and present additional details of our analytical calculations in the Appendices.

II BCS model and quench protocol

II.1 Pseudospin formalism for equilibrium two-band superconductors

We start from the reduced BCS Hamiltonian Bardeen et al. (1957) for two-band superconductors

[TABLE]

where $\alpha,\,\beta\in\left\{1,\,2\right\}$ are the band indices, $\varepsilon_{\mathbf{k},\alpha}$ is the electronic dispersion near the Fermi level in band $\alpha$ (including the chemical potential), and $V_{\alpha\beta}$ is the effective pairing interaction between band $\alpha$ and band $\beta$. Although not important in the following, one may assume parabolic dispersions, $\varepsilon_{\mathbf{k},\alpha}=\mathbf{k}^{2}/2m_{\alpha}-\mu$. The interaction constants $V_{\alpha\beta}$ are positive (negative) if the interaction is repulsive (attractive). In multi-band systems, different bands develop different values of the superconducting gap, depending on the values of the intra-band interactions, $V_{11}$ and $V_{22}$, and the inter-band interactions, $V_{12}$ and $V_{21}$ (see Fig. 2 (A)) as well as the density of states of the two bands at the Fermi level, $\mathcal{N}_{\alpha}$. We assume that the two bands have the same intra-band electronic interactions such that $V_{11}=V_{22}\equiv U$; by definition, $V_{12}=V_{21}\equiv V$. Due to the different density of states $\mathcal{N}_{1}\neq\mathcal{N}_{2}$, electrons in different bands experience different effective interaction strengths. The BCS gap equation is therefore band-dependent:

[TABLE]

Going from summation over momenta to integrations over energy using the density of states, we write the equilibrium BCS gap equations explicitly in matrix form in the band-space.

[TABLE]

where $\Lambda$ is a high-energy cutoff and

[TABLE]

with $\eta=\mathcal{N}_{2}/\mathcal{N}_{1}$ being the ratio of the density of states of the two bands, $E_{\alpha}=\sqrt{\varepsilon^{2}+\Delta_{\alpha}^{2}}$ is the Bogoliubov quasiparticle dispersion in band $\alpha$ and $T$ is the temperature of the system. In the following, we restrict our analysis to the $T=0$ ground state as the initial pre-quench state of the system. We have also defined the dimensionless inter-band interaction coupling constant $v=V\mathcal{N}_{1}$, and the dimensionless ratio $r=-U/V$ between intra-band and inter-band interactions. Here, we include the minus sign in the definition, as we will assume that $U<0$ is negative, corresponding to attractive intra-band interaction.

Note that the ratio of the density of states in the two bands, $\eta=\mathcal{N}_{2}/\mathcal{N}_{1}$, determines the relative sizes of the superconducting gaps of the two bands. If the two bands have the same density of states near the Fermi energy, i.e. $\eta=1$, the matrix $\hat{\gamma}$ becomes symmetric. Therefore, the gap equations are solved by $\Delta_{1}=-\Delta_{2}$ for repulsive inter-band interaction ($v>0$), corresponding to $s^{+-}$ pairing, and $\Delta_{1}=\Delta_{2}$ for attractive inter-band interaction ($v<0$), corresponding to $s^{++}$ pairing. In this paper, we will focus on the case with $\eta\neq 1$, in which case the amplitude of the two gaps is different in equilibrium $|\Delta_{1}|\neq|\Delta_{2}|$ and the multi-band nature of the system has a pronounced imprint on the non-equilibrium dynamics of the superconducting gap.

It is convenient to use the pseudospin formalism Anderson (1958b) to study the non-equilibrium dynamics of the superconducting state. In the mean-field approach, which is exact in the BCS regime we consider here, the BCS Hamiltonian can be described by pseudospins exposed to an effective magnetic field:

[TABLE]

where $\mathbf{B}_{\mathbf{k},\alpha}=2\left(\Delta_{\alpha}^{{}^{\prime}},-\Delta_{\alpha}^{{}^{\prime\prime}},-\varepsilon_{\mathbf{k},\alpha}\right)$ and

[TABLE]

The constant term contributes to the condensation energy, which will be ignored because it is not relevant to the dynamics out of equilibrium. The mapping between pseudo-spins and electronic pair operators is summarized in Fig. 2(B). The anti-commutation relation between the electronic operators ensures the spin commutation relation between $\mathbf{\hat{S}}_{\mathbf{k},\alpha}$. Notice that despite the simple form of the pseudospin Hamiltonian, the effective magnetic field is self-consistently determined by the pseudospins collectively via the gap equation:

[TABLE]

where $S_{\mathbf{k},\alpha}^{-}=\left\langle\hat{S}_{\mathbf{k},\alpha}^{-}\right\rangle=\left\langle c_{-\mathbf{k},\downarrow,\alpha}c_{\mathbf{k},\uparrow,\alpha}\right\rangle$.

In equilibrium, the pseudospins are parallel to the effective magnetic field. It is convenient to work in a gauge where both the gaps are real. Then the expectation values of the pseudo-spins at temperature $T$ are given by

[TABLE]

Note that the length of the pseudospins in equilibrium is determined by the Fermi-Dirac distribution, $n_{\text{F}}$, of the Bogoliubov quasiparticles, i.e. $\left|\mathbf{S}_{\mathbf{k},\alpha}\right|=\frac{1}{2}-n_{\text{F}}$. As mentioned above, we will focus hereafter on initial pre-quench states at zero temperature ($T=0$).

II.2 Equations of motion for the pseudospins

We consider the situation where the system is driven out of equilibrium by a sudden quench of the pairing interaction. Specifically, we focus on a sudden change of the inter-band coupling $v_{i}\rightarrow v_{f}$ while keeping the ratios between intra- and inter-band interactions, $r=U/V$, and between the densities of states, $\eta$, unchanged, i.e., $r_{i}=r_{f}$ and $\eta_{i}=\eta_{f}$. The subscript $i$ and $f$ denote the initial and final values of the respective dimensionless constants. Note that this requires quenching both intra- and inter-band interactions $U$ and $V$ in such a way to keep their ratio $r$ fixed. We focus on these quench protocols to constrain the parameter space. Generally, one can also consider quenches of $r$, however, this is expected to not lead to qualitative changes to the non-equilibrium dynamics, as it corresponds to a different way to prepare the initial conditions.

If the two bands have different densities of states, i.e. $\eta\neq 1$, the quench dynamics is intrinsically different from single-band systems. In the pseudospin formalism, the superconducting gap determines the intrinsic frequency of the pseudospin precession. Therefore, once the two bands have different densities of states, they develop different values of the gap, leading to two distinct intrinsic frequencies. In addition, the gap also serves as the effective magnetic field that drives the precession. Through the inter-band interaction, each band experiences an oscillating magnetic field with the intrinsic frequency of the other band. Hence, the dephasing of the pseudospin oscillations in multi-band systems is fundamentally different from single-band systems. The dynamics is described by two sets of equations of motion for the two bands, which are derived from Eq. (5) in terms of expectation values of the pseudospins operators,

[TABLE]

which are similar to the one-band case, but now with an extra band index $\alpha$. More importantly, the pseudospin dynamics in the two bands are coupled via the gap equations with a time-dependent inter-band coupling strength $v\left(t\right)=v_{i}\theta\left(-t\right)+v_{f}\theta\left(t\right)$:

[TABLE]

The equations of motion for the pseudospins, Eq. (11) and the time-dependent gap equation, Eq. (12), determine the post-quench gap dynamics of two-band superconductors.

In the following, we first solve these equations numerically and describe our results. Then, we analytically find the long-term asymptotic behavior of the gap oscillations using Laplace transforms. We develop a generalization of the well-known procedure pioneered by Volkov and Kogan in Ref. Volkov and Kogan, 1974 (see also Refs. Yuzbashyan and Dzero (2006); Yuzbashyan et al. (2015)). By expanding around the long-time non-equilibrium pseudo-spin steady state, instead of the final equilibrium state, we are able to not only determine the power-law decay of the gap oscillations, but also the steady-state non-equilibrium gap values $\Delta_{\alpha,\infty}$. We also explicitly show how our solution approaches the known single-band result as $\eta\rightarrow 1$.

III Numerical results for the post-quench gap dynamics

We solve the equations of motion (11), together with the gap equation (12), numerically using the Runge-Kutta method. We focus on the weak-quench limit, to later compare with our analytical expansion. Results for two different ratios of initial and final inter-gap couplings $v_{i}/v_{f}=0.95$ and $0.9$ (with fixed $v_{f}=0.2$) are shown in Fig. 3. The other parameters are kept fixed: $r_{i}=r_{f}=0$, $\eta=0.8$, $T_{i}=0$. In equilibrium, this corresponds to the following gap ratios $\Delta_{1,i}/\Delta_{2,i}=-0.8852$ for $v_{i}=0.19$ , $\Delta_{1,i}/\Delta_{2,i}=-0.8857$ for $v_{i}=0.18$ and $\Delta_{1,f}/\Delta_{2,f}=-0.8847$ for $v_{f}=0.2$. The figure contains both the time traces of the gap oscillations as well as their Fourier transforms.

There are two important qualitative features that emerge in the two-band case: first, the gap oscillations are characterized by two frequencies, corresponding to the steady-state values $\Delta_{1,\infty}$ and $\Delta_{2,\infty}$. This leads to pronounced beating when these two frequencies are sufficiently close to each other. This phenomenon has been described previously in numerical studies of two-band (multi-band) superconductors exposed to terahertz laser pulses Akbari et al. (2013); Krull et al. (2016); Zachmann et al. (2013) Second, the algebraic decay of the gap oscillations ($\propto t^{-\alpha}$) occurs more rapidly than in the single-band case. We numerically determine the exponent to be $\alpha_{\text{2-band}}=3/2$ as opposed to $\alpha_{\text{1-band}}=1/2$.

This behavior seems insensitive to the actual value of $r$. In Fig. 4, we compare the behavior of $\Delta_{1}(t)$ for the cases in which $r=0.5$ and $r=0$. The other parameters used were $v_{i}=0.19$ and $v_{f}=0.2$. We note that an exponent of $\alpha=3/2$ also emerges if one considers deep quenches into the Bose-Einstein condensate (BEC) regime in a three-dimensional system Gurarie (2009); Yuzbashyan et al. (2015).

IV Long-time asymptotic gap dynamics

In order to gain more insights on the transient dynamics of the superconducting gap in two-band systems, it is instructive to have analytic solutions for the superconducting gap evolution. The gap dynamics in single-band conventional superconductors with isotropic gap structures can be solved exactly due to the integrability of the BCS model Richardson and Sherman (1964); Dukelsky et al. (2004); Yuzbashyan et al. (2005b, c, a, 2006); Barankov and Levitov (2006); Yuzbashyan et al. (2015). The two-band BCS model doubles the number of degrees of freedom compared to the single-band model. Due to the coupling between the two distinct bands, the integrals of motion that were constructed previously for the single-band BCS model Yuzbashyan et al. (2005c, 2015) do not commute between the two bands, except in the symmetric case $\eta=1$. In the single-band case, it was determined that there are three different “phases” depending on the strength of the quench $\Delta_{i}/\Delta_{f}$: in phase I, corresponding to $\Delta_{i}/\Delta_{f}>\mathrm{e}^{\pi/2}$, the gap asymptotically approaches zero in an exponential fashion; in phase II, for $\mathrm{e}^{-\pi/2}<\Delta_{i}/\Delta_{f}<\mathrm{e}^{\pi/2}$, the gap shows damped $t^{-1/2}$ oscillations around one asymptotic value; and in phase III, which takes place for $\Delta_{i}/\Delta_{f}<\mathrm{e}^{-\pi/2}$, the gap shows persistent oscillations between two asymptotic values.

Whether the two-band BCS model is integrable or not is beyond the scope of this work. Given the difficulties in finding the integrals of motion of the two-band case, in this section we employ instead a perturbative method to extract the long-time asymptotic dynamics of the superconducting gap in phase II, where the gap shows damped oscillations. This is precisely the behavior found numerically for weak quenches, shown in Fig. 3. In particular, the method we develop here is a modified version of the one pioneered by Volkov and Kogan in Ref. Volkov and Kogan, 1974, which allows us to also analytically determine the steady-state gap values $\Delta_{\alpha,\infty}$.

For convenience, we briefly review our notation scheme: subscripts $i$ and $f$ denote the thermal equilibrium value before ($i$) and after ($f$) the quench. The subscript $\infty$ denotes the long-time asymptotic steady-state value of the gap. For example, $\Delta_{\alpha,i}$ ($\Delta_{\alpha,f}$) is the equilibrium value of gap $\alpha$ before (after) the quench, and $\Delta_{\alpha,\infty}$ is its long-time asymptotic steady-state value following the time evolution governed by the BCS Hamiltonian. We note that our analysis is restricted to weak quenches, resulting in the system being in phase II, where the gap experiences Volkov-Kogan-like behavior.

IV.1 Linearized equations of motion

To analytically describe the post-quench gap dynamics at long times, we generalize the method used first by Volkov and Kogan in Ref. Volkov and Kogan, 1974. Instead of expanding around the final equilibrium state $S_{\alpha,f}^{i}$ and $\Delta_{\alpha,f}$, however, we expand around the long-time non-equilibrium steady-state values $S_{\alpha,\infty}^{i}$ and $\Delta_{\alpha,\infty}$. Importantly, these steady-state values will be determined self-consistently in our calculation using Laplace’s final value theorem. We thus assume that in the long-time limit the superconducting gaps reach their long-time asymptotic values $\Delta_{\alpha,\infty}$. Specifically, we expand the equations of motion and the gap equations around the asymptotic steady-state values

[TABLE]

where, from the stationary condition of the equations of motion, $S_{\alpha,\infty}^{\pm}=S_{\alpha,\infty}^{x}$, $S_{\alpha,\infty}^{y}=0$, $\varepsilon S_{\alpha,\infty}^{x}=-\Delta_{\alpha,\infty}S_{\alpha,\infty}^{z}$. Note that $f_{\alpha}$ describes pairing amplitude fluctuations and $g_{\alpha}$ describes density fluctuations. The deviation of the gap from its long-time asymptotic value is denoted by $\delta_{\alpha}$, which is determined by the pairing-amplitude fluctuations $f_{\alpha}$ via the gap equation:

[TABLE]

where $\gamma_{\alpha\beta}$ is given in Eq. 4. As we will show below, because $f^{\prime\prime}_{\alpha}$ is an odd function of $\varepsilon$, $\delta_{\alpha}$ is real, as long as we choose the initial equilibrium gaps of the two bands $\Delta_{\alpha,i}$ to be real. With this in mind, we linearize the equations of motion by inserting Eqs. (13) into Eq. (11) to obtain

[TABLE]

where $f_{\alpha}=f_{\alpha}^{\prime}+if_{\alpha}^{\prime\prime}$ and the notation $\dot{f}\equiv\frac{df}{dt}$ is used. Note that, as anticipated, $f_{\alpha}^{\prime\prime}$ remains an odd function of $\varepsilon$ for all times, since $S_{z,\infty}$ and $g_{\alpha}$ are odd while $f_{\alpha}^{\prime}$ is even. As a result, the gap remains real for all times. The fact that the phases of the gaps are constants of motion follows directly from the particle-hole symmetry of the BCS Hamiltonian Barankov and Levitov (2006). Therefore, the relative phase of the two gaps is also a constant of motion and the Leggett (relative phase) mode, which would in any case be overdamped in the regime we study here of inter-band pairing interaction only, is not excited in our quench protocol. In order to excite it, one must break the particle-hole symmetry of the BCS Hamiltonian, for example, by external perturbations as in the pump-probe setups Krull et al. (2016).

The linearized equations of motion faithfully describe the long-time dynamics since at the long-time limit, the deviations from the asymptotic values are small, i.e. $\left(g_{\alpha},f_{\alpha},\delta_{\alpha}\right)\ll\left(S_{\alpha,\infty}^{z},S_{\alpha,\infty}^{-},\Delta_{\alpha,\infty}\right)$. To have a better description of the gap dynamics over a wider time range, we focus on relatively weak quenches where $v_{f}/v_{i}$ is close to 1. In this case, the oscillations around $\Delta_{\alpha,\infty}$ are small already at earlier times, allowing for a better comparison between numerics and analytics. Such weak quench regime is also the most relevant to experiments, where excess heating is suppressed.

Since we are interested in $\delta_{\alpha}$, which is only related to $f_{\alpha}$, see Eq. 14, we can further simplify the above equations by eliminating $g_{\alpha}$ to find

[TABLE]

where $E_{\alpha,\infty}^{2}=\varepsilon^{2}+\Delta_{\alpha,\infty}^{2}$. Eqs. 16a and 16b describe the dynamics of the imaginary and real parts of the pairing amplitude fluctuations, respectively, which determine the time evolution of the gap.

IV.2 Solution in Laplace space

To solve the differential equations (16a) and (16b), it is useful to perform a Laplace transformation $y\left(s\right)=\int_{0}^{\infty}y\left(t\right)e^{-st}dt$. We find the the following algebraic equations:

[TABLE]

Here, $s$ is the complex frequency in the Laplace domain and the subscript [math] indicates an initial condition, i.e. $f_{\alpha,0}\equiv f_{\alpha}\left(\varepsilon,t=0^{+}\right)$, $\delta_{\alpha,0}\equiv\delta_{\alpha}\left(t=0^{+}\right)$, etc. Physically, Eqs. (17a) and (17b) describe the phase and amplitude dynamics of the gap, respectively.

Since $\delta_{\alpha}$ and $f_{\alpha}$ are related through the gap equation 14, it is convenient to integrate both sides of the above equations over $\varepsilon$. Then, Eq. (17a) is trivially satisfied, since $S_{\alpha,\infty}^{z}$ is an odd function of $\varepsilon$, by virtue of Eq. 10c, $f_{\alpha,0}^{{}^{\prime\prime}}=0$ by construction, and $\dot{f}_{\alpha,0}^{\prime\prime}$ is an odd function of $\varepsilon$, by virtue of the second equation of 15c. Consequently, we are left with a single equation for $f_{\alpha}^{\prime}$ and $\delta_{\alpha}$, which are related through the gap equation 14.

Expressing $f$ in terms of $\delta$, and recasting Eq. (17b) in matrix form, the deviations of the superconducting gaps from their asymptotic values, $\delta_{\alpha}$, are given by:

[TABLE]

where the hat (arrow) denote a matrix (vector) in band space. Here, we defined:

[TABLE]

where $\mathbb{I}$ is the identity matrix in band space and the following notation is used:

[TABLE]

For convenience, we write $\hat{\Phi}_{\alpha\beta}^{\infty}\left(s\right)\equiv\mathbb{I}_{\alpha\beta}\Phi_{\alpha}^{\infty}$ and define:

[TABLE]

The function $\vec{I}\left(s\right)$ on the right-hand side is given by (detailed derivation in Appendix A)

[TABLE]

with:

[TABLE]

The solution for $\vec{\delta}(s)$ in Laplace space is then simply given by

[TABLE]

It is clear that without inter-band interaction, $V=0$, $\hat{\mathcal{M}}_{\alpha\beta}$ becomes a diagonal matrix, since $\hat{\gamma}_{\alpha\beta}$ in Eq. 4 is diagonal. As a result, Eq. 18 becomes diagonal in band space as well, and the two-band model reduces to two independent one-band models.

In the following subsections, we will extract the dynamics of the gaps in the long-time limit from their analytic behaviors in Laplace space. These are determined by the functions $\hat{\Phi}^{\infty}\left(s\right)$ and $\vec{I}\left(s\right)$, as they are the only $s$-dependent functions in Eq. (18). Their $s$-dependence comes from the two functions $\Phi_{\alpha}^{\infty}\left(s\right)$ and $\Phi_{\alpha}^{i}\left(s\right)$ defined above.

The function $\Phi_{\alpha}^{i}\left(s\right)$ is straightforward to calculate since the initial pseudospin configuration is given by the equilibrium value of the gap at $T=0$, i.e. $S_{\alpha,i}^{x}/\Delta_{\alpha,i}=\frac{1}{2\sqrt{\varepsilon^{2}+\Delta_{\alpha,i}^{2}}}$. Inserting this initial pseudo-spin state into Eq. (24), this can be brought to the form

[TABLE]

where we defined the dimensionless ratio $\tilde{\Delta}_{\alpha,i}=\Delta_{\alpha,i}/\Delta_{\alpha,\infty}$ and the function

[TABLE]

To find an explicit expression for $\Phi_{\alpha}^{\infty}\left(s\right)$, given by Eq. 22, we first need to compute the function $S_{\alpha,\infty}^{x}/\Delta_{\alpha,\infty}$. The gap equation (see Eq. (9)), which is satisfied regardless of whether the system is in thermal equilibrium or not, restricts the expectation value of this quantity to:

[TABLE]

As we discussed above, the non-zero inter-band interactions render the matrix $\hat{\mathcal{M}}$ off-diagonal and make the two-band model fundamentally different than the single-band case. While a generic discussion of arbitrary inter- and intra-band interactions is possible, the analysis is simplified considerably by focusing on the case of inter-band repulsion only, i.e. $r=0$. Indeed, our numerical results discussed in Fig. 4 show that the general behavior of the two-band problem is the same for $r=0$ and $r\neq 0$. Setting $r=0$ in Eq. (4) yields an off-diagonal matrix \hat{\gamma}=\left(\begin{array}[]{cc}0&-\eta\\ -1&0\end{array}\right). As result, the equation above becomes:

[TABLE]

Note that this ratio involves the pseudospin of band $\alpha$ and the gap of the other band $\bar{\alpha}$, where $\bar{\alpha}=1(2)$ for $\alpha=2(1)$. To proceed, we note that, in equilibrium, the same relationship holds between the ratios of the pseudospin and the gap:

[TABLE]

The difference is that, in equilibrium, from Eq. 10c, we know precisely the expression for $S_{\alpha,f}^{x}$:

[TABLE]

Based on this similarity, we propose the following ansatz:

[TABLE]

where $\tilde{\Delta}_{\alpha,f}=\Delta_{\alpha,f}/\Delta_{\alpha,\infty}$ is defined analogously to $\tilde{\Delta}_{\alpha,i}$. Clearly, this ansatz satisfies the constraint 29. For $r\neq 0$, the constraint will likely have a more complicated form; thus, for the sake of clarity, we focus on the case $r=0.$ We will verify the validity of this ansatz later by an explicit comparison to numerical calculations and by comparison with the exact solution of the single-band case. For now, we proceed with this ansatz and perform the energy integration in the expression of $\Phi_{\alpha}^{\infty}\left(s\right)$. We obtain:

[TABLE]

IV.3 Asymptotic gap values

In this subsection, we show how to extract the long-time asymptotic steady-state gap values $\Delta_{\alpha,\infty}$ self-consistently. To set the stage, and validate the ansatz proposed in the previous subsection, we first present the calculation for the single-band case, comparing the perturbative solution with the exact one.

IV.3.1 Asymptotic gap for the single-band model

In the single-band BCS model with attractive pairing interaction $u\equiv U\mathcal{N}$, a quench suddenly changes the pairing interaction $u_{i}\rightarrow u_{f}$. It is convenient to use $\Delta_{i}/\Delta_{f}$ as the quench parameter, where $\Delta_{i}$ ($\Delta_{f}$) is the equilibrium value of the gap with pairing interaction $u_{i}$ ($u_{f}$). We employ the same linearization scheme for the single-band model as above in Eqs. (13a)-(15c) for the two-band case, and expand around the long-time asymptotic values, $S_{\infty}^{\alpha}$ and $\Delta_{\infty}$. The equation for the gap deviation $\delta$ in Laplace space, Eq. 18, becomes in the single-band case:

[TABLE]

where

[TABLE]

Here, $S_{i}^{x}/\Delta_{i}=1/(2E_{i})$ with $E_{i}=\sqrt{\varepsilon^{2}+\Delta_{i}^{2}}$ is given by its value in the initial $T=0$ ground state state prior to the quench. The ratio $S_{\infty}^{x}/\Delta_{\infty}$, according to our ansatz (32), becomes in the single-band case:

[TABLE]

This ansatz can be recast in an alternative way as an ansatz for the non-equilibrium distribution function. From the definition of $S_{f}^{x}$, Eq. 10c, we have:

[TABLE]

where we defined the equilibrium distribution function $n_{0}(\varepsilon)=\tanh(\sqrt{\varepsilon^{2}+\Delta_{f}^{2}}/\left(2T\right)).$ From the gap equation, it follows that $\left\langle\frac{S_{f}^{x}}{\Delta_{f}}\right\rangle=1$. Analogously, we can express $S_{\infty}^{x}$ in terms of the non-equilibrium quasiparticle distribution function $n_{\text{eff}}\left(\varepsilon\right)$:

[TABLE]

Because the gap equation has to be satisfied also in non-equilibrium, it follows that:

[TABLE]

The ansatz 36 thus can be recast as an ansatz for the non-equilibrium distribution function:

[TABLE]

Having obtained an explicit expression for $S_{\infty}^{x}/\Delta_{\infty}$, we can derive analytic expressions for $\Phi_{i}\left(s\right)$ and $\Phi_{\infty}\left(s\right)$:

[TABLE]

To find the long-time asymptotic value of the gap, we use the self-consistency condition that $\lim_{t\rightarrow\infty}\Delta\left(t\right)=\Delta_{\infty}$, or equivalently, $\lim_{t\rightarrow\infty}\delta\left(t\right)=0$. Using the final value theorem in Laplace space, this condition becomes

[TABLE]

Using Eq. (34) and inserting the explicit expressions from Eq. (41), we find that the asymptotic value of the gap $\Delta_{\infty}$ must satisfy

[TABLE]

It is straightforward to show that this equation is identical to the one that emerges in the exact solution of the single-band BCS gap dynamics using the method of the Lax vector Barankov and Levitov (2006); Yuzbashyan et al. (2006). In Fig. 5, we compare the results from both methods, which match perfectly. Interestingly, in phase III (persistent oscillations), our method gives the average value of the gap. Of course, our method formally breaks down in this phase, because the Laplace final value theorem ceases to hold for an oscillatory long-time solution.

This comparison validates the ansatz 32 for the single-band case, giving us confidence to apply it to the two-band case as well. Note that the perfect agreement with the exact solution does not necessarily imply that the non-equilibrium distribution function 40 is also exact .

IV.3.2 Asymptotic gap for the two-band model

We now perform the same calculation for the two-band model with pure inter-band repulsion ($r=0$). Using Eqs.(26) and (33), we obtain the following expression for $\delta_{\alpha}\left(s\right)$ from Eq. 18:

[TABLE]

where, for convenience of notation, we introduced $\eta_{1}=1$ and $\eta_{2}\equiv\eta$, $I_{\alpha}\left(s\right)$ is given by Eq. (23), and:

[TABLE]

To find the asymptotic long-time value of the gaps $\Delta_{\alpha,\infty}$, we employ once again the final value theorem in Laplace space, Eq. 42. We numerically solve for $\Delta_{1,\infty}$ and $\Delta_{2,\infty}$ for a given quench protocol, $v_{i}\rightarrow v_{f}$, or equivalently $\Delta_{1,i}\rightarrow\Delta_{1,f}$. As shown in Fig. 6, we find that, in the case of pure inter-band interactions ($r=0$), the ratios between the asymptotic and final equilibrium gaps $\Delta_{\alpha,\infty}/\Delta_{\alpha,f}$ are, to a very good approximation (i.e. with a numerical deviation of less than $0.01\%$), equal for both bands, i.e. $\tilde{\Delta}_{1,f}=\tilde{\Delta}_{2,f}$. They are also identical to the single-band ratio if we adjust the definition of the quench amplitude accordingly, such that the $x$-axis corresponds to $\Delta_{i}/\Delta_{f}$ in the single-band case and to $\Delta_{1,i}/\Delta_{1,f}$ in the two-band case.

Using the result obtained here that $\tilde{\Delta}_{1,f}=\tilde{\Delta}_{2,f}$, the pre-factor of Eq. 33 becomes 1. Thus, both $\Phi_{\alpha}^{\infty}\left(s\right)$ and $\Phi_{\alpha}^{i}\left(s\right)$ have the same functional dependence: $\Phi_{\alpha}^{\infty}\left(s\right)=\Upsilon\left(\tilde{\Delta}_{\alpha,f},\frac{s}{2\Delta_{\alpha,\infty}}\right)$, $\Phi_{\alpha}^{i}\left(s\right)=\Upsilon\left(\tilde{\Delta}_{\alpha,i},\frac{s}{2\Delta_{\alpha,\infty}}\right)$.

IV.4 Damped gap oscillations in the long-time limit

The long-time behavior of the gap in the time domain $\Delta(t)$ can be obtained by applying the inverse Laplace transformation to Eq. (44). In order to perform the inverse Laplace transformation, we first need to study the analytical behavior of the solution in Laplace space and find its poles and branch cuts. They are determined by the analytic properties of the function $\Upsilon(\Delta,x)$, defined in Eq. (27) and repeated here for convenience:

[TABLE]

The reason why only the analytical properties of $\Upsilon(\Delta,x)$ matter is because we can express both $\Phi_{\alpha}^{i}(s)$ and $\Phi_{\alpha}^{\infty}(s)$ in terms of this function:

[TABLE]

where $z=\frac{s}{2\Delta_{1,\infty}}$, $\tilde{\Delta}_{\alpha,i/f}=\frac{\Delta_{\alpha,i/f}}{\Delta_{\alpha,\infty}}$, and $\kappa=\frac{\Delta_{1,\infty}}{\Delta_{2,\infty}}$. For concreteness, in this section we consider the gap with $\alpha=1$ to be the one that is asymptotically smaller, implying that $\left|\kappa\right|<1$. But note that our results can be straightforwardly applied also to the case $\left|\kappa\right|>1$.

The function $\Upsilon\left(\Delta,\,z\right)$ has two branch cuts, one between $\left(-i\infty,\ -i\right)$ and another one between $\left(i,\ i\infty\right)$. The function is analytic elsewhere. Applying the Cauchy’s residue theorem (see Appendix C and Fig. 9 for details), we convert the Bromwich integral into four integrals along the sides of the two branch cuts. Note that we have already eliminated the pole at the origin by imposing the final value theorem in Section IV.3.1. In addition, we also use the following properties of the function $\Upsilon$:

[TABLE]

As a result, the inverse Laplace transformation is given by the following integral:

[TABLE]

where $z\delta_{\alpha}\left(z\right)$ is given by

[TABLE]

with

[TABLE]

In the long-time limit, where $2\Delta_{1,\infty}t\gg 1$, the integrand of Eq. (49) is highly oscillatory. Only singular behaviors of $\text{Im}\left[z\delta_{\alpha}\left(z\right)\right]$ will therefore make a contribution to the long-time dynamics of the superconducting gap. Indeed, $\text{Im}\left[z\delta_{\alpha}\left(z\right)\right]$ has two branch points along $z\in\left[i,\,i\infty\right)$: one is located at $z=i$ and the other one is located at $z=i/\left|\kappa\right|$. We expand $\text{Im}\left[z\delta_{\alpha}\left(z\right)\right]$ near these two branch points, i.e. $z=i+i\epsilon$ and $z=i/\left|\kappa\right|\pm i\epsilon$, and find that both exhibit $\sqrt{\epsilon}$ behavior (details shown in Appendix (B)). This is sharply distinct from the single-band case, where only one branch point is present along $z\in\left[i,\,i\infty\right)$. More importantly, the asymptotic behavior in the vicinity of the branch point in the single-band case is $1/\sqrt{\epsilon}$ rather than $\sqrt{\epsilon}$. The two cases are plotted and compared in Fig. 7. The $1/\sqrt{\epsilon}$ behavior leads to a $t^{-1/2}$ decay of the gap oscillation amplitude at long times in the single-band case Volkov and Kogan (1974). In contrast, the $\sqrt{\epsilon}$ behavior in Laplace space leads to a faster $t^{-3/2}$ decay in the two-band model

[TABLE]

for $2\Delta t\gg 1$ (details are shown in Appendix C). The damping of the gap oscillations thus occurs faster for two-band superconductivity.

To find the full long-time expressions of the gap, including prefactors and oscillatory factors, we perform a careful asymptotic analysis of $\text{Im}\left[z\delta_{\alpha}\left(z\right)\right]$. The final result for the long-time gap oscillations reads

[TABLE]

where the pre-factors $\mathcal{A}_{\alpha}$, $\mathcal{B}_{\alpha}$ and $\mathcal{C}_{\alpha}$ are calculated from the asymptotic analysis and explicitly shown in Appendix B. The gap oscillation frequencies are determined by the asymptotic values of the gaps in the two different bands $\Delta_{\alpha,\infty}$. As discussed in the previous sections, the asymptotic values of the gaps are determined by the quench amplitude $\Delta_{\alpha,i}/\Delta_{\alpha,f}$ and the ratio of the density of states $\eta$ between the two bands. In general, they will also depend on $r=-U/V$, which we have set to zero for simplicity here. The same holds for the prefactors of the sinusoidal oscillations.

In Fig. 8, we compare our analytical results to the numerical solution of the equations of motion for two different weak quench amplitudes in phase II. We find an excellent quantitative agreement between the two, which also justifies our analytical ansatz a posteriori.

We finish this section by commenting on how our solution gives the known single-band result in the limit where the ratio between the two densities of states approaches one, $\eta\rightarrow 1$. In this limit, the gaps have the same asymptotic magnitude, i.e. $\Delta_{1,\infty}=\left|\Delta_{2,\infty}\right|$. The equilibrium gaps also have the same magnitude, leading to $\Upsilon\left(\tilde{\Delta}_{\alpha,f},z\right)=\Upsilon\left(\tilde{\Delta}_{\bar{\alpha},f},z\right)=\Upsilon\left(\tilde{\Delta}_{1,f},z\right)$. As a result, Eq. 50 becomes

[TABLE]

where $D\left(z\right)=\Upsilon^{2}(\tilde{\Delta}_{1,f},z)-2\Upsilon(\tilde{\Delta}_{1,f},z)$. Further simplification of the above equation gives:

[TABLE]

In writing this last equation, we used the fact that $\Upsilon\left(\tilde{\Delta}_{\alpha,f},z\right)=\Upsilon\left(\tilde{\Delta}_{1,f},z\right)$. This is the same expression as the solution of the single-band case in Laplace-space, Eq. 34. Using the asymptotic behavior of $\Upsilon\left(\tilde{\Delta}_{\alpha},iy\right)$ near the branch point $y\rightarrow 1$ (details shown in Appendix B, see Eq. 62), we arrive at the following asymptotic behavior:

[TABLE]

By applying the inverse Laplace transformation, we find that the gap dynamics is characterized by oscillations with frequency $2\Delta_{\infty}$ and $t^{-1/2}$ damping:

[TABLE]

V Conclusions

In this paper, we developed a generalization of the Volkov-Kogan Laplace-space analysis for the post-quench dynamics of $s$-wave BCS superconductors in the collisionless regime Volkov and Kogan (1974), and applied it to interaction quenches of two-band BCS superconductors. We showed that the two-band case is fundamentally different from the single-band case. Not only do the gap oscillations display beating associated with the two different gap values on the two bands, but they also display a faster $t^{-3/2}$ power-law damping, as opposed to the $t^{-1/2}$ damping of the single-band case. For weak quenches, our analytical results agree very well with the numerical results in the long-time limit, demonstrating that the gap dynamics of multi-band systems cannot be simply decomposed into the sum of the gap dynamics of single-band systems. Formally, this new power-law decay can be understood as arising from the “splitting” of the relevant branch point in Laplace space in two, as shown in Fig. (7). As a result, one expects the same $t^{-3/2}$ behavior to take place even when the number of bands is larger than $2$. From a more physical perspective, the stronger damping in the two-band case arises because the Cooper-pairs dephasing involves states from both bands due to the inter-band coupling. Such a dephasing is thusintrinsic to multi-band systems and independent on the quench amplitude.

From a methodological viewpoint, our analysis is distinct from the one introduced by Volkov and Kogan Volkov and Kogan (1974) (see also the more recent works by Yuzbashyan and co-workers Yuzbashyan and Dzero (2006); Yuzbashyan et al. (2015)), because we linearize the equations of motion around the asymptotic long-time pseudo-spin states as opposed to the final equilibrium states. This allows us to self-consistently determine the asymptotic long-time steady-state values of the gaps over the full range of quench amplitudes in phase II (and phase I, where the steady-state gaps vanish). We explicitly showed that the self-consistent equation for the steady-state gap in the single-band case agrees with the exact expression derived within the Lax vector analysis Yuzbashyan et al. (2005a); Barankov and Levitov (2006); Yuzbashyan et al. (2006). Like in the two-band model we investigate here, our method can be very useful in cases where an exact solution is not (yet) available, for example, to investigate quenches towards more exotic fully gapped pairing states such as $s+is$ or $s+id$. Other interesting future directions are to include a finite intra-band pairing interaction $r\neq 0$, competing electronic order parameters such as spin-density waves Schütt et al. (2018), or generalize and apply our Laplace method to study quenches in superconductors with a nodal gap structure such as those with $d$-wave symmetryPeronaci et al. (2015).

Acknowledgements.

This work was supported by the by U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award DE-SC0012336. The authors acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota, where the numerical calculations were performed. P.P.O. acknowledges support from Iowa State University Startup Funds. T.C. also acknowledges the support from the Doctoral Dissertation Fellowship awarded by the University of Minnesota.

Appendix A Initial conditions for the interaction quench

The system is at equilibrium before the interaction quench. For systems with only inter-band repulsion, the superconducting gap is given by

[TABLE]

where $v_{i}=V_{i}\mathcal{N}_{1}$ is the dimensionless inter-band repulsion, and $\eta=\frac{\mathcal{N}_{2}}{\mathcal{N}_{1}}$ is the ratio between the density of states near the Fermi level of the two bands. The pseudospins are

[TABLE]

After the interaction quench, the inter-band repulsion is suddenly changed to a different value, $v_{f}$. The initial conditions of the post-quench dynamics of the superconducting gaps are thus given by replacing the inter-band repulsion with its post-quench value $v_{f}$.

[TABLE]

Substituting in the linearized equations 13a and 15c, the initial conditions on the pseudospin deviations $f_{\alpha}$ become

[TABLE]

We recall that $f_{\alpha,0}^{\prime\prime}$ and $\dot{f}_{\alpha,0}^{\prime\prime}$ are related to the dynamics of the superconducting gap in Laplace space via $I_{\alpha}\left(s\right)=\left\langle\frac{2\varepsilon\left[sf_{\alpha,0}^{\prime}+\dot{f}_{\alpha,0}^{\prime\prime}\right]}{s^{2}+4E_{\alpha,\infty}^{2}}\right\rangle$, which yields Eq. 23.

Appendix B Asymptotic analysis of the superconducting gap in Laplace space

In this appendix, we analyze the asymptotic behavior of the gap in Laplace space near the branch points. From Eq. 50, there are 7 terms that determine the analytic behavior of the gap. The branch points all come from the function $\Upsilon\left(\Delta,z\right)$, which opens branch cuts at at $\left(-i\infty,\ -i\right)$ and $\left(i,\ i\infty\right)$, as shown in Fig. 9. Let $z=iy$, then, around $y=1$, we have

[TABLE]

where $\epsilon$ is an infinitesimal positive number.

We use the asymptotic behavior of $\Upsilon\left(\Delta,y\right)$ to expand all the terms in Eq. 50, and obtain the following results:

[TABLE]

Appendix C Inverse Laplace transformation and useful integrals

The inverse Laplace transformation is given by the Bromwich integral:

[TABLE]

where $\sigma$ is a real number that is larger than the real parts of all the singularities of $Y\left(s\right)$.

All the asymptotic behaviors of the gap in Laplace space are square root like. Consequently, transforming back to real time domain leads to a $t^{-3/2}$ decay.

[TABLE]

for $2\Delta_{1,\infty}t\gg 1$. Similarly, $\int_{\frac{1}{\left|\kappa\right|}}^{\infty}\frac{\sqrt{\left|\kappa\right|y-1}}{y}\cos\left(2\Delta_{1,\infty}yt\right)dy\simeq-\frac{\sqrt{\pi}\sin\left(2\left|\Delta_{2,\infty}\right|t+\frac{\pi}{4}\right)}{2\left(2\left|\Delta_{2,\infty}\right|t\right)^{3/2}}$, and

[TABLE]

Appendix D Analytical expressions for the gap dynamics

We wrote the long-time asymptotic expressions of the gap oscillations in Eqs. 53a and 53b. In this appendix, we provide the explicit expressions for the prefactors $\mathcal{A}_{\alpha}$, $\mathcal{B}_{\alpha}$ and $\mathcal{C}_{\alpha}$ that appear in the two equations.

[TABLE]

[TABLE]