# Theory of nonequilibrium dynamics of multiband superconductors

## Abstract

We study the nonequilibrium dynamics of multiband BCS superconductors subjected to ultrashort pump pulses. Using density-matrix theory, the time evolution of the Bogoliubov quasiparticle densities and the superconducting order parameters are computed as a function of pump pulse frequency, duration, and intensity. Focusing on two-band superconductors, we consider two different model systems. The first one, relevant for iron-based superconductors, describes two-band superconductors with a repulsive interband interaction which is much larger than the intraband pairing terms. The second model, relevant for MgB, deals with the opposite limit where the intraband interactions are dominant and the interband pair scattering is weak but attractive. For ultrashort pump pulses, both of these models exhibit a nonadiabatic behavior which is characterized by oscillations of the superconducting order parameters. We find that for nonvanishing , the superconducting gap on each band exhibits two oscillatory frequencies which are determined by the long-time asymptotic values of the gaps. The relative strength of these two frequency components depends sensitively on the magnitude of the interband interaction .

###### pacs:

74.40.Gh###### pacs:

78.47.J-###### pacs:

78.20.BhNonequilibrium processes in superconductivity Femtosecond techniques; in spectroscopy of solid state dynamics Optical properties; theory and models

The discovery of iron-based superconductors[1] and MgB[2], where superconductivity is characterized by more than one order parameter, has lead to a renewed interest in multiband superconductors. The existence of multigap superconductivity was predicted more than fifty years ago in the pioneering works by Moskalenko[3] and Suhl et al.[4], where it was argued that a band-dependent interaction potential can give rise to different order parameters on different Fermi surface sheets. The presence of multiple superconducting order parameters has many interesting physical consequences, most of which are related to the relative amplitudes and phase differences between the gaps on different Fermi surfaces.

In recent years, numerous studies of the nonequilibrium dynamics of multiband superconductors have been performed using femtosecond time-resolved spectroscopy[5, 6, 7, 8, 9, 10]. The relaxation kinetics measured in these experiments gives important information on the electronic band structure, on the electron-phonon coupling strengths, as well as on the symmetry of the superconducting order parameters. For example, time-resolved measurements on the iron pnictide superconductor BaKFeAs[10] have revealed two distinct quasiparticle relaxation processes: a fast one, whose decay rate increases linearly with excitation density, and a slow one, which is independent of excitation density. This behavior has been attributed to the multigap structure of the superconductor. Furthermore, from a careful analysis of the temperature dependence of the decay rates, the authors of Ref. [10] have been able to deduce the pairing symmetries of the superconducting state.

In this paper, we shall not be concerned with the relaxation dynamics, but rather with the nonequilibrium evolution of multiband superconductors at times shorter than the relaxation time, i.e., with the nonadiabatic behavior of these system. It has been shown that the nonequilibrium response at these short time scales exhibits interesting properties, such as order parameter oscillations[11, 12, 13, 14]. In a pioneering work, B. Mansart et al.[15] have been able to detect these coherent oscillations of the order parameter in a high-T cuprate superconductor via transient broad-band reflectivity measurements. The recently discovered multiband superconductors MgB and iron pnictide superconductors, which all have a relatively large superconducting gap, offer an excellent opportunity to analyze this nonadiabatic dynamics also in multiband systems, where the interplay of the intraband pairing interactions and interband pair scattering terms play an important role. One goal of the present work, is to provide guidance for experimental studies of these interesting nonadiabatic response properties in multigap superconductors.

In the folowing, we employ the density-matrix theory to study the temporal evolution of Bogoliubov quasiparticle distributions in multiband superconductors after excitation by a short pump pulse. We focus on the nonadiabtic regime, i.e., the regime where the pump pulse duration is much shorter than the dynamical time scales of the superconducting order parameters. In this case, the pump pulse creates rapid oscillations in the quasiparticle densities. The essential physics of multiband superconductors can be captured by a minimal two-band model with Hamiltonian , with and

(1) | |||||

where is the electron (hole) annihilation operator with band index , spin , and momentum . The bare dispersions of the two bands, and , are assumed to be quadratic with , where and represent the effective electron (hole) mass and Fermi energy of band , respectively. The intraband pairing interactions are denoted by and , while the interband pair scattering is represented by . Assuming that the splitting of the bands is much larger than the superconducting energy scales, we neglect all interband pairing terms, i.e., we set for all . Furthermore, we restrict ourselves to momentum-independent pairing interactions with , , and . The intraband pair couplings are assumed to be attractive, i.e., . The phase of is taken to be either or , depending on whether the interband pair scattering is attractive or repulsive.

The dynamics of superconductors at times shorter than the quasiparticle energy relaxation time can be fully described within mean-field BCS theory[12, 11]. We therefore perform a mean-field decoupling of the quartic interactions in Eq. (1), which yields the BCS Hamiltonian

(2) |

with the mean-field gaps

(3) |

The sums in Eqs. (2) and (Theory of nonequilibrium dynamics of multiband superconductors) are over the set of momentum vectors with and the cut-off frequency. The quasiparticle spectrum of is given by , where , with . In the equilibrium, the order parameters and are determined by the self-consistency equation

(4) |

where , with . The gap equation (4) has a nontrivial solution whenever the determinant of the corresponding secular matrix vanishes, i.e.,

(5) |

From this condition, the transition temperature can be obtained.

The two-band superconductor (2) is perturbed by a short Gaussian pump pulse, which in the Coulomb gauge is described by a transverse vector potential , with amplitude , photon wave vector , and full width at half maximum . The optical pump pulse is centered at and has a central frequency . The coupling of the vector potential to the two-band superconductor is given by

(6) | |||

By means of the density-matrix formalism, we numerically compute the coherent response of the above model Hamiltonian after excitation by a short pump pulse. To this end, we need to derive equations of motion for the Bogoliubov quasiparticle densities. It is advantageous to perform this calculation in the basis in which , Eq. (2), at , i.e., in the initial state, takes diagonal form. Hence, we consider the following Bogoliubov transformation

(7) |

where the Bogoliubov coefficients and are time-independent, with and , where . The phase factors of the superconducting gaps in the initial state are chosen to be either or , depending on whether in the initial state is attractive or repulsive. Substituting Eq. (7) into Eq. (2) yields

(8) | |||||

where , and . Note that here, as before, all terms except are time-independent. In particular, we have assumed that has negligible time-dependence. Combining Eqs. (Theory of nonequilibrium dynamics of multiband superconductors) and (7) we find that the superconducting order parameters at time can be expressed as

(9) | |||||

Hence, the temporal evolution of the gaps is fully determined by the time-dependence of the Bogoliubov quasiparticle densities , , , and . As shown in the Appendix, it is straightforward to derive a closed set of equations of motion for these quantities, see Eq. (0.2).

We numerically solve the equations of motion for the Bogoliubov quasiparticle densities (see Appendix) to determine the time evolution of the order parameter amplitudes . As is evident from Eq. (0.2), off-diagonal elements in the quasiparticle density matrices, such as, e.g., , are of order . Thus, for sufficiently small , the off-diagonal entries decrease rapidly as increases. Therefore we set all off-diagonal elements with to zero, which substantially reduces the computational costs. For the numerical simulations we choose the BCS ground state at zero temperature as the initial state. In what follows, we consider two model systems, one with weak attractive , relevant for MgB, and the other with strong repulsive , relevant for iron-based superconductors.

We start by discussing the model appropriate for the multiband superconductor MgB. The electronic structure of MgB consists of two sets of bands: -bands forming quasicylindrical hole pockets and -bands giving rise to tubular networks of electronlike Fermi surface sheets. Superconductivity in MgB arises from couplings between the electrons and the E-phonon mode. These couplings lead to attractive intraband interactions, which are stronger in the -bands than in the -bands. Since the interband pair scattering is weak, this gives rise to two gaps with different magnitudes on the - and -bands.

We describe the band structure of MgB in a somewhat oversimplified manner by two parabolic bands[16], i.e, a hole band, with Fermi velocity m/s and effective mass , and an electron band, with m/s and . Here, denotes the free electron mass. The intraband pairing interactions on the - and -bands are chosen to be meV and meV, respectively, whereas the attractive interband coupling is taken to be much smaller than either or , i.e., we set . With these material parameters, the superconducting gaps in the initial state can be computed from the self-consistency equation (4). Assuming that the cut-off energy is approximately equal to the E-phonon energy, i.e, meV, we find that in the initial state the gaps on the - and -bands are given by meV and meV, respectively. The central energy of the optical pump pulse is chosen to be either in between and , or equal to twice the gap on the -band, i.e., meV.

Let us first discuss the most interesting case, namely the situation where the pump photon energy lies in between and . In Figs. 1(a) and 1(b) we plot the calculated time evolution of and , respectively, for meV, , and three different pump pulse lenghts , 1.0, and ps. In accordance with previous results on single band superconductors[12, 17, 14, 13], we observe two distinct dynamical regimes: a nonadiabatic regime for pulse widths much smaller than the dynamical time-scales of the superconductor and an adiabatic regime for . In the nonadiabatic regime the Bogoliubov quasiparticle densities build up coherently while the system is out of equilibrium. This gives rise first to a monotonic increase and then to fast oscillations in the quasiparticle occupations. Similarly, the superconducting gaps first decrease monotonically and then start to oscillate rapidly (solid green curves in Fig. 1). Since our model system does not contain any relaxation terms, such as interactions or couplings to an external bath, the oscillations in the quasiparticle occupations are undamped. In contrast, the order parameter oscillations decay as . That is, due to decoherence effects, the superconducting gap approaches a constant value as time tends to infinity. To be more precise, we find that to a good approximation the time evolution of is described by

(10) |

where and are constants that depend on the initial conditions. In other words, the superconducting gaps oscillate with an amplitude decaying as and two frequencies, which are determined by the long-time asymptotic gap values [see Figs. 1(a) and 1(b)]. Interestingly, the relative strength of the two frequency components and depends sensitively on the magnitude of . That is, the amplitudes and increase with increasing interband pair scattering . For , on the other hand, the amplitudes and vanish, and hence Eq. (10) reduces to the well-known result for single band models, where the gap oscillates with a single frequency[12, 17, 14, 13, 18].

A nonvanishing interband pair scattering not only affects the relative strength of the oscillation frequencies but also strongly influences the long-time asymptotic gap values and hence the frequencies of the order parameter oscillations. This is illustrated in Figs. 2(a) and 2(b), which show the dependence of on the integrated pump pulse intensity for meV and both for zero and nonzero . Comparing Figs. 2(a) and 2(b), we observe that in the presence of nonzero interband pair scattering the time evolutions of the two gaps are coupled together. The change in due to interband interactions seems to be roughly proportional to V, with [cf. Eq. (9)]. That is, the larger gap influences the smaller gap more strongly than vice versa. With increasing integrated laser intensity both gaps are suppressed simultaneously, until for sufficiently large they vanish at the same value of . For zero , in contrast, the order parameters can be suppressed independently. For example, we find that in the nonadiabatic regime the larger gap vanishes at a smaller integrated intensity than the smaller gap [solid green curve in Fig. 2(b)]. We observe that in general the asymptotic gap values decrease linearly at small , but deviate from this behavior at larger integrated intensities. With increasing , the curves corresponding to short pump pulses ( ps) show a downward bend, while those with longer pump pulses ( and ps) flatten due to Pauli blocking.

Let us also briefly discuss the case where the pump pulse energy is equal to meV. In general, the behavior of the gaps in this case is rather similar to the previously discussed case, where meV. As before, we find that for pump pulse widths both gap amplitudes decrease monotonically from their initial equilibrium values to their final values , whereas for the gaps approach their asymptotic values in an oscillatory fashion according to Eq. (10). However, as can be seen from Figs. 2(c) and 2(d), the dependence of on the integrated intensity is quite different than in Figs. 2(a) and 2(b). Interestingly, we find that a pump pulse with central energy depletes the superconducting condensate more efficiently on the -band than on the -band. In particular, for the gap on the -band is almost unaffected by the pump pulse, even for large , whereas the asymptotic gap value on the -band decreases steadily with increasing . For nonzero the behavior is similar, although here the gap on the -band is slightly suppressed, which is due to the coupling with the gap on the -band.

Secondly, we consider a model appropriate for the iron pnictide superconductors[19]. Many features of iron-based superconductors can be captured within a simple two-band model with quasinested holelike and electronlike pockets centered at the and points in the Brillouin zone, respectively [20]. Due to the quasinesting of the Fermi surfaces, the repulsive interband pair scattering in such a description is much larger than the intraband pairing potentials. The latter are assumed to be overscreened due to the nesting effects and become attractive for low energies relevant for superconductivity. Therefore, in what follows, we model the band structure of iron-based superconductors by two parabolic bands, a hole band with effective mass and Fermi velocity m/s, and an electron band with and m/s[21]. For convenience we shift the location of the electron Fermi surface from the point to the point, i.e., we set , where . The attractive intraband pairing interactions on the electron and hole Fermi surfaces are chosen to be meV and meV, while the repulsive interband pair scattering is set to meV. With these material parameters and assuming that the cut-off energy is determined by the energy of the spin fluctuations in the paramagnetic state, meV[22], we find that the equilibrium gaps are given by meV. That is, the gaps on the hole and electron bands are equal in magnitude but opposite in phase, which is roughly in line with experimental findings in some of the hole-doped iron-based superconductors [21]. The two-band superconductor is driven out of equilibrium by a pump pulse with energy meV, which is of the same order but slightly larger than twice the gap amplitudes in the initial state.

Based on this model for pnictide superconductors, we compute the temporal evolution of the superconducting order parameters on the hole and electron bands. In Fig. 3(a), the time dependence of is shown for a pump pulse with length ps, i.e., for the nonadiabatic regime. We observe that both gaps exhibit almost the same behavior: they first decrease monotonically and then start to oscillate with different frequencies producing a beatinglike pattern in the amplitudes. Due to the smaller Fermi velocity and the larger effective mass of the hole band, the gap on the hole band is slightly less suppressed than the gap on the electron band . Because of the nonzero interband coupling , both gaps oscillate with two frequencies, and , which differ very little, i.e., . This in turn leads to a pronounced beating phenomenon, as can be seen in Fig. 3(a).

The dependence of the asymptotic gap values on the integrated pump pulse intensity is plotted in Fig. 3(b). Due to the large interband coupling the gaps on both bands decrease almost synchronously with increasing . However, the asymptotic value of the gap on the electron band is always slightly smaller than the asymptotic value of the gap on the hole band . As discussed above, this is because of the larger effective mass and the smaller Fermi velocity on the hole band. In general, we find that the relative difference between and gradually increases with increasing pump pulse intensity. Correspondingly, the difference between the two oscillation frequencies and increases, and hence the beating phenomenon becomes less pronounced, as the integrated pump pulse intensity is increased.

In this work we analyzed the nonequilibrium dynamics of two-band superconductors after excitation by a short opitcal pump pulse. We considered two model systems: one with dominant intraband pairing and weak attractive interband pair scattering, and one where the repuslive intrerband interactions are much larger than the intraband pairing potentials. The former model is relevant for MgB superconductors, whereas the latter one is appropriate for iron-based superconductors. For both of these model systems we numerically computed the time evolution of the order parameters as a function of pump pulse duration and integrated pump pulse intensity. Our main observation is that the ratio between the gaps for asymptotically large times depends sensitively on the interband Cooper-pairing strength and differs from its equilibrium value. This allows in principle to use pump-probe experiments for direct identification of the interband Cooper-pair scattering strength. In addition, sufficiently short pump pulses create fast oscillations in the gap amplitudes and the quasiparticle densities of these two-band superconductors. We found that for nonzero interband pair scattering these oscillations are characterized by two frequencies which are determined by the long-time asymptotic values of the gaps on the two bands (see Fig. 1). We showed that the relative strength of the two frequency components sensitively depends on the magnitude of the interband interaction . When the gaps on the two bands are of similar magnitude (which is the case, for example, in iron-based superconductors) the relative difference between the two oscillation frequencies is small, and hence the gaps oscillate with a beatinglike pattern in the amplitudes.

## 0.1 Acknowledgments

We thank A. Avella, U. Bovensiepen, B. Kamble, N. Hasselmann, and G. Uhrig for useful discussions. The work of A.A. and I.E. is supported by the Merkur Foundation.

## 0.2 Appendix

In this Appendix, we present one of the equations of motion (as an example) which determine the time evolution of the quasiparticle densities , , , and . Here, as in the main text, denotes the band index. Since and are related by Hermitian conjugation, only six out of the eight quasiparticle densities need to be evaluated numerically. By use of Heisenberg’s equation of motion, the following equation for can be obtained

where we have introduced the short-hand notation and , with the coherence factors and . The functions and are defined in the main text. In deriving Eqs. (0.2), we used the relation .

### References

- Y. Kamihara, T. Watanabe, M. Hirano, and H. Hosono, J. Am. Chem. Soc. 130, 3296 (2008).
- J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani, and J. Akimitsu, Nature (London) 410, 63 (2001).
- V.A. Moskalenko, Fiz. Met. Metalloved. 8, 503 (1959) [Phys. Met. Metalog. 8, 25 (1959)].
- H. Suhl, B. T. Matthias, and L. R. Walker, Phys. Rev. Lett. 3, 552 (1959).
- J. Demsar, R.D. Averitt, A.J. Taylor, V.V. Kabanov, W. N. Kang, H. J. Kim, E. M. Choi, and S. I. Lee, Phys. Rev. Lett. 91, 267002 (2003).
- T. Mertelj, V. V. Kabanov, C. Gadermaier, N. D. Zhigadlo, S. Katrych, J. Karpinski, and D. Mihailovic Phys. Rev. Lett. 102, 117002 (2009); T. Mertelj, P. Kusar, V.V. Kabanov, L. Stojchevska, N.D. Zhigadlo, S. Katrych, Z. Bukowski, S. Weyeneth, J. Karpinski, D. Mihailovic, Phys. Rev. B 81, 224504 (2010).
- B. Mansart, D. Boschetto, A. Savoia, F. Rullier-Albenque, A. Forget, D. Colson, A. Rousse, and M. Marsi, Phys. Rev. B 80, 172504 (2009).
- B. Mansart, D. Boschetto, A. Savoia, F. Rullier-Albenque, F. Bouquet, E. Papalazarou, A. Forget, D. Colson, A. Rousse, and M. Marsi, Phys. Rev. B 82, 024513 (2010).
- L. Rettig, R. Cortes, S. Thirupathaiah, P. Gegenwart, H.S. Jeevan, M. Wolf, J. Fink, and U. Bovensiepen, e-print arXiv:1110.3938 (unpiblished).
- D. H. Torchinsky, G. F. Chen, J. L. Luo, N. L. Wang, and N. Gedik, Phys. Rev. Lett. 105, 027005 (2010); D.H Torchinsky, J. W. McIver, D. Hsieh, G.F. Chen, J.L. Luo, N.L. Wang, and N. Gedik, Phys. Rev. B 84, 104518 (2011)
- A. F. Volkov and Sh. M. Kogan, Zh. Eksp. Teor. Fiz. 65, 2039 (1973) [Sov. Phys. JETP 38, 1018 (1974)]
- E. A. Yuzbashyan, B. L. Altshuler, V. B. Kuznetsov, and V. Z. Enolskii, J. Phys. A 38, 7831 (2005).
- E. A. Yuzbashyan, O. Tsyplyatyev, and B. L. Altshuler, Phys. Rev. Lett. 96, 97005 (2006).
- T. Papenkort, V. M. Axt, and T. Kuhn, Phys. Rev. B 76, 224522 (2007).
- B. Mansart, J. Lorenzana, M. Scarongella, M. Chergui, F. Carbone, e-print arXiv:1112.0737v2 (unpiblished).
- J.M. An and W.E. Pickett, Phys. Rev. Lett. 86, 4366 (2001); Y. Kong, O. V. Dolgov, O. Jepsen, and O. K. Andersen, Phys. Rev. B 64, 020501(R) (2001); D. Varshney, and M Nagar, Supercond. Sci. Technol. 20 930 (2007).
- A. P. Schnyder, D. Manske, and A. Avella, Phys. Rev. B 84, 214513 (2011).
- T. Papenkort, T. Kuhn, and V. M. Axt, Phys. Rev. B 78, 132505 (2008)
- P. J. Hirschfeld, M. M. Korshunov, I. I. Mazin, Rep. Prog. Phys. 74, 124508 (2011).
- A. V. Chubukov, D. V. Efremov, and I. Eremin, Phys. Rev. B 78, 134512 (2008).
- H. Ding, K. Nakayama, P. Richard, S. Souma, T. Sato, T. Takahashi, M. Neupane, Y.-M. Xu, Z.-H. Pan, A. V. Fedorov, Z. Wang, X. Dai, Z. Fang, G.F. Chen, J.L. Luo, and N.L. Wang, J. Phys.: Condens. Matter 23 135701 (2011).
- D.S. Inosov, J.T. Park, P. Bourges, D.L. Sun, Y. Sidis, A. Schneidewind, K. Hradil, D. Haug, C.T. Lin, B. Keimer, and V. Hinkov, Nature Phys. 6, 178 (2010).