Nonlinear optical response of collective modes in multiband superconductors assisted by nonmagnetic impurities

TL;DR

This paper theoretically studies how nonmagnetic impurities enhance the nonlinear optical response of Higgs modes in multiband superconductors, showing a significant effect on Higgs but not Leggett modes, with implications for third-harmonic generation.

Contribution

It demonstrates that nonmagnetic impurities selectively enhance the nonlinear optical response of Higgs modes in multiband superconductors, using an extended density matrix approach.

Findings

Nonmagnetic impurities significantly enhance Higgs mode nonlinear response.

Leggett mode response remains unaffected by impurity scattering.

Quasiparticle excitations contribute subdominantly to third-harmonic generation.

Abstract

In multiband superconductors, multiple collective modes exist associated with the multiple order parameters. Oscillations of the amplitude and the relative phase of the order parameters are called Higgs and Leggett modes, respectively. Recently, it has been suggested that nonmagnetic impurity scattering would enhance nonlinear coupling between the Higgs mode and an electromagnetic wave with a frequency located in the superconducting gap region, while its effect on the Leggett mode is still unresolved. Here, we theoretically investigated the nonlinear optical response of multiband Bardeen-Cooper-Schrieffer-type superconductors in the presence of nonmagnetic impurities with a density matrix approach extending the Mattis-Bardeen model of linear response. We found that the drastic enhancement of nonlinear optical response due to the nonmagnetic impurity scattering occurs only for the Higgs modes and not for the Leggett mode. As a result, both the light-induced dynamics of the superconducting gaps and the resulting third-harmonic generation are dominated by the Higgs modes. We also examined the role of quasiparticle excitations to find that they give the subdominant contribution to the third-harmonic generation.

Tables

Table 1. Table 1: Order of the light-induced collective modes, in the unit of the ponderomotive energy e 2 ​ 𝐀 2 / 2 ​ m superscript 𝑒 2 superscript 𝐀 2 2 𝑚 e^{2}\mathbf{A}^{2}/2m . Arrows connects the clean limit ( γ ≪ 2 ​ Δ much-less-than 𝛾 2 Δ \gamma\ll 2\Delta ) on the left side and the dirty limit ( γ ≫ 2 ​ Δ much-greater-than 𝛾 2 Δ \gamma\gg 2\Delta ) on the right side. For a reference, square brackets show the case when the energy dispersion relation exhibits a nonparabolicity, calculated in the clean limit Tsuji2015 . The asterisked part is neglected in the main text.

Channel	Para ($𝐩\cdot 𝐀$)	Dia (${𝐀}^2$)
HM ($\delta {\Delta }^{\prime }$)	${ϵ}_F\gamma /{\Delta }^2\to {ϵ}_F/\gamma$	0 [$\Delta /{ϵ}_F$]
LM ($\delta {\Delta }^{\prime \prime }$)	$\gamma /\Delta \to \Delta /\gamma (\ast )$	1

Table 2. Table 2: Order of the third-order current 𝐣 ​ ( t ) | 3 evaluated-at 𝐣 𝑡 3 \mathbf{j}(t)|_{3} in the unit of ( e 4 ​ n / m 2 ​ ϵ F ) ​ A 3 superscript 𝑒 4 𝑛 superscript 𝑚 2 subscript italic-ϵ F superscript 𝐴 3 (e^{4}n/m^{2}\epsilon_{\mathrm{F}})A^{3} . Arrows connects the clean limit ( γ ≪ 2 ​ Δ much-less-than 𝛾 2 Δ \gamma\ll 2\Delta ) on the left side and the dirty limit ( γ ≫ 2 ​ Δ much-greater-than 𝛾 2 Δ \gamma\gg 2\Delta ) on the right side. For reference, square brackets show the case when the energy dispersion relation exhibits a nonparabolicity, calculated in the clean limit Tsuji2015 ; Murotani2017 . The asterisked part is neglected in the main text. The combination of QP and Dia is specially called charge-density fluctuations.

Channel	Para ($𝐩\cdot 𝐀$)	Dia (${𝐀}^2$)
HM	${\left({ϵ}_F\gamma /{\Delta }^2\right)}^2\to {\left({ϵ}_F/\gamma \right)}^2$	0 [${\left(\Delta /{ϵ}_F\right)}^2$]
QP	${\left({ϵ}_F\gamma /{\Delta }^2\right)}^2\to {\left({ϵ}_F/\gamma \right)}^2$	0 [1]
LM	${\left(\gamma /\Delta \right)}^2\to {\left(\Delta /\gamma \right)}^2(\ast )$	1

Equations

$${\mathcal{H}}_0=\sum _{i\mathbf{k}\sigma }{ϵ}_{i\mathbf{k}}{c}_{i\mathbf{k}\sigma }^{†}{c}_{i\mathbf{k}\sigma },$$

$${ϵ}_{i\mathbf{k}}=s_i\left(\frac{{\mathbf{k}}^2}{2m_i}-{ϵ}_{\mathrm{F}i}\right),s_i=\{\begin{array}{ll}+1& \mathrm{electron}\mathrm{band}\\ -1& \mathrm{hole}\mathrm{band}\end{array}$$

$${\mathcal{H}}_{\mathrm{int}}$$

$${\mathcal{H}}_1=-\sum _{i\mathbf{k}{\mathbf{k}}^{\prime }\sigma }\left({\mathbf{J}}_{i\mathbf{k}{\mathbf{k}}^{\prime }}\cdot \text{A}\right)c_{i\mathbf{k}\sigma }^{†}{c}_{i{\mathbf{k}}^{\prime }\sigma }.$$

$${\mathcal{H}}_2$$

$$\mathcal{H}={\mathcal{H}}_0+{\mathcal{H}}_1+{\mathcal{H}}_2+{\mathcal{H}}_{\mathrm{int}}.$$

$${⟨{\left|\mathbf{e}\cdot {\mathbf{J}}_{i\mathbf{k}{\mathbf{k}}^{\prime }}\right|}^2⟩}_{\mathrm{Av}}$$

$$\simeq \frac{(e{v}_{\mathrm{F}i}{)}^2}{3N_i(0)}{W}_i({ϵ}_{i\mathbf{k}},{ϵ}_{i{\mathbf{k}}^{\prime }}),$$

$$W_i(ϵ,{ϵ}^{\prime })$$

$$\gamma \ll {ϵ}_{\mathrm{F}}.$$

$${\Delta }_i$$

$$2{\Delta }^{\mathrm{eq}}\ll {ϵ}_{\mathrm{F}}.$$

$$\{\begin{array}{ll}\gamma \ll 2{\Delta }^{\mathrm{eq}}& \mathrm{clean}\mathrm{limit},\\ \gamma \gg 2{\Delta }^{\mathrm{eq}}& \mathrm{dirty}\mathrm{limit}.\end{array}$$

$${\rho }_{i\mathbf{k}{\mathbf{k}}^{\prime }}^{ab}=⟨{\psi }_{i\mathbf{k}}^{a†}{\psi }_{i{\mathbf{k}}^{\prime }}^b⟩(a,b=1,2),$$

$${\rho }_{i\mathbf{k}{\mathbf{k}}^{\prime }}^{ab}$$

$$=-⟨\frac{\delta \mathcal{H}}{\delta \mathbf{A}}⟩={\text{j}}_{\mathrm{P}}+{\text{j}}_{\mathrm{D}},$$

$$\displaystyle\mathbf{j}_{\mathrm{P}}(t)\Big{|}_{1}=\mathbf{e}\sum_{i}\frac{e^{2}n_{i}}{m_{i}}\iint d\epsilon~{}d\epsilon^{\prime}~{}W_{i}(\epsilon,\epsilon^{\prime})$$

$$\times \left[l_i(ϵ,{ϵ}^{\prime }{)}^2\mathrm{Re}{F}_i^{11}(ϵ,{ϵ}^{\prime })+p_i(ϵ,{ϵ}^{\prime }{)}^2\mathrm{Re}{F}_i^{21}(ϵ,{ϵ}^{\prime })\right],$$

$$\left[i\frac{\partial }{\partial t}-(E^{\prime }-E)\right]F_i^{11}(ϵ,{ϵ}^{\prime })$$

$$\left[i\frac{\partial }{\partial t}-(E^{\prime }+E)\right]F_i^{21}(ϵ,{ϵ}^{\prime })$$

$$F_i^{22}(ϵ,{ϵ}^{\prime })=F_i^{11}(ϵ,{ϵ}^{\prime }{)}^{\ast },F_i^{12}$$

$$\displaystyle\mathbf{j}_{\mathrm{D}}(t)\Big{|}_{1}$$

$$\sigma (\omega )$$

$$A(t)=A_0\exp \left(-\frac{t^2}{2{\tau }^2}\right)\cos \Omega t,$$

$${\sigma }_{\mathrm{Drude}}(\omega )=\sum _i\frac{i{e}^2{n}_i}{m_i(\omega +i{\gamma }_i)},$$

$$\displaystyle\mathbf{j}(t)\Big{|}_{2}=0.$$

$$\displaystyle\delta\Delta_{i}^{\prime}\Big{|}_{2}$$

$$+\frac{1}{2}[u_j(ϵ{)}^2-v_j(ϵ{)}^2][r_j^{21}(ϵ)+r_j^{12}(ϵ)]\},$$

$$\displaystyle\delta\Delta_{i}^{\prime\prime}\Big{|}_{2}$$

$$\displaystyle r_{i}^{ab}(\epsilon_{i\mathbf{k}})=\int\frac{d\Omega_{\mathbf{k}}}{4\pi}\rho_{i\mathbf{k}\mathbf{k}}^{ab}\Big{|}_{2}.$$

$$i\frac{\partial }{\partial t}{r}_i^{11}(ϵ)$$

$$r_i^{22}(ϵ)$$

$$S_i^{11}(ϵ)$$

$$-p_i(ϵ,{ϵ}^{\prime }{)}^2\mathrm{Im}{F}_i^{21}(ϵ,{ϵ}^{\prime })]W_i(ϵ,{ϵ}^{\prime }).$$

$$r_i^{11}(ϵ)=r_i^{11}(-ϵ).$$

Full text

I Introduction

Collective modes in superconductors have recently gained great interest owing to the development of nonlinear terahertz spectroscopy Matsunaga2017-1 . Because a superconducting order parameter is a complex quantity in general, one superconducting degree of freedom accommodates one amplitude mode Volkov1973 and one phase mode Nambu1960 ; Goldstone1961 ; Goldstone1962 . In superconductors, the Anderson-Higgs mechanism elevates the energy of the phase mode to the plasma frequency far larger than the superconducting gap Anderson1963 ; Higgs1964-1 ; Englert1964 ; Higgs1964-2 ; Guralnik1964 . As a result, the amplitude mode is left stable in the low-energy region and is specially called a Higgs mode Higgs1964-2 ; Varma2002 ; Pekker2015 . The Higgs mode was initially identified in a superconductor with a coexisting charge-density wave order by Raman spectroscopy Sooryakumar1980 ; Littlewood1981 ; Measson2014 . Nonlinear terahertz spectroscopy has enabled the observation of the Higgs mode in conventional superconductors as well, in the form of a free or forced oscillation and the resulting third-harmonic generation (THG) Matsunaga2013 ; Matsunaga2014 ; Matsunaga2017-2 . This technique has been also applied to unconventional superconductors to investigate elementary excitations and paring mechanisms Katsumi2018 ; Chu2019 .

Multiband superconductors can accommodate a variety of collective modes because of multiple degrees of freedom. For example, two interacting phase modes are transformed into a high-energy plasma oscillation and a low-energy gapped mode, called a Leggett mode Leggett1966 ; Sharapov2002 ; Burnell2010 . The Leggett mode has been observed in a two-gap superconductor MgB2 by Raman spectroscopy Blumberg2007 ; Klein2010 ; Cea2016-2 . Multiple Higgs modes should also interact with each other through the interband paring interaction Murotani2017 , though it is not yet experimentally observed. Even more interestingly than the single-band cases, nonlinear terahertz spectroscopy offers a unique opportunity to investigate interaction among these collective modes, and the origin of paring in multiband superconductors Krull2016 . Recently, a terahertz pump-probe measurement on MgB2 has been reported Giorgianni2019 , where the pump-induced oscillation of the transmitted probe electric field is attributed to the Leggett mode. However, the theoretical analysis presented there neglects the effect of impurity scattering, which may underestimate the contribution from the Higgs mode.

In single-band superconductors, it has been theoretically shown that optical transitions mediated by nonmagnetic impurity scattering significantly enhances nonlinear optical response of the Higgs mode Jujo2015 ; Jujo2018 ; Silaev2019 . Typical multiband superconductors, such as MgB2 and iron-based superconductors, actually often exhibit characteristics of impurity scattering in the energy region around the superconducting gaps even in the linear response Ortolani2008 ; Pimenov2013 . Therefore, it is indispensable to take into account the effects of impurities beyond the clean-limit analysis Tsuji2015 ; Cea2016-1 ; Murotani2017 ; Cea2018 ; Giorgianni2019 to consider nonlinear optical response of multiband superconductors. For that purpose, we extend the Mattis-Bardeen (MB) model for linear response of single-band Bardeen-Cooper-Schrieffer (BCS)-type superconductors Mattis1958 ; Zimmermann1991 ; Berlinsky1993 to a multiband system and nonlinear regime. We found that optical response of the Leggett mode is insensitive to nonmagnetic impurity scattering. Consequently, magnitude of the light-induced oscillation of the Leggett mode is far smaller than that of the Higgs mode assisted by impurities. We also revealed that the Higgs mode contributes to THG dominating the contributions from the Leggett mode and quasiparticle excitations.

This paper is organized as follows. Section II formulates the Hamiltonian which takes into account both impurity scattering and light-matter interaction, and introduces a density matrix to describe light-induced collective modes and nonlinear optical response. Section III analyzes the linear response of a two-band superconductor. Section IV considers light-induced collective modes in the presence of impurity scattering. First, we confirm the effect of impurities in a single-band superconductor which was previously tested by Green’s function method Jujo2015 ; Jujo2018 ; Silaev2019 . We then examine the excitation of the Higgs and Leggett modes in a two-band superconductor for several excitation conditions. Section V discusses THG, decomposing contributions from the Higgs, Leggett modes and quasiparticles. Finally, Sec. VI summarizes the obtained results.

II Density matrix formulation

II.1 Hamiltonian

We extend the Mattis-Bardeen (MB) model Mattis1958 of light-matter interaction in superconductors to a multiband system. First, the noninteracting part of the Hamiltonian is written as

[TABLE]

where $c_{i\mathbf{k}\sigma}^{\dagger}$ ($c_{i\mathbf{k}\sigma}$) creates (annihilates) an electron with the crystal momentum $\mathbf{k}$ and the spin $\sigma$ in an energy band labeled by $i$. We assume the inversion symmetry of the system requiring $\epsilon_{i\mathbf{k}}=\epsilon_{i(-\mathbf{k})}$. In the following, we restrict ourselves to the simplest case where every band is well described by a parabolic dispersion relation,

[TABLE]

with $m_{i}(>0)$ and $\epsilon_{\mathrm{F}i}(>0)$ being the effective mass and the Fermi energy, respectively. The band extremum, set at the $\Gamma$ point ($\mathbf{k}=0$) in Eq. (2) for simplicity, can be moved to any point in the Brillouin zone, as long as the overall inversion symmetry is maintained.

To be exact, $\mathbf{k}$ cannot be regarded as the crystal momentum in the presence of impurities that break the translational symmetry of the crystal. Nevertheless, it serves as a quantum number labeling the true energy eigenstates. When the concentration of impurities is not too high, the energy eigenvalues should be hardly modified, which allows us to use Eq. (2) even for dirty systems.

For later convenience, the Fermi momentum $k_{\mathrm{F}i}=\sqrt{2m_{i}\epsilon_{\mathrm{F}i}}$ and the Fermi velocity $v_{\mathrm{F}i}=k_{\mathrm{F}i}/m_{i}$ are introduced. Note that we have chosen $\hbar=1$. We assume band-dependent quantities (e.g., $\epsilon_{\mathrm{F}i}$ and $k_{\mathrm{F}i}$) to be in the same order for every band. The respective order will be represented by symbols without the band index (e.g, $\epsilon_{\mathrm{F}}$ and $k_{\mathrm{F}}$).

Next, the pairing interaction is given by

[TABLE]

where $U_{ij}$ is the pairing potential between $i$th and $j$th bands Suhl1959 . Throughout the paper, the volume of the system is set to unity.

In single-band superconductors, Anderson’s theorem guarantees robustness of the paring potential (3) against nonmagnetic impurities Anderson1959 . In multiband systems, however, Anderson’s theorem generally fails because of interband scattering Stanev2014 . As a result, $U_{ij}$ can depend on density of impurities and concomitantly on momenta $\mathbf{k}$ and $\mathbf{k}^{\prime}$. It then should modify the transition temperature and the magnitude of superconducting gaps Golubov1997 . However, it has been shown that both quantities only slightly depend on the concentration of impurities in MgB2 Mazin2002 ; Ortolani2008 . In FeSe being another typical multiband superconductor, superconductivity is more robust against impurities than expected from theory considering interband impurity scattering Urata2016 . Therefore, we neglect this effect and adopt Eq. (3) in the following. Physically, this assumption requires that the impurities cannot scatter electrons over the separation between different Fermi surfaces.

According to the above reasoning, we neglect interband impurity scattering also in optical transitions. Then, only intraband transitions suffice for consideration of low-energy optical response, described by an interaction Hamiltonian

[TABLE]

Here, $\mathbf{J}_{i\mathbf{k}\mathbf{k}^{\prime}}$ is the matrix element of the current operator $e\hat{\mathbf{p}}/m_{0}$ with the momentum operator $\hat{\mathbf{p}}$, the charge $e(<0)$ and the mass $m_{0}$ of a free electron, and A is vector potential of the external light field assumed to be spatially homogeneous upon dipole approximation. We also have to take into account the nonlinear coupling to the light field,

[TABLE]

to preserve the gauge invariance of the model. Finally, the total Hamiltonian is given by

[TABLE]

$\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ describe the paramagnetic and diamagnetic coupling to the light field, respectively. In an ideal crystal (the clean limit explained below), the former coupling is known to be negligible Matsunaga2014 ; Tsuji2015 . In dirty systems, however, the paramagnetic coupling is essential. For example, $\mathcal{H}_{1}$ is indispensable to describe linear response of dirty single-band BCS superconductors Mattis1958 ; Zimmermann1991 ; Berlinsky1993 . It has been predicted that this coupling would dominate also the nonlinear optical response Jujo2015 ; Jujo2018 ; Yu2017 ; Yang2018a ; Yang2018b ; Silaev2019 . In terms of the Green’s function method, this is because impurity scattering produces nonvanishing Feynman diagrams that vanish in the clean limit Tsuji2016 . We adopt the MB model rather than employing Green’s functions to examine effects of the paramagnetic coupling assisted by the impurity scattering.

The essence of the MB model Mattis1958 lies in assuming the transition matrix element to be well approximated by a Lorentzian function,

[TABLE]

where $\mathbf{e}$ is the polarization vector of light ($\mathbf{A}(t)=A(t)\mathbf{e}$), $N_{i}(0)=m_{i}k_{\mathrm{F}i}/2\pi^{2}$ is the density of states per spin on the Fermi surface, and $\gamma_{i}$ is the impurity scattering rate. See Appendix A for derivation. Here, it is assumed that

[TABLE]

Condition (9) requires the impurity scattering to occur only within a thin shell around the Fermi surfaces.

When the considered energy region is far smaller than $\gamma$ (the dirty limit), the transition matrix element can be approximated by an energy-independent constant. Such an approach has been used to calculate transient optical conductivity of photoexcited superconductors Chou2017 ; Kennes2017 .

II.2 Energy scales

The interaction Hamiltonian $\mathcal{H}_{\mathrm{int}}$ is treated in the mean-field approximation. As a result, the superconducting gap function

[TABLE]

is defined and regarded as the multicomponent order parameter. At the thermal equilibrium, $2\Delta_{i}^{\mathrm{eq}}$ is self-consistently determined and serves as a gap in the excitation spectrum (see Appendix B). As in usual BCS theory, we assume

[TABLE]

Now, we have three characteristic energy scales:

(A) The Fermi energy $\epsilon_{\mathrm{F}}$.

(B) The scattering rate $\gamma$.

(C) The superconducting gap 2$\Delta^{\mathrm{eq}}$.

We have already assumed two inequalities, namely Eq. (9) between (A) and (B) and Eq. (11) between (A) and (C). The remaining relationship between (B) and (C) concerns the distinction between “clean” and “dirty” systems:

[TABLE]

In the clean limit, we can neglect the impurity scattering in the frequency region around $2\Delta^{\mathrm{eq}}$ that we are interested in. On the other hand, the impurity scattering plays a significant role in mediating nonlinear interaction between light and collective modes in the dirty limit Jujo2015 ; Jujo2018 . We also call more general cases “dirty,” when $\gamma$ and $2\Delta^{\mathrm{eq}}$ are in the same order.

II.3 Equation of motion

To consider linear and nonlinear optical response, we introduce the density matrix for Bogoliubov quasiparticles,

[TABLE]

where the two-component spinor $\psi_{i\mathbf{k}}$ is defined by Eq. (67) in Appendix B. We solve the corresponding equation of motion, i.e., Eq. (76) in Appendix B, in a perturbative manner with respect to the external field $\mathbf{A}(t)$:

[TABLE]

where the additional subscript denotes the order of $\mathbf{A}$.

Motion of the gap function can be calculated through Eq. (10). Here, one has to pay attention to consistency of the formulation. The variation $\delta\Delta_{i}(t)=\Delta_{i}(t)-\Delta_{i}^{\mathrm{eq}}$ formally acts as an external field; the explicit perturbation Hamiltonian is given by Eq. (107) in Appendix B. The resulting motion of quasiparticles in turn affects $\delta\Delta_{i}(t)$ itself through the definition (10). This feedback process induces the collective modes and must be taken into account in numerical calculations.

To consider optical response itself, we need to calculate electric current density

[TABLE]

with the paramagnetic component $\textbf{j}_{\mathrm{P}}=-\langle\delta\mathcal{H}_{1}/\delta\mathbf{A}\rangle$ and the diamagnetic component $\textbf{j}_{\mathrm{D}}=-\langle\delta\mathcal{H}_{2}/\delta\mathbf{A}\rangle$. Their explicit expressions are given by Eqs. (117), (118) in Appendix B.

III Linear response

First, we consider the linear response. It is known that the collective modes do not respond to light in the linear response regime in the absence of any other external fields. As a result, we can neglect $\delta\mathcal{H}_{\mathrm{int}}$ and consider only $\mathcal{H}_{1}$ as the perturbation Hamiltonian. We find that the paramagnetic component is given by

[TABLE]

where $n_{i}=k_{\mathrm{F}i}^{3}/3\pi^{2}$ is the density of carriers (either electrons or holes), $l_{i}(\epsilon_{i\mathbf{k}},\epsilon_{i\mathbf{k}^{\prime}})=l_{i\mathbf{k}\mathbf{k}^{\prime}}$ and $p_{i}(\epsilon_{i\mathbf{k}},\epsilon_{i\mathbf{k}^{\prime}})=p_{i\mathbf{k}\mathbf{k}^{\prime}}$ are coherence factors defined by Eq. (116), and the function $F_{i}^{ab}(\epsilon,\epsilon^{\prime})$ follows

[TABLE]

with $E=E_{i\mathbf{k}}$, $f^{\prime}=f_{i\mathbf{k}^{\prime}}$, etc. On the other hand, the diamagnetic component is given by

[TABLE]

For derivation, see Appendix C.

After the summation $\mathbf{j}(t)|_{1}=\mathbf{j}_{\mathrm{P}}(t)|_{1}+\mathbf{j}_{\mathrm{D}}(t)|_{1}$, we can obtain the optical conductivity

[TABLE]

where $j(\omega)$ is the Fourier transform of $\textbf{e}\cdot\mathbf{j}(t)|_{1}$ and $E(\omega)$ is that of the electric field $E(t)=-\partial A(t)/\partial t$. The real part $\sigma_{1}(\omega)$ corresponds to absorption of light while the imaginary part $\sigma_{2}(\omega)$ corresponds to refraction.

For simplicity, we considered a two-band system with the following parameters:

(Band 1) $s_{1}=+1$, $m_{1}=1.0$, $\epsilon_{\mathrm{F}1}=500$, $\gamma_{1}=10$.

(Band 2) $s_{2}=-1$, $m_{2}=1.2$, $\epsilon_{\mathrm{F}2}=300$, $\gamma_{2}=10$.

(Paring interaction) $U_{11}=0.08$, $U_{22}=0.18$, $U_{12}=0.05$, $\omega_{\mathrm{D}}=10$.

The chosen parameters qualitatively simulate MgB2 with electron- and hole-like bands with $\epsilon_{\mathrm{F}}\sim 1$ eV, $\gamma\sim 0.1$ eV, and $2\Delta\sim 0.01$ eV Ortolani2008 . To avoid formal failure of the weak-coupling approximation in our theory, we chose relatively small values for the paring potential such that $N_{2}(0)U_{22}\simeq 0.3$ while the reported values are relatively large, e.g., $N_{2}(0)U_{22}\simeq 1$ Blumberg2007 . We expect that the results presented below will be qualitatively the same even for strong couplings.

Temperature dependence of the superconducting gaps obtained from the gap equation (69) is shown in Fig. 1 (a). Both bands establish superconductivity ($2\Delta_{i}\neq 0$ with the superscript “eq” omitted) at a common critical temperature $T_{\mathrm{c}}=0.54$. They take the maximum values at the absolute zero, $2\Delta_{1}(0)=0.97$ and $2\Delta_{2}(0)=2.1$.

We introduced a monocycle electric field described by

[TABLE]

with $\tau=2\pi\times 0.1$, $\Omega=0$. The electric field waveform is plotted in Fig. 1 (b). We solved Eqs. (17) and (18) in the time domain to calculate the electric current. The simulated results are plotted in Fig. 1 (c) for $T=0.6>T_{\mathrm{c}}$ (dotted) and $T=0.02<T_{\mathrm{c}}$ (solid). While the former closely follows the electric field waveform, the latter shows a characteristic oscillatory structure. We calculated the optical conductivity (21) and show it in Fig. 1 (d) for $T=0.6$ and in (e) for $T=0.02$. The solid and open circles are the real and imaginary parts, respectively. Above $T_{\mathrm{c}}$, the optical conductivity is structureless in the plotted frequency range as seen in (d). It can be shown that the response exactly coincides with the Drude model Mattis1958 ,

[TABLE]

which is shown in Fig. (d) as solid lines. Small discrepancy between the calculated result and the Drude model in $\sigma_{2}(\omega)$ arises from the cutoff introduced in numerical integration.

On the other hand, below $T_{\mathrm{c}}$ (Fig. 1 (e)), the imaginary part (open circles) diverges toward $\omega\to 0$ being an indicator of superconductivity. Even more remarkably, the real part (solid circles) shows a double-gapped structure. The absorption edges coincide with the superconducting gaps $2\Delta_{i}$ indicated by arrows. The dashed line depicts the contribution from the band with the smaller gap $2\Delta_{1}$ as a guide to the eye. It can be shown that the response exactly coincides with the sum of MB conductivity explicitly given in Refs. Zimmermann1991 ; Berlinsky1993 over all bands. MB model conductivity is plotted as the solid lines in Fig. 1 (e), displaying a good agreement with the simulation. This ensures the validity of our time-domain calculation.

IV Nonlinear excitation of collective modes

Let us proceed to the next order. The relevant equation of motion is given in Appendix D. It is well known that an isotropic system exhibits no even-order nonlinearities, so that

[TABLE]

Therefore, we can concentrate on $\delta\Delta_{i}(t)|_{2}$. Assuming the particle-hole symmetry, the real and imaginary parts are rewritten as

[TABLE]

respectively, where $u_{i}(\epsilon_{i\mathbf{k}})=u_{i\mathbf{k}}$, $v_{i}(\epsilon_{i\mathbf{k}})=v_{i\mathbf{k}}$, and

[TABLE]

Thus we can drop the dependence of the density matrix on the angle of $\mathbf{k}$. The induced motion of the density matrix can be also decomposed into the quasiparticle, Higgs mode, and Leggett mode.

IV.0.1 Quasiparticles

The diagonal components $r_{i}^{11}(\epsilon)$ and $r_{i}^{22}(\epsilon)$ correspond to the quasiparticle excitation. They follow

[TABLE]

where

[TABLE]

Because $S_{i}^{11}(\epsilon)=S_{i}^{11}(-\epsilon)$, we have

[TABLE]

IV.0.2 Higgs mode

The non-diagonal components $r_{i}^{21}(\epsilon)=r_{i}^{12}(\epsilon)^{*}$ correspond to the collective modes. To separate the Higgs mode and the Leggett mode, we further decompose them into the odd and even parts,

[TABLE]

Each of them satisfies

[TABLE]

Among them, the odd component corresponds to the Higgs mode. This can be easily understood because only the odd component contributes to the right-hand side of Eq. (25). Note that $r_{i}^{11}(\epsilon)$ and $r_{i}^{22}(\epsilon)$ also contribute, but because their motion is determined by a closed equation of motion (28), quasiparticle excitations only trigger the Higgs mode without any feedback effect. The odd component $r_{i}^{21,\mathrm{odd}}(\epsilon)$ follows

[TABLE]

where

[TABLE]

The source term $S_{i}^{21}(\epsilon)$ arises from the paramagnetic coupling with light mediated by the impurity scattering. Note that Eq. (35) has to be solved keeping Eq. (25) to be always satisfied. This self-consistency condition induces the Higgs mode resonance.

For the incident frequency in the order of $2\Delta$, the amplitude of the induced Higgs mode is estimated as follows (see Appendix E):

[TABLE]

This estimation predicts the most efficient excitation of the Higgs mode at $\gamma\sim 2\Delta$, which will be confirmed numerically. Equation (37) arises from only the paramagnetic coupling (abbreviated as “Para”) because the diamagnetic coupling (to be abbreviated as “Dia”) vanishes in Eq. (35). However, it is known that the latter can also induce the Higgs mode when the energy dispersion relation exhibits a nonparabolicity in the clean limit Tsuji2015 . Because the diamagnetic term (4) of the Hamiltonian is less sensitive to impurity scattering than the paramagnetic term (5), it is reasonable to assume that the correction by a small nonparabolicity does not depend on $\gamma$:

[TABLE]

Here, we have used that $U\sim\epsilon_{\mathrm{F}}$ which is valid for phonon-mediated interactions Coleman2015 . Because the diamagnetic contribution (38) is smaller than paramagnetic one (37), the paramagnetic coupling is more important in optical excitation of the Higgs mode. Result of the order estimation is summarized in the second row of Table 1.

The excitation mechanism of the Higgs mode can be understood as follows. Because this mode has an even parity in momentum space, dipole-allowed intermediate states are necessary to excite it through a two-photon process via the paramagnetic coupling. In the clean limit, there is no dipole-allowed excitation because all the spectral weight concentrates on the zero frequency. However, presence of impurity scattering produces dipole-allowed excitations as exemplified by Fig. 1 (e), which then enhance the two-photon excitation of the Higgs mode Jujo2015 .

It has been suggested that retardation of the phonon-mediated interaction also enhances the nonlinear optical response of the Higgs mode Tsuji2016 . Even in this case, the enhancement arises from the paramagnetic coupling, termed there “resonant coupling” from an analogy with Raman scattering. It is an interesting problem to compare it with the effect of impurity scattering quantitatively, which is outside the scope of this paper. However, retarded interaction predicts comparable contributions by the Higgs mode induced via the paramagnetic or “resonant” coupling and by quasiparticles excited via the diamagnetic or “nonresonant” coupling in THG. This contrasts with our discussion in the next section which predicts dominance of the former contribution. Therefore, it is possible that impurity scattering is more important in enhancing the coupling between the Higgs mode and light.

IV.0.3 Leggett mode

The even component $r_{i}^{21,\mathrm{even}}(\epsilon)$ corresponds to the Leggett mode. It follows

[TABLE]

The paramagnetic coupling vanishes here. This equation has to be solved consistently with Eq. (26), which induces the Leggett mode. Its amplitude is estimated as

[TABLE]

This result is corrected by a particle-hole asymmetry which modifies energy integrals through the finite slope of the density of states on the Fermi surface. The resulting correction is in the following order:

[TABLE]

This correction is not essential, because it does not exceed Eq. (40). Therefore, we will neglect it in the following. Result of the order estimation is summarized in the third row of Table 1.

The diamagnetic contribution (40) has been explained in terms of the coupling between the Leggett mode and the potential difference between Fermi surfaces induced by the diamagnetic coupling Murotani2017 . Along this line, one can infer from the above result that quasiparticle excitation caused by the paramagnetic coupling induces no potential difference in the presence of particle-hole symmetry. This resembles the behavior of chemical potential being independent of temperature when carriers exhibit a constant density of states. For a system slightly lacking the particle-hole symmetry, the paramagnetic contribution (41) can be viewed as a two-photon excitation similar to that of the Higgs mode.

IV.1 Single-band case

Before discussing the multiband case, let us verify the effect of impurity scattering on the Higgs mode in a single-band superconductor, previously investigated by the Green’s function method Jujo2015 . We used $N(0)U=0.27$ and $\omega_{\mathrm{D}}=10$, which leads to $2\Delta=1$ at $T=0$. The real part of the optical conductivity $\sigma_{1}(\omega)$ is plotted in Fig. 2 (a) for $\gamma/2\Delta$=0.1, 0.2, 1, 3, and 5, as blue, green, yellow, orange, and red lines, respectively. We introduced a multicycle pulse described by vector potential (22) with $\tau=2\pi\times 1$, $\Omega=0.5$. Its waveform $E(t)$ is shown in Fig. 2 (b) while the power spectrum $|E(\omega)|^{2}$ is plotted as the dotted line in (a). The latter is located inside the gap $2\Delta$, indicating no optical excitation of quasiparticles. The simulated dynamics of $\delta\Delta^{\prime}(t)$ is plotted in Fig. 2 (c) with the same parameters and colors as in (a). All curves show a clear oscillation with the doubled frequency $2\Omega$. Free oscillation of the Higgs mode remains after excitation ($2\Delta t/2\pi>6$), because the incident pulse satisfies the resonance condition for the Higgs mode, $2\Omega=2\Delta$ Murotani2017 ; Jujo2018 . The amplitude of oscillation, however, depends on $\gamma$. We plotted the maximum of $|\delta\Delta^{\prime}(t)|$ as a function of $\gamma$ as circles in Fig. 2 (d), which takes the maximum value at $\gamma\sim 2\Delta$ consistently with the above order estimation. In Fig. 2 (d), we also plotted the maximum value of $\sigma_{1}(\omega)$, that clearly correlates with the amplitude of the Higgs mode. This result indicates that virtual excitation of optically active intermediate states concerns excitation of the Higgs mode. All these properties successfully reproduce the results obtained by the Green’s function method Jujo2015 .

IV.2 Two-band case

Now, we return to the two-band system considered in the last section. We concentrate on $T=0.02$ well below $T_{\mathrm{c}}$. First, we examine a nonadiabatic excitation in which the electric field varies faster than the superconducting response time $2\pi/2\Delta$. The monocycle pulse used in the previous section fits this purpose. In Fig. 3 (a), the electric field waveform is shown again. Its power spectrum plotted in Fig. 3 (d) exhibits a broad bandwidth covering both gaps indicated by dotted lines. Figure 3 (b) shows the induced dynamics of $\delta\Delta_{1,2}^{\prime}(t)$. Oscillations after excitation correspond to the Higgs mode. Both variables asymptotically approach negative values because excitation of Bogoliubov quasiparticles by a broadband pulse suppresses the superconducting order. Figure 3 (e) shows the power spectra of $\delta\Delta_{1,2}^{\prime}(t)$. Both of them exhibit peaks at $\omega=2\Delta_{1}$ and $2\Delta_{2}$, being the resonance energies of the two Higgs modes Murotani2017 . Due to the interband interaction $U_{12}$, they mutually interact so that both gaps oscillate with both frequencies. On the other hand, Fig. 3 (c) shows dynamics of the phase difference

[TABLE]

which displays a damped oscillation of the Leggett mode. Its power spectrum is shown in Fig. 3 (f). Because relatively large interband interaction $U_{12}$ pushes up the resonance energy above the lower gap $2\Delta_{1}$, the mode acquires finite lifetime coming from decay into the quasiparticle continuum. These oscillations of the gaps may be detectable in pump-probe experiments Matsunaga2013 . However, having seen that the Leggett mode is far smaller than the Higgs mode, $\delta\Delta^{\prime\prime}(t)\ll\delta\Delta^{\prime}(t)$, the observation of the Leggett mode is expected to be difficult.

Next, we consider the excitation by multicycle pulses. In Fig. 4 (a), we plot the power spectra of pulses with $\Omega=0.25$, 0.5, 0.75, 1, 1.5, and 3, from left to right. Another parameter $\tau$ is set to $\tau=2\pi\times 1$ for the former five and $\tau=2\pi\times 0.3$ for the last one.

Figures 4 (b1)-(b3) show $E(t)$, $\delta\Delta_{1,2}^{\prime}(t)$, and $\delta\varphi(t)$, respectively, for $\Omega=0.25$, $\tau=2\pi\times 1$. In this case, photon energy $\Omega<2\Delta_{1}$ is insufficient to excite Bogoliubov quasiparticles. Also, it deviates from the resonance condition for the Higgs mode ($2\Omega\simeq 2\Delta_{1,2}$). As a result, the gap function varies only through virtual excitation of quasiparticles and approaches 0 right after the electric field vanishes.

Figures 4 (c1)-(c3) show the case with $\Omega=0.5$. Photon energy is again insufficient to excite quasiparticles ($\Omega<2\Delta_{1}$), but now satisfies the resonance condition for the lower-energy Higgs mode, $2\Omega\simeq 2\Delta_{1}$. As a result, a small free oscillation remains in both gaps after illumination.

Figures 4 (d1)-(d3) show the case with $\Omega=0.75$. Still, quasiparticles are not excited because $\Omega<2\Delta_{1}$, and the resonance condition $2\Omega\simeq 2\Delta_{1,2}$ is not met again. Correspondingly, the gap functions rapidly approach 0 after forced oscillation. In fact, $2\Omega=1.5$ is close to the resonance condition for the Leggett mode, because its energy is about 1.3 in the present model Murotani2017 . However, Fig. 4 (d3) shows no remarkable structure because it is damped.

Figures 4 (e1)-(e3) show the case with $\Omega=1$. Now, the power spectrum overlaps with the onset of conductivity ($2\Delta_{1}$) as seen in Fig. 4 (a), so that quasiparticles are excited. In addition, the resonance condition for the larger-energy Higgs mode $2\Omega\simeq 2\Delta_{2}$ is satisfied, so that the Higgs mode is also excited. As a result, $\delta\Delta_{1,2}^{\prime}(t)$ oscillates even after the excitation because of the induced Higgs mode, and approaches a negative value because of the excited quasiparticles.

Figures 4 (f1)-(f3) show the case with $\Omega=1.5$. Only quasiparticles are excited because the photon energy is apart from the resonance conditions for collective modes. Correspondingly, the gap functions show a forced oscillation only under the electric field, and quickly approach the final values.

Finally, we examine the impulsive stimulated Raman scattering of low-energy light by collective modes. To this end, we considered a pulse with $\Omega=3$, $\tau=2\pi\times 0.3$. The rightmost curve in Fig. 4 (a) plots the corresponding power spectrum, showing a bandwidth comparable with the smaller gap and thus sufficient to excite the lower-energy Higgs mode with a Raman process. The simulated result is shown in Figs. 4 (g1)-(g3). Now the oscillation of the incident pulse is so fast that the real parts $\delta\Delta_{1,2}^{\prime}(t)$ cannot follow it. As a result, they almost monotonically decrease within the pulse duration. After excitation ($t/2\pi>4.5$), however, a small oscillation remains, which reflects the Higgs mode induced by the Raman process. But its oscillation is blurred by the large gap reduction by quasiparticle excitations. On the other hand, oscillation of the Leggett mode is more easily seen in Fig. 4 (c).

In ordinary Raman experiments, the incident photon energy exceeds $\gamma$ ($\sim 0.1~{}\mathrm{eV}$), which makes the intraband transitions assisted by the impurity scattering irrelevant. Then, electronic Raman scattering will be dominated by interband transitions Klein1984 ; Klein2010 . This consideration may explain the reason why the Higgs mode in MgB2 has not been observed through the spontaneous Raman scattering Blumberg2007 .

To be exact, relaxation processes will modify light-induced dynamics of the gap functions. For example, strong electron-phonon interaction gives rise to damping of the Higgs mode at high temperatures Murakami2016 while carrier-carrier scattering redistributes the photoexcited quasiparticles. However, it is reasonable to expect that the calculated results will not be significantly altered, especially just under illumination of the optical pulse. Thus we do not consider relaxation processes other than the impurity scattering already taken into account.

V Third-harmonic generation

Let us now consider THG. We start from a brief discussion of its origin. Figure 5 schematically shows the processes that concern THG in BCS superconductors. First, Fig. 5 (a) gives the simplest one in which a Bogoliubov quasiparticle absorbs three photons and emits one photon. No collective mode concerns this process. However, quasiparticles can emit not only photons but also collective modes. In this point of view, Eqs. (25) and (26) can be regarded as the solution of a wave equation for $\delta\Delta$ with the source $r^{ab}(\epsilon)$. The emitted $\delta\Delta$ then formally acts as an external field in Eq. (107) and thus is absorbed by quasiparticles again. Such a sequence induces the collective modes as depicted in Fig. 5 (b). As a result, the collective modes can also contribute to THG through diagrams shown in Fig. 5 (c).

Below, we assume that an isotropic impurity scattering mixes the states close to the Fermi energy almost equally. Then, direction of a quantum number $\mathbf{k}$ can be neglected, which enables us to approximate $|\mathbf{e}\cdot\mathbf{J}_{i\mathbf{k}\mathbf{k}^{\prime}}|^{2}\simeq\langle|\mathbf{e}\cdot\mathbf{J}_{i\mathbf{k}\mathbf{k}^{\prime}}|^{2}\rangle_{\mathrm{Av}}$. Then, it follows that

[TABLE]

In fact, the right-hand side gives the minimum among the possible values of the left-hand side. The possible maximum is $9/5$ times larger, which is not far from the above equation. Therefore, we expect that the exact form of the transition matrix elements does not matter significantly to the final results.

The third-order current density is given by

[TABLE]

where

[TABLE]

and

[TABLE]

The function $R_{i}^{ab}(\epsilon,\epsilon^{\prime})$ follows the equation of motion

[TABLE]

under the particle-hole symmetry. Only the Higgs mode and quasipaticle excitations contribute to the paramagnetic component. An order estimation gives

[TABLE]

Here, HM and QP are abbreviations of the Higgs mode and quasiparticles, respectively.

When a small particle-hole asymmetry is taken into account, the Leggett mode (LM) also contributes to $\mathbf{j}_{\mathrm{P}}(t)|_{3}$. However, its order is small compared to the others:

[TABLE]

Therefore, this contirbution is negligible.

We turn to the diamagnetic component. For a parabolic and isotropic dispersion relation, only $r_{i}^{21,\mathrm{even}}(\epsilon)$ contributes to Eq. (46). In other words, only the Leggett mode contributes to the diamagnetic component. This is understood in terms of Eq. (20), which can be extended to nonlinear current by allowing carrier density $n_{i}$ to vary. While the Leggett mode changes $n_{i}$ through the interband Josephson coupling $U_{12}$, the Higgs mode and quasiparticles do not. Therefore, the latter two do not contribute to the diamagnetic component of the third-order current. Its order is then estimated as

[TABLE]

In the clean limit, however, it is known that a nonparabolicity of the dispersion relation enables the Higgs mode and quasiparticles to induce nonzero $\mathbf{j}_{\mathrm{D}}(t)$ Tsuji2015 ; Murotani2017 . The corresponding correction is given by

[TABLE]

Again, it is reasonable to assume the validity of this order estimation even for dirty cases, because of the insensitivity of the interaction Hamiltonian (5) against the impurity scattering. The order estimation is summarized in Table 2.

As already mentioned, only the diamagnetic coupling contributes to THG in the clean limit. And the phase degree of freedom is negligible in single-band systems. As a result, the combination of quasiparticles and the diamagnetic coupling (called charge-density fluctuations) gives the dominant origin of THG Cea2016-1 , unless one includes the retardation effect beyond the BCS mean-field treatment Tsuji2016 . The Leggett mode also contributes in the same order in multiband systems Murotani2017 . In dirty systems, however, Eq. (54) is negligible compared to Eq. (51). Therefore, in most superconductors exhibiting a dirty nature ($2\Delta<\gamma$), the dominant origin of THG will be the Higgs mode and quasiparticles excited by the paramagnetic coupling.

The above order estimation cannot reveal the relative importance of the Higgs mode and quasiparticles. To examine it, we considered a two-band system with the same parameters as in Sec. III. For simplicity, we concentrated on $T=0.02$ well below $T_{\mathrm{c}}=0.54$ and on a multicycle pulse (22) with $\Omega=1$ and $\tau=2\pi\times 1$, which was used also in Fig. 4 (e). The chosen frequency satisfies the resonance condition $2\Omega\simeq 2\Delta_{2}$ both for the Higgs mode with the larger energy and for the quasiparticle excitations with the larger gap.

The electric field waveform is plotted in the top panel of Fig. 6 (a) with the induced third-order current $j^{(3)}(t)=\mathbf{e}\cdot\mathbf{j}(t)|_{3}$ in the next panel. We can decompose the latter into contributions from the Higgs mode, quasiparticles, and the Leggett mode, which are shown in the subsequent three panels in order. In our approximations, the diamagnetic component $j^{(3)}_{\mathrm{L}}(t)=\mathbf{e}\cdot\mathbf{j}_{\mathrm{D}}(t)|_{3}$ corresponds to the Leggett mode (green). This is about 5000 times smaller than the others, being consistent with the order estimation predicting a $(\epsilon_{\mathrm{F}}/\gamma)^{2}\sim 10^{3}$ times smaller contribution for the chosen parameters. The contribution from quasiparticles $j^{(3)}_{\mathrm{Q}}(t)=\mathbf{e}\cdot\mathbf{j}_{\mathrm{P}}(t)|_{3,\delta\Delta^{\prime}=0}$ is obtained by neglecting the self-consistency condition (25), i.e., artificially putting $\delta\Delta^{\prime}=0$ (yellow). The remaining part $j^{(3)}_{\mathrm{H}}(t)=\mathbf{e}\cdot[\mathbf{j}_{\mathrm{P}}(t)|_{3}-\mathbf{j}_{\mathrm{P}}(t)|_{3,\delta\Delta^{\prime}=0}]$ arises from the Higgs mode (orange).

All these contributions consist of components with the frequency $\Omega+\Omega+\Omega=3\Omega$ (THG) and with $\Omega-\Omega+\Omega=\Omega$ (two-photon absorption). First, we examine the former. In Fig. 6 (b), we plot the power spectrum of the induced current $j^{(3)}(t)$ for $2\leq\omega\leq 4$. The topmost (red) curve corresponds to the total third-order current, peaked at the third-harmonic frequency $\omega=3\Omega=3$. The lower three lines are given by considering only the Higgs mode, only quasiparticles, and only the Leggett mode, from top to bottom, respectively. As seen easily, the total third harmonic is dominated by the Higgs mode (orange), while the contribution from quasiparticles is about one order of magnitude smaller (yellow). This behavior is similar to the conductivity spectrum of a single-band superconductor with a dc supercurrent featured by a modification by the Higgs mode larger than by quasiparticles Moor2017 ; Nakamura2018 . Because the contribution from the Leggett mode is so small that it is originally outside the plotted region, the corresponding curve is multiplied by $10^{5}$ (green). This numerical simulation thus reveals that the Higgs mode plays the dominant role in THG,

Finally, we mention the two-photon absorption. We introduce a nonlinear absorption spectrum

[TABLE]

This function gives the nonlinear correction to the absorption spectrum effectively felt by the incident electric field itself. We plotted the calculated spectra in Fig. 6 (c). Because Eq. (55) condenses all sum- and difference-frequency generation processes into dependence on a single frequency $\omega$, an unphysical upturn at the edges of the pump bandwidth appears. However, it does not matter because the absorbed energy is proportional to $\operatorname{Re}\sigma^{(3)}(\omega)$ times $|E(\omega)|^{2}$. In Fig. 6 (c), it can be seen that contributions from the Higgs mode (orange) and quasiparticles (yellow) are comparable, while that from the Leggett mode (green) is small (it is shown after multiplication by 5000). This observation confirms the result of Sec. IV which found that two-photon excitation of the Higgs mode can occur more efficiently than the Leggett mode.

VI Conclusion

In summary, we investigated the nonlinear optical response of weak-coupling multiband superconductors containing nonmagnetic impurities with a density matrix approach. We found that impurity scattering enhances the light-induced Higgs mode through the paramagnetic coupling with the gauge field, while the Leggett mode is left hardly affected. Consequently, light-induced non-equilibrium dynamics of superconducting gaps is dominated by the Higgs mode for a low-energy excitation. We also studied THG in dirty multiband superconductors, revealing the dominant contribution from the Higgs mode. Contribution from quasiparticle excitations is smaller, and Leggett mode will be negligible in this phenomenon.

It is an interesting problem to quantitatively compare the enhancement of the Higgs mode by nonmagnetic impurities with another enhancement by retarded interaction Tsuji2016 . In addition, to the best of our knowledge, it is not known how the interband impurity scattering affects the resonance structure of collective modes. Use of Green’s function method, which can take into account the retardation effect of phonon-mediated interaction Murakami2016 ; Tsuji2016 , interband impurity scattering Sung1967 , and also magnetic impurities Jujo2018 etc., may pave the way for more detailed understanding of non-equilibrium properties of superconductors.

Acknowledgements.

We wish to thank D. Manske and K. Tomita for fruitful discussions. R.S. acknowledges partial support by JSPS KAKENHI Grant Nos. 18H05324 and 15H02102. Y.M. is supported by JSPS Research Fellowship for Young Scientists.

Appendix A Derivation of transition matrix element

In this section, we derive Eq. (7) according to the method of MB Mattis1958 . MB assumed a correlation function in a form of

[TABLE]

where $\phi_{\mathbf{k}}(\mathbf{r})$ is an eigenfunction of the noninteracting Hamiltonian including impurity potential, $\langle~{}\rangle_{\mathrm{Av}}$ denotes an average over the angle of $\mathbf{k}$, $\mathbf{R}=\mathbf{r}-\mathbf{r}^{\prime}$, and $l$ is the mean free path. We assume the above equation for each band independently. From definition, the transition matrix element is given by

[TABLE]

for $i=x,y,z$. Averaging over directions of $\mathbf{k}$ and $\mathbf{k}^{\prime}$ yields

[TABLE]

Note that we have set the volume of the system as unity. Using Eq. (2.13) in MB, this integral is reduced to

[TABLE]

for $k,k^{\prime}\simeq k_{\mathrm{F}}$. Substituting $\epsilon_{\mathbf{k}}\simeq v_{\mathrm{F}}(k-k_{\mathrm{F}})$ into Eq. (59), we obtain Eq. (7) with $\gamma=v_{\mathrm{F}}/l$.

Appendix B Quasiparticle representation

In this section, we construct a density matrix method for Bogoliubov quasiparticles. First, in the thermal equilibrium where the mean-field $\Delta_{i}$ is constant, the Hamiltonian can be diagonalized in a form of

[TABLE]

where we have performed Bogoliubov transformation

[TABLE]

with

[TABLE]

Here, the superscript “eq” labels the equilibrium values. Assuming Fermi statistics of Bogoliubov quasiparticles, Eq. (10) is explicitly rewritten as

[TABLE]

where $\omega_{\mathrm{D}}$ is the Debye frequency. Equation (69) gives the gap equation for a multiband BCS superconductor. We will consider $U_{ij}$ to be a real number, which allows us to choose real values for $\Delta_{i}^{\mathrm{eq}}$ and $v_{i\mathbf{k}}$ at the equilibrium.

Next, we present the equation of motion for the density matrix (13). For that purpose, we construct a 4-component vector

[TABLE]

When the total Hamiltonian is expressed as

[TABLE]

the equation of motion for the 4-component vector is given by

[TABLE]

with

[TABLE]

The equilibrium values are given by

[TABLE]

In addition to the equilibrium Hamiltonian (62), there are three non-equilibrium terms given below. First, when the gap function $\Delta_{i}(t)\equiv\Delta_{i}^{\mathrm{eq}}+\delta\Delta_{i}(t)$ is in motion, $\mathcal{H}_{\mathrm{int}}$ produces a perturbation Hamiltonian

[TABLE]

where $\delta\Delta_{i}^{\prime}=\operatorname{Re}\delta\Delta_{i}$, $\delta\Delta_{i}^{\prime\prime}=\operatorname{Im}\delta\Delta_{i}$. Using the density matrix, the gap function (10) is rewritten as

[TABLE]

Motion of the real and imaginary parts corresponds to the Higgs and Leggett modes, respectively. The latter can be viewed as follows. When $\Delta^{\prime\prime}_{j}$ varies from 0, it induces motion of $\rho_{j\mathbf{k}\mathbf{k}}^{21}$ and $\rho_{j\mathbf{k}\mathbf{k}}^{12}$ due to the coupling Hamiltonian (107). It then changes the value of $\Delta_{i}^{\prime\prime}$ with $i\neq j$ through Eq. (109). Such an interaction between imaginary parts of different bands induces the normal mode of the phase difference, i.e., the Leggett mode.

On the other hand, the light-matter interaction is rewritten as

[TABLE]

where

[TABLE]

are coherence factors Schrieffer1964 . We have used $\mathbf{J}_{i(-\mathbf{k}^{\prime})(-\mathbf{k})}=-\mathbf{J}_{i\mathbf{k}\mathbf{k}^{\prime}}$ which characterizes the case II interaction in the classification by BCS BCS1957 . Using the density matrix, the paramagnetic and diamagnetic components of electric current density are given by

[TABLE]

Appendix C Derivation of linear response

In this section, we derive Eq. (16). The corresponding equation of motion to be solved is

[TABLE]

where

[TABLE]

Equation (119) is solved in a form of

[TABLE]

where the function $F_{i}^{ab}(\epsilon,\epsilon^{\prime})$ follows Eqs. (17)-(19). We can confirm that the gap function does not respond linearly to vector potential, i.e.,

[TABLE]

by substituting Eq. (139) into Eqs. (108) and (109). Because we have considered an isotropic system, the induced current is parallel to $\mathbf{A}$. As a result, the paramagnetic component of the current becomes

[TABLE]

Now, we replace the summation over $\mathbf{k}$ by integration over energy, following

[TABLE]

where $\Omega_{\mathbf{k}}$ denotes the solid angle of $\mathbf{k}$. Extension of the integration interval to $(-\infty,\infty)$ is justified because large $|\epsilon_{i\mathbf{k}}|$ contributes only negligibly to the integral. With the above replacement, we obtain

[TABLE]

Substitution of Eq. (7) into Eq. (143) yields Eq. (16) in the main text.

As for the diamagnetic component, it is more convenient to go back to the interaction Hamiltonian (5) than to use expression (118). From Eq. (5), we obtain

[TABLE]

For an electron band ($s_{i}=+1$), this equation is reduced to Eq. (20). Because Eq. (20) does not depend on the sign of $e$, it is also applicable to a hole band. Strictly speaking, one has to take into account the nonparabolicity in the dispersion relation far below the Fermi surface to derive Eq. (20) for a hole band. This follows from a more general expression of the diamagnetic component of intraband current,

[TABLE]

where the band and spin indices are omitted for simplicity. $n_{\mathbf{k}}$ is the occupation number of an electron in a pure Bloch state (and thus is not $\langle c_{\mathbf{k}}^{\dagger}c_{\mathbf{k}}\rangle$ in the main text). The second line follows from $\nabla_{\mathbf{k}}\epsilon_{\mathbf{k}}=0$ at the edges of the Brillouin zone. For a hole band, we have assumed Eq. (2) with $s=-1$ in the region where $1-n_{\mathbf{k}}\neq 0$. Then, Eq. (20) follows also for a hole band.

In numerical calculations, we introduce a cutoff in the energy integrals in Eq. (16). Then, the real part of the optical conductivity $\sigma_{1}(\omega)$ is incorrectly given at a large photon energy $\omega$ beyond the cutoff. This leads to incorrect evaluation of the imaginary part $\sigma_{2}(\omega)$ even for a low photon energy because of the Kramers-Kronig relation. To improve the evaluation, it is convenient to follow MB’s method Mattis1958 , rewriting

[TABLE]

where $f(\epsilon)=1/(e^{\beta\epsilon}+1)$. Equation (146) is identical to Eq. (20) as long as the integration interval ranges infinity. Even when a cutoff is introduced, this expression guarantees $\sigma_{2}(\omega\to 0)=0$ for the normal state above the critical temperature.

Appendix D Second-order equation of motion

In this section, we give the equation of motion which describes the second-order response of the system. The relevant equation of motion is

[TABLE]

The second-order Hamiltonian can be decomposed into contributions from the diamagnetic coupling (D), the Higgs mode (H), and the Leggett mode (L) as

[TABLE]

First, the diamagnetic contribution is

[TABLE]

Second, the Higgs mode contribution is

[TABLE]

Finally, the Leggett mode contribution is

[TABLE]

A straightforward calculation transforms the equation of motion (147) into Eqs. (28), (35), and (39), for the angle-averaged density matrix $r_{i}^{ab}(\epsilon)$.

Appendix E Order estimation

In this section, we outline the order estimation of light-induced quantities. We exemplify Eq. (37) for $\delta\Delta^{\prime}$ and Eq. (41) for $\delta\Delta^{\prime\prime}$. For simplicity, the absolute zero will be assumed below so that $F^{11}=0$. When the incident photon energy is in the order of $\Delta$, we can approximate $i\partial/\partial t\sim\Delta$ for the order estimation. This transforms Eqs. (18), (28), (35) into

[TABLE]

where we have also used $E\sim\Delta$, $u^{2}\sim 1$, etc. In the last, we have assumed $\delta\Delta^{\prime}\sim S^{21}$ to be confirmed later. We now evaluate the order of $S^{ab}$ in Eqs. (30), (36). In the clean limit, $W(\epsilon,\epsilon^{\prime})$ in the integrand restricts the integral region to $|\epsilon-\epsilon^{\prime}|\lesssim\gamma\ll 2\Delta$. In this case, we have

[TABLE]

which yields

[TABLE]

The arrow follows from Eq. (25) putting $UN(0)\sim 1$ and $\int d\epsilon\sim\Delta$. The above result confirms $\delta\Delta^{\prime}\sim S^{21}$ and gives the upper part of Eq. (37). On the other hand, in the dirty limit, $W(\epsilon,\epsilon^{\prime})\sim\gamma^{-1}$ is almost constant in the relevant energy scale. We can then estimate

[TABLE]

so that

[TABLE]

which gives the lower part of Eq. (37).

To estimate the magnitude of the Leggett mode excited by the paramagnetic coupling, we have to take particle-hole asymmetry into account. For example, particle-hole asymmetry induces a small even component of $S^{21}$ which should appear in the right-hand side of Eq. (39) as a source term. The order of the even component is about $N^{\prime}(0)\Delta/N(0)\sim\Delta/\epsilon_{\mathrm{F}}$ times that of the odd component. This observation enables us to conclude that the magnitude of $\delta\Delta^{\prime\prime}$ induced by the paramagnetic coupling is $\epsilon_{\mathrm{F}}/\Delta$ times smaller than that of $\delta\Delta^{\prime}$ estimated above. Then Eq. (41) follows.

Order of the third-order current can be estimated in the same way.