Time-resolved optical conductivity and Higgs oscillations in two-band dirty superconductors

Recent studies have emphasized the importance of impurity scattering for the optical Higgs response of superconductors. In the dirty limit, an additional paramagnetic coupling of light to the superconducting condensate arises which drastically enhances excitation. So far, most work concentrated on the periodic driving with light, where the third-harmonic generation response of the Higgs mode was shown to be enhanced. In this work, we additionally calculate the time-resolved optical conductivity of single- and two-band superconductors in a two-pulse quench-probe setup, where we find good agreement with existing experimental results. We use the Mattis-Bardeen approach to incorporate impurity scattering and calculate explicitly the time-evolution of the system. Calculations are performed both in a diagrammatic picture derived from an effective action formalism and within a time-dependent density matrix formalism.


I. INTRODUCTION
When a continuous symmetry is spontaneously broken, collective excitations emerge. In the case of a superconductor, where the complex order parameter ∆e iϕ acquires a finite value below a critical temperature T C , two bosonic modes appear: the massive Higgs mode and a massless Goldstone mode [1,2]. They may be seen as amplitude δ∆ and phase δθ fluctuations of the complex order parameter in the Mexican hat-shaped free energy potential.
In the presence of Coulomb interaction, the Goldstone mode is shifted to optical frequencies by means of the Anderson-Higgs mechanism while the Higgs mode remains a stable gapped excitation in the Terahertz regime [3].
In a two-band superconductor, two gapped Higgs modes and two phase modes exist. While the global phase fluctuation occurs again only at energies close to the plasma frequency in the presence of long-range Coulomb interactions, the relative phase fluctuation, quantized as the Leggett mode, persists as a gapped excitation at low energies [4].
Experimental observation of Higgs and Leggett collective modes is difficult. Since these fields are scalar quantities, no linear coupling to the electromagnetic field exists [2]. Thus, there are no direct experimental signatures in linear response. Only with a coexisting CDW order, a signature of the Higgs mode is visible in Raman spectroscopy [5][6][7][8][9]. As a consequence, experiments need to be performed in the non-linear regime. Here, the challenge is twofold: intense light sources are required but experiments also have to be performed on energy scales mostly within the superconducting gap such that optical excitation of quasiparticles does not deplete the condensate.
Recent developments in ultrafast Terahertz spectroscopy have caused a surge in interest to study collective excitations in non-equilibrium superconductors both in theory [10][11][12][13][14][15] and experiment [16][17][18][19][20][21][22][23], where first exper-imental signatures of the Higgs mode have been observed in various materials. The main excitation scheme so far consist of two approaches. First, samples are illuminated in a pump-probe setup where an excitation of the Higgs mode by a single-cycle THz pump acting as a quench appears as an oscillation of the probe signal as a function of pump-probe delay [16]. In a second type of experiments, the Higgs mode is resonantly driven by an intense multi-cycle pulse that yields an electrical field component of three times the pump frequency in the reflected or transmitted beam [17,20,22].
The fact that characteristics of the Higgs mode in superconductors are observable in experiments is not selfevident. Early theoretical calculations in the clean limit predicted extremely weak experimental signatures that relied on breaking of the particle-hole symmetry. Therefore, the first observations [17] of the third-harmonic response generated by the Higgs mode was doubted [24] as it should be overlaid by much stronger charge fluctuations. Only recently, the role of impurities has been appreciated as it drastically enhances the coupling of light to the Higgs mode due to an additional paramagnetic coupling absent in the clean limit [25][26][27][28]. This coupling becomes the dominant contribution even for small disorder. It was further shown that impurity scattering yields qualitatively different behavior in the polarization dependence of the driving pulses [28].
While previous studies on impurities concentrated mostly on the excitation scheme with periodic driving, in this work, we additionally explore the excitation with a two-pulse quench-probe scheme. We consider both oneand two-band superconductors where the bands can be in different impurity regimes.
We also calculate the individual contributions of quasiparticles, Leggett mode and Higgs mode to the thirdharmonic generation response. Our results support the findings of a recent work, where the third-harmonic response in the two-band superconductor MgB 2 shows a resonance only for the lower gap [22]. This can be under-stood from the fact that the upper band is either in the clean limit or that the Fermi surface is very small.
We incorporate the effect of impurities in our model using the Mattis-Bardeen approximation [29]. This approach constitutes an excellent description for many conventional superconductors [30]. To calculate the timeresolved optical conductivity, we extend the densitymatrix approach of [25] to a two-pulse excitation scheme. Here, the short first pulse acts as a quench, while the second probe pulse with variable time-delay probes the dynamics of the system. In addition, we consider a diagrammatic approach derived from an effective action formalism, where the Mattis-Bardeen ansatz is incorporated by an effective finite momentum interaction vertex. This diagrammatic approach is equivalent to the density matrix formalism but allows to understand the involved processes in more detail.
This article is organized as follows. In Sec. II we formulate the model (a) in terms of a diagrammatic expansion of an effective action and (b) in terms of a density matrix equation of motion approach that was previously established by Murotani and Shimano [25]. The two formulations are equivalent. We then discuss results in the case of a single-band superconductor in a pump-probe scenario in Sec. III. In Sec. IV we study in detail the case of a two-band superconductor motivated by material parameters of MgB 2 . Here we focus on both pump-probe and third harmonic generation (THG) experiments. We summarize all results in Sec. V.

A. Hamiltonian
We consider the BCS multiband Hamiltonian where ik = s i k 2 /2m i − Fi is the parabolic dispersion of the i-th band with Fermi-energy Fi and electron mass m i . The factor s i = ± determines electron-or hole-like character of the respective band. At the mean-field level the interacting term is decoupled in the pairing channel, where order parameters ∆ i are self-consistently determined by the BCS gap equation ∆ i = jk U ij c j−k↓ c jk↑ [31]. The order parameters of different bands are mixed by off-diagonal terms in the coupling matrix U ij . In the present work, we parametrize gap-mixing by a parameter v and define For given ∆ i and v we can find U 11 and U 22 such that the gap equation is satisfied.
To model an experimental probe with a laser pulse, we introduce a time-dependent vector potential A(t) = A(t) e with polarization vector e by means of minimal coupling, where J ikk = ik| epi mi |ik are intraband transition matrix elements of the current operator. The two terms in H 1 corresponds to the paramagnetic and diamagnetc coupling of the laser field, respectively. The full Hamiltonian is given by H = H 0 + H 1 .

B. Impurity scattering
In a clean system momentum conservation yields J ikk ∼ δ kk , or J ikk ∼ δ k,k ±q if a photon wavevector q is considered. In disordered systems, translational invariance is broken, so that transitions between states of different momenta are allowed. Here, we adopt the approach of Murotani and Shimano [25] and model the effects of impurities within the Mattis-Bardeen (MB) approximation [29]. Explicitly, impurities enter through the approximation with Fermi velocity v Fi , density of states at the Fermi level N i (0) and impurity scattering rate γ i . A derivation of this matrix element is given in Ref. [25]. We see that impurity scattering broadens the δ kkdistribution into a Lorentzian of width γ i centered at zero momentum transfer. The bandstructure defined by H 0 remains unaffected in this approximation. Instead of broadening the momentum resolution of the bandstructure, one may view impurities as effectively broadening the momentum of the photon.

C. Effective Action
We first present a perturbative solution of above Hamiltonian by a path-integral formalism in imaginary time τ [32][33][34][35]. The full problem is formally captured by the partition function Z = D(c † c)e −S with the euclidean action k, ω n +Ω The coupling matrix J is defined in Eq. (A16). (e) Paramagnetic coupling of Higgs modes with susceptibility χ σ 0 σ 0 σ 1 . (f) Paramagnetic coupling of phase modes with χ σ 3 σ 3 . Other couplings at Gaussian level vanish in the presence of particle-hole symmetry.
As detailed in Appendix A, we decouple the interacting part of H in the paring channel, introducing collective fields ∆ i (ω n ) exp (iθ(ω n )). ∆ i and θ i describe amplitude and phase fluctuations, respectively, of the superconducting condensate. In the present work we restrict ourselves to collective fluctuations in time only, i.e. we focus on k = 0 excitations of Higgs and phase fields.
Performing the fermionic path integral results in an effective action S[∆ i , θ i , A] in bosonic and classical EM fields (see Eq. (A12)), where now Z = We only keep terms quadratic in collective fields ∆ i , θ i and to fourth order in A. The resulting terms are diagrammatically presented in Fig. 1 and Fig. 2 and their integral expressions are derived in Appendix A.
The diagrammatic representation contains Higgs fields ∆ i (ω) (blue-dashed lines), phase fields (green-dotted lines), and EM fields (wavy lines). Paramagnetic coupling to the photon field corresponds to vertices with a single photon field line, implying the factor A(ω). Diamagnetic vertices with two photon field lines contribute the term A 2 (ω) = dω A(ω − ω )A(ω ). Only paramagnetic vertices introduce external momentum. Solid black k, ω n +Ω k' ,ω n k,ω n k, ω n +Ω+Ω' k' , ω n +Ω  lines correspond to mean-field Nambu Green's functions and loops imply a trace over Nambu indices, frequencies, and momenta. Figs. 1 and 2 are a complete representation of all terms in the quadratic action in the presence of particle-hole symmetry and impurities in the MB approximation.
In the clean limit paramagnetic photon lines no longer carry momentum and, as a consequence, diagrams 1(e) and 2(a) vanish. The inclusion of paramagnetic diagrams with vertices J ikk determined by the MB model is the main difference of the diagramatic formalism from other literature [24,36].
Absence of diagram 1(e) in the clean limit implies that the Higgs mode does not couple to light without impurities. However, when a non-parabolicity of the bandstructure is taken into account, a diamagnetic coupling to the Higgs mode arises, yielding an additional, non-vanishing diagram [12,24,25].
We note that paramagnetic and diamagnetic terms do not mix in the present model. Consequently, the partition function factors into two contributions Z = Z para Z dia .
Since only the paramagnetic part is affected by impurities, and since Z para does not contain phase contributions, we conclude that only the Higgs mode and quasiparticles are sensitive to impurity scattering in the MB approximation.
The path integrals over ∆ i , θ i can be performed exactly at Gaussian level. This is equivalent to an RPA renormalization of the quasiparticle terms diagrammatically represented in Fig. 3 where the dashed and dotted lines correspond to coupling matrices U ij /2 and Josephson coupling matrices J −1 ij , respectively. One is left with S[A(ω)], explicitly given in Eq. (A28). A functional derivative with respect to A(ω) gives the current

D. Density matrix equations of motion
We solve for the time dynamics of above Hamiltonian using a density matrix approach. To this end, we define the density matrix ρ = |ψ 0 ψ 0 |, or, in the basis of Bogoliubov-de Gennes, we have The time dependence of ρ is given by Heisenberg's equation of motion, where H is the operator H in the BdG basis. We are interested in computing the dynamics of the current j = − δH δA = j P + j D , consisting of a paramagnetic and diamagnetic contribution, as well as the dynamics of the superconducting order parameter To apply the MB substitution, we further expand above equations of motion in orders of A(t). To account for effects of a THG response, we consider terms up to third order. As detailed in Appendix D, the current only has odd order components j = j 0 + j 3 + . . . and the gap contains even contributions of A, ∆ = ∆ 0 + δ∆ 2 + . . . .
Finally, we exploit the rotational invariance of our model and perform the integral over angular degrees of freedom explicitly. Thus, by replacing all momentum summations by an integral k → N i (0) d ik dΩ k 4π , we effectively reduce the model to a one-dimensional system. Note that rotational invariance of our continuum model neglects polarization dependence of observable quantities.
We are left to compute the equations of motion of the first order quasiparticle expectation values, ρ ikk 1 , and the angle-averaged quantities We solve them numerically using a Runge-Kutta solver on a discretized energy grid |ki| of up to 10 3 points in the interval [−ω D , ω D ]. A detailed derivation and explicit presentation of the full equations of motion is given in Appendix D.

A. Optical conductivity
We begin by computing the optical conductivity in linear response, This can be done in either of two ways. First, by implementing a time-dependent density matrix simulation with pulse A(t). The numerically evaluated current j(t) 1 and the pulse are then Fast-Fourier transformed and Eq. (18) is evaluated. Here, one needs to choose a pulse of sufficient ω-bandwidth such that the region of interest is covered. The second way involves the functional derivative of the diagrams in Fig. 2(a,b) according to Eq. (9). At T = 0 one obtains the expression for the real part where E = √ ∆ 2 + 2 , W ( , ) is the Lorentzian of Eq. (6), N the density of states at the Fermi surface, and η is an infinitesimal positive constant.
We can understand the analytical structure of σ (ω) by inspecting the susceptibility χ σ0σ0 ( , , ω). For ω < 2∆ it vanishes exactly. For ω > 2∆ its structure is exemplary shown in Fig. 4. We observe two straight spectral lines at = ±ω + . These features can be understood in the picture of a particle-hole or hole-particle excitation process, illustrated in Fig. 4(a). χ σ0σ0 ( , , ω) has nonzero spectral weight at given , if an occupied state at can be excited into a state at by a photon of frequency ω. Multiplication of the integrand in Eq. (19) with W ( − ) enforces momentum conservation.
In this picture it is easy to see that the total spectral weight χ σ0σ0 (ω) = d d χ σ0σ0 ( , , ω) should be approximately proportional to Θ(ω − 2∆)(ω − 2∆), where Θ is the Heaviside function. Since W ( − ) is constant along contours = ±ω + , we find the simple analytical approximation that holds for ω 2∆ in the dirty limit γ ∆. In Fig. 5 we plot numerically evaluated real and imaginary parts σ (ω), σ (ω) of the optical conductivity for various impurity concentrations and temperatures. σ shows a clear conductivity gap below 2∆. In the clean limit, a pronounced coherence peak is observed around 2∆, reflecting the additional density of states amassed above the quasiparticle gap. The conductivity peak grows and shifts to higher ω as γ is increased. It then broadens into the characteristic dome shape frequently observed in experiment [16,29,38]. In the T → 0 limit, the conductivity is expected to show a condensate δ-peak at ω = 0 which is not numerically resolveable. Instead, we observe a buildup of spectral weight around ω = 0 as the condensate peak is broadened at finite temperatures. The imaginary part σ follows a 1/ω power law as expected for a superconducting state.
The linear response optical conductivity contains information of the bandstructure only and is unaffected by collective modes. This can be inferred from the diagrammatic description where all terms in the RPA renormalization of diagram Fig. 2(a) containing k = 0 collective fluctuations vanish exactly. To reveal the presence of collective modes, we turn to the dynamics of the superconducting order parameter and the non-linear current j 3 and additionally model realistic THz pulses in a pump-probe setting.

B. Excitation of Higgs mode
We choose the electromagnetic pulse form A(t) = A 0 exp −(t − t ) 2 /2τ 2 cos Ωt with coefficients to match the reported data of Ref. [16]. The resulting waveform is shown in Fig. 6(a).
A characteristic property of a pump pulse is its pulse length τ compared to the natural timescale of the superconductor 1/∆. For τ 1/∆ the superconductor is quenched, while it is adiabatically driven in the opposite limit of τ 1/∆.
The different behavior in the two limits can be intuitively understood within the diagrammatic picture. Here, the pulse induced change of the order parameter δ∆(ω) is given by the diagram in Fig. 7(a) which has the integral expression Presence of a collective Higgs mode translates into a peak of the kernel K(ω) = (χ σ1σ1 (ω) + 2/U ) −1 at the characteristic mode energy ω H = 2∆. Excitation of the collective mode, however, is only possible if energy conservation is satisfied, i.e. if A(ω )A(−ω H − ω ) is finite for some ω . Higgs oscillations are therefore expected when the Fourier transform of the squared vector potential A 2 (ω) = dω A(ω − ω )A(ω ) overlaps with the mode-energy ω H . The double-peaked structure of A(ω) is shown in Fig. 6(c). The first peak, centered at ω = 0, corresponds to a difference frequency generation process (DFG), while the second peak at ω = 2Ω corresponds to a sum frequency generation process (SFG). The resonance frequency of the Higgs mode, ω H , is illustrated by a vertical line. Remaining terms in Eq. (21) describe the coupling to light in presence of impurities and ensure momentum conservation in a virtual two step excitation process.
Let us now consider two limiting cases of the optical pulse width. For ∆τ 1, the frequency spectrum of A 2 (ω) is very broad. The response of δ∆(ω) is then dominated by the sharp resonance peak of K(ω) giving rise to pronounced 2∆-oscillations of the superconducting gap in the time domain. Since the DFG peak is guaranteed to overlap with the Higgs resonance, these oscillations will always be present, independent of the frequency of H ij ω n +ω ω n ω n +ω+Ω' the optical pulse. The SFG process only contributes if the pulse frequency lies in the vicinity of Ω ≈ ∆.
In the transient limit, ∆τ 1, the spectrum of δ∆(ω) is finite only for a narrow region around 2Ω. In the timedomain, the gap shows forced 2Ω-oscillations which are resonantly enhanced for 2Ω ≈ 2∆.
Following Matsunaga [16], we choose a pulse with ∆τ = 0.68, closest to the quench scenario, and perform simulations within the density-matrix formulation. The order parameter responds to the THz pulse by a marked drop followed by damped oscillations around a new asymptotic value ∆ ∞ = ∆(t → ∞) of frequency 2∆ = 0.6 THz as displayed in Figs. 6(d-e). The drop of the equilibrium gap is captured by the ω = 0 component of δ∆. Evaluating Eq. (21) for ω = 0, one finds that χ σ0σ0σ1 (ω = 0, ω , k, k ) is finite only for ω > 2∆, similar to the discussion in Sec. III A. Consequently, δ∆(0) is non-zero only if |A(ω)| 2 overlaps with the quasiparticle continuum, which is illustrated in Fig. 6(b). In physical terms, depletion of the superconducting order parameter is a consequence of quasiparticle excitation by A(ω).
Both the oscillation amplitude and ∆ ∞ show a strong dependence on the impurity scattering rate and are peaked at γ ≈ ∆ as shown in Fig. 6(f). This is a consequence of momentum conservation. For γ → 0, Higgs oscillations vanish exactly.
We note that order parameter dynamics are expected to show oscillations of frequency 2∆ ∞ and not, as in our case, 2∆(t = 0) [11,39]. 2∆ ∞ oscillations have also been observed in experiment [16]. The discrepancy can be attributed to the expansion in powers of the pump field A(t) performed in the time-dependent density matrix formalism. If contributions to δ∆ beyond the second order are considered, the oscillation frequency of the order parameter should correctly reflect the non-equilibrium value 2∆ ∞ . Strictly speaking, our model is fully valid only in the limit of small pump fields where the difference between ∆ and ∆ ∞ is negligible.

C. Pump-probe spectroscopy
Higher orders of the optical conductivity include contributions of collective modes that smooth out the absorption edge and add spectral weight inside the conductivity gap. Here, we calculate the non-linear contribution, in a pump-probe setting of the time-dependent densitymatrix formalism. To this end, we pump the system with an intense pulse of fluence A 0 = 0.5 × 10 −8 J s C −1 m −1 and, after a delay δt pp , apply a weak probe pulse. Following experimental schemes [40], we perform two simulations. First, we simulate both a pump and a probe pulse to compute j pp . In a second simulation we apply the pump only, obtaining j p . We then compute the optical conductivity from the difference in currents j = j pp − j p . This ensures that resilient contributions of the pump do not affect the optical conductivity. Figs. 8(a-b) show the real and imaginary part of the optical conductivity σ(ω, δt pp ) as a function of frequency and pump-probe delay. It can be seen that the third-order contribution j 3 adds spectral weight to the conductivity below absorption gap. The conductivity shows clear oscillations in δt pp , as emphasized in Fig. 8(c) where the dome-shaped envelope has been subtracted and the remaining signal was normalized for each ω. A Fourier transform of these oscillations, shown in Fig. 8(e), indicates that the oscillation frequency matches the resonance frequency of the Higgs mode 2∆.
Our results show that signatures of the Higgs mode are measurable in the pump-probe response of the optical conductivity. Yet, to excite the Higgs mode, impurities are crucial. We find that the calculated time-resolved optical response of a single-band superconductor in the dirty-limit is in good agreement with the experimentally measured response [16].

IV. MULTI-BAND SUPERCONDUCTIVITY
Motivated by the good agreement of the theory with experimental data for a single-band superconductor, we now turn to the case of a two-band superconductor. For concreteness, we focus on the superconducting state of MgB 2 . We model the π-and σ-bands believed to be responsible for superconductivity by choosing material parameters ∆ π = 3 meV, ∆ σ = 7 meV, F,π = 2.9 eV, F,σ = 0.7 eV, m π = 0.85m e , m σ = 1.38m e , ω D = 50 meV, s π = 1, s σ = −1 [41].
Convincing evidence for the two-band character of MgB 2 has been found in tunneling measurements [42,43] and ARPES [44]. However, optical linear response probes have only revealed signatures of a superconducting gap in the π-band [22,45]. A recent work [22] on third harmonic generation suggests strong evidence of a collective Higgs resonance in the π-band, but no collective response in the σ-band was observed.

A. Optical conductivity
The linear response optical conductivity of multi-band superconductors is additively composed of contributions from the two bands, σ(ω) = σ π + σ σ , where the bandspecific conductivities are determined by a straightforward generalization of Eq. (19). Figures 9(a-b) show optical conductivities for various different combinations of band impurity concentrations.
Experimental measurements of the optical conductivity of MgB 2 below T C show a clear absorption gap below 2∆ π and a dome shaped onset above 2∆ π . A second onset at ω = 2∆ σ has so far not been observed. Our simulations reproduce these findings in two different parameter regimes: in the dirty-clean limit, where only the first gap contributes to σ(ω), and in the dirty-dirty limit shown in Fig. 9(c-d). Latter case only shows a weak onset of the σ-gap which may be unnoticeable with experimental uncertainties. The reason of the subdominant contribution of the second gap lies in the small Fermi surface of the σ-band. Explicitly, this can be seen from the prefactor v Fi N i in Eq. (19). For our choice of parameters, which include a high estimate of Fσ , this yields a suppression of the σ-gap conductivity by a factor v Fπ N π /v Fσ N σ = 6.6. For a more conservative estimate of Fσ , the suppression should be even more pronounced.

B. Collective modes
Pulse induced changes of the two order parameters ∆ i with i = π, σ in the two-band case are given by where and where susceptibilities χ σ0σ0σ1 i , χ σ1σ1 i are listed in Appendix A. The gaps exhibit two resonances which are determined by the Higgs propagator. In Fig. 10 we show a logarithmic false-color plot of the quantity | det H| −1 , responsible for any divergence, as a function of frequency ω and interband coupling strength v. As expected the two resonance energies are at 2∆ π and 2∆ σ , illustrated by solid green horizontal lines. Resonances are sharp at small v but decrease and broaden in the strong interband coupling regime.
Energy conservation in Eq. (23) is established by the factor A(ω )A(−ω − ω ). Oscillation of the gaps is therefore only possible for a finite overlap of A 2 (ω) with the resonance frequencies. The matrix structure of H ij further implies that both gaps will oscillate with all excited modes at finite v. Dynamics of the phase modes θ i in the frequency domain are determined by Due to the Anderson-Higgs mechanism only the dynamics of the phase difference δϕ = δθ π − δθ σ is physical. Inserting Eq. (25) yields the expression Solid and dashed orange lines in Fig. 10 trace the maximum and full width at half max (FWHM) of δϕ(ω)/A 2 (ω). Red diamonds are the dominant oscillation frequency of the phase evaluated by computing δ∆ i in a time-dependent density matrix formulation for a broadband optical pulse. The two methods show excellent agreement. At small coupling the phase exhibits completely undamped oscillations due to the absence of decay channels. The Leggett frequency ω L increases for stronger coupling. Once its energy reaches the quasiparticle threshold it is increasingly damped and the resonance broadens. The present results reproduce the findings of Refs. [36,46] which were obtained in the clean limit. This should come at no surprise since impurities do not change the frequency of the collective resonance within the MB approach and additionally the Leggett mode only couples diamagnetically to electromagnetic fields.

C. Pump-probe simulations
We proceed to model the pump-probe response of a two-band superconductor. Analogous to the single-band case we consider non-linear contributions to the optical conductivity and pump the system with an intense pulse. After some time delay δt pp , the optical conductivity is probed in the linear response regime by a weak probe pulse.
In Fig. 11 we adopt the dirty-dirty limit with γ π = 100 meV and γ σ = 50 meV as a potential description of MgB 2 and select various pump pulses shown in the leftmost panels. Gray and dark gray areas illustrate the onset of the quasiparticle continuum of the two bands. Lower panels show A 2 (ω) where Higgs resonance frequencies are marked by gray vertical lines. The second column shows the gap dynamics δ∆ i (t) following the pump pulse. The third column shows the real part of the time-resolved nonlinear optical conductivity σ (ω, δt pp ). Last two columns plot isolated and normalized conductivity oscillations, obtained by subtraction of the constant dome shaped background, as well as their Fourier transform along the δt pp axis.
The first pump has a broad frequency spectrum such that it overlaps with both Higgs resonances. Following the excitation, both gaps oscillate with both frequencies.
The overlap of A(ω) with the quasiparticle continuum induces a small drop of δ∆ . The optical conductivity shows oscillations in the pump-probe delay δt pp with mostly 2∆ π and a small 2∆ σ component. We attribute the subdominance of the δ∆ σ -contribution to the small σ-band Fermi surface.
For a narrowband pulse centered at ω = ∆ π (second row), we observe 2∆ π oscillations only. Here, the pulse A(ω) does not overlap with the quasiparticle continuum. As a result, the gap oscillates around its equilibrium value ∆ ∞ = ∆.
When the narrowband pulse is centered around the second Higgs resonance at ω = 2∆ σ (third row), the gap performs 2∆ σ oscillations only. However, the response is weak and numerically hard to resolve in the optical conductivity.
We note that simulations presented in Fig. 11 show no signatures of the Leggett mode. This is because Higgs and quasiparticle contributions dominate the optical response even at small disorder.

D. Third harmonic generation
Finally, we simulate the non-linear response of a multiband superconductor in a THG setup within the timedependent density matrix framework. We model a realistic multi-cycle pulse of frequency Ω, exemplary shown in Fig. 12, and compute the third order current j(t) 3 . The Fourier transform of j(t) 3 reveals a 3Ω third harmonic (TH) component next to the original first harmonic (FH) peak.
We adopt the dirty-dirty band description of MgB 2 with γ π = 100 meV, γ σ = 50 meV and choose two different interband coupling strengths, v = 0.05 and v = 0.4. Then, we sweep temperature to investigate the resonant behaviour of the TH component. We consider three pulses of frequencies Ω = 0.5, 0.6, 0.7 THz and expect the TH component to be resonantly enhanced when 2Ω = 2∆ i . Figs. 13(a-b) show the temperature dependence of the BCS gap. Horizontal lines mark pulse frequencies Ω used in independent simulations. Resonance conditions are satisfied at intersections with a gap. In the second row, Figs. 13(c-d), the amplitude of the TH peak is found as a function of temperature. In the weak coupling case, v = 0.05, the THG signal for the lower two frequencies exhibits a pronounced peak at the resonance condition for the lower gap. The THG signal peak of the largest frequency is less pronounced, as this frequency is almost equal to the lower gap for a range of temperatures. We also observe much smaller peaks at temperatures where pulses are in resonance with the larger σ-gap.
In the strong coupling case, v = 0.4, we no longer observe a peak-like resonance for the lower π-band gap. This can be understood as a result of broadening of the Higgs resonance at large v, further discussed in Appendix C. The σ-gap still induces a sharp resonance peak, albeit small in comparison to the low-temperature signal.
Panels (e),(f) of Fig. 13 decompose the THG signal for the Ω = 0.5 THz pulse into contributions from the Higgs mode, quasiparticles, and Leggett mode. The Leggett mode contribution is found numerically by considering only the diamagnetic component of the current j D 3 . The quasiparticle contribution is found by forcing δ∆ i = 0 when solving the equations of motion, removing the selfconsistency condition that induces collective modes. In both the weak coupling and large-v case the THG response is dominated by the Higgs mode. The relative contribution of quasiparticles increases in the strong interband coupling regime. The Leggett contribution is vanishingly small. The present results are interesting when compared to the experimental findings of Ref. [22]. Our results affirm the claim that the THG response is mainly attributed to the Higgs resonance of the π-band. The small contribution of the the σ-band Higgs mode and the Leggett mode in our simulation is consistent with the experiment where no signatures of the Leggett or second Higgs mode were observed. We have further computed the THG response in the dirty-clean limit where we found nearly identical results, apart from the absence of the small σ resonance peak at temperatures close to T C .
The failure of our theory to produce resonance peaks of the σ-Higgs mode at large v suggests that the MB approximation might not correctly describe the THG response in the strong coupling limit as assumed for MgB 2 [36,47]. A recent study has found that incorporating impurities beyond Mattis-Bardeen as random onsite-energies in a lattice model shows a stronger contribution of quasiparticles [28]. This, however, is beyond the scope of this paper and will be explored in future investigations.

V. CONCLUSION
We have calculated the time-resolved optical response of dirty multiband superconductors. We have incorporated impurity scattering within the Mattis-Bardeen approximation that effectively broadens the photon momentum distribution of the optical pulse. This approach is known to accurately describe superconductors [16,48], yet deviations in the strong disorder regime are possible [49][50][51]. The response was calculated within two different frameworks. First, the time-evolution of the system after an excitation with a light pulse was calculated explicitly using a time-dependent density matrix formalism. Here, the Mattis-Bardeen ansatz enters through a replacement of the matrix element of the current operator with a Lorentzian-shaped momentum transfer distribution [25]. In the second approach, we calculated the relevant susceptibilities in a diagrammatic formalism derived from an effective action approach. Here, impurities enter through a paramagnetic electromagnetic coupling vertex that carries external momentum. As a consequence, additional diagrams arise that usually vanish in the clean limit. While both approaches yield equivalent results, the diagrammatic approach allows to understand the relevant processes in more detail and is numerically more efficient in certain cases.
In accordance with previous literature [24,25,28], we find that the collective Higgs response is drastically enhanced even for small impurity concentrations. The Leggett mode is unaffected by impurity scattering and hence becomes subdominant. This may change slightly when realistic, non-parabolic bandstructures and weak violations of particle-hole symmetry are taken into account. An interesting further question is the inclusion of Coulomb interactions within the MB approach.
As a first result, we calculated the dynamics of superconducting order paramater of a single-band superconductor in the dirty limit excited by a short THz quench pulse. Using a second probe-pulse after a variable time delay, we further computed the time-resolved optical conductivity. Both quantities show oscillations with the Higgs frequency ω H = 2∆. The optical response is in good agreement to measurements on NbN [16].
Extending the model to two bands, we studied pumpprobe optical conductivities of the model as an effective low-energy description of MgB 2 for various impurity limits of the π-and σ-band. We found that experimental results are well reproduced either when both bands are dirty, or when only the lower π-band is dirty. Here, the collective contribution to the non-linear optical response is always dominated by the amplitude mode of ∆ π .
Finally, we presented the third-harmonic generation (THG) response of the two-band model. Interestingly, results obtained in the weak interband coupling regime seem to match available experimental data, showing a pronounced THG resonance mostly due to the π-band Higgs mode [22]. However, our theory shows deviations in the strong interband coupling case, believed to be representative of MgB 2 , where no obvious π-Higgs-resonance is present.
In summary, we have presented simulations of timeresolved optical conductivities and modelled THG resonance experiments within the Mattis-Bardeen approximation using both a diagrammatic approach and a timedependent density matrix formalism for single-band and two-band superconductors. As studies of collective excitations in superconductors with THz spectroscopy become more and more common, it is important to understand the correct excitation scheme in the presence of impurities. The Mattis-Bardeen approximation, as shown in this work for either a density-matrix formalism or a diagrammatic calculation, allows to describe the effects of impurities in a simple way. This will help further studies in achieving more realistic models of experimental results. The problem is stated with the partition function We decouple the interacting term in the pairing channel via the Hubbard Stratonovich transformation The bosonic field ∆ is complex, i.e. it permits amplitude and phase fluctuations. We decompose it into real fields and additional express fluctuation with respect to the meanfield saddlepoint, . Note that we are neglecting spatial fluctuations of the collective fields, i.e. ∆ i (τ ), θ i (τ ) depend on time only. The action is where Integrating out the Fermions gives where the trace is performed over time, momenta, and Nambu indices, but not over band-indices i. To separate amplitude ∆(τ ) and phase θ(τ ) fields, we introduce a unitary transformation V = exp(iθ(τ )σ 3 /2) [52], We splitG −1 i into a meanfield part, G −1 0,i , and all remaining contributions Σ i . In frequency space, this gives Note that phase fluctuations θ live in the σ 3 channel, i.e. the charge channel.
We expand S at Gaussian level. To compute currents j 1 and j 3 we additionally keep terms up to fourth order in the classical field A. In expanding the trace, we use The quadratic action is given by the terms which are explicitly Here, λ = 8∆1∆2v U22−v 2 U11 . The Higgs and Leggett terms, Eqs. (A14-A17), are diagrammatically shown in Fig. 1. The quasiparticle terms, Eqs. (A18-A21) are shown in Fig. 2. The susceptibilities, given by the fermionic bubbles, are Additional diagrams, that vanish due to particle-hole symmetry, are listed in Fig. 14. We proceed to integrate out all collective fields. This gives We note that phase and fourth order diamagnetic quasiparticle term combine to give the Leggett contributioñ Above equations, obtained by Gaussian integration, have the diagrammatic representation of an RPA summation shown in Fig. 3. For the Higgs propagator this can be seen by expanding where X ij = χ σ1σ1 i δ ij corresponds to fermionic bubbles and U ij /2 corresponds to the to dashed lines. The case of the Leggett mode is analogous. The currents can now, after analytic continuation of all external frequencies, be computed by a functional derivative of the action, ( 1 cm 1 ) density matr. eff. action Figure 15. Real and imaginary part of optical conductivity computed in the time-dependent density matrix formalism (blue lines) and from diagrams Fig. 2(a,b) in the effective action approach. There is perfect agreement between the two methods.
In Fourier space this results in Latter equality follows in full generality from the chain rule of functional derivation, where we have denoted the Fourier transformÃ(ω) = FT [A](ω) by a tilde. The Mattis-Bardeen approximation enters by replacing according to Eq. (5). For the fourth-order paramagnetic quasiparticle contribution, Eq. (A21), we follow Ref. [25] and further approximate kk k Appendix B: First order currents and optical conductivity The paramagnetic first order current j P 1 is represented by the diagram in Fig. 2(a) and explicitly given by a functional derivative of Eq. (A20). After analytical continuation and MB substitution one arrives at The diamagnetic first order current j P 1 reads where we have used that δ δA(ω) A 2 (0) = δ δA(ω) dω A(−ω )A(ω ) = 2A(−ω). Note that the k-sum does not vanish away from the Fermi surface and therefore strongly depends on the numerical cutoff. Here, we follow Murotani [25] and regularize the integral as with and Fermi function f ( ) and the band specific carrier density n i = k 3 Fi /3π 2 . For THG experiments, we are interested in the non-linear current j(3Ω) evaluated at ω = 3Ω where Ω is the dominant frequency of the optical pulse A(Ω): A diagrammatic representation of Eq. (C1) is shown in Fig. 16. The field A(ω) with respect to which the functional derivative is performed is colored red. All four choices are equivalent. The functional derivative forces the external frequency of the field A to be 3Ω. In principal one now needs to integrate over all remaining external frequencies, while satisfying energy conservation. This can be numerically challenging. Here, we focus instead on the case of a monochromatic field A(t) = A 0 cos Ωt, A(ω) = A0 2 (δ(ω − Ω) + δ(ω + Ω)) where external fields possess two discrete frequencies ±Ω. Then energy conservation dictates all remaining external legs to carry frequency −Ω. Note that the energy flow through collective Higgs or Leggett propagators is 2Ω, i.e. THG probes the optical kernel at twice the driving frequency as expected for a non-linear process. Fig. 17 shows magnitude and phase of the Higgs contribution to the THG current j H (3Ω) as a function of Ω and T for two interband couplings v = 0.05, 0.4. Panels (a,b,e,f) correspond to the limit of a dirty π-band and a clean σ-band, whereas the remaining panels are computed for two dirty bands. Both cases are possible descriptions of MgB 2 . Yellow spectral lines map out the Higgs resonance that follow 2∆ π , 2∆ σ . In all cases the π-resonance is dominant, although the relative σ-contribution is enhanced in the dirty-dirty limit and for strong v. Increased interband coupling v decreases and broadens the overall Higgs response.
The Higgs resonance is sharp at small v, but much broader in the v = 0.4 case. Therefore, slices along the T -axis for a given drive frequency Ω do not exhibit a pronounced resonance peak. The observation of a resonance peak in Ref. [22] when experimentally sweeping the temperature would be therefore suggestive of a small v coupling in MgB 2 . This is in disagreement to Refs. [19,47] that experimentally determined a large v based on evidence of the Leggett mode above 2∆ π .
Lower panels in Fig. 17 show a phase jump of π in the THG current across the first Higgs resonance along the Ω direction that is most pronounced at low temperatures. The phase also shows features of the σ-Higgs resonance, albeit  less clearly. Approaching the resonance along the T axis does not yield a phase behavior that is consistently simple to interpret. These results are to be contrasted to the clean case where one expects a phase jump of π/2. Fig. 18 shows the amplitude and phase response of the Leggett THG signal for three different coupling strengths v = 0.02, 0.2, 0.5. The overall contribution is about three magnitudes smaller than the Higgs contribution and therefore negligible. At large coupling, the Leggett resonance is very broad but sharpens at high temperatures. This observation was first reported in Ref. [46]. The phase shows a clear π/2-jump across the resonance for all temperatures below T C .
The quasiparticle contribution is shown in Fig. 19 for v = 0.05 in different impurity cases. Here clean refers to γ = 0.01 meV and dirty specifies γ = 100 meV. Results at different v are nearly identical since the only v-dependent quantity in Eq. (A21) is the superconducting order parameter at finite T . For all impurity concentrations, the quasiparticle THG signal is peaked at the onset of the quasiparticle continuum of the π-band. The signal is about one order of magnitude smaller than the Higgs contribution in the small v case. For v = 0.4, the quasiparticle signal remains nearly identical but the Higgs signal increases, so that the Higgs contribution is only slightly larger. In all but the dirty-dirty case, the quasiparticle signal has a large contribution for small Ω and large T .
The bottom row of Fig. 19 shows the phase of the non-linear THG signal. In the dirty-clean and dirty-dirty cases we observe a clear phase jump of π/2 across a resonance. where and ∆ eq i is the equilibrium value of the gap. We are free to choose the initial phase so that u ik and v ik are both real. Next, we construct the density matrix ρ in the Bogoliubov quasiparticle basis We can now rewrite the Hamiltonian in terms of the quasiparticles and their expectation values In equilibrium, where ∆ i = ∆ eq i , this reduces to (D8) In non-equilibrium, we can write the time-dependent gap as ∆ i (t) = ∆ eq i + δ∆ i (t). Furthermore we can decompose δ∆ i (t) = δ∆ i (t) + iδ∆ i (t) into real and imaginary parts. Then, we find The coupling to the pulse becomes where Notice the total Hamiltonian can be expressed as Next, we decompose the gap into real and imaginary parts ∆ i = ∆ i + i∆ i . After the transformation, the gap equation becomes Finally, we write the paramagnetic and diamagnetic current densities in terms of the quasiparticles (D16) the density matrix in vector form, Heisenberg's equation of motion is stated as where we have defined the following two matrices: Note that we work in natural units where = 1. In order to solve the equation of motion with the impurity scattering replacement, we proceed using a perturbative approach with respect to A. This involves expanding all relevant quantities in powers of A. For instance, we expand ρ ikk as ρ ikk = ρ ikk 0 + ρ ikk 1 + ρ ikk 2 + . . . , where ρ ikk 1 is proportional to A, ρ ikk 2 is proportional to A 2 , etc. We now solve each order separately. The zeroth-order components are simply the equilibrium values, where the quasiparticles are occupied according to Fermi statistics We can also use these values to calculate the equilibrium value of the gap, exchanging the sum in the gap equation with an integral over the energy: Now we proceed with the first order. Considering only the terms proportional to A, the equation of motion becomes where We solve the equation by writing it in terms of new functions F ab ikk . The first order expression of ρ ikk becomes where F ab i is defined according to We have introduced a simplified notation = ik , E = E ik , F ab ikk = F ab i ( , ) etc. The collective modes do not couple linearly to light, and correspondingly δ∆ i 1 = 0. This can be confirmed by using Eq. (D26) in the gap equation. Hence, we only examine the response of the current. By substituting Eq. (D26) into Eq. (D15), changing the momentum sums into integrals over energy, and making the Mattis-Bardeen replacement (5), we find We can also derive the induced diamagnetic current, This derivation is carefully discussed in [25]. We now consider the second order solution. Keeping only terms proportional to A, we find that the off-diagonal terms in the equation of motion vanish, and as a result (D18) becomes We can decompose H 1 ikk 2 and H 2 ikk 2 into contributions from the diamagnetic quasiparticle current, Higgs mode, and Leggett mode as follows.
where the diamagnetic quasiparticle current contribution is the Higgs contribution is To simplify the equation of motion, we introduce a new angle-averaged quantity, r ab i ( ), defined by The equations of motion can then be written as follows The terms r 11 i and r 22 i correspond to the quasiparticle excitations. The remaining r ab i terms correspond to the collective modes. We break up the remaining terms into odd and even components, which are responsible for the Higgs and Leggett modes respectively.
The equations of motion for these terms are [i∂ t − 2E] r 21,even After exchanging momentum sums into energy integrals, the gap equations written in terms of the angle averaged quasiparticles become The equations of motion are solved numerically, and must be solved self-consistently with the gap equations at each time step. This condition induces the collective modes. Finally, to consider pump-probe simulations and THG, we must go to third order. First, the diamagnetic third order current can be directly calculated from the angle-averaged quantities r ab i .
To find the paramagnetic third order current, we start from the third order equation of motion i∂ t ρ ikk 3 = H We proceed by computing the equation explicitly for ρ 11 ikk . Here, we do not consider any contributions from the Leggett mode (i.e. terms involving δ∆ i ) or the A 2 part of the EM field. These terms vanish because of particle hole symmetry. The remaining contributions, consisting of quasiparticles and Higgs mode, are [i∂ t − (E ik − E ik )] ρ 11 ikk 3 = J ik k · e l ik k A ρ 11 ik k 2 − ρ 11 ikk 2 + δ∆ i 2 The next step is to insert this into the expression for the paramagnetic current (D15). Replacing summation over k with integrals over energy, we find j P 3 = e ikk e · J ikk l ikk ρ 11 ikk 3 + ρ 22 ikk 3 + p ikk ρ 21 ikk 3 − ρ 12 occur. The differential equation for these quantities is Now we make the final approximation |J ik k · e| 2 ρ ab ikk 2 Av ≈ |J ik k · e| 2 Av ρ ab ikk Av 2 .
With this, the differential equation becomes ikk Av 3 |J ikk · e| 2 Av = l i ( , ) A ρ 11 ik k Av 2 − ρ 11 ikk Av 2 + δ∆ i 2 Note that |J ik k · e| 2 Av = |J ikk · e| 2 Av . Defining and noting we rewrite this as which precisely gives Eq. 48 in [25]. Similar equations can be derived for the remaining R ab i . Rewriting the expression for the paramagnetic current j P in terms of the R ab i ( , ), and then making the Mattis-Bardeen impurity replacement (5)