# Contribution of Collective Excitations to Third Harmonic Generation in
Multiband Superconductors: the case of MgB2

Jacopo Fiore  Department of Physics and ISC-CNR, “Sapienza” University of
Rome, P.le A. Moro 5, 00185 Rome, Italy  Mattia Udina  Department of Physics
and ISC-CNR, “Sapienza” University of Rome, P.le A. Moro 5, 00185 Rome, Italy
Marco Marciani  Department of Physics and ISC-CNR, “Sapienza” University of
Rome, P.le A. Moro 5, 00185 Rome, Italy  Goetz Seibold  Institut für Physik,
BTU Cottbus-Senftenberg, PBox 101344, 03013 Cottbus, Germany  Lara Benfatto
Department of Physics and ISC-CNR, “Sapienza” University of Rome, P.le A. Moro
5, 00185 Rome, Italy

###### Abstract

Multiband superconductors can host collective excitations with marked
differences with respect to their single-band counterpart. We first study the
spectrum of collective amplitude fluctuations in a clean two-bands
superconductor, showing that the spectral weight of the Higgs mode rapidly
deviates from the naive extension of the single band case as the interband
coupling is turned on. These results are then used to critically analyze the
non-linear optical response in MgB2, providing an explanation for the
apparently contradictory results of recent experiments, pointing towards a
selective relevance either of the Leggett mode or of the amplitude
fluctuations at twice the lower gap. By using exact numerical simulations and
realistic estimate of disorder we compute the relative contribution of the
quasiparticle, amplitude and phase fluctuations to the non-linear optical
response. We show that at low pumping frequency only the resonance at twice
the smaller gap emerges, as due to the BCS response, while the Leggett mode
dominates only in a narrow range of higher pumping frequencies matching its
low-temperature value. Our findings provide a fresh perspective on the
potential of non-linear THz spectroscopy to detect collective modes in other
multiband systems, as e.g. iron-based superconductors.

##  I Introduction

Multiband superconductivity has been reported for a variety of interesting
systems, and it has been linked to both conventional (phonon-mediated) and
unconventional pairing mechanisms. To the former category belongs the
benchmark case of MgB2[1], that has been historically the first well-studied
case of multiband superconductors, but also more recent examples as
superconducting (SC) interfaces between insulating oxides[2, 3] or layered
dichalcogenides[4]. To the latter category belong the widely studied families
of iron-based superconductors, where pairing originates from electronic
correlations, as possibly due to the exchange of spin fluctuations among
quasi-nested hole-like and electron-like bands[5, 6, 7]. Even though rather
often metals display multiple bands at the Fermi level, the characteristic
feature of multiband superconductors that is common to all the above-mentioned
cases is the presence of markedly different SC gaps $\Delta_{\alpha}$ on the
various sheets of the Fermi surface, labelled here with a band index $\alpha$.
This fact can lead to a series of spectroscopic signatures rather different
from the single-band case, since the ratio $\Delta_{\alpha}/T_{c}$ between the
gap values and the critical temperature $T_{c}$ can deviate from the standard,
single-band BCS case, even in the presence of weak-coupling pairing[8, 9].
This effect manifests in all the probes sensitive to thermal activation of
single-particle excitations, as e.g. in the specific heat, and in the
spectroscopic probes sensitive to the particle-hole excitations above the
$2\Delta_{\alpha}$ optical gap in each band. In both cases, the observed
behavior differ from the single-band case, and very often only the lower
optical gap can be resolved in the optical absorption, since this is the only
one in the dirty limit [10, 11, 12, 13, 14].

In addition to the above modifications of the quasiparticle excitations,
multiband superconductors host also a spectrum of collective fluctuations much
richer than in the single-band case[15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26]. Indeed, by considering for the sake of discussion the two-bands case,
each order parameter admits both amplitude and phase fluctuations below
$T_{c}$. In the phase sector, besides the Goldstone mode of the broken $U(1)$
symmetry, associated with the fluctuations of the overall
$\theta_{1}+\theta_{2}$ phase, a second so-called Leggett mode connected to
relative phase fluctuations $\theta_{1}-\theta_{2}$ emerges[15]. While the
Goldstone mode is dynamically pushed to the plasma energy by its coupling to
charge fluctuations, the mass $\omega_{L}$ of the Leggett mode is only
controlled by the pairing strength[15, 16]. For weak interband coupling it can
occur below the lower optical gap, and it can even soften completely in
systems with three bands at the verge of a time-reversal-symmetry-breaking
transition[17, 18, 19, 20]. For what concerns the amplitude fluctuations, the
general behavior is rather different from the single-band case, where the
spectral function of the so-called Higgs mode has a well-defined resonance at
twice the gap value. As it has been already noticed in previous works[22, 23,
25, 26], for multiband superconductors amplitude fluctuations in each band
carry information on the optical gaps of both the bands.

This rich phenomenology, which already challenges the interpretation of
experiments with conventional spectroscopies, makes the understanding of high-
field THz spectroscopies rather complex[21, 23, 25, 26]. An example is
provided by recent findings in MgB2, where both time-resolved protocols[24]
and third-harmonic generation in transmission[27] have been used to
investigate the non-linear optical response triggered by strong THz light
pulses. As it has been widely discussed within the context of single-band
models[28, 29, 30, 31, 32, 33] the main essence of these experiments is to
trigger a non-linear current $j^{NL}$ which scales with the third power of the
THz field $E$, $j^{NL}\sim\chi^{(3)}E^{3}$, where $\chi^{(3)}$ is the non-
linear optical kernel and time convolution has been omitted for simplicity.
Whenever the $\chi^{(3)}$ admits a marked maximum (resonance) at a given
frequency $\omega_{res}$ one can observe, under specific experimental
conditions, oscillations in the differential probe signal at $\omega_{res}$ as
a function of the pump-probe time-delay $t_{pp}$ in time-resolved experiments,
or enhanced THG for a multicycle light pulse centred at approximately
$\omega_{res}/2$ [34, 35]. In general, $\omega_{res}$ can be identified with
any collective excitation of the system non-linearly coupled to the light. In
a superconductor this is represented either by the threshold for the BCS
particle-hole continuum at the optical gap $2\Delta$, or by the collective
amplitude and phase modes. As a consequence, in a multiband superconductor one
could expect in principle multiple signatures: quasiparticles or Higgs
fluctuations at $2\Delta_{\alpha}$ of each band, and the Leggett mode at
$\omega_{L}$. In MgB2 $\omega_{L}$ lies in between the lower $2\Delta_{\pi}$
and the higher $2\Delta_{\sigma}$ optical gap, as it has been proven by
spontaneous Raman measurements[36]. On the other hand, non-linear optics
reported possible resonances in $\chi^{(3)}$ at either $\omega_{L}$[24] or
$2\Delta_{\pi}$[27], with the two features appearing mutually exclusive in the
two experimental set-ups, challenging their interpretation.

Figure 1: Sketch of the various contributions to THG stemming from the
spectrum of the collective fluctuations in MgB2. Dotted lines denote the
temperature dependence of the two order parameters and of the Leggett mode,
while shaded areas mark the typical frequency broadening of the corresponding
resonances in the non-linear optical kernel. Around $\Delta_{\sigma}(T)$ both
the BCS and Higgs fluctuations have a sizeable weight, while around
$\Delta_{\pi}(T)$ BCS fluctuations dominate. The dashed line marks the
temperatures $T^{*}_{\pi}$ and $T^{*}_{\sigma}$ such that
$\Delta_{\pi}(T^{*}_{\pi})/T^{*}_{\pi}=\Delta_{\sigma}(T^{*}_{\sigma})/T^{*}_{\sigma}\approx
2$, that is the condition at which the BCS/Higgs resonances are smeared out
and damped by thermal fluctuations. The dashed-dotted lines denote the pumping
frequencies of the experiment by Giorgianni et al. [24] (red) and by Kovalev
et al. [27] (blue).

In the present manuscript we provide a detailed analysis of the non-linear
response of multiband superconductors, with the twofold aim from one side to
clarify the nature of the collective fluctuations, in particular the Higgs
mode, and from the other side to establish the experimental conditions for
their detection via non-linear optics, that has been recently applied both to
MgB2[24, 27] and iron-based superconductors[37, 38, 39]. For what concerns the
amplitude fluctuations, we will show that the redistribution of spectral
weight between the two peaks at $2\Delta_{\alpha}$ in the spectral functions
depends crucially on the structure of the coupling matrix. This leads in
particular to a pile up of the spectral weight at twice the largest gap for
MgB2 and to a less straightforward redistribution for iron-based
superconductors, crucially depending on the system parameters. In the latter
case, an analogous effect has been observed for the Leggett mode, whose
signature moves at twice the largest case as the SC coupling becomes
predominantly interband[21]. Once established the nature of the collective
fluctuations of each order parameter, we will critically analyze their
manifestation in the non-linear optical response for the specific case of
MgB2. Indeed, as it has been widely discussed in the single-band case[31, 32,
33], the relative contribution of both quasiparticle and Higgs-mode
fluctuations to the $\chi^{(3)}$ response at $2\Delta_{\alpha}$ critically
depend on the disorder level of the sample, with the former ones dominating at
weak and intermediate disorder and the latter one dominating at very strong
disorder. As a consequence, no definitive conclusion can be drawn without an
accurate quantitative estimate of the disorder level of the sample, that can
be obtained by a direct inspection of linear spectroscopy. We will use the
experimentally estimated values $\gamma_{\alpha}/(2\Delta_{\alpha})$ in each
band to extract the contribution of each mode to the non-linear kernel, with
$\gamma_{\alpha}$ being the scattering rate for carriers in the $\alpha$ band.
In particular we extend the clean-limit calculations of Ref. [21, 24] to an
effective model accounting for the disorder-mediated coupling of BCS,
amplitude and phase fluctuations to light, in order to establish the dominant
contributions to the non-linear response as a function of the temperature and
the frequency of the light pulse, as summarized in Fig. 1. We then show that
the apparent discrepancy among the two experiments of Ref. [24] and [27] is
purely linked to the different experimental conditions. In particular the
strong signatures at $2\Delta_{\pi}$ observed in Ref. [27], where the pumping
frequency is of the order of $\Delta_{\pi}(T\approx 0)$, should be ascribed to
BCS quasiparticle excitation in the $\pi$ band. In contrast, the measurements
of Ref. [24] are performed with a larger frequency, matching one half of the
low-temperature value of the Leggett mode, making it the most relevant
contribution to the non-linear response, see Fig. 1. Finally, we show how
multiband MgB2 is the ideal system to estimate the effects of the linear-
response screening in the SC state to correctly match the theoretical
predictions with the experimental results. This issue has been recently raised
by the analysis of experiments on $d$-wave cuprate superconductors[40], but
seems not relevant for conventional $s$-wave systems as NbN. Here we show that
in the multiband case, due to the presence of a low $\Delta_{\alpha}/T_{c}$
ratio in one band, screening effects must be included to capture the correct
temperature dependence of the THG, in analogy with what is seen in the case of
a nodal single-band superconductor.

The plan of the manuscript is the following. In Sec. II we study the spectrum
of the amplitude fluctuations in a two-band superconductor, and we show
representative results for the case of MgB2 and iron-based superconductors. In
Sec. III we introduce the formalism to study the non-linear optical response
of MgB2 in the clean or disordered case. in Sec. IV we define a proper
phenomenological model to account for disorder effect, where the couplings
between light and electronic fluctuations are computed numerically via an
exact procedure. In Sec. V we present the theoretical results and we compare
them with the existing experimental data. Sec. VI contains the final
discussion and conclusions. The Appendices contain additional material.
Appendix A presents a detailed study of the experimental optical conductivity
of MgB2, aimed at extracting realistic estimates for the ratio between
scattering rate and SC gap in each band. Appendix B contains additional
details on the normalization procedure for the light-matter coupling constants
while Appendix C explains the effect of short-range interactions on the value
of the Leggett-mode frequency.

##  II Amplitude fluctuations: from the intraband to the interband dominated
case.

As we mentioned in the introduction, the recent experimental advances in the
use of strong THz pulses to probe SC systems has been responsible for an
increasing theoretical interest on the nature of collective fluctuations of
the SC order parameter. In particular the observation of an increased non-
linear optical response below $T_{c}$ calls for a deeper understanding of the
nature of collective excitations, and their coupling to the electric field.
Indeed, as discussed in details e.g. in Ref. [34, 33], in order to understand
the experiments one has to address two separate problems: (i) the nature of
the collective excitations, encoded in their spectral function, and (ii) their
coupling to light. This is the scheme developed in detail to study more
conventional collective excitations like phonons[41, 42]: their contribution
to the non-linear response depends on their nature (oscillation frequency and
damping) and on their coupling to light, dictated by symmetry (with a
distinction among infrared or Raman-active phonons). In the case of single-
band superconductors the ongoing discussion on the contribution of the Higgs
mode to THG has been focused on the second aspect: indeed, while it is well
understood that the Higgs mode peaks at twice the SC gap in a single-band
supercondutor, its non-linear coupling to light has been shown to depend on
disorder, being negligible in the clean limit[29] and sizeable in the strong-
disorder limit[31, 32, 33]. However, as we shall see in the present section,
for a multiband superconductor also the nature of the Higgs fluctuations is
not trivial, being strongly affected by the nature of the interband coupling.
To address this issue, we introduce a two-bands ($\alpha,\beta=1,2$) BCS
hamiltonian, in close analogy with previous works [22, 21, 23, 24, 25, 26],

$H=\sum_{\alpha\mathbf{k}\sigma}\xi_{\alpha\mathbf{k}}c_{\alpha\mathbf{k}\sigma}^{\dagger}c_{\alpha\mathbf{k}\sigma}-\sum_{\alpha\beta\mathbf{q}}g_{\alpha\beta}\Psi_{\alpha\mathbf{q}}^{\dagger}\Psi_{\beta\mathbf{q}},$
(1)

where $\Psi_{\alpha\mathbf{q}}^{\dagger}$ is written explicitly as

$\Psi_{\alpha\mathbf{q}}^{\dagger}=\sum_{\mathbf{k}}c_{\alpha\mathbf{k}+\mathbf{q}/2\,\uparrow}^{\dagger}\,c_{\alpha-\mathbf{k}+\mathbf{q}/2\,\downarrow}^{\dagger}.$
(2)

The kinetic term contains
$\xi_{\alpha\mathbf{k}}=\epsilon_{\alpha\mathbf{k}}-\mu$, where
$\epsilon_{\alpha\mathbf{k}}$ is the band dispersion and $\mu$ is the chemical
potential. We consider a continuum model with free-electron like dispersions
of the form

$\epsilon_{\alpha\mathbf{k}}=s_{\alpha}\Bigl{(}\frac{k^{2}}{2m_{\alpha}}-\epsilon_{\alpha
0}\Bigr{)},$ (3)

where $\epsilon_{\alpha 0}$ is the distance of the bottom ($s_{\alpha}=1$,
electron-like band) or top ($s_{\alpha}=-1$, hole-like band) of the band from
the chemical potential in the normal state and $m_{\alpha}$ is the band mass.
Such an approximation is good both for MgB2 [43] and pnictides [44]. The
interacting term contains the symmetric matrix $g_{\alpha\beta}$ that includes
an effective $s$-wave interaction among Cooper pairs belonging to the same
band (diagonal entries) or different band (off-diagonal entries). When the
determinant $\lvert g\rvert$ of the interaction matrix is positive (negative)
the pairing is said to be intraband (interband) dominated.

To study the collective fluctuations we derive the effective action for the
bosonic pairing fields, by following a standard procedure outlined in several
manuscripts, see e.g. Ref. [16, 19, 21, 45] for the multiband case. We perform
a canonical Hubbard-Stratonovich (HS) transformation to decouple the
interaction along the superconducting channel of each band by introducing the
complex fields $\Delta_{\alpha}$. The resulting action is then quadratic in
the fermionic variables that can be integrated out, obtaining the effective
action for the HS fields. We separate then the mean-field contribution from
the fluctuating part truncated at Gaussian order in the HS fields.
Minimization of the mean-field contribution leads to the equations for the gap
amplitudes that read

$(\Pi-g^{-1})\Delta=0,$ (4)

where $\Delta=(\Delta_{1},\Delta_{2})$ can be chosen real due to the $U(1)$
symmetry of the action and being $\Pi-g^{-1}$ a symmetric matrix. The term
$\Pi_{\alpha\beta}=\Pi_{\alpha}\delta_{\alpha\beta}$ is the so-called Cooper’s
bubble reading in terms of
$E_{\alpha\mathbf{k}}=\sqrt{\xi_{\alpha\mathbf{k}}^{2}+\Delta_{\alpha}^{2}}$

$\Pi_{\alpha}=\sum_{\mathbf{k}}\frac{\tanh{(\beta
E_{\alpha\mathbf{k}}/2)}}{2E_{\alpha\mathbf{k}}}.$ (5)

The remaining action for the uniform ($\lvert\mathbf{q}\rvert=0$) gaussian
fluctuations of the mean-field gaps can be written as

$S_{\text{fluc}}\sim\operatorname{Re}{\delta\Delta}^{T}M_{A}\operatorname{Re}{\delta\Delta}+\operatorname{Im}{\delta\Delta}^{T}M_{P}\operatorname{Im}{\delta\Delta}.$
(6)

The first (second) contribution involves fluctuations $\delta\Delta$ of the
real (imaginary) part of the order parameter with respect to the mean-field
gaps. Since we chose a real mean-field gap, the first term corresponds to
amplitude fluctuations, whose propagator is given by $M_{A}^{-1}$. Being then
$\operatorname{Im}\delta{\Delta}=(\Delta_{1}\delta\theta_{1},\Delta_{2}\delta\theta_{2})$,
the second contribution describes the fluctuations of the phase
$\theta_{\alpha}$ in each band, governed by the propagator $M_{P}^{-1}$. We
notice that in this continuum model the amplitude and phase fluctuations are
decoupled at $\lvert\mathbf{q}\rvert=0$[21].

We then study the behavior of the amplitude fluctuations of a generic two-
bands superconductor as a function of the interband pairing $g_{12}$. At
gaussian level these are expressed as

$\langle\lvert\delta\Delta_{\alpha}\rvert\lvert\delta\Delta_{\beta}\rvert\rangle=(M_{A}^{-1})_{\alpha\beta},$
(7)

where the matrix $M_{A}$ reads,

$M_{A}=\begin{pmatrix}\chi_{1}^{\sigma_{1}\sigma_{1}}+2g_{11}^{-1}&2g_{12}^{-1}\\\
2g_{12}^{-1}&\chi_{2}^{\sigma_{1}\sigma_{1}}+2g_{22}^{-1}\\\ \end{pmatrix}.$
(8)

The fermionic bubble $\chi_{\alpha}^{\sigma_{1}\sigma_{1}}$ reads explicitly,
introducing the fermionic Matsubara frequencies $i\omega_{n}$,

$\chi_{\alpha}^{\sigma_{1}\sigma_{1}}(i\omega_{n})=(4\Delta_{\alpha}^{2}-(i\omega_{n})^{2})F_{\alpha}(i\omega_{n})-2\Pi_{\alpha},$
(9)

where the function $F_{\alpha}(i\omega_{n})$ is defined as

$F_{\alpha}(i\omega_{n})=\sum_{\mathbf{k}}\frac{\tanh{(\beta
E_{\alpha\mathbf{k}}/2)}}{E_{\alpha\mathbf{k}}(4E_{\alpha\mathbf{k}}^{2}-(i\omega_{n})^{2})}.$
(10)

To put in evidence the main difference between the single and the two-bands
cases, let us write the determinant of the matrix $M_{A}$

$\lvert
M_{A}\rvert=(4\Delta_{1}^{2}-\omega^{2})F_{1}(\omega)(4\Delta_{2}^{2}-\omega^{2})F_{2}(\omega)+\frac{2g_{12}}{\lvert
g\rvert}\frac{\Delta_{2}}{\Delta_{1}}(4\Delta_{2}^{2}-\omega^{2})F_{2}(\omega)+\frac{2g_{12}}{\lvert
g\rvert}\frac{\Delta_{1}}{\Delta_{2}}(4\Delta_{1}^{2}-\omega^{2})F_{1}(\omega).$
(11)

A zero of the determinant signals the presence of a well-defined mode of the
system. Since $(4\Delta_{\alpha}^{2}-\omega^{2})F_{\alpha}(\omega)=0$ at
$\omega=2\Delta_{\alpha}$, it is possible to get a zero in (11) if either the
two bands are identical or the interband coupling is zero, otherwise the
determinant is always nonzero and we only find a finite minimum at twice the
amplitude of the two gaps.

In the latter case the spectra of the amplitude fluctuations in (7) presents
two maxima and not two divergences. This is a first striking difference with
respect to the single-band case, where the Higgs-mode fluctuations emerge out
of a proper divergence of the spectral function at $2\Delta$, at least at weak
and intermediate coupling [46, 29]. Most importantly, at finite interband
coupling both $\langle\lvert\delta\Delta_{1}\rvert^{2}\rangle$ and
$\langle\lvert\delta\Delta_{2}\rvert^{2}\rangle$ present a finite spectral
weight at both $\omega=2\Delta_{1}$ and $\omega=2\Delta_{2}$. To assess their
relative weight as a function of interband coupling, we introduce the
parameter $\lvert\eta\rvert<1$ [21] such that

$\eta=\frac{g_{12}-\sqrt{g_{11}g_{22}}}{g_{12}+\sqrt{g_{11}g_{22}}}.$ (12)

For $-1<\eta<0$ ($0<\eta<1$) we are in the intraband (interband) dominated
regime. Notice that $\eta=-1$ corresponds to a purely intraband interaction,
i.e. to a diagonal coupling matrix. Instead for $\eta=1$ the product of the
diagonal terms is zero. Together with $\eta$, we fix the ratio
$\omega_{c}/\Delta_{1}$, where $\omega_{c}$ is the usual bosonic energy scale
of the pairing, the ratio $\Delta_{2}/\Delta_{1}$ of the two gap amplitudes,
and the ratio of the DOS at Fermi surface $N_{2}/N_{1}$. The entries of the
coupling matrix will then be consistently determined in order to satisfy the
gap equations (4).

We analyzed two cases in order to get an insight into the spectra of the
amplitude fluctuations in MgB2, see left panels of Fig. 2, and in a hole doped
compound of the parent pnictide BaFe2As2 (BKFA), see right panels of Fig. 2.
In the first case we use $N_{2}/N_{1}=0.75$ and $\Delta_{2}/\Delta_{1}=3$ [47,
1], where label $1$ refers to the $\pi$ band and $2$ to the $\sigma$ one. In
the second case we set $N_{2}/N_{1}=0.35$ and $\Delta_{2}/\Delta_{1}=2$ [9,
44], with label $1$ assigned to the $\beta$ band (hole like, at $\Gamma$) and
$2$ to the $\gamma$ one (electron like, at M).

We plot the values of $\langle\lvert\delta\Delta_{1}\rvert^{2}\rangle$ and
$\langle\lvert\delta\Delta_{2}\rvert^{2}\rangle$ for different values of
$\eta$. As we will see, in our phenomenological model for the non-linear
response in the disordered case, the
$\langle\lvert\delta\Delta_{\alpha}\rvert^{2}\rangle$ fluctuations associated
with the most disordered band will be the most coupled to light and thus will
dominate the Higgs sector. Therefore we avoid to report explicitly the
behavior of
$\langle\lvert\delta\Delta_{1}\rvert\lvert\delta\Delta_{2}\rvert\rangle$.

For both sets of band parameters, at $\eta=-1$ we have two well-defined
divergences at $\omega=2\Delta_{\alpha}$ for
$\langle\lvert\delta\Delta_{\alpha}\rvert^{2}\rangle$. We notice that as soon
as the interband coupling is activated, the spectral weight of
$\langle\lvert\delta\Delta_{1}\rvert^{2}\rangle$ at $\omega=2\Delta_{2}$ stays
constant and the one at $\omega=2\Delta_{1}$ gradually decreases, while the
complementary behavior is observed for
$\langle\lvert\delta\Delta_{2}\rvert^{2}\rangle$. The situation for $\eta=1$
is different in the two cases of MgB2 and pnictides. In MgB2, the spectral
weight at $\omega=2\Delta_{2}$ is largely predominant over the one at
$\omega=2\Delta_{1}$ for both channels, while in BKFA there is a comparable
weight at both $\omega=2\Delta_{1}$ and $\omega=2\Delta_{2}$.

Figure 2: The modulus of the amplitude fluctuations for the smaller (a) (c)
and bigger (b) (d) gaps at fixed gap and DOS ratio $\Delta_{2}/\Delta_{1}=3$,
$N_{2}/N_{1}=0.75$ (as appropriate for MgB2), (a) (b), and
$\Delta_{2}/\Delta_{1}=2$, $N_{2}/N_{1}=0.35$ (as appropriate for BaFe2As2),
(c) (d), as a function of interband coupling (blue tones for intraband and red
tones for interband). Higgs fluctuations in MgB2 (e) and in 122 pnictides (f)
for realistic values of the interband coupling. Values of
$\langle\lvert\delta\Delta_{1}\rvert^{2}\rangle$ are magnified $\times 4$ for
clarity.

To further clarify this behavior, we plot the ratio of the peaks at
$\omega=2\Delta_{1}$ and $\omega=2\Delta_{2}$ for both channels for the
extreme case $\eta=1$, as a function of the DOS and gap anisotropy. In panel
(a) of Fig. 3 we plot

$r_{1}=\frac{\langle\lvert\delta\Delta_{1}\rvert^{2}\rangle(\omega=2\Delta_{1})}{\langle\lvert\delta\Delta_{1}\rvert^{2}\rangle(\omega=2\Delta_{2})},$
(13)

while in panel (b) we show

$r_{2}=\frac{\langle\lvert\delta\Delta_{2}\rvert^{2}\rangle(\omega=2\Delta_{2})}{\langle\lvert\delta\Delta_{2}\rvert^{2}\rangle(\omega=2\Delta_{1})}.$
(14)

We start by observing that in the case $\eta=1$ we do not necessarily have a
completely off diagonal ($g_{11}=g_{22}=0$) coupling matrix. As can be seen by
explicitly solving the gap equations, this happens only when the gap and DOS
ratios satisfy

$\frac{N_{2}}{N_{1}}=\Bigl{(}\frac{\Delta_{1}}{\Delta_{2}}\Bigr{)}^{2}\frac{\Pi_{1}}{\Pi_{2}}\approx\Bigl{(}\frac{\Delta_{1}}{\Delta_{2}}\Bigr{)}^{2},$
(15)

where the approximation is valid in the limit
$\omega_{c}\gg\Delta_{1},\Delta_{2}$. The solution of this equation is plot as
a solid black line in Fig. 3. Changing the DOS ratio away from this condition
at fixed gap ratio, or the converse, implies to increase one among $g_{11}$ or
$g_{22}$ in order to keep the gap equation satisfied, while the other remains
zero. Thus, one order parameter can now sustain amplitude fluctuations due to
both interband and intraband scattering channels, while the other only through
interband scattering, resulting in a predominant Higgs weight at twice the
order parameter that has more scattering channels available. This is reflected
in the peak ratios $r_{1}$ and $r_{2}$, that are respectively bigger than one
roughly below ($g_{11}\neq 0$, $g_{22}=0$) and above ($g_{11}=0$, $g_{22}\neq
0$) the black line.

In summary, the naive expectation that in a multiband superconductor the
amplitude fluctuations in each band $\alpha$ have a sharp maximum at the
optical gap $2\Delta_{\alpha}$ for the same band is in general not correct. In
general, the absence of a zero in the determinant (11) makes the maxima of the
spectral function for amplitude fluctuations in a multiband system much
smoother than in a single-band case. In addition, signatures at one or the
other optical gap can be strongly suppressed by means of the interband
pairing.

In panels (e) and (f) of Fig. 2 we show the spectral function of MgB2 and 122
compounds, computed for realistic values of $\eta$. In the MgB2 case, which is
of interest in the following sections, we derive a suitable coupling matrix in
order to reproduce the behavior in temperature of the $\pi$ gap of the sample
used by Kovalev et al. [27]. We fix the gap amplitude at zero temperature at
$2\Delta_{\pi}^{0}=1.1$ THz and the critical temperature at $T_{c}=36$ K as
measured by the experiment. The Debye frequency is set at $\omega_{c}=55$ meV,
i.e. the average phonon frequency from ab initio calculations, and we take the
calculated ratio between density of states per spin at the Fermi level,
$N_{\sigma}/N_{\pi}=0.73$. We fix as well the level of interband coupling [21]
$\eta=-0.46$ desumed from the same ab initio calculations [1]. We can solve
the gap equations at $T=0$ and $T=T_{c}$ to obtain the adimensional coupling
matrix elements $\tilde{g}_{\alpha\beta}=g_{\alpha\beta}N_{\beta}$, namely
$\tilde{g}_{\pi\pi}=0.13$, $\tilde{g}_{\sigma\sigma}=0.29$,
$\tilde{g}_{\pi\sigma}=0.06$ and $\tilde{g}_{\sigma\pi}=0.08$, a similar
procedure being used in Ref. [48]. The obtained values are comparable with the
usual coupling matrices used for MgB2 [47, 1] when the renormalization due to
the Coulomb pseudopotential is taken into account in the static limit [10,
49], reading in the two-bands case
$\tilde{g}_{\alpha\beta}\to({\tilde{g}_{\alpha\beta}-\mu_{\alpha\beta})/(1+\sum_{\beta}\tilde{g}_{\alpha\beta}})$.

As one can see in Fig. 2, while for pnictides the amplitude fluctuations in
each band retain a reasonable amount of spectral weight on both optical gaps,
see panel (f), for MgB2 the Higgs fluctuations in both bands have most of the
spectral weight at the largest optical gap, i.e. at $2\Delta_{\sigma}$, see
panel (e). This result must be taken in mind while discussing the Higgs-mode
contribution in THz experiments, as we will do in the next section.

Figure 3: The values $r_{1}$, equation (13), in panel (a) and $r_{2}$,
equation (14), in panel (b) as a function of DOS anisotropy and gap ratio for
dominated interband scattering $\eta=1$. Values greater than one have been
saturated in order to better appreciate the region of anomalous behavior.
Solid black line is the solution of gap equations for $g_{11}=g_{22}=0$. Blue
dots are the position of the appropriate gap and DOS ratio for BKFA.

##  III Non-linear optical response in MgB2

Once clarified the nature of the amplitude fluctuations in a multiband
superconductor, we want now to provide an estimate of their contribution to
the non-linear optical response, and their relative importance as compared to
BCS particle-hole excitations and to the Leggett mode, whose spectrum has been
studied previously by several authors[36, 16, 21]. We will focus on the case
of MgB2, whose THz response at strong fields has been recently measured[24,
27], by providing a phenomenological estimate of the contributions of the
different superconducting fluctuations induced by the paramagnetic coupling of
light triggered by disorder.

Let us first recall the main idea behind the computation of the non-linear
optical response, as it appears both in pump-probe experiments and in
measurements of THG in transmission. As mentioned in the introduction and
detailed in several manuscripts, we need to compute the non-linear current
$j^{NL}\sim\chi^{(3)}E^{3}$. By using the effective action approach mentioned
before, one can compute [29, 45, 34, 26] $j^{NL}\sim\delta S^{(4)}/\delta
A\sim AKA^{2}$ as derivative of the fourth-order action $S^{(4)}[A]\sim
A_{i}A_{j}K_{ij;kl}A_{k}A_{l}$ expressed in terms of the gauge field
$\mathbf{A}$, that is minimally coupled to electrons in the original
Hamiltonian (1). After integrating out both the fermions and the collective
modes one obtains the various contributions to the non-linear optical kernel
$K$. The bare processes involve direct excitations of particle-hole processes
by light, and will be denoted as BCS response in what follows, while vertex
corrections represent the contribution of the collective modes, like Higgs and
Leggett fluctuations. Finally, the intensity of the THG scales as
$I_{THG}=K(2\Omega_{p})$, where $\Omega_{P}$ is the central frequency of the
THz light pulse, and the factor of two denotes the sum-frequency process
behind non-linear driving of the system. We refer for further details to
previous work [29, 21, 34, 33].

###  III.1 Clean system

For the clean case the only BCS contribution is due to diamagnetic
processes[29], originating from a term like $\sim
A^{2}s_{\alpha}e^{2}/(2m_{\alpha})c^{\dagger}_{\alpha\mathbf{k}\sigma}c_{\alpha\mathbf{k}\sigma}$
in the microscopic Hamiltonian, corresponding to the usual diamagnetic
contribution to electronic current. $K^{BCS}_{dia}$ in the two-bands model is
simply given by the sum of two single band-like expressions, reading

$K^{BCS}_{dia}=\sum_{\alpha}\Bigl{(}\frac{s_{\alpha}e^{2}}{2m_{\alpha}}\Bigr{)}^{2}\chi_{\alpha}^{\sigma_{3}\sigma_{3}},$
(16)

where
$\chi_{\alpha}^{\sigma_{3}\sigma_{3}}=4\Delta_{\alpha}^{2}F_{\alpha}(\omega)$.
The corresponding bubble is shown in panel (a) of Fig. 4. The enhancement of
THG response when $2\omega=2\Delta_{\alpha}$, that we interpret as particle-
hole excitations of BCS quasiparticles, is determined by the divergence of the
$\chi_{\alpha}^{\sigma_{3}\sigma_{3}}$ bubble, the prefactor in round brackets
in Eq. (16) being only a measure of the strength of diamagnetic coupling with
the two bands, since it includes the mediating charge $e$ and the inverse band
mass $1/m_{\alpha}$.

The Higgs contribution $K^{H}_{dia}$ in the clean two-bands case is given by a
more involved expression reading

$K^{H}_{dia}=\sum_{\alpha\beta}\Bigl{(}\frac{s_{\alpha}e^{2}}{2m_{\alpha}}\Bigr{)}\Bigl{(}\frac{s_{\beta}e^{2}}{2m_{\beta}}\Bigr{)}\chi_{\alpha}^{\sigma_{3}\sigma_{1}}\chi_{\beta}^{\sigma_{3}\sigma_{1}}(M_{A}^{-1})_{\alpha\beta}.$
(17)

We notice that if the interband coupling is zero we recover also here the sum
of two single band-like amplitude fluctuations, as shown e.g. in Ref. [45].
The fermionic bubble $\chi_{\alpha}^{\sigma_{3}\sigma_{1}}$ reads

$\chi_{\alpha}^{\sigma_{3}\sigma_{1}}=2\Delta_{\alpha}\sum_{\mathbf{k}}\frac{\xi_{\alpha\mathbf{k}}\tanh{(\beta
E_{\alpha\mathbf{k}}/2)}}{E_{\alpha\mathbf{k}}(4E_{\alpha\mathbf{k}}^{2}-(\omega+i0^{+})^{2})}.$
(18)

As it has been discussed in Ref. [29], such a bubble vanishes in the particle-
hole symmetric cases since the integral (18) scales as
$N_{\alpha}\int_{-\infty}^{\infty}d\xi G_{\alpha}(\xi)$, with
$G_{\alpha}(\xi)$ odd function of $\xi$ and $N_{\alpha}$ DOS at the Fermi
level of band $\alpha$. The coupling becomes however finite if one can account
for a finite particle-hole asymmetry of the band, by retaining a linear term
in the expression for the DOS, namely $N_{\alpha}(\xi)\approx
N_{\alpha}+\xi/(2\epsilon_{0\alpha})$, so that

$\chi_{\alpha}^{\sigma_{3}\sigma_{1}}=\frac{\Delta_{\alpha}}{4\epsilon_{0\alpha}}\Bigl{(}2g_{\alpha\alpha}^{-1}+2g_{12}^{-1}\frac{\Delta_{\bar{\alpha}}}{\Delta_{\alpha}}-(4\Delta_{\alpha}^{2}-\omega^{2})F_{\alpha}(\omega)\Bigr{)},$
(19)

with $\bar{\alpha}=2,1$ if $\alpha=1,2$. Even including such an asymmetry, for
MgB2 the $\chi^{\sigma_{3}\sigma_{1}}$ remains extremely small, since usually
$\Delta_{\alpha}/\epsilon_{0\alpha}\ll 1$,. The situation could be different
for pnictides, where a similar prefactor $\Delta_{\alpha}/\epsilon_{0\alpha}$
is indeed obtained by considering a constant DOS but accounting for the fact
that the bottom/top of the electron/hole like band lies inside the bosonic
frequency range $\sim\omega_{c}$. In that case the integration runs over an
asymmetric shell above and below the Fermi surface and Eq. (19) is finite.
However, whenever the $\chi^{\sigma_{3}\sigma_{1}}$ bubble is not zero one
should also retain the coupling of Higgs to the density/phase sector, that
again suppresses this effect [29, 45]. As we shall see below, in what follow
we will consider as main coupling of the Higgs to light the one induced by
paramagnetic processes, so the structure of Eq. (19) will only be used to
rescale such a coupling as a function of temperature.

While the Higgs-mode contribution is negligible in the clean limit, the
coupling to phase fluctuations is crucial in order to guarantee gauge-
invariance of the theory, especially for a system with an approximated
parabolic band structure as the one considered here[29, 21, 50]. For the clean
multiband case this has been done in Ref. [21, 50], leading to a contribution
to the non-linear kernel completely analogous to the Higgs case (17), see
panel (b) of Fig. 4, where now we should consider coupling with the phase
fluctuations
$\chi_{\alpha}^{\sigma_{3}\sigma_{1}}\to\chi_{\alpha}^{\sigma_{3}\sigma_{2}}$
and $(M_{A}^{-1})_{\alpha\beta}\to(M_{P}^{-1})_{\alpha\beta}$, with

$M_{P}(i\omega_{n})=\begin{pmatrix}\chi_{1}^{\sigma_{2}\sigma_{2}}(i\omega_{n})+2g_{11}^{-1}&2g_{12}^{-1}\\\
2g_{12}^{-1}&\chi_{2}^{\sigma_{2}\sigma_{2}}(i\omega_{n})+2g_{22}^{-1}\\\
\end{pmatrix}.$ (20)

The resulting expression can be rearranged as

$\tilde{K}_{dia}^{L}=-K_{dia}^{BCS}+K_{dia}^{L}.$ (21)

The first term cancels the BCS fluctuations, so that in the clean case one is
left only with the Leggett-mode contribution $K_{dia}^{L}$, which reads

$K_{dia}^{L}=\Bigl{(}\frac{s_{\sigma}e^{2}}{2m_{\sigma}}-\frac{s_{\pi}e^{2}}{2m_{\pi}}\Bigr{)}^{2}\frac{\lambda}{L},$
(22)

where $\lambda=8g_{\sigma\pi}^{-1}\Delta_{\sigma}\Delta_{\pi}$. The zeros of
the denominator

$L=\omega^{2}-\lambda\frac{\chi_{1}^{\sigma_{3}\sigma_{3}}+\chi_{2}^{\sigma_{3}\sigma_{3}}}{\chi_{1}^{\sigma_{3}\sigma_{3}}\chi_{2}^{\sigma_{3}\sigma_{3}}},$
(23)

determine the dispersion of the Leggett mode. Notice that the Leggett-mode
response only survives when the bands have different masses and/or opposite
(hole-like vs electron-like) character, so that the term in round bracket of
Eq. (22) is not zero. In contrast, when the two bands have same mass and same
character the non-linear optical kernel is proportional to the density
response of the system $K\sim\langle\rho\rho\rangle$, that must vanish at
$\lvert\mathbf{q}\rvert=0$ and finite $\omega$ because of charge
conservation[21, 50]. Since the BCS approximation by itself is not gauge
invariant this result can only be achieved by integrating out the phase mode,
as encoded indeed in the above Eq. (21). In summary, in the clean case where
only diamagnetic-like coupling to light matters, only the Leggett mode
contributes to the non-linear response, in agreement also with the
experimental observation in spontaneous Raman experiments in the symmetric
channel[36], whose response function reduces in the effective-mass
approximation to the same density-like correlation function considered here.

###  III.2 Disordered system

In presence of disorder the non-linear kernel acquires a finite contribution
also from processes mediated by paramagnetic-like matter-light coupling terms,
originating from terms like $\sim
Aev_{\alpha\mathbf{k}}c^{\dagger}_{\alpha\mathbf{k}\sigma}c_{\alpha\mathbf{k}\sigma}$
in the microscopic Hamiltonian, with $v_{\alpha\mathbf{k}}$ band velocity, as
shown by the diagrams (c) and (d) in Fig. 4. In particular, the diagram (c)
represents the BCS response, while diagram (d) represents the correction due
to amplitude or phase fluctuations. Recently, an exact numerical calculations
of all these processes in the presence of disorder has been presented for the
single-band case[33]. The relative weight of the two is mainly determined by
the level of disorder of the band, expressed as a function of the adimensional
ratio $\gamma/(2\Delta)$ of the quasiparticle scattering rate over the gap
amplitude. Since the extension of these results to the multiband case is
rather challenging, we will construct here a phenomenological model where the
frequency dependence of the various contributions is borrowed from the clean-
limit calculation, while the coupling of the various modes to light in
estimated numerically by the exact results. Such an approach is motivated by
the observation that for the disorder level of MgB2 the simulations of Ref.
[33] have proven that the overall frequency dependence of the third harmonic
current induced by both BCS and Higgs fluctuations has a rather similar
structure as in the clean case, the only relevant difference being the overall
spectral weight of the response.

Let us start from the BCS contribution. We will retain the structure in
frequency of the $\sigma_{3}\sigma_{3}$ bubble to model the spectrum of the
BCS fluctuations, and we will consider the following phenomenological
expression for the paramagnetic contributions activated by disorder

$K^{BCS}_{para}=\sum_{\alpha}w_{\alpha}^{BCS}\chi_{\alpha}^{\sigma_{3}\sigma_{3}},$
(24)

where $w_{\alpha}^{BCS}$ is a suitable frequency independent weight depending
on the level of disorder and on the dispersion of the band. In principle we
could also consider the diamagnetic contribution that is still present in a
disordered system, but it has been shown that already at small disorder level
($\gamma/(2\Delta)\approx 0.1$) it is overwhelmed by paramagnetic processes
[33]. For what concerns the Higgs modes we will consider the analogous of
(24), i.e.

$K^{H}_{para}=\sum_{\alpha\beta}w_{\alpha}^{H}w_{\beta}^{H}(M_{A}^{-1})_{\alpha\beta},$
(25)

by introducing again proper Higgs-light couplings $w_{\alpha}^{H}$. Finally,
the paramagnetic coupling of the Leggett mode with light can be rewritten in
terms of a phenomenological weight $w^{L}$, so that

$K_{para}^{L}=\frac{w^{L}}{L}.$ (26)

The overall non-linear response will then be given by

$K=K^{BCS}_{para}+K^{H}_{para}+K_{para}^{L},$ (27)

with suitable effective couplings $w_{\alpha}^{BCS},w_{\alpha}^{H}$ and
$w^{L}$ obtained by a microscopic calculation with disorder, as detailed in
the next Section.

Figure 4: The bubbles involved in our phenomenological model of THG in MgB2.
Wavy lines denote the gauge field $\mathbf{A}$, that does not count in the
expressions for $K(\omega)$ reported in the main text, while continuous lines
are fermionic Nambu-Green functions in the superconducting state. The green
vertices in panels (a) and (b) denote diamagnetic coupling of light with
fermions and contain an overall factor $e^{2}/m_{\alpha}$, while red vertices
in panels (c) and (d) denote paramagnetic processes with an overall factor
$ev_{\alpha}$. Blurred fermionic bubbles in panels (c) and (d) indicate that
in presence of disorder we should take into account the diagram resummation
associated with scattering by impurities, as done in [33], that is included in
our model through the weights $w_{\alpha}$. The blue vertices in panels (b)
and (d) represent coupling with collective modes, specifically the Higgs and
the Leggett, dashed lines.

##  IV Effective couplings to light in the presence of disorder

We determine the values of $w_{\alpha}^{BCS}$, $w_{\alpha}^{H}$ and $w^{L}$ by
considering the results by Seibold et al. [33] which treat the effect of
disorder exactly including every possible vertex and self-energy correction.
The starting point is an attractive Hubbard hamiltonian on a square lattice
with local disorder

$H=-\sum_{ij\sigma}t_{ij}c_{i\sigma}^{\dagger}c_{j\sigma}-U\sum_{i}n_{i\uparrow}n_{i\downarrow}+\sum_{i\sigma}V_{i}n_{i\sigma},$
(28)

where $n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$. The first term is the
kinetic term with next-neighbor hopping $t_{ij}=t$, the second is an
attractive interacting term and the last is a random impurity potential
$-V_{0}\leq V_{i}\leq V_{0}$ coupling to the local density, with $V_{0}$
determining the strength of disorder. The time-evolution of the mean-field
density matrix is then computed in presence of a monochromatic vector
potential $A(t)$, which allows one to calculate the third order non-linear
current $j$ averaging over a certain number of disorder configurations. We
refer to the original work [33] for further details.

To estimate the effective couplings we compute the non-linear current at zero
frequency, considering that its behavior is representative of the disorder
dependence of the non-linear response at all frequencies. This quantity is
shown in panel (a) of Fig. 5, where we report the results of Ref. [33] for the
BCS response along with its corrections due to amplitude and phase
fluctuations, as a function of $\gamma/(2\Delta)$. In Ref. [33] it was
considered a relatively large value of the SC coupling, $U/t=2$. Here, to
properly model the case of MgB2, we performed the same calculations for a
smaller SC coupling $U/t=0.8$ and the disorder level desumed from the
experimental conductivity. In general, for a given disorder level, the
absolute value of the non-linear current depends both on the band structure
and on the value of $\Delta$, being larger for smaller $U$ at weak disorder.
Nevertheless, we suggest that the ratio of the relative contributions of BCS,
Higgs and phase response depends mostly on the disorder level and not on the
SC coupling. This can be seen in panel (b) of Fig. 5, where the results at
$U/t=0.8$ and $U/t=2$ are compatible within the error bars, once renormalized
to have a similar scale. This observation makes it reasonable to assume that
the relative strength of the amplitude (phase) couplings with respect to the
BCS one $w^{H}/w^{BCS}$ ($w^{P}/w^{BCS}$) are indeed quite universal, and only
depend on the disorder level appropriate for each band,
$\gamma_{\alpha}/(2\Delta_{\alpha})$.

Figure 5: In panel (a) we report the non-linear current at zero frequency
computed by Seibold et al. in Ref. [33] for the three separated contributions:
BCS (blue line), Higgs (red line) and phase (green line). These results are
used to extract the phenomenological weights for the coupling of light to
quasiparticles ($w^{BCS}$), to the Higgs modes ($w^{H}$) and to the phase
modes ($w^{P}$) as a function of disorder. In panel (b) we make a comparison
between the zero frequency third order currents allowing only BCS fluctuations
(blue line) and all collective modes (red line) in the $U/t=2$ and $U/t=0.8$
cases. The $U/t=0.8$ values are all normalized by the same factor in order to
have a similar vertical scale.

We consider a final correction accounting for the different Fermi velocities
$v_{F\alpha}$ of the two bands of MgB2, that are of paramount importance in
determining the visibility of the contributions from the different order
parameters. This can be intuitively justified by realizing that each
paramagnetic vertex of the diagrams depicted in Fig. 4 carries a factor
$ev_{F\alpha}$. This factor is indeed foreseen in analytical approximations at
both strong and small disorder levels [31, 26] for BCS and Higgs fluctuations.
For the BCS contribution we then take

$w_{\alpha}^{BCS}=e^{4}v_{F\alpha}^{4}\tilde{w}^{BCS}\Bigl{(}\frac{\gamma_{\alpha}}{2\Delta_{\alpha}}\Bigr{)},$
(29)

where $\tilde{w}^{BCS}$ is the value of the zero frequency non linear current
reported in Fig. 5 at the disorder level $\gamma_{\alpha}/(2\Delta_{\alpha})$.
For the Higgs contribution we use instead

$w_{\alpha}^{H}=e^{2}v_{F\alpha}^{2}\zeta_{\alpha}\sqrt{\tilde{w}^{H}\Bigl{(}\frac{\gamma_{\alpha}}{2\Delta_{\alpha}}\Bigr{)}},$
(30)

the square root taking into account the presence of two fermionic bubbles in
panel (d) of Fig. 4. The $\zeta_{\alpha}$ factor is introduced to account for
the softening of the light-Higgs coupling by approaching the critical
temperature, and it corresponds to the numerator of (19), that is indeed
proportional to $\Delta_{\alpha}$. We notice that including or not the
frequency dependent part of (19) makes no detectable difference. For the
Leggett contribution, mimicking (22) we take

$w^{L}=\Bigl{[}e^{2}v_{F\sigma}^{2}\sqrt{\tilde{w}^{P}\Bigl{(}\frac{\gamma_{\sigma}}{2\Delta_{\sigma}}\Bigr{)}}-e^{2}v_{F\pi}^{2}\sqrt{\tilde{w}^{P}\Bigl{(}\frac{\gamma_{\pi}}{2\Delta_{\pi}}\Bigr{)}}\Bigr{]}^{2}\lambda,$
(31)

where also here the $\lambda$ factor, being proportional to
$\Delta_{\sigma}\Delta_{\pi}$, softens the coupling by approaching the
critical temperature. We stress that in the recent analysis of Ref. [25, 26]
impurity scattering has been introduced at the level of the Mattis-Bardeen
formalism, that corresponds to neglect vertex corrections due to disorder. In
this limit the only coupling to the phase is diamagnetic, and then subleading
in the presence of disorder. However, as previous work both in the linear[51]
and non-linear[33] regime has shown, the paramagnetic coupling of light to the
phase induced by an exact treatment of disorder is crucial at low energy. In
the present multiband case our phenomenological model takes into account also
this latter type of processes, which are of paramount importance to account
for the Leggett-mode contribution to $K$.

##  V Results

We show in Fig. 6 the results of $I_{THG}$ from a monochromatic field at
frequency $\omega$ as a function of temperature $T$, corresponding to the
modulus squared of the non-linear kernel defined before, $I_{THG}\propto\lvert
K(2\omega)\rvert^{2}$ [29]. We focus first on the separate Higgs and BCS
contributions, in order to distinguish their effects at the two optical gaps.

To highlight the importance of the disorder-mediated coupling of the modes to
light we also show in panels (a) and (b) the fictious response obtained by
using $w_{\alpha}^{H}=w_{\alpha}^{BCS}=1$, i.e. by assuming that the coupling
of light is the same for both fluctuation channels in both bands, in order to
understand the structure of their spectra in the clean limit. We stress two
crucial aspects. First, due to the different absolute values of the order
parameters, for a fixed pumping frequency $\omega$ the resonant condition
$\omega\approx\Delta_{\alpha}(T)$ will in most cases occur with only one band,
due to the thermal softening of the signal as the temperature increases and
the two gap values approach each other. This is due to the presence of the
$\sim\tanh(\beta\Delta_{\alpha}/2)$ factor in the numerator of the function
$F_{\alpha}(\omega)$, see (10), that softens the resonance in the region
$\Delta_{\alpha}/T\lesssim 2$, see dashed line in Fig. 1. Secondly, we notice
in panel (b) that at the level of interband coupling of MgB2, the spectral
weight of amplitude fluctuations is concentrated on the bigger order
parameter, with a very faint prominence of the peak in correspondence of the
$\pi$ gap.

In panels (c) and (d) we report the results using the appropriate weights due
to disorder. We stress that in order to have a realistic estimate of the ratio
$\gamma_{\alpha}/(2\Delta_{\alpha})$ in each band we performed a careful
analysis of the absorption in linear response, as detailed in Appendix D,
obtaining that $\gamma_{\sigma}=4.5$ meV and $\gamma_{\pi}=6.2$ meV. These
values should be contrasted with the values quoted e.g. in Ref. [26], where
the unrealistic value $\gamma_{\pi}=100$ meV has been used, that would be far
beyond the limit where the concept itself of quasiparticle is valid. For both
BCS and Higgs channels the visibility of the $\pi$ gap is enhanced with
respect to the $\sigma$ one because the Fermi velocity of the $\pi$ gap is
larger, the ratio $v_{F\pi}/v_{F\sigma}\approx 1.85$ being extrapolated from
ab initio calculations [52]. We notice as well that at least at the smaller
optical gap the Higgs fluctuations are subdominant with respect to the BCS
ones, see the intensity scale for panel (d). In addition, even if disorder
enhances the Higgs contribution in the lower $\pi$ band, as we noticed before
the Higgs fluctuations of the $\pi$ order parameter are mostly located at the
larger optical gap, see panel (e) of Fig. 3 and panel (b) of Fig. 6. On the
other hand BCS fluctuations in the $\pi$ bands have a maximum at
$2\Delta_{\pi}$, so that once disorder triggers the BCS paramagnetic response
in the lower $\pi$ band one recovers a larger signal at the lower optical gap,
see panel (c).

Figure 6: Comparison between the BCS (a) (c) and Higgs (b) (d) contributions
to the THG signal as a function of temperature and pump frequency. In panels
(a) and (b) the clean diamagnetic response is reported assuming equal band
mass, so that the relative visibility of the two bands in the Higgs channel is
only determined by the level of interband coupling, see Fig. 2. The prefactor
$\Delta/\epsilon_{0}$ in panel (b) is set to $10^{-3}$. In panels (c) and (d)
the dominant paramagnetic response activated by disorder is reported. The
stronger visibility of the fluctuations in $\pi$ band is determined by its
stronger disorder and larger Fermi velocity.

We then performed simulations of the THG intensity including all the
contributions of BCS, Higgs and Leggett fluctuations, whose results are shown
in Fig. 7. We selected two pump frequencies in order to capture the
composition of the THG signal in the two experimental situations of Giorgianni
et al. [24] and Kovalev et al. [27]. In the first case we are in resonance
with the Leggett mode at $T=0$, where the contribution of the relative phase
fluctuations is maximized, see the red horizontal cut in panel (a). As
expected, the leading contribution for $\Omega_{P}=1.25$ THz is by far given
by the Leggett excitation, while the relatively small bump observed at
$T\approx 0.8T_{c}$ is given by the resonance with the larger $\sigma$ order
parameter, as shown in panel (b). Interestingly, at that point the weight of
Higgs and BCS fluctuations is comparable. In the second case, the pump
frequency $\Omega_{P}=0.5$ THz corresponds to the blue horizontal cut shown in
panel (c). In this situation the resonance condition for THG occurs with the
smaller $\pi$ order parameter at temperatures where the corresponding BCS part
of the non-linear kernel still has a sizable spectral weight, while the
resonance with the larger $\sigma$ order parameter occurs for temperatures
$T\approx T_{c}$, where the non-linear kernel is already small due to the
thermal effect discussed before. In addition, at this pumping frequency the
Higgs contribution is very small since the spectral weight of amplitude
fluctuations in the $\pi$ band move to the largest optical gap, as explained
before, being far from resonance condition in the experiments of Ref. [27].
Also the Leggett gives a vanishingly small contribution, that is not
surprising: if the pumping frequency is much lower than half the $T\approx 0$
value of the Leggett mode, the resonance condition
$\omega\approx\omega_{L}(T)/2$ can only occur for $T\approx T_{c}$, where the
Leggett-mode contribution to the non-linear optical kernel is negligible. This
is also the condition studied in previous theoretical work[25], whose
conclusions about the irrelevance of the Leggett are only due to the frequency
mismatch between the pumping frequency and the Leggett mode and to the lack of
inclusion of paramagnetic coupling of light to the phase. In contrast, in the
experimental situation of Ref. [24] the set up has been designed ad hoc to
match the condition $2\Omega_{P}\approx\omega_{L}(T\approx 0)$ where the
Leggett-mode contribution is maximized.

Figure 7: Contribution to the THG signal coming from BCS (blue), Higgs
(orange) and Leggett (green) fluctuations for pump frequencies
$\Omega_{P}=1.25$ THz (b) used by Giorgianni et al. [24] and $\Omega_{P}=0.5$
THz (c) used by Kovalev et al. [27]. Even thought the $I_{THG}$ has been
computed by summing all contributions, in panels (a) and (c) we show the THG
response coming from the dominant one alone, along with the temperature
dependence of the gaps (white dashed lines) and the corresponding pumping
frequency (horizontal solid lines).

At last, to make a closer comparison with the experimental results of Ref.
[27] we account for two additional effects. First, we observe that in order to
reproduce the experimental data, reported in panel (d) of Fig. 8, the internal
screening of the applied field must be considered. This issue has been
recently discussed in Ref. [53] within the context of THG measurements in
cuprate superconductors. In these systems the $d$-wave symmetry of the order
parameter makes the sub-gap optical absorption finite even at very low
temperature, leading to a significant temperature dependence of the internal
screening. This fact affects the comparison between the experimentally
measured THG and the theoretical computation of the non-linear current, that
in principle accounts for the response to the internal field, i.e.
$j^{NL}(3\Omega_{P})\sim K(2\Omega_{P})\bar{A}_{int}^{3}$. To account for the
difference between the external field, whose spectral components one can
easily simulate, and the internal one, that accounts for the response of the
material, the authors of Ref. [53] suggested to normalize the measured THG to
the third power of the measured transmitted first harmonics. For the moment
being we use the same approximation and we compute $j^{NL}(3\Omega_{P})\sim
K(2\Omega_{P})\bar{A}_{ext}^{3}t^{3}$, where $t(\Omega_{P})=\lvert
E_{t}(\Omega_{P})\rvert/\lvert E_{ext}(\Omega_{P})\rvert$ is the temperature
dependent transmission ratio measured on the same MgB2 sample in Ref. [27].

We show in panel (b) of Fig. 8 the effect of the internal screening for the
$I_{THG}^{BCS}$ from the $\pi$ band obtained for a monocromatic incident field
at frequencies $\omega_{1}=0.3$ THz and $\omega_{3}=0.5$ THz. For the smaller
pump frequency, when the screening effect is neglected (grey dotted line), we
observe a small peak at the temperature $T^{*}\approx 0.8T_{c}$ such that
$\omega_{1}=\Delta_{\pi}(T^{*})$, followed by an increase of the THG signal at
lower temperatures. This is again due to the smallness of the
$\sim\tanh{(\beta^{*}\Delta_{\pi}/2)}$ factor at the resonance and is
ultimately linked to the two-bands nature of MgB2, for which
$\beta_{c}\Delta_{\pi}(T=0)\approx 0.73$ is much smaller than the single band
value $\beta_{c}\Delta(T=0)\approx 1.76$, for which this effect is usually not
observed. To reconcile this behavior with the experimental curve, it is then
essential to consider the action of screening, the transmission of the
material decreasing rapidly near $T_{c}$ and then being roughly constant down
to zero temperature, see grey full line also in same panel. At the larger pump
frequency, where correspondingly $T^{*}\ll T_{c}$, the effect of the screening
does not influence the shape of the peak but just contributes to renormalize
its overall intensity.

Finally, we compute the THG intensity by simulating also a realistic
narrowband multicycle pulse with form
$A(t)=A_{0}e^{-(t/\tau)^{2}}\cos{(\Omega_{P}t)}$, where the parameter $\tau$,
akin to the temporal duration, is fitted to match the temporal duration of the
experimental pulses used by Kovalev et al. [27]. The results obtained by
including all the fluctuations channel and the internal screening is reported
in panel (c) of Fig. 8. We observe an excellent agreement for the grey and
blue curves ($\omega_{1}$ and $\omega_{3}$ frequencies) while the red one
($\omega_{2}$) presents an experimentally larger peak than what we found, due
to the fact that for $\omega_{2}$ the resonance condition occurs in a range of
temperatures where the largest relative deviations occur between the
theoretical gap values and the experimental ones, see panel (a).

Figure 8: (a) The theoretical results for the temperature dependence of the
$\pi$ (dashed line) and $\sigma$ (dotted line) gaps are compared to the
experimental data (full circles) obtained by Kovalev et al. [27]. The
experimental pump frequencies $\omega_{1}=0.3$ THz (grey line),
$\omega_{2}=0.4$ THz (red line) and $\omega_{3}=0.5$ THz (blue line) are also
reported. (b) BCS contribution to the THG signal from the $\pi$ band for a
monocromatic field at the two pump frequencies $\omega_{1}$ and $\omega_{3}$
taking into account screening of the applied field (full lines) and neglecting
it (dotted lines). (c) Simulated temperature dependence of the THG signal
including all fluctuation channels and a realistic pump field. (d)
Experimental results of Ref. [27] for the same pumping frequencies. The curves
at each frequency are normalized to the corresponding maximum THG signal in
temperature for a direct comparison.

##  VI Discussion and conclusions

In the present manuscript we investigated in a qualitative and quantitative
way the non-linear optical response in multiband superconductors, showing that
THz spectroscopy represents an excellent tool to investigate the collective
modes in this category of superconductors. A first issue, addressed in Section
II, concerns the nature of the Higgs modes in a multiband superconductor. In
general, while in a single-band superconductor both the quasiparticle BCS
continuum and the Higgs mode have a sharp resonance at the optical gap, i.e.
at twice the value of the SC order parameter, in the multiband case one can
have multiple amplitude fluctuations with different intensity at the two
optical gaps. In general, the amplitude fluctuations in a given band $\alpha$
do not necessarily have their maximum at the optical gap $2\Delta_{\alpha}$ of
the same band, due to a subtle interplay between the gap anisotropy and the
interband pairing strength. For the specific case of MgB2 amplitude
fluctuations in both bands have larger spectral weight at $2\Delta_{\sigma}$,
i.e. at the larger optical gap, see panel (e) of Fig. 2. This result reflects
also on the non-linear optical kernel, where a crucial ingredient is the
correct estimate of the coupling of collective excitations to light, in
particular in the presence of disorder where paramagnetic-like processes
should be considered, as demonstrated by extensive work in the context of
single-band superconductors[31, 32, 33]. Indeed, despite the fact that the
lower band is in the relatively large disorder limit, due to the large
estimated ratio $\gamma_{\pi}/(2\Delta_{\pi})$ for the samples used in Ref.
[27], the Higgs fluctuations of the $\pi$ band have a large spectral weight at
$2\Delta_{\sigma}$, giving overall a small contribution of the Higgs at
$2\Delta_{\pi}$, see panel (d) of Fig. 6. On the other hand, the relatively
large $\gamma_{\pi}/(2\Delta_{\pi})$ ratio strongly enhances BCS fluctuations
in the $\pi$ band, whose spectral weight piles up at $2\Delta_{\pi}$ in full
analogy with the single-band case. These results fully account for the
experimental results of Ref. [27], where the system is driven with three
pumping frequencies matching the resonant condition in correspondence of the
lower gap $\Delta_{\pi}$, see panels (c) and (d) of Fig. 8. On the other hand,
under this driving conditions the Higgs and Leggett mode give a negligible
contribution to the non-linear response.

For what concerns the results of Ref. [24], in that case the experiment has
been designed with a driving frequency matching half the value of the Leggett
mode around zero temperature, namely in the temperature range where the
Leggett resonance (that is usually rather broad since it is located above the
lower optical gap) is reasonably well defined. In this range of driving
frequencies the largest contribution to the non-linear response is due to the
Leggett mode, see panels (a) and (b) of Fig. 7, confirming then the original
interpretation provided in Ref. [24]. However, it should be emphasized that in
contrast to Ref. [24] the theoretical estimate of the phase-fluctuations
contribution to the non-linear kernel has been based here on the finite
paramagnetic coupling of the light to the phase mode, triggered by disorder.
This coupling emerges in the exact treatment of disorder implemented in Ref.
[33], and confirmed here by calculations at smaller SC coupling, while it is
absent in the Mattis-Bardeen like treatment of disorder considered in previous
theoretical work[25, 26]. Including only the diamagnetic coupling to the
phase, as done in Ref. [25], surely underestimates the Leggett-mode
contribution to the non-linear optical kernel. In addition, if one analyzes
the THG for a pumping frequency matching the largest gap value,
$\omega\approx\Delta_{\sigma}(T\approx 0)$, as done in Ref. [25], it is not
surprising that the Leggett mode gives a small contribution. This is the
direct counterpart on the band $\sigma$ of the results shown in panel (d) of
Fig. 7 on the band $\pi$, where the Leggett-mode contribution is largely
subdominant. Nonetheless, we stress that the overall scale of the THG due to
the Leggett mode appears smaller than the one due to the BCS excitations, as
one can see by comparing panels (a) and (c) of Fig. 7. To test this prediction
experimentally it would be interesting to perform THG measurements on the same
sample for a wider range of pumping frequencies, in order to measure the
absolute variation of the $I_{THG}$ as different collective modes are excited
by the light pulse.

Finally, envisioning possible applications of this analysis, we would like to
comment on recent experiments in electron-doped[37] or hole-doped[39] 122
multiband iron-based superconductors. So far, the non-linear response in these
systems has been probed via time-domain protocols either with multidimensional
THz pulses[37] or via THz pump-optical probe response[39]. In the former case,
the authors observed marked oscillations of the transmitted signal as a
function of the pump-probe delay $t_{pp}$ at a frequency corresponding to
twice the smaller gap. This result[34, 35] implies a marked resonance of the
non-linear kernel at the lowest optical gap, as also suggested by THG
experiments in FeSe0.5Te0.5[38]. However, to really establish the predominance
of the Higgs over the BCS response one would need again to estimate the role
of disorder for each band, that was not explicitly provided in Ref. [38].
Disorder effects were neglected in the theoretical analysis of Ref. [37], but
the effects of the applied electromagnetic field were included within a non-
perturbative scheme, so they could access a regime not captured by the
theoretical scheme implemented in the present work. The relative importance of
the two effects is an interesting open question that deserves further
investigation. Finally, in Ref. [39] the authors measured the time-resolved
signal for a pumping frequency much smaller than both optical gaps, so that
one probes in time the square of the applied electric field weighted with
$K(\omega\approx 0)$, as theoretically expected by computing the convolution
of the kernel with the pump profile[54, 34, 35]. The results confirm the
increase of the non-linear response in the SC state, with the additional
observation of a larger $B_{1g}$ signal when superconductivity coexists with
nematic order. Also in this case, the claim of a relevance of the Higgs is not
supported by realistic theoretical calculations: indeed, the authors of Ref.
[39] do not consider the role of disorder, that is crucial not only to enhance
all the paramagnetic processes, including the BCS ones, but also to change the
symmetry of the non-linear kernel, that is the main focus of Ref. [39]. As a
consequence, one cannot expect, as claimed in Ref. [39], that the clean-limit
calculation reproduces the same polarization dependence of the dirty-limit
results. More specifically, detailed recent calculations in Ref. [35] have
shown, for the case of cuprates, that the BCS contribution due to paramagnetic
processes activated by disorder has a rather different polarization dependence
of the diamagnetic one, and both should be considered to properly account for
the measured polarization dependence of the signal. We finally observe that
the calculation itself of the Higgs-mode contribution presented in Ref. [39]
has been done assuming only intraband pairing, that is in striking contrast
with the predominance of interband pairing in these systems[5, 6, 7]. As we
explained in Sec. II, to correctly account for the amplitude-mode spectrum in
a multiband superconductor it is crucial to use a realistic pairing matrix,
and in particular the interband pairing strength has a crucial role in mixing
the amplitude fluctuations in the two bands. In addition, in the presence of
disorder nematicity has a peculiar effect in mixing the fully symmetric A1g
and the B1g channel already in the non-superconducting state[55], as evidenced
by measurements of the spontaneous Raman spectroscopy[56]. At present, a
striking outcome of Ref. [39] is that the case of pnictides evidences how THz
light pulses can access non-linear optical processes complementary to the ones
visible in ordinary Raman measurements. This is probably the most interesting
perspective, which motivates future work aimed at establishing the selection
rules enconding how light can access the different interaction mechanisms at
play in such complex systems.

###### Acknowledgements.

We acknowledge useful discussions with C. Castellani and J. Lorenzana. This
work has been supported by PRIN 2017 No. 2017Z8TS5B, and by Sapienza
University via Grant No. RM11916B56802AFE and RM120172A8CC7CC7.

##  Appendix A Study of the normal state conductivity of MgB2.

In order to determine the contribution of the various superconducting modes to
the observed THG signal, it is necessary to assess the disorder level in the
$\pi$ and $\sigma$ bands of the MgB2 13 nm thin film analyzed by Kovalev et
al. [27]. Here we focus mainly on the available experimental data for the
normal state conductivity.

The impurity scattering rates $\gamma_{\sigma}$ and $\gamma_{\pi}$ of MgB2
have been widely debated in literature. A survey of selected results obtained
by modeling linear response functions with one or two Drude terms, plus
Lorentz oscillators where included, is reported in Table 1. We considered both
transmission protocols [57, 58, 59] of THz light through various films and
reflectivity measurements [60, 61, 62, 63, 11, 64] in a wide spectral range on
single crystals, polycrystals and films. While transmission protocols are
limited to a few THz region, in the thin film approximation they give
immediate access to the complex conductivity without using Kramers-Kronig
relations. Conversely, broad spectrum reflectivity protocols require to use
the latter to derive the real and imaginary part of linear response functions,
but can be combined with ellipsometry, as done in [61, 11], to obtain more
reliable extrapolations even at low frequency. Thus they can both provide a
reasonable benchmark of the impurity scattering rates. Other experimental
techniques used to infer the impurity scattering rates are Raman spectroscopy
[12] and magnetoresistance [65, 66, 67], which usually indicate smaller
scattering rates than those coming from linear response functions.

Ref. | Type | Exp. |  $\omega_{p,\sigma}$ (eV) |  $\omega_{p,\pi}$ (eV) |  $\gamma_{\sigma}$ (meV) |  $\gamma_{\pi}$ (meV) | Fit
---|---|---|---|---|---|---|---
[60] | SC | R (300 K) | 2.6 | 4.2 | 37 | 87 | DD+L
[61] | SC | R (300 K) | $4.14^{*}$ | 4.72 | 12.4 | 85.6 | DD+L
[62] | 1000 nm F | R (50 K) | 2.54 | 1.27 | 6 | 248 | DD
[57] | 100 nm F | T (45 K) | - | 4.10 | - | 12 | SD
[58] | 50 nm F | T (40 K) | $3.74^{+0.44}_{-0.24}$ | - | $99^{+25}_{-13}$ | - | SD**
[63] | SC | R (300 K) | $4.14^{*}$ | $5.89^{*}$ | 25 | 248 | DD+L
[59] | 100 nm F | T (40 K) | 1.5 | - | 37 | - | SD**
[11] | PC | R (300 K) | 1.39 | 4.94 | 30 | 1160 | DD+L
[64] | 400 nm F | T (45 K) | $1.69\pm 0.02$ | - | $9\pm 1$ | - | SD**
Table 1: Survey of scattering rates listed in literature. The type labels are
single crystal (SC), film (F) and polycrystal (P); the experiment labels are
reflectance (R) and transmission (T) with temperature of acquisition of the
spectra in brackets; the fit labels are single Drude (SD), double Drude (DD)
and Lorentz (L) when present. Values marked with * have been fixed by the
authors of the experiment to theoretical ab initio values from [68]. Single
Drude fittings marked with ** do not attribute the observed spectral weight
specifically to the $\sigma$ or $\pi$ band, thus the results are arbitrarily
listed in the $\sigma$ column.

The spectral range analyzed by Kovalev et al. [27] is restricted between 0.5
THz and 2.5 THz, where we expect Drude-like behavior with a constant
scattering rate to be dominant. Indeed, the large values of scattering rates
obtained in most of the works listed in Table 1 arise from the attempt to use
a constant $\gamma$ in a range of frequency where inelastic effects are
relevant. Since here we need to estimate only elastic impurity effects, it is
imperative to restrict the analysis to a small range of frequency, as the one
provided by the equlibrium measurements of Ref. [27]. Due to the presence of
two electronic bands crossing the Fermi surface, in principle we expect both
of them to contribute to the normal state conductivity.

A first indication on the disorder level of the film can then be extrapolated
directly from the experimental data. We can interpolate the frequency
$\omega_{0}=1.34$ THz at which the real ($\sigma_{1}$) and imaginary
($\sigma_{2}$) part of the conductivity intersect and the corresponding value
of the real conductivity $\sigma_{1}(\omega_{0})=3.2\times
10^{5}\,\Omega^{-1}\text{cm}^{-1}$. In the case of a single Drude component
$\omega_{0}$ corresponds to the scattering rate $\gamma$ and the zero
frequency conductivity is $\sigma_{1}(0)=2\sigma_{1}(\omega_{0})$. In our case
of two Drude components we can show that irrespectively on the specific
scattering rates and plasma frequencies of the two bands we have
$\sigma_{1}(0)>2\sigma_{1}(\omega_{0})$, thus the zero frequency conductivity
of this film is at least $\sigma_{1}(0)=6.4\times
10^{5}\,\Omega^{-1}\text{cm}^{-1}$. This value is higher than other previously
reported extrapolations from linear response functions [62, 57, 63], thus
indicating an extremely clean film.

It is now imperative to distinguish the $\pi$ and $\sigma$ contributions to
the conductivity. We evaluate two different scenarios according to the ratio
$\alpha=\gamma_{\sigma}/\gamma_{\pi}$. In pristine MgB2 samples both
theoretical [68] and all the previously cited experimental studies support the
hypothesis that the scattering rate in the $\sigma$ band is smaller than the
one in the $\pi$ band, thus we can assume $\alpha<1$.

In the first scenario we consider $\gamma_{\sigma}\ll\gamma_{\pi}$, as
implicitly assumed by Kovalev et al. [27] and reported also by Hwang et al.
[62] and Kuz’menko et al. [11], the behavior in frequency is thus completely
determined by the $\sigma$ carriers as in a single band case. If we perform a
fit of the experimental conductivity with a single band Drude model, the
obtained best fit value $\omega_{p,\sigma}=5.1$ eV is notably higher than the
theoretical in-plane plasma frequency of the $\sigma$ band
$\omega_{p,\sigma}=4.14$ eV [68]. We stress that including an overall static
dielectric constant, as done in [60, 61, 63, 11], will renormalize the plasma
frequency to a still higher value, incompatible with the theoretical
prediction.

To reinforce the rebuttal of this hypothesis, we report a previous study by
Zhang et al. [69] on epitaxial MgB2 thin films of thicknesses between 6 nm and
40 nm grown by the same hybrid physical-chemical vapor deposition (HPCVD)
technique used by Kovalev et al.. The study concluded that the conductivity of
those samples is mainly determined by the $\pi$ band plasma. This result has
been established by analyzing the behavior of DC resistivity with temperature
of the various films and is compatible with the clean sample scenario
described with ab initio simulations [68]. The observed low resistivity at low
temperatures and its mild growth with temperature suggest indeed that the
scattering rates of $\sigma$ and $\pi$ bands are comparable, thus ruling out
the $\gamma_{\sigma}\ll\gamma_{\pi}$ scenario.

In the second scenario we then consider $\gamma_{\sigma}$ and $\gamma_{\pi}$
comparable. We contemporarily fit the real and imaginary part of the
conductivity with the sum of two Drude peaks with different plasma frequencies
and impurity scattering rates. Due to the very narrow frequency range accessed
by the experiment, we decide to fix the ratio
$\omega_{p,\pi}/\omega_{p,\sigma}=1.42$ according to ab initio calculations
[68] and the ratio $\alpha$ of the scattering rates. We fix the value
$\omega_{0}=1.34$ THz at which $\sigma_{1}$ and $\sigma_{2}$ intersect and we
plot the expected position of the maximum of the imaginary part as a function
of $\alpha$ in the inset of Fig. 9. We notice that for $\alpha\lesssim 0.5$
the predicted position of the maximum is incompatible with the experimental
observation, that occurs in close proximity to or even for a frequency larger
than $\omega_{0}$. We then fix $\alpha=N_{\pi}/N_{\sigma}$ following the
prescription of the Born approximation when we assume similar intraband
scattering matrix elements of the impurity potential between $\sigma$ and
$\pi$ Bloch states. The result of the fitting procedure is shown in the main
panel of Fig. 9. We obtain the impurity scattering rates $\gamma_{\sigma}=4.5$
meV and $\gamma_{\pi}=6.2$ meV and the plasma frequencies
$\omega_{p,\sigma}=3$ eV and $\omega_{p,\pi}=4.3$ eV.

We notice that both values of scattering rates are smaller than what is
reported in Table 1, with the notable exception of the result by Jin et al.
[57], that was indeed obtained by a fit in a similar frequency range. As we
mentioned before, this is however consistent with the expectation that elastic
effects only dominate at low frequencies, so that pursuing a Drude fit with a
constant scattering rate in a range of frequencies where inelastic effects
dominate leads to an unphysical enhancement of $\gamma$. For what concerns the
plasma frequencies, we obtain a smaller value than the theoretically predicted
one. The values can be reconciled by introducing a dielectric constant as
recalled before or by using an extended Drude formalism with frequency
dependent scattering rates [61]. Inclusion of inelastic effects leads in
general to a larger scattering rate with a consequent increase of the spectral
weight. This compensates the small underestimate of the plasma frequencies
that we have obtained from such a low-frequency fit.

Figure 9: Main panel: the fit (full lines) of the real (blue) and imaginary
(red) part of the experimental conductivity (full circles) from Kovalev et al.
[27]. Inset: the expected position of the maximum $\omega_{max}$ of the
imaginary part of the conductivity as a function of
$\gamma_{\sigma}/\gamma_{\pi}$ (full line) with respect to intersection
$\omega_{0}$ between the real and imaginary part of the conductivity (dashed
line).

##  Appendix B Normalizations

There is a certain degree of arbitrariness in defining the weights of our
phenomenological model of a two-bands disordered system. To clarify the
procedure that we have adopted we start from the single-band description. The
paramagnetic BCS contribution at $\omega=0$ and $T=0$ reads

$K^{BCS}_{para}=w^{BCS}\chi^{\sigma_{3}\sigma_{3}}(\omega,T=0)\equiv\tilde{w}^{BCS},$
(32)

since at $\omega=0$ and $T=0$ we want the result to coincide with the non-
linear current calculated in [33], namely $\tilde{w}^{BCS}$ in the main text.
This means that the weight should be normalized as

$w^{BCS}=\frac{\tilde{w}^{BCS}}{\chi^{\sigma_{3}\sigma_{3}}(\omega,T=0)}.$
(33)

In presence of two bands we want to generalize the previous result. Then
considering again that at $\omega=0$ and $T=0$ we have

$K^{BCS}_{para}=\sum_{\alpha}w_{\alpha}^{BCS}\chi_{\alpha}^{\sigma_{3}\sigma_{3}}(\omega,T=0),$
(34)

we mimick Eq. (33) by taking

$w_{\alpha}^{BCS}=\frac{e^{4}v_{F\alpha}^{4}\tilde{w}_{\alpha}^{BCS}}{\sum_{\beta}\chi_{\beta}^{\sigma_{3}\sigma_{3}}(\omega,T=0)},$
(35)

that corresponds to take an average weight of the $w_{\alpha}^{BCS}$ by the
$\omega,T=0$ values of the $\chi^{\sigma_{3}\sigma_{3}}$ clean bubble. A
similar procedure is adopted for $w_{\alpha}^{H}$

$w_{\alpha}^{H}=\frac{e^{2}v_{F\alpha}^{2}\zeta_{\alpha}\sqrt{\tilde{w}^{H}}}{\sum_{\beta\gamma}\zeta_{\beta}\zeta_{\gamma}(M_{A}^{-1})_{\beta\gamma}(\omega,T=0)},$
(36)

and $w^{L}$

$w^{L}=\frac{(e^{2}v_{F\sigma}^{2}\sqrt{\tilde{w}^{P}_{\sigma}}-e^{2}v_{F\pi}^{2}\sqrt{\tilde{w}^{P}_{\pi}})^{2}\lambda}{\lambda
L(\omega,T=0)}.$ (37)

##  Appendix C Short-range interaction

Plugging the calculated ab initio coupling matrix into expression (23), the
predicted Leggett-mode frequency is not compatible with both Raman [36] and
pump-probe experiments [24]. A better agreement is obtained in Ref. [36]
considering an additional repulsive short-range interaction to the hamiltonian

$H\to
H+\sum_{\alpha\beta\mathbf{q}}U_{\alpha\beta}\Phi_{\alpha\mathbf{q}}^{\dagger}\Phi_{\beta\mathbf{q}},$
(38)

where $\Phi_{\alpha\mathbf{q}}^{\dagger}$ is written explicitly as

$\Phi_{\alpha\mathbf{q}}^{\dagger}=\sum_{\mathbf{k}\sigma}c_{\alpha\mathbf{k}+\mathbf{q}/2\sigma}^{\dagger}c_{\alpha\mathbf{k}-\mathbf{q}/2\sigma},$
(39)

and $U$ a suitable two by two matrix. It is possible to decouple this
interaction in the particle-hole channel by introducing the Hubbard-
Stratonovich fields $\rho_{\alpha}$, associated with the density fluctuations
of the two bands. An analogous procedure has been adopted in previous work
[21] to treat the repulsive Coulomb interaction in the two-bands case. As it
has been already pointed out in that context, the phase and density
fluctuations are then coupled and cannot be integrated out separately. The
action (6) is then modified as follows

$S_{\text{fluc}}\sim{\lvert\delta\Delta\rvert}^{T}M_{A}\lvert\delta\Delta\rvert+\begin{pmatrix}\delta\theta&\delta\rho\end{pmatrix}M_{P,D}\begin{pmatrix}\delta\theta\\\
\delta\rho\end{pmatrix},$ (40)

having defined the vector $\delta\rho=(\delta\rho_{1},\delta\rho_{2})$. The
matrix of coupled phase and density fluctuations is a straightforward
generalization of equation (B8) of Ref. [21]

$M_{P,D}=\begin{pmatrix}(\frac{i\omega_{n}}{2})^{2}\chi_{1}-\frac{\lambda}{4}&\frac{\lambda}{4}&\frac{i\omega_{n}}{2}\chi_{1}&0\\\
\frac{\lambda}{4}&(\frac{i\omega_{n}}{2})^{2}\chi_{2}-\frac{\lambda}{4}&0&\frac{i\omega_{n}}{2}\chi_{2}\\\
-\frac{i\omega_{n}}{2}\chi_{1}&0&-\chi_{1}+2U_{11}^{-1}&2U_{12}^{-1}\\\
0&-\frac{i\omega_{n}}{2}\chi_{2}&2U_{12}^{-1}&-\chi_{2}+2U_{22}^{-1}\\\
\end{pmatrix}.$ (41)

Here we denoted $\chi_{\alpha}=\chi^{\sigma_{3}\sigma_{3}}_{\alpha}$. With
lengthy calculations it can be shown that the non-zero solution of $\lvert
M_{P,D}\rvert=0$ reads

$\omega^{2}=\lambda\frac{(\chi_{1}+\chi_{2})}{\chi_{1}\chi_{2}}-\lambda\frac{U_{11}+U_{22}-2U_{12}}{2},$
(42)

that being $\lambda<0$ in the case of MgB2 shifts the frequency of the Leggett
mode to a higher value, as already stated in Ref. [36]. There is anyway no
physical reason why one should retain $U=g$, i.e. using the same entries for
the interaction respectively in the particle-hole and particle-particle
channels. In our model we used than the last term of the previous equation as
a fitting parameter in order to match the experimental $\omega_{L}\approx 2.5$
THz.

## References

  * Liu  _et al._ [2001] A. Y. Liu, I. I. Mazin, and J. Kortus, Beyond Eliashberg superconductivity in ${\mathrm{MgB}}_{2}$: Anharmonicity, two-phonon scattering, and multiple gaps, Phys. Rev. Lett. 87, 087005 (2001).
  * Biscaras  _et al._ [2012] J. Biscaras, N. Bergeal, S. Hurand, C. Grossetête, A. Rastogi, R. C. Budhani, D. LeBoeuf, C. Proust, and J. Lesueur, Two-dimensional superconducting phase in ${\mathrm{LaTiO}}_{3}/{\mathrm{SrTiO}}_{3}$ heterostructures induced by high-mobility carrier doping, Phys. Rev. Lett. 108, 247004 (2012).
  * Singh  _et al._ [2019] G. Singh, A. Jouan, G. Herranz, M. Scigaj, F. Sánchez, L. Benfatto, S. Caprara, M. Grilli, G. Saiz, F. Couëdo, C. Feuillet-Palma, J. Lesueur, and N. Bergeal, Gap suppression at a Lifshitz transition in a multi-condensate superconductor, Nature Materials 18, 948 (2019).
  * Saito  _et al._ [2016] Y. Saito, Y. Nakamura, M. S. Bahramy, Y. Kohama, J. Ye, Y. Kasahara, Y. Nakagawa, M. Onga, M. Tokunaga, T. Nojima, Y. Yanase, and Y. Iwasa, Superconductivity protected by spin-valley locking in ion-gated ${\mathrm{MoS}}_{2}$, Nature Physics 12, 144 (2016).
  * Fernandes  _et al._ [2022] R. M. Fernandes, A. I. Coldea, H. Ding, I. R. Fisher, P. J. Hirschfeld, and G. Kotliar, Iron pnictides and chalcogenides: a new paradigm for superconductivity, Nature 601, 35 (2022).
  * Mazin  _et al._ [2008] I. I. Mazin, D. J. Singh, M. D. Johannes, and M. H. Du, Unconventional superconductivity with a sign reversal in the order parameter of ${\mathrm{LaFeAsO}}_{1-x}{\mathrm{F}}_{x}$, Phys. Rev. Lett. 101, 057003 (2008).
  * Chubukov [2012] A. Chubukov, Pairing mechanism in Fe-based superconductors, Annual Review of Condensed Matter Physics 3, 57 (2012).
  * Dolgov  _et al._ [2009] O. V. Dolgov, I. I. Mazin, D. Parker, and A. A. Golubov, Interband superconductivity: Contrasts between Bardeen-Cooper-Schrieffer and Eliashberg theories, Phys. Rev. B 79, 060502 (2009).
  * Benfatto  _et al._ [2009] L. Benfatto, E. Cappelluti, and C. Castellani, Spectroscopic and thermodynamic properties in a four-band model for pnictides, Phys. Rev. B 80, 214522 (2009).
  * Ummarino [2013] G. Ummarino, Eliashberg Theory, in  _Emergent Phenomena in Correlated Matter_ , Schriften des Forschungszentrums Jülich. Reihe modeling and simulation, Vol. 3, edited by E. Pavarini, E. Koch, and U. Schollwöck (Forschungszentrum Jülich GmbH Zentralbibliothek, Verlag, Jülich, 2013).
  * Kuz’menko  _et al._ [2002] A. Kuz’menko, F. Mena, H. Molegraaf, D. van der Marel, B. Gorshunov, M. Dressel, I. Mazin, J. Kortus, O. Dolgov, T. Muranaka, and J. Akimitsu, Manifestation of multiband optical properties of ${\mathrm{MgB}}_{2}$, Solid State Communications 121, 479 (2002).
  * Quilty  _et al._ [2003] J. W. Quilty, S. Lee, S. Tajima, and A. Yamanaka, $c$-axis Raman scattering spectra of ${\mathrm{M}\mathrm{g}\mathrm{B}}_{2}$: Observation of a dirty-limit gap in the $\pi$ bands, Phys. Rev. Lett. 90, 207006 (2003).
  * Tu  _et al._ [2010] J. J. Tu, J. Li, W. Liu, A. Punnoose, Y. Gong, Y. H. Ren, L. J. Li, G. H. Cao, Z. A. Xu, and C. C. Homes, Optical properties of the iron arsenic superconductor ${\mathrm{BaFe}}_{1.85}{\mathrm{Co}}_{0.15}{\mathrm{As}}_{2}$, Phys. Rev. B 82, 174509 (2010).
  * Charnukha  _et al._ [2011] A. Charnukha, O. V. Dolgov, A. A. Golubov, Y. Matiks, D. L. Sun, C. T. Lin, B. Keimer, and A. V. Boris, Eliashberg approach to infrared anomalies induced by the superconducting state of ${\mathrm{Ba}}_{0.68}{\mathrm{K}}_{0.32}{\mathrm{Fe}}_{2}{\mathrm{As}}_{2}$ single crystals, Phys. Rev. B 84, 174511 (2011).
  * Leggett [1966] A. J. Leggett, Number-phase fluctuations in two-band superconductors, Progress of Theoretical Physics 36, 901 (1966).
  * Sharapov  _et al._ [2002] S. G. Sharapov, V. P. Gusynin, and H. Beck, Effective action approach to the Leggett’s mode in two-band superconductors, The European Physical Journal B - Condensed Matter and Complex Systems 30, 45 (2002).
  * Ota  _et al._ [2011] Y. Ota, M. Machida, T. Koyama, and H. Aoki, Collective modes in multiband superfluids and superconductors: Multiple dynamical classes, Phys. Rev. B 83, 060507 (2011).
  * Carlström  _et al._ [2011] J. Carlström, J. Garaud, and E. Babaev, Length scales, collective modes, and type-1.5 regimes in three-band superconductors, Phys. Rev. B 84, 134518 (2011).
  * Marciani  _et al._ [2013] M. Marciani, L. Fanfarillo, C. Castellani, and L. Benfatto, Leggett modes in iron-based superconductors as a probe of time-reversal symmetry breaking, Phys. Rev. B 88, 214508 (2013).
  * Maiti and Chubukov [2013] S. Maiti and A. V. Chubukov, $s+is$ state with broken time-reversal symmetry in Fe-based superconductors, Phys. Rev. B 87, 144511 (2013).
  * Cea and Benfatto [2016] T. Cea and L. Benfatto, Signature of the Leggett mode in the ${A}_{1g}$ Raman response: From ${\text{MgB}}_{2}$ to iron-based superconductors, Phys. Rev. B 94, 064512 (2016).
  * Krull  _et al._ [2016] H. Krull, N. Bittner, G. S. Uhrig, D. Manske, and A. P. Schnyder, Coupling of Higgs and Leggett modes in non-equilibrium superconductors, Nature Communications 7, 11921 (2016).
  * Murotani  _et al._ [2017] Y. Murotani, N. Tsuji, and H. Aoki, Theory of light-induced resonances with collective Higgs and Leggett modes in multiband superconductors, Phys. Rev. B 95, 104503 (2017).
  * Giorgianni  _et al._ [2019] F. Giorgianni, T. Cea, C. Vicario, C. P. Hauri, W. K. Withanage, X. Xi, and L. Benfatto, Leggett mode controlled by light pulses, Nature Physics 15, 341 (2019).
  * Murotani and Shimano [2019] Y. Murotani and R. Shimano, Nonlinear optical response of collective modes in multiband superconductors assisted by nonmagnetic impurities, Phys. Rev. B 99, 224510 (2019).
  * Haenel  _et al._ [2021] R. Haenel, P. Froese, D. Manske, and L. Schwarz, Time-resolved optical conductivity and Higgs oscillations in two-band dirty superconductors, Phys. Rev. B 104, 134504 (2021).
  * Kovalev  _et al._ [2021] S. Kovalev, T. Dong, L.-Y. Shi, C. Reinhoffer, T.-Q. Xu, H.-Z. Wang, Y. Wang, Z.-Z. Gan, S. Germanskiy, J.-C. Deinert, I. Ilyakov, P. H. M. van Loosdrecht, D. Wu, N.-L. Wang, J. Demsar, and Z. Wang, Band-selective third-harmonic generation in superconducting ${\mathrm{MgB}}_{2}$: Possible evidence for the Higgs amplitude mode in the dirty limit, Phys. Rev. B 104, L140505 (2021).
  * Tsuji and Aoki [2015] N. Tsuji and H. Aoki, Theory of Anderson pseudospin resonance with Higgs mode in superconductors, Phys. Rev. B 92, 064508 (2015).
  * Cea  _et al._ [2016] T. Cea, C. Castellani, and L. Benfatto, Nonlinear optical effects and third-harmonic generation in superconductors: Cooper pairs versus Higgs mode contribution, Phys. Rev. B 93, 180507 (2016).
  * Tsuji  _et al._ [2016] N. Tsuji, Y. Murakami, and H. Aoki, Nonlinear light-Higgs coupling in superconductors beyond BCS: Effects of the retarded phonon-mediated interaction, Phys. Rev. B 94, 224519 (2016).
  * Silaev [2019] M. Silaev, Nonlinear electromagnetic response and Higgs-mode excitation in BCS superconductors with impurities, Phys. Rev. B 99, 224511 (2019).
  * Tsuji and Nomura [2020] N. Tsuji and Y. Nomura, Higgs-mode resonance in third harmonic generation in NbN superconductors: Multiband electron-phonon coupling, impurity scattering, and polarization-angle dependence, Phys. Rev. Research 2, 043029 (2020).
  * Seibold  _et al._ [2021] G. Seibold, M. Udina, C. Castellani, and L. Benfatto, Third harmonic generation from collective modes in disordered superconductors, Phys. Rev. B 103, 014512 (2021).
  * Udina  _et al._ [2019] M. Udina, T. Cea, and L. Benfatto, Theory of coherent-oscillations generation in terahertz pump-probe spectroscopy: From phonons to electronic collective modes, Phys. Rev. B 100, 165131 (2019).
  * Udina  _et al._ [2022] M. Udina, J. Fiore, T. Cea, C. Castellani, G. Seibold, and L. Benfatto, Thz non-linear optical response in cuprates: predominance of the BCS response over the Higgs mode, Faraday Discuss. , (2022).
  * Blumberg  _et al._ [2007] G. Blumberg, A. Mialitsin, B. S. Dennis, M. V. Klein, N. D. Zhigadlo, and J. Karpinski, Observation of Leggett’s collective mode in a multiband ${\mathrm{MgB}}_{2}$ superconductor, Phys. Rev. Lett. 99, 227002 (2007).
  * Vaswani  _et al._ [2021] C. Vaswani, J. H. Kang, M. Mootz, L. Luo, X. Yang, C. Sundahl, D. Cheng, C. Huang, R. H. J. Kim, Z. Liu, Y. G. Collantes, E. E. Hellstrom, I. E. Perakis, C. B. Eom, and J. Wang, Light quantum control of persisting Higgs modes in iron-based superconductors, Nature Communications 12, 258 (2021).
  * Isoyama  _et al._ [2021] K. Isoyama, N. Yoshikawa, K. Katsumi, J. Wong, N. Shikama, Y. Sakishita, F. Nabeshima, A. Maeda, and R. Shimano, Light-induced enhancement of superconductivity in iron-based superconductor ${\mathrm{Fe}}{\mathrm{Se}}_{0.5}{\mathrm{Te}}_{0.5}$, Communications Physics 4, 160 (2021).
  * Grasset  _et al._ [2022] R. Grasset, K. Katsumi, P. Massat, H.-H. Wen, X.-H. Chen, Y. Gallais, and R. Shimano, Terahertz pulse-driven collective mode in the nematic superconducting state of ${\mathrm{Ba}}_{1-x}{\mathrm{K}}_{x}{\mathrm{Fe}}_{2}{\mathrm{As}}_{2}$, npj Quantum Materials 7, 4 (2022).
  * Chu  _et al._ [2020a] H. Chu, M.-J. Kim, K. Katsumi, S. Kovalev, R. D. Dawson, L. Schwarz, N. Yoshikawa, G. Kim, D. Putzky, Z. Z. Li, H. Raffy, S. Germanskiy, J.-C. Deinert, N. Awari, I. Ilyakov, B. Green, M. Chen, M. Bawatna, G. Cristiani, G. Logvenov, Y. Gallais, A. V. Boris, B. Keimer, A. P. Schnyder, D. Manske, M. Gensch, Z. Wang, R. Shimano, and S. Kaiser, Phase-resolved higgs response in superconducting cuprates, Nature Communications 11, 1793 (2020a).
  * Merlin [1997] R. Merlin, Generating coherent THz phonons with light pulses, Solid State Communications 102, 207 (1997), highlights in Condensed Matter Physics and Materials Science.
  * Juraschek and Maehrlein [2018] D. M. Juraschek and S. F. Maehrlein, Sum-frequency ionic Raman scattering, Phys. Rev. B 97, 174302 (2018).
  * Kortus  _et al._ [2001] J. Kortus, I. I. Mazin, K. D. Belashchenko, V. P. Antropov, and L. L. Boyer, Superconductivity of metallic boron in ${\mathrm{MgB}}_{2}$, Phys. Rev. Lett. 86, 4656 (2001).
  * Mazin and Schmalian [2009] I. Mazin and J. Schmalian, Pairing symmetry and pairing state in ferropnictides: Theoretical overview, Physica C: Superconductivity 469, 614 (2009).
  * Cea  _et al._ [2018] T. Cea, P. Barone, C. Castellani, and L. Benfatto, Polarization dependence of the third-harmonic generation in multiband superconductors, Phys. Rev. B 97, 094516 (2018).
  * Cea  _et al._ [2015] T. Cea, C. Castellani, G. Seibold, and L. Benfatto, Nonrelativistic dynamics of the amplitude (Higgs) mode in superconductors, Phys. Rev. Lett. 115, 157002 (2015).
  * Golubov  _et al._ [2002] A. A. Golubov, J. Kortus, O. V. Dolgov, O. Jepsen, Y. Kong, O. K. Andersen, B. J. Gibson, K. Ahn, and R. K. Kremer, Specific heat of $\mathrm{Mg}{\mathrm{B}}_{2}$ in a one- and a two-band model from first-principles calculations, Journal of Physics: Condensed Matter 14, 1353 (2002).
  * Jin  _et al._ [2003] B. B. Jin, T. Dahm, A. I. Gubin, E.-M. Choi, H. J. Kim, S.-I. Lee, W. N. Kang, and N. Klein, Anomalous coherence peak in the microwave conductivity of $c$-axis oriented ${\mathrm{M}\mathrm{g}\mathrm{B}}_{2}$ thin films, Phys. Rev. Lett. 91, 127006 (2003).
  * Ummarino  _et al._ [2004] G. Ummarino, R. Gonnelli, S. Massidda, and A. Bianconi, Two-band Eliashberg equations and the experimental $T_{\mathrm{c}}$ of the diboride ${\mathrm{Mg}}_{1-x}{\mathrm{Al}}_{x}{\mathrm{B}}_{2}$, Physica C: Superconductivity 407, 121 (2004).
  * Maiti  _et al._ [2017] S. Maiti, A. V. Chubukov, and P. J. Hirschfeld, Conservation laws, vertex corrections, and screening in Raman spectroscopy, Phys. Rev. B 96, 014503 (2017).
  * Seibold  _et al._ [2017] G. Seibold, L. Benfatto, and C. Castellani, Application of the Mattis-Bardeen theory in strongly disordered superconductors, Phys. Rev. B 96, 144507 (2017).
  * Brinkman  _et al._ [2002] A. Brinkman, A. A. Golubov, H. Rogalla, O. V. Dolgov, J. Kortus, Y. Kong, O. Jepsen, and O. K. Andersen, Multiband model for tunneling in ${\mathrm{MgB}}_{2}$ junctions, Phys. Rev. B 65, 180517 (2002).
  * Chu  _et al._ [2020b] H. Chu, M.-J. Kim, K. Katsumi, S. Kovalev, R. D. Dawson, L. Schwarz, N. Yoshikawa, G. Kim, D. Putzky, Z. Z. Li, H. Raffy, S. Germanskiy, J.-C. Deinert, N. Awari, I. Ilyakov, B. Green, M. Chen, M. Bawatna, G. Cristiani, G. Logvenov, Y. Gallais, A. V. Boris, B. Keimer, A. P. Schnyder, D. Manske, M. Gensch, Z. Wang, R. Shimano, and S. Kaiser, Phase-resolved Higgs response in superconducting cuprates, Nature Communications 11, 1793 (2020b).
  * Katsumi  _et al._ [2018] K. Katsumi, N. Tsuji, Y. I. Hamada, R. Matsunaga, J. Schneeloch, R. D. Zhong, G. D. Gu, H. Aoki, Y. Gallais, and R. Shimano, Higgs mode in the $d$-wave superconductor $\mathrm{Bi_{2}Sr_{2}{CaCu}_{2}O_{8+x}}$ driven by an intense terahertz pulse, Phys. Rev. Lett. 120, 117001 (2018).
  * Udina  _et al._ [2020] M. Udina, M. Grilli, L. Benfatto, and A. V. Chubukov, Raman response in the nematic phase of FeSe, Phys. Rev. Lett. 124, 197602 (2020).
  * Zhang  _et al._ [2021] W. Zhang, S. Wu, S. Kasahara, T. Shibauchi, Y. Matsuda, and G. Blumberg, Quadrupolar charge dynamics in the nonmagnetic ${\mathrm{Fe}}{\mathrm{Se}}_{1-x}{\mathrm{S}}_{x}$ superconductors, Proceedings of the National Academy of Sciences 118, e2020585118 (2021).
  * Jin  _et al._ [2005] B. B. Jin, P. Kuzel, F. Kadlec, T. Dahm, J. M. Redwing, A. V. Pogrebnyakov, X. X. Xi, and N. Klein, Terahertz surface impedance of epitaxial ${\mathrm{M}\mathrm{g}\mathrm{B}}_{2}$ thin film, Applied Physics Letters 87, 092503 (2005).
  * Jung  _et al._ [2002] J. H. Jung, K. W. Kim, H. J. Lee, M. W. Kim, T. W. Noh, W. N. Kang, H.-J. Kim, E.-M. Choi, C. U. Jung, and S.-I. Lee, Far-infrared transmission studies of a c-axis-oriented superconducting ${\mathrm{MgB}}_{2}$ thin film, Phys. Rev. B 65, 052413 (2002).
  * Kaindl  _et al._ [2001] R. A. Kaindl, M. A. Carnahan, J. Orenstein, D. S. Chemla, H. M. Christen, H.-Y. Zhai, M. Paranthaman, and D. H. Lowndes, Far-infrared optical conductivity gap in superconducting ${\mathrm{MgB}}_{2}$ films, Phys. Rev. Lett. 88, 027003 (2001).
  * Di Castro  _et al._ [2006] D. Di Castro, M. Ortolani, E. Cappelluti, U. Schade, N. D. Zhigadlo, and J. Karpinski, Infrared properties of ${\mathrm{Mg}}_{1-x}{\mathrm{Al}}_{x}{({\mathrm{B}}_{1-y}{\mathrm{C}}_{y})}_{2}$ single crystals in the normal and superconducting state, Phys. Rev. B 73, 174509 (2006).
  * Guritanu  _et al._ [2006] V. Guritanu, A. B. Kuzmenko, D. van der Marel, S. M. Kazakov, N. D. Zhigadlo, and J. Karpinski, Anisotropic optical conductivity and two colors of ${\mathrm{MgB}}_{2}$, Phys. Rev. B 73, 104509 (2006).
  * Seo  _et al._ [2017] Y.-S. Seo, J. H. Lee, W. N. Kang, and J. Hwang, Revisiting optical properties of ${\mathrm{MgB}}_{2}$ with a high-quality sample prepared by a HPCVD method, Scientific Reports 7, 8977 (2017).
  * Kakeshita  _et al._ [2006] T. Kakeshita, S. Lee, and S. Tajima, Anisotropic Drude response and the effect of anisotropic C substitution in $\mathrm{Mg}({\mathrm{B}}_{1-x}{\mathrm{C}}_{x}{)}_{2}$, Phys. Rev. Lett. 97, 037002 (2006).
  * Tu  _et al._ [2001] J. J. Tu, G. L. Carr, V. Perebeinos, C. C. Homes, M. Strongin, P. B. Allen, W. N. Kang, E.-M. Choi, H.-J. Kim, and S.-I. Lee, Optical properties of c-axis oriented superconducting ${\mathrm{MgB}}_{2}$ films, Phys. Rev. Lett. 87, 277001 (2001).
  * Pallecchi  _et al._ [2005] I. Pallecchi, V. Ferrando, E. G. D’Agliano, D. Marré, M. Monni, M. Putti, C. Tarantini, F. Gatti, H. U. Aebersold, E. Lehmann, X. X. Xi, E. G. Haanappel, and C. Ferdeghini, Magnetoresistivity as a probe of disorder in the $\pi$ and $\sigma$ bands of $\mathrm{Mg}{\mathrm{B}}_{2}$, Phys. Rev. B 72, 184512 (2005).
  * Pallecchi  _et al._ [2006] I. Pallecchi, V. Ferrando, E. G. D’Agliano, D. Marré, M. Monni, M. Putti, C. Tarantini, F. Gatti, H. U. Aebersold, E. Lehmann, X. X. Xi, E. G. Haanappel, and C. Ferdeghini, Erratum: Magnetoresistivity as a probe of disorder in the $\pi$ and $\sigma$ bands of ${\mathrm{MgB}}_{2}$ [Phys. Rev. B 72, 184512 (2005)], Phys. Rev. B 73, 029901 (2006).
  * Yang  _et al._ [2008] H. Yang, Y. Liu, C. Zhuang, J. Shi, Y. Yao, S. Massidda, M. Monni, Y. Jia, X. Xi, Q. Li, Z.-K. Liu, Q. Feng, and H.-H. Wen, Fully band-resolved scattering rate in ${\mathrm{MgB}}_{2}$ revealed by the nonlinear hall effect and magnetoresistance measurements, Phys. Rev. Lett. 101, 067001 (2008).
  * Mazin  _et al._ [2002] I. I. Mazin, O. K. Andersen, O. Jepsen, O. V. Dolgov, J. Kortus, A. A. Golubov, A. B. Kuz’menko, and D. van der Marel, Superconductivity in ${\mathrm{M}\mathrm{g}\mathrm{B}}_{2}$: Clean or dirty?, Phys. Rev. Lett. 89, 107002 (2002).
  * Zhang  _et al._ [2013] C. Zhang, Y. Wang, D. Wang, Y. Zhang, Z.-H. Liu, Q.-R. Feng, and Z.-Z. Gan, Suppression of superconductivity in epitaxial ${\mathrm{M}\mathrm{g}\mathrm{B}}_{2}$ ultrathin films, Journal of Applied Physics 114, 023903 (2013).