Influence of scattering on optical response of superconductivity

TL;DR

This paper analytically examines how scattering affects the optical properties of superconductors, revealing a crossover in absorption and a phase shift in Higgs mode responses, with implications for experimental detection.

Contribution

It provides a detailed analytical study of scattering effects on superconducting optical responses, including linear and second-order regimes, using gauge-invariant kinetic equations.

Findings

Optical absorption exhibits a crossover at ω=2|Δ| due to scattering.

Finite low-temperature absorption occurs for ω<2|Δ|, contrasting previous theories.

Scattering induces a phase shift in Higgs mode responses, with a π-jump at ω=|Δ|.

Abstract

By using the gauge-invariant kinetic equation, we analytically investigate the influence of the scattering on the optical properties of superconductors in the normal-skin-effect region. Both linear and second-order responses are studied under a multi-cycle terahertz pulse. In the linear regime, we reveal that the optical absorption $\sigma_{1s}(\omega)$, which origins from the scattering, exhibits a crossover point at $\omega=2|\Delta|$. Particularly, it is further shown that when $\omega<2|\Delta|$, $\sigma_{1s}(\omega)$ from the scattering always exhibits a finite value even at low temperature, in contrast to the vanishing $\sigma_{1s}(\omega)$ in the anomalous-skin-effect region as the Mattis-Bardeen theory revealed. In the second-order regime, responses of the Higgs mode during and after the optical pulse are studied. During the pulse, we show that the scattering causes a phase shift in the second-order response of the Higgs mode. Particularly, this phase shift exhibits a significant $\pi$-jump at $\omega=|\Delta|$, which provides a very clear feature for the experimental detection. After the pulse, by studying the damping of the Higgs-mode excitation, we reveal a relaxation mechanism from the elastic scattering, which shows a monotonic enhancement with the increase of the impurity density.

Tables

Table 1. Table 1: The used parameters in our calculations. With the specific values of Δ 00 subscript Δ 00 \Delta_{00} and ω D subscript 𝜔 𝐷 \omega_{D} , the effective electron-electron attractive potential g 𝑔 g is determined by Eq. ( 24 ) at T = 0 𝑇 0 T=0~{} K.

${\Delta }_{00}$	$1.268$meV	${\omega }_D$	$15.856$meV
$E_F$	$700$meV	$e_0$	$8$meV
${\left(e{E}_0/i\omega \right)}^2/m$	${10}^{-4}{\Delta }_{00}$	$A_0^{⟂}$	$0$

Equations

$$\displaystyle\partial_{t}\rho^{c}_{\bf k}\!+\!i\Big{[}(\xi_{k}\!+\!e\phi\!+\!\mu_{H}+\!\mu_{F})\tau_{3}\!+{\hat{\Delta}({R})},\rho^{c}_{\bf k}\Big{]}\!$$

$$\displaystyle\mbox{}\!+\!i\Big{[}\frac{e^{2}A^{2}}{2m}\tau_{3},\rho^{c}_{\bf k}\Big{]}\!-\!i\Big{[}\frac{1}{8m}\tau_{3},{\bm{\nabla}}^{2}_{\bf R}\rho^{c}_{\bf k}\Big{]}\!+\frac{1}{2}\!\Big{\{}\frac{\bf k}{m}\tau_{3},{\bm{\nabla}}_{\bf R}\rho^{c}_{\bf k}\Big{\}}$$

$$\text{mbox}+\frac{1}{2}\left\{e\mathbf{E}{\tau }_3-({\nabla }_{\mathbf{R}}-2ie\mathbf{A}{\tau }_3)\hat{\Delta }(R),{\partial }_{\mathbf{k}}{\rho }_{\mathbf{k}}^c\right\}$$

$$\displaystyle\mbox{}\!-\!\frac{i}{8}\Big{[}({\bm{\nabla}_{\bf R}}\!-\!2ie{\bf A}\tau_{3})({\bm{\nabla}_{\bf R}}\!-\!2ie{\bf A}\tau_{3}){\hat{\Delta}}(R),\partial_{\bf k}\partial_{\bf k}\rho^{c}_{\bf k}\Big{]}$$

$$\displaystyle\mbox{}\!-\!e\Big{[}\frac{2{\bf A}\cdot{\bm{\nabla}}_{\bf R}\!+\!{\bm{\nabla}}_{\bf R}\cdot{\bf A}}{4m}\tau_{3},\tau_{3}\rho^{c}_{\bf k}\Big{]}\!=\!\partial_{t}\rho^{c}_{\bf k}\Big{|}_{\rm sc}.$$

$$\Delta (R)$$

$${\Delta }^{\ast }(R)$$

$${\mu }_F(R)$$

$${\mu }_H(R)$$

$$e\mathbf{E}$$

$$en$$

$$\mathbf{j}$$

$$e{A}_{\mu }$$

$$\theta (R)$$

$${\partial }_{t}e\delta n+{\nabla }_{\mathbf{R}}\cdot \mathbf{j}=0,$$

$${\partial }_t{\rho }_{\mathbf{k}}+i[({\xi }_k+{\mu }_{\mathrm{eff}}){\tau }_3+|\Delta |{\tau }_1,{\rho }_{\mathbf{k}}]+\frac{i}{8}[|\Delta |{\tau }_1,({\mathbf{p}}_{\mathbf{s}}\cdot {\partial }_{\mathbf{k}}{)}^2{\rho }_{\mathbf{k}}]$$

$$\displaystyle\mbox{}+\frac{1}{2}\left\{e{\bf E}\tau_{3}\!+\!{\bf p}_{s}|\Delta|\tau_{2},\partial_{\bf k}\rho_{\bf k}\right\}\!=\!\partial_{t}\rho_{\bf k}\Big{|}_{\rm sc},$$

$${\mathbf{p}}_s$$

$${\mu }_{\mathrm{eff}}$$

$$g{\sum _{\mathbf{k}}}^{\prime }{\rho }_{\mathbf{k}1}$$

$$g{\sum _{\mathbf{k}}}^{\prime }{\rho }_{\mathbf{k}2}$$

$$\partial_{t}\rho_{\bf k}\Big{|}_{\rm sc}=-[S_{\bf k}(>,<)-S_{\bf k}(<,>)+h.c.],$$

$$S_{\mathbf{k}}(>,<)$$

$${\rho }_{\mathbf{k}}^{>/<}(t^{\prime })=e^{-i(t^{\prime }-t)H_{\mathbf{k}}}{\rho }_{\mathbf{k}}^{>/<}(t)e^{i(t^{\prime }-t)H_{\mathbf{k}}},$$

$${\rho }_{\mathbf{k}}={\rho }_{\mathbf{k}}^0+{\rho }_{\mathbf{k}}^{\omega }{e}^{i\omega t}+{\rho }_{\mathbf{k}}^{2\omega }{e}^{2i\omega t},$$

$${\rho }_{\mathbf{k}}^0=\frac{1}{2}-\frac{1-2f(E_k)}{2}\left(\frac{{\xi }_k}{E_k}{\tau }_3+\frac{{\Delta }_0}{E_k}{\tau }_1\right).$$

$$\theta$$

$$|\Delta |$$

$${\Delta }_0=-g{\sum _{\mathbf{k}}}^{\prime }{\rho }_{\mathbf{k}1}^0=g{\sum _{\mathbf{k}}}^{\prime }\left[{\Delta }_{0}g(E_k)\right],$$

$${\partial }_t{\rho }_{\mathbf{k}}{|}_{\mathrm{sc}}^{n\omega }$$

$$Y_{\mathbf{k}{\mathbf{k}}^{\prime }}^i(n\omega )$$

$$N_{\mathbf{k}{\mathbf{k}}^{\prime }}^i(n\omega )$$

$$i\omega {\rho }_{\mathbf{k}}^{\omega }+i[{\xi }_k{\tau }_3+{\Delta }_0{\tau }_1,{\rho }_{\mathbf{k}}^{\omega }]+(e{\mathbf{E}}_0\cdot {\partial }_{\mathbf{k}}){\rho }_{\mathbf{k}3}^0{\tau }_0={\partial }_t{\rho }_{\mathbf{k}}{|}_{\mathrm{sc}}^{\omega }.$$

$${\rho }_{\mathbf{k}0}^{\omega }={\rho }_{\mathbf{k}0}^{\omega }{|}_{\mathrm{cl}}-\frac{n_i\pi }{i\omega }\sum _{{\mathbf{k}}^{\prime }}|V_{\mathbf{k}-{\mathbf{k}}^{\prime }}{|}^2{Y}_{\mathbf{k}{\mathbf{k}}^{\prime }}^0(\omega )({\rho }_{\mathbf{k}0}^{\omega }-{\rho }_{{\mathbf{k}}^{\prime }0}^{\omega }),$$

$${\rho }_{\mathbf{k}0}^{\omega }\approx \frac{e{\mathbf{E}}_0\cdot {\mathbf{k}}_F}{im\omega }l(E_k)+\frac{e{\mathbf{E}}_0\cdot {\mathbf{k}}_F}{m{\omega }^2}{\eta }_{\mathbf{k}},$$

$$\frac{{\sigma }_{1s}(\omega )}{{\sigma }_{1n}(\omega )}$$

$${\sigma }_{2s}(\omega )$$

$$2i\omega {\rho }_{\mathbf{k}}^{2\omega }+i[{\xi }_k{\tau }_3+{\Delta }_0{\tau }_1,{\rho }_{\mathbf{k}}^{2\omega }]+i[{\mu }_{\mathrm{eff}}^{2\omega }{\tau }_3+\delta |\Delta {|}^{2\omega }{\tau }_1,{\rho }_{\mathbf{k}}^0]$$

Full text

I Introduction

In the past few decades, the optical properties of the superconducting states have attracted much attention in both linear and nonlinear regimes. The linear response is focused on the behavior of the optical conductivity L1 ; L2 ; L3 ; L4 ; L5 ; L6 ; L7 ; L8 ; NSL0 ; NSL1 ; NSL2 ; NSL3 ; NSL4 ; NSL5 ; NSL6 ; NSL8 ; NSL9 ; NSL10 ; NSL11 ; NSL12 , which was first discussed by Mattis and Bardeen (MB) within the framework of the Kubo current-current correlation approach in the anomalous-skin-effect region MB ; MBo . In this region, the excited current at one space point, depends not only on the electric field at that point but also on the ones nearby. This non-local effect dominates in systems with a small skin depth $\delta$ in comparison with the mean free path $l$, as usually the case in thin-film superconductors or clean type-I superconductors, whereas the scattering effect in this circumstance is marginal. The MB theory suggests that the optical absorption at zero temperature is realized by breaking the Cooper pairs into quasiparticles when the optical frequency $\omega$ is larger than twice the superconducting gap amplitude $|\Delta|$ MB . Thus, the real part of the optical conductivity $\sigma_{1s}(\omega)$ vanishes at $T=0~{}$K when $\omega<2|\Delta|$ but becomes finite above $2|\Delta|$, leading to a crossover point at $2|\Delta|$. At finite temperature, an additional quasiparticle contribution appears below $2|\Delta|$. This theory so far has successfully described the observed data in the anomalous-skin-effect region, as experiments in In L1 , Pb L2 ; L6 ; L7 , Al L8 , thin-film Nb L5 and NbN L3 ; L4 ; NSL0 superconductors demonstrated.

The counterpart of the anomalous-skin-effect region is known as the normal-skin-effect one NSL0 ; NS ($l<\delta$) where the dirty type-II superconductors lie in and the scattering effect becomes important. The optical absorption in the normal-skin-effect region, as experiments in dirty Nb NSL0 ; NSL2 , MgB2 NSL3 ; NSL4 , NbTiN NSL5 ; NSL6 , NbN NSL10 ; NSL8 , MoN NSL11 and Al NSL9 ; NSL12 superconductors, always exhibits a finite $\sigma_{1s}(\omega)$ even at low temperature for $\omega<2|\Delta|$, in contrast to the vanishing $\sigma_{1s}(\omega)$ in the anomalous-skin-effect region. Moreover, with the decrease of $\omega$ in terahertz (THz) regime from $\omega\gg 2|\Delta|$, the observed $\sigma_{1s}(\omega)$ first decreases at $\omega>2|\Delta|$ and then shows an upturn below $2|\Delta|$, leading to a crossover point at $2|\Delta|$. Although the experimental observations are very convincing, theories in the normal-skin-effect region where the scattering effect dominates, are still in progress. The difficulty within the Kubo formalism comes from the inevitable calculation of the vertex correction due to the scattering, which becomes hard to tackle in superconductors NS ; G1 . Whereas the Eilenberger equation is restricted by the normalization condition Eilen ; Ba7 ; Ba8 ; Ba20 , and is also hard to handle for calculation of the scattering. So far, to fit the experimental data, the MB theory derived from the anomalous region is excessively used NSL0 ; NSL2 ; NSL4 ; NSL5 ; NSL6 ; NSL8 ; NSL9 ; NSL10 ; NSL11 ; NSL12 . Nevertheless, such an unphysical fit underestimates $\sigma_{1s}(\omega)$ below $2|\Delta|$ particularly at low temperature where the quasiparticle contribution from MB theory is too small to count for finite experimental result NSL0 ; NSL2 ; NSL4 ; NSL5 ; NSL6 ; NSL8 ; NSL9 ; NSL10 ; NSL11 ; NSL12 . To explain the residual $\sigma_{1s}(\omega)$, several works NSL8 ; NSL9 ; NSL12 considered the influences of the collective gapful Higgs Am1 ; Am2 ; Am3 ; Apm ; Am8 and gapless Nambu-Goldstone Apm ; Am8 ; gi0 ; AK ; Gm1 ; Gm2 ; Ba0 ; pm0 ; pm1 ; pm2 ; pm3 ; gi1 ; Ba9 ; Ba10 (NG) modes which describe the amplitude and phase fluctuations of the order parameter respectively. However, the Higgs mode is charge neutral and does not manifest itself in the linear regime Am1 ; Apm ; Am8 unless under the dc supercurrent injection DCSI . The linear response of the NG mode does not occur either due to its coupling with the long-range Coulomb interaction AK ; Apm ; pm0 ; gi1 ; Ba0 ; Am1 ; gi0 which causes the original gapless energy lifted up to the plasma frequency as a result of Anderson-Higgs mechanism AHM . Therefore, a detailed study capable of clarifying the scattering effect is necessary.

As for the non-linear regime, it was recently realized that through the intense THz pulse, one can excite the oscillation of the superfluid density in the second-order response, which is attributed to the excitation of the Higgs mode NL1 ; NL2 ; NL3 ; NL4 ; NL5 ; NL6 ; NL7 . The most convincing evidence comes from the observed resonance at $2\omega=2|\Delta|$ NL2 ; NL3 ; NL4 , in consistency with the energy spectrum of the Higgs mode Am1 ; Am3 ; Apm . After the THz pulse, a fast damping of this oscillation is observed, and then, a suppressed gap is further observed as a consequence of the thermal effect NL1 ; NL2 ; NL3 . Theory in the literature for these findings is based on Bloch NL3 ; NL4 ; NL5 ; NL6 ; NL7 ; B1 ; B2 ; B3 ; B4 ; B5 ; B6 or Liouville Liou1 ; Liou2 ; Liou3 ; Liou4 equation derived in the Anderson pseudospin picture As . The vector potential ${\bf A}$ naturally involves in this description as a second-order term, which pumps up the fluctuation of the order parameter (pump effect). Nevertheless, the microscopic scattering is absent in the literature. In order to describe the observed damping after the optical pulse, the phenomenological relaxation time is further introduced into the Anderson pseudospin picture NL6 ; NL7 . Very recently, this whole set of approach is challenged. On one hand, this approach with no drive effect fails in the linear regime to give the optical current. On the other hand, symmetry analysis from the Anderson pseudospin picture implies that the pump effect excites the NG mode rather than the observed Higgs mode symmetry . Besides these deficiencies, without the microscopic origin, the introduced phenomenological relaxation mechanism is not exact and convincing.

Very recently, by using the equal-time non-equilibrium $\tau_{0}$-Green function, the gauge-invariant kinetic equation (GIKE) of superconductivity with the microscopic scattering is developed in our previous papers GOBE1 ; GOBE2 ; GOBE3 ; GOBE4 ; GOBE5 . We have proved that the retained gauge invariance in this theory directly leads to the charge conservation in the electromagnetic response GOBE5 , in consistency with Nambu’s conclusion that the gauge invariance in superconductors is equivalent to the charge conservation gi0 . In fact, neither the Bloch NL3 ; NL4 ; NL5 ; NL6 ; NL7 ; B1 ; B2 ; B3 ; B4 ; B5 ; B6 nor Liouville Liou1 ; Liou2 ; Liou3 ; Liou4 equation mentioned above are gauge invariant under the gauge transformation in superconductors gi0 . In contrast, in the GIKE, thanks to the gauge invariance, both pump and drive effects mentioned above are kept GOBE1 ; GOBE2 ; GOBE3 ; GOBE4 ; GOBE5 . Moreover, both superfluid and normal-fluid dynamics are involved in the GIKE GOBE4 ; GOBE5 , beyond the previous Boltzmann equation of superconductors with only the quasiparticle physics retained Ba3 ; Bol ; Ba5 .

Consequently, the well-known clean-limit results such as the Ginzburg-Landau equation and Meissner supercurrent in the magnetic response as well as the optical current captured by the two-fluid model can be directly derived from the GIKE GOBE4 . Particularly, we show that the normal fluid is present only when the excited superconducting velocity $v_{s}$ is larger than a threshold GOBE4 . Moreover, the linear responses of the collective modes from the GIKE also agree with the well-known results in the literature GOBE5 . Whereas the second-order response from the GIKE exhibits interesting physics. On one hand, a finite second-order response of the Higgs mode, attributed solely to the drive effect rather than the widely considered pump effect, is revealed GOBE5 , in contrast to the above theory from Anderson pseudospin picture NL3 ; NL4 ; NL5 ; NL6 ; NL7 ; B1 ; B2 ; B3 ; B4 ; B5 ; B6 ; Liou1 ; Liou2 ; Liou3 ; Liou4 . On the other hand, a finite second-order response of the NG mode, survived from the Anderson-Higgs mechanism, is predicted as a consequence of charge conservation. An experimental scheme for this response is further proposed GOBE5 . Actually, thanks to the equal-time scheme, the microscopic scattering in superconductors, which is hard to deal with in the literature as mentioned above, becomes easy to handle within the GIKE approach. Thus, rich physics from the scattering can be expected. Particularly, at low frequency (i.e., large $v_{s}$), we have analytically shown that due to the scattering, there exists viscous superfluid besides the non-viscous one GOBE4 . Then, together with the normal fluid, a three-fluid model is proposed GOBE4 .

In this work, by extending the previous scattering terms in Ref. GOBE4, into the THz regime via carefully implementing the Markovian approximation, we further apply the GIKE to investigate the influence of the scattering on the optical properties of superconductors in the normal-skin-effect region ($l<\delta$). Both linear and second-order responses are analytically studied under a multi-cycle THz pulse. In the linear regime, we show that the optical absorption $\sigma_{1s}(\omega)$, induced by the scattering, always exhibits a finite value even at low temperature when $\omega<2|\Delta|$, in contrast to the vanishing $\sigma_{1s}(\omega)$ in the anomalous-skin-effect region as MB theory revealed MB . Moreover, with the decrease of the optical frequency from $\omega\gg 2|\Delta|$, $\sigma_{1s}(\omega)$ first increases and then drops abruptly around $2|\Delta|$. By further decreasing $\omega$ below $2|\Delta|$, an upturn of $\sigma_{1s}(\omega)$ is observed, leading to a crossover point at $2|\Delta|$. In the second-order regime, responses of the Higgs mode during and after the optical pulse are revealed. During the pulse, it is found that the scattering causes a phase shift in the optical response of the Higgs mode. Particularly, this phase shift exhibits a significant $\pi$-jump at $\omega=|\Delta|$, which provides a very clear feature for the experimental detection. After the pulse, the damping of the Higgs-mode excitation is studied. In this situation, we reveal a relaxation mechanism due to the elastic scattering, which shows a monotonic enhancement with the increase of the impurity density.

This paper is organized as follows. We first present the GIKE of superconductivity in Sec. II. Then, we perform the analytic analysis of the influence from the scattering on the optical properties of superconductors in Sec. III. We summarize in Sec. IV.

II MODEL

In this section, we first introduce the complete GIKE. Then, we present a simplified GIKE to study the optical response of superconductors in the normal-skin-effect region. The microscopic scattering terms of the non-magnetic impurity scattering are also addressed in this section.

II.1 GIKE

The GIKE of the $s$-wave BCS superconductors, which is developed in our previous papers GOBE4 ; GOBE5 , reads:

[TABLE]

Here, $[~{},~{}]$ and $\{~{},~{}\}$ represent the commutator and anti-commutator, respectively; $\xi_{k}=\frac{k^{2}}{2m}-\mu$ with $m$ and $\mu$ being the effective mass and chemical potential; ${\hat{\Delta}(R)}=\Delta(R)\tau_{+}+\Delta^{*}(R)\tau_{-}$; $R=(t,{\bf R})$ stands for the center-of-mass coordinate; $\tau_{i}$ are the Pauli matrices in the particle-hole space; ${\phi}$ and ${\bf A}$ denote the scalar and vector potentials, respectively; $\rho^{c}_{\bf k}$ is the density matrix in the Nambu space; on the right-hand side of Eq. (1), the scattering term \partial_{t}\rho^{c}_{\bf k}\Big{|}_{\rm sc} is added for the completeness, whose explicit expression is shown in Sec. II.3.

The superconducting order parameter $\Delta$, Fock field $\mu_{F}$ and Hartree field $\mu_{H}$ in Eq. (1) are written as

[TABLE]

where $\delta{n(R)}$ represents the density fluctuation; $V_{\bf R-R^{\prime}}$ denotes the Coulomb potential whose Fourier component $V_{\bf q}=e^{2}/(q^{2}\epsilon_{0})$; $g$ stands for the effective electron-electron attractive potential in the BCS theory BCS . $\sum^{\prime}_{\bf k}$ here and hereafter represents the summation restricted in the spherical shell ($|\xi_{k}|<\omega_{D}$) with $\omega_{D}$ being the Debye frequency BCS .

The effective electric field ${\bf E}$ in Eq. (1), as a gauge-invariant measurable quantity, is given by

[TABLE]

The gauge-invariant density $n$ and current ${\bf j}$ read GOBE4 :

[TABLE]

We emphasize that Eq. (1) is gauge-invariant under the gauge transformation first revealed by Nambu gi0 :

[TABLE]

where the four vectors $A_{\mu}=(\phi,{\bf A})$ and $\partial_{\mu}=(\partial_{t},-{\bm{\nabla}}_{\bf R})$; $\theta$ denotes the phase of the superconducting order parameter. Thanks to the retained gauge invariance, the charge conservation:

[TABLE]

is naturally satisfied during the electromagnetic response as we proved in our latest work GOBE5 . This agrees with the Nambu’s conclusion via the Ward’s identity that the gauge invariance in the superconducting states is equivalent to the charge conservation gi0 . Moreover, due to the gauge invariance, both the pump [third term in Eq. (1)] and drive [sixth and seventh terms in Eq. (1)] effects mentioned in the introduction are kept.

II.2 Simplified GIKE in normal-skin-effect region

In this part, we present a simplified GIKE in the normal-skin-effect region. We first choose a specific gauge by transforming Eq. (1) under the gauge transformation $\rho_{\bf k}(R)={e^{-i\tau_{3}\theta(R)/2}}\rho^{c}_{\bf k}(R)e^{i\tau_{3}\theta(R)/2}$. Then, under a spatially uniform (i.e., long-wave-limit) optical field in the normal-skin-effect region, the spatial gradient terms in the kinetic equation can be neglected. Consequently, Eq. (1) becomes:

[TABLE]

with the gauge-invariant superconducting momentum ${\bf p}_{s}$ and effective field $\mu_{\rm eff}$ written as

[TABLE]

Moreover, by expanding the density matrix as $\rho_{\bf k}=\sum^{3}_{i=0}\rho_{{\bf k}i}\tau_{i}$, the gap equations [Eqs. (2) and (3)] correspondingly read:

[TABLE]

As shown in our latest work GOBE5 , Eq. (15) gives the gap equation, from which one can self-consistently obtain the Higgs mode. The NG mode can be self-consistently determined by Eq. (16). Moreover, under the uniform optical response, one finds that ${\bm{\nabla}}_{\bf R}\cdot{\bf j}=0$. Therefore, as a consequence of the charge conservation [Eq. (11)], the density fluctuation $\delta{n}$ and hence both the Hartree $\mu_{H}$ and Fock $\mu_{F}$ fields vanish.

II.3 Microscopic Scattering

We next present the scattering terms $\partial_{t}\rho_{\bf k}|_{\rm sc}$ in Eq. (12) which are derived based on the generalized Kadanoff-Baym ansatz spintronic ; DS1 ; DS2 ; GKB . Considering the fact that the electron-phonon scattering is weak at low temperature, we mainly consider the electron-impurity scattering. The specific impurity scattering terms read (detailed derivation can be found in Refs. DS1, ; spintronic, ; DS2, ):

[TABLE]

with

[TABLE]

Here, $\rho^{<}_{\bf k}=\rho_{\bf k}$ and $\rho^{>}_{\bf k}=1-\rho_{\bf k}$; $H_{\bf k}=\xi_{{\bf k}+m{\bf v}_{s}\tau_{3}}\tau_{3}+|\Delta|\tau_{1}$ denotes the BCS Hamiltonian in the presence of the superconducting velocity ${\bf v}_{s}$; $n_{i}$ is the impurity density; $U_{\bf kk^{\prime}}=V_{\bf k-k^{\prime}}\tau_{3}$ stands for the electron-impurity interaction in the Nambu space. This scattering term [Eq. (17)] is non-Markovian.

It is well established in semiconductor optics DS2 and spintronics spintronic that the clean-limit solution of the corresponding kinetic (i.e., Liouville) equation:

[TABLE]

is substituted into the scattering terms as the Markovian approximation to obtain the conventional energy conservation in the scattering. In our previous works GOBE1 ; GOBE2 ; GOBE3 ; GOBE4 , we also take such approach in Eq. (18) to derive the scattering in superconductors. In the present work, this approach is sublated in the presence of the multi-cycle THz optical field, since the free coherent oscillation in this circumstance does not hold, i.e., Eq. (19) is no longer the clean-limit solution of the GIKE in superconductors [Eq. (12)]. In fact, as shown in the next section, during the multi-cycle THz pulse, the response of the density matrix is forced to oscillate with the multiples of the optical frequency.

III Analytic Analysis

In this section, by solving the simplified GIKE [Eq. (12)] in the normal-skin-effect region, we analytically investigate the scattering effect in the optical response of superconductors under multi-cycle THz pulse. In this circumstance, analytic analyses for two extreme cases: during and after the pulse, are performed to carefully handle the Markovian approximation in order to turn the non-Markovian scattering in Eq. (17) into the Markovian one. The multi-cycle THz pulse, as applied in recent experiments NL7 , possesses a stable phase as well as a narrow frequency bandwidth. Consequently, during the optical pulse, the system is under a periodic drive scheme at a well-defined frequency, similar to the case under a continuous waveform field. In this situation, the response of the superconductivity is forced to oscillate with the multiples of the optical frequency. Whereas after the optical pulse, the system is free from the optical field, and the study in this situation reveals the relaxation mechanism of the optically excited non-equilibrium states.

III.1 Forced oscillation

During the multi-cycle THz pulse, by assuming the electromagnetic potential $\phi=\phi_{0}({\bf R})e^{i\omega{t}}$ and ${\bf A}={\bf A}_{0}e^{i\omega{t}}$, the density matrix $\rho_{\bf k}$ reads:

[TABLE]

with the equilibrium-state density matrix $\rho^{0}_{\bf k}$ given by GOBE1 ; GOBE4 ; GOBE5

[TABLE]

Here, $\rho^{\omega(2\omega)}_{\bf k}$ denotes the linear (second-order) response of the density matrix; $E_{k}=\sqrt{\xi_{k}^{2}+\Delta_{0}^{2}}$; $f(x)$ represents the Fermi-distribution function.

Correspondingly, the responses of the phase $\theta$ and amplitude $|\Delta|$ of the superconducting order parameter are written as

[TABLE]

From Eqs. (15) and (21), with $g(E_{k})=\frac{1-2f(E_{k})}{2E_{k}}$, the equilibrium-state order parameter $\Delta_{0}$ is determined by

[TABLE]

which is exactly the gap equation in the BCS theory BCS . Moreover, as shown in our latest work GOBE5 , the Higgs mode dose not manifest itself in the linear regime ($\delta|\Delta|^{\omega}=0$). The linear response of the NG mode from the GIKE GOBE5 , due to its coupling to the long-range Coulomb interaction, does not effectively occur either (i.e., $\mu_{\rm eff}^{\omega}=0$ and ${\bf p}^{\omega}_{s}=-2e{\bf A}_{0}^{\perp}$ with ${\bf A}_{0}^{\perp}$ being the physical transverse vector potential) as a result of the Anderson-Higgs mechanism AHM , in agreement with the previous works in the literature AK ; Apm ; pm0 ; gi1 ; Ba0 ; Am1 ; gi0 .

Furthermore, it is noted that in the presence of the multi-cycle THz pulse, the response of the density matrix [Eq. (20)], as the solution of the kinetic equation, is forced to oscillate with the multiples of the optical frequency, rather than the free coherent oscillation mentioned above. Then, substituting this forced oscillation [Eq. (20)] into the scattering term [Eq. (17)], the $n$-th order of the scattering during the optical pulse can be obtained (refer to Appendix A):

[TABLE]

with

[TABLE]

Here, $\eta=\pm$; the projection operators $\Gamma^{\pm}_{k}$ are written as $\Gamma^{\pm}_{k}=U^{\dagger}_{k}Q^{\pm}U_{k}$ with $Q^{\pm}=(1\pm\tau_{3})/{2}$ and $U_{k}=u_{k}\tau_{0}-v_{k}\tau_{+}+v_{k}\tau_{-}$ being the unitary transformation matrix from the particle space to the quasiparticle one. $u_{k}=\sqrt{1/2+\xi_{k}/(2E_{k})}$ and $v_{k}=\sqrt{1/2-\xi_{k}/(2E_{k})}$; $E_{\bf k}^{\pm}={\bf k}\cdot{\bf v}_{s}\pm{E_{k}}$ denotes the tilted quasiparticle energy. It is noted that at low frequency $\omega{\ll}\Delta_{0}$, the scattering term in Eq. (25) recovers the one in our previous work where we propose the three-fluid model as mentioned in the introduction GOBE4 . In the present work for the optical properties, we focus on the THz regime where $\omega{\sim}\Delta_{0}$. Moreover, considering a weak and fast-oscillating optical field, the tilt in quasiparticle energy (i.e., Doppler shift ${\bf k}\cdot{\bf v}_{s}$), related to electromagnetic field GOBE4 , can be neglected (i.e., $E_{\bf k}^{\pm}=\pm{E_{k}}$).

Then, as seen from Eq. (25), due to the forced oscillation of the density matrix by the influence of the multi-cycle THz pulse, the optical frequency $\omega$ is involved in $\delta(E_{\bf k^{\prime}}^{\eta_{1}}+n\omega-E_{\bf k}^{\eta_{2}})$ (i.e., the energy conservation of the scattering). Consequently, besides the intraband scattering ($\eta_{1}=\eta_{2}$), the interband scattering channel ($\eta_{1}=-\eta_{2}$) is opened.

III.1.1 Linear response: optical conductivity

We first investigate the optical conductivity in the linear regime. The linear order of the GIKE [Eq. (12)] reads:

[TABLE]

From above equation, it is noted that only the $\tau_{0}$ component of $\rho^{\omega}_{{\bf k}}$ is optically excited:

[TABLE]

and the other components of $\rho^{\omega}_{{\bf k}}$ are zero, in consistency with the above mentioned vanishing $\delta|\Delta|^{\omega}$ [Eq. (15)] and $\mu_{\rm eff}^{\omega}$ [Eq. (16)]. Here, $\rho^{\omega}_{{\bf k}0}|_{\rm cl}=\frac{e{\bf E}_{0}\cdot{\bf k}}{im\omega}l(E_{k})$ is the clean-limit solution with $l(E_{k})=\partial_{\xi_{k}}[\xi_{k}g(E_{k})]=-\frac{\Delta_{0}^{2}}{E_{k}}\partial_{E_{k}}g(E_{k})-\partial_{E_{k}}f(E_{k})$ consisting of superfluid [$-\frac{\Delta_{0}^{2}}{E_{k}}\partial_{E_{k}}g(E_{k})$] and quasiparticle [$-\partial_{E_{k}}f(E_{k})$] contributions, exactly same as the one in our previous works GOBE4 ; GOBE5 . The second term on the right-hand side of Eq. (29) comes from the scattering.

The exact analytic solution of $\rho^{\omega}_{{\bf k}0}$ from Eq. (29) is difficult in the presence of the scattering. Nevertheless, at the relatively weak scattering (i.e., $\xi<l$ with $\xi$ being the coherence length), after the first-order iteration by substituting $\rho^{\omega}_{{\bf k}0}|_{\rm cl}$ into the scattering term [second term on the right-hand side of Eq. (29)], $\rho^{\omega}_{{\bf k}0}$ can be directly solved:

[TABLE]

with $\eta_{\bf k}=n_{i}\pi\sum_{\bf k^{\prime}}|V_{\bf k-k^{\prime}}|^{2}Y^{0}_{\bf kk^{\prime}}(\omega)[l(E_{k})-\cos\theta_{\bf kk^{\prime}}l(E_{k^{\prime}})]$.

Then, substituting the solved $\rho^{\omega}_{{\bf k}0}$ into Eq. (8), the optical conductivity in the superconducting states $\sigma_{s}(\omega)=\sigma_{1s}(\omega)+i\sigma_{2s}(\omega)$ is obtained (refer to Appendix B):

[TABLE]

where $\theta(x)$ is the step function; $\sigma_{1n}(\omega)=\frac{ne^{2}}{m\omega^{2}\tau_{p}}$ with $\frac{1}{\tau_{p}}=\Gamma_{0}-\Gamma_{1}$ exactly being the momentum relaxation rate in normal metals and $\Gamma_{i}=2n_{i}{\pi}D\int\frac{d\Omega_{\bf k^{\prime}}}{4\pi}|V_{\bf k_{F}-k_{F}^{\prime}}|^{2}\cos{i}\theta_{\bf kk^{\prime}}$. $D$ is the density of states. It is noted that the first term in $\sigma_{2s}(\omega)$ recovers the clean-limit one in the superfluid as revealed in our previous work GOBE4 .

Firstly, we point out that the obtained optical conductivity from the GIKE [Eqs. (31) and (32)] becomes $\sigma_{n}=\frac{ne^{2}}{m\omega^{2}\tau_{p}}+\frac{ne^{2}}{im\omega}$ in the normal states at $T>T_{c}$ with $\Delta_{0}=0$ (refer to Appendix C), exactly recovering the one in normal metals as the Drude model or conventional Boltzmann equation revealed. To the best of our knowledge, so far there is no theory of the optical conductivity in the literature that can rigorously recover the conductivity in normal metals from $T<T_{c}$ to $T>T_{c}$, due to the difficulty in calculating the vertex correction in superconductors NS ; G1 ; vertex . The GIKE here actually provides an efficient approach to deal with the scattering.

We then discuss the frequency dependence of the optical absorption $\sigma_{1s}(\omega)$ in the superconducting states. In Eq. (31), the first term originates from the intraband scattering. Whereas the second one comes from the interband scattering, leading to the step function. The frequency dependence of $\sigma_{1s}(\omega)$ is plotted in Fig. 1. As seen from the figure, $\sigma_{1s}(\omega)$ shows a significant crossover at $\omega=2\Delta_{0}(T)$, which comes from the step function (i.e., opened interband-scattering channel) for $\omega>2\Delta_{0}$ in Eq. (31). Secondly, at $T=0~{}$K with the finite superfluid contribution $l(E_{k})=\frac{\Delta_{0}^{2}}{2E_{k}^{3}}$ in Eq. (31), one finds that $\sigma_{1s}(\omega)$, shown by the solid curve in Fig. 1, always exhibits a finite value even when $\omega<2\Delta_{0}$, in sharp contrast to the vanishing $\sigma_{1s}(\omega)$ in the anomalous-skin-effect region as MB theory revealed. Moreover, as shown in Fig. 1, with the decrease of $\omega$ from $\omega\gg 2\Delta_{0}$, $\sigma_{1s}(\omega)$ first increases and then drops abruptly around $2\Delta_{0}$. By further decreasing $\omega$ below $2\Delta_{0}$, due to the fast increase of $\sigma_{1n}(\omega)$ in Eq. (31), a significant upturn of $\sigma_{1s}(\omega)$ is observed.

Results in the dirty limit ($l<\xi$) require a full numerical calculation of Eq. (29) and go beyond the analytic analysis. Nevertheless, from Eq. (29), thanks to the finite value of superfluid contribution in $l(E_{k})$ at $T=0~{}$K, the finite $\sigma_{1s}(\omega)$ at low temperature when $\omega<2\Delta_{0}$ is unlikely changed even in the dirty limit. In addition, due to the existence of $\delta(\omega-E_{k}+E_{k^{\prime}})$ in $Y^{0}_{\bf kk^{\prime}}(\omega)$ [Eq. (26)], the crossover point at $\omega=2\Delta_{0}(T)$ can also be obtained in the dirty limit. These two points, by the full numerical calculation of Eq. (29) in the dirty limit, are justified (refer to Appendix D), in qualitative agreement with the experimental findings NSL1 ; NSL2 ; NSL10 ; NSL3 ; NSL4 ; NSL0 ; NSL11 ; NSL5 ; NSL6 ; NSL8 ; NSL9 ; NSL12 .

As mentioned in the introduction, the MB theory derived from anomalous-skin-effect region is excessively used in the literature NSL1 ; NSL2 ; NSL4 ; NSL5 ; NSL6 ; NSL8 ; NSL9 ; NSL10 ; NSL11 ; NSL12 to fit the experimental data in the normal-skin-effect region where the scattering effect dominates. Nevertheless, as shown by the dotted curve in Fig. 1, at low temperature, when $\omega<2\Delta_{0}$, $\sigma_{1s}(\omega)$ from MB theory derived at the anomalous-skin-effect region, which comes from the quasiparticle contribution MB , is too small in comparison with the finite experimental observation NSL1 ; NSL2 ; NSL4 ; NSL5 ; NSL6 ; NSL8 ; NSL9 ; NSL10 ; NSL11 ; NSL12 . Therefore, such an unphysical fit underestimates the upturn of $\sigma_{1s}(\omega)$ below $2\Delta_{0}$ particularly at low temperature, and hence, is incapable of capturing the experimental findings NSL1 ; NSL2 ; NSL3 ; NSL4 ; NSL5 ; NSL6 ; NSL8 ; NSL9 ; NSL10 ; NSL11 ; NSL12 .

III.1.2 Second-order response: excitation of Higgs mode

We next investigate the second-order response of the Higgs mode. The second-order GIKE is written as

[TABLE]

from which $\rho^{2\omega}_{\bf k}$ can be analytically solved at the relatively weak scattering.

Substituting the solved $\rho^{2\omega}_{{\bf k}1}$ into Eq. (15), the second-order response of the Higgs mode can be self-consistently derived (refer to Appendix E):

[TABLE]

where

[TABLE]

with $d(E)=\frac{\partial_{E}l(E)}{E}$, $o(E)=\frac{\sqrt{E^{2}-\Delta^{2}_{0}}}{E^{2}-\omega^{2}}$ and functional function

[TABLE]

It is noted that in the absence of the scattering (i.e., $\Gamma_{0}=0$), Eq. (34) exactly reduces to the clean-limit one revealed in our latest work GOBE5 .

As seen from Eq. (34), $\gamma_{H}$ from the scattering causes the broadening of the Higgs-mode spectrum whereas $s_{H}$ represents the second-order optical absorption through the scattering. The existences of $s_{H}$ and $\gamma_{H}$ result in an imaginary part in the second-order response of the Higgs mode, and hence, lead to a phase shift in this response. The magnitude $A^{2\omega}_{H}(\omega)$ and phase shift $\phi(\omega)$ of the second-order response of the Higgs mode $\delta|\Delta|^{2\omega}=A_{H}^{2\omega}e^{i\phi(\omega)}$ are plotted in Fig. 2(a) and (b), respectively. As seen from Fig. 2(a), the magnitude of the second-order response of the Higgs mode exhibits a resonant peak at $2\omega=2\Delta_{0}(T)$, in consistency with the experimental observation NL2 ; NL3 ; NL4 . The phase shift $\phi$ of this second-order response [Fig. 2(b)] exhibits a $\pi$-jump at $\omega=\Delta_{0}(T)$. This is natural since from Eq. (34), the real part of $\delta|\Delta|^{2\omega}$ at the weak scattering is proportional to $(\omega^{2}-\Delta^{2}_{0})^{-1}$ whereas the imaginary one is proportional to $(\omega^{2}-\Delta^{2}_{0})^{-2}$, leading to $\tan\phi\propto(\omega^{2}-\Delta^{2}_{0})^{-1}$.

III.2 Free decay

In the previous subsection, we have investigated the response of the superconducting states during the optical pulse. In this part, we focus on the situation of the temporal evolution of the optically excited collective modes after the optical pulse.

III.2.1 Simplified model

The GIKE after the optical pulse is written as

[TABLE]

The density matrix is given by

[TABLE]

where $\delta\rho_{\bf k}$ denotes the part deviated from the equilibrium state due to the optical excitation. The fluctuations of the amplitude (i.e., $\delta|\Delta|$) and phase (i.e., $\mu_{\rm eff}$) of the order parameter can be obtained from Eqs. (15) and (16), respectively.

It is noted that in Eq. (39), the second term on the left-hand side causes the coherent oscillation of the density matrix whereas the one on the right-hand side provides the scattering. In this circumstance, as established in the semiconductor optics DS2 and spintronics spintronic , Eq. (19) as a clean-limit solution of Eq. (39), can be safely used into Eq. (18) as the Markov approximation to further derive the scattering terms. Then, the scattering which becomes free from the influence from the optical frequency, is given by

[TABLE]

where $w_{kk^{\prime}}=\sum_{\eta}\Gamma^{\eta}_{k}\tau_{3}\Gamma^{\eta}_{k^{\prime}}=w^{1}_{kk^{\prime}}\tau_{1}+w^{3}_{kk^{\prime}}\tau_{3}$ with $w^{1}_{kk^{\prime}}=\frac{\Delta_{0}(\xi_{k^{\prime}}+\xi_{k})}{2E_{k}E_{k^{\prime}}}$ and $w^{3}_{kk^{\prime}}=u_{k^{\prime}}^{2}u_{k}^{2}+v^{2}_{k^{\prime}}v^{2}_{k}-2u_{k}v_{k}u_{k^{\prime}}v_{k^{\prime}}$.

Since only the isotropic part of the density matrix in the momentum space survives the summation in Eqs. (15) and (16), i.e., contributes to the calculations of the amplitude and phase of the order parameter, we neglect the anisotropic part in $\delta\rho_{\bf k}$. Then, considering the fact $w_{kk^{\prime}}|_{\xi_{k}=-\xi_{k^{\prime}}}=0$, the scattering term in Eq. (41) is simplified after the summation of ${\bf k^{\prime}}$, and the GIKE becomes

[TABLE]

Particularly, it is pointed out that Eq. (42) in the Anderson pseudospin picture As is written as

[TABLE]

where ${\bf b}_{k}=(\Delta_{0}+\delta\Delta,0,\xi_{k}+\mu_{\rm eff})$ and ${\bf s}_{k}=(\rho_{k1},\rho_{k2},\rho_{k3})$ denote the Anderson pseudo field and spin, respectively; ${\bf a}_{1}=(0,1,0)$ and ${\bf a}_{2}=(\xi_{k}/E_{k},0,-\Delta_{0}/E_{k})$ are two transverse directions to the equilibrium-state pseudo field ${\bf b}^{0}_{k}$. It is noted that in Eq. (43), the second term on the left-hand side of the equation causes the coherent precession of the Anderson pseudospin, exactly same as the one in the previous works NL3 ; NL4 ; NL5 ; NL6 ; NL7 ; B1 ; B2 ; B3 ; B4 ; B5 ; B6 . The terms on the right-hand side come from the scattering, which provide the relaxation of the non-equilibrium states. Particularly, since $s^{x}_{k}$ and $s^{y}_{k}$ contribute to the calculations of the Higgs [Eq. (15)] and NG [Eq. (16)] modes separately, one immediately finds that the first term on the right-hand side of Eq. (43) provides the damping of the excited Higgs mode whereas the second term causes the damping of the NG mode.

We point out that in the present work, the relaxation terms on the right-hand side of Eq. (43), exactly come from the microscopic scattering, differing from and going beyond the previous phenomenological relaxation in the Anderson pseudospin picture mentioned in the introduction NL6 ; NL7 . In fact, the previous phenomenological relaxation mechanism, without the microscopic origin, is not exact and convincing. Specifically, in Ref. NL6, , in analogy with the real spin precession, the longitudinal and transverse relaxation processes, which describe the damping of the components of $\delta{\bf s}_{k}$ along and perpendicular ${\bf b}^{0}_{k}$, are introduced into the Anderson pseudospin picture through the phenomenological relaxation time. Nevertheless, one finds that the longitudinal component of the pseudospin $\delta{\bf s}_{k}\cdot{\bf b}^{0}_{k}/E_{k}=(\Delta_{0}\delta\rho_{k1}+\xi_{k}\delta\rho_{k3})/E_{k}=\delta\rho^{q}_{k3}$. Since the diagonal $\delta\rho^{q}_{k3}$ is related to the quasiparticle distribution, the longitudinal relaxation process [i.e., terms like ($\delta{\bf s}_{k}\cdot{\bf b}^{0}_{k}$)] directly describes the damping of the quasiparticles in which only the inelastic scattering contributes and the elastic scattering makes no contribution at all. Hence, in superconductors, considering the weak inelastic electron-phonon scattering at low temperature, the longitudinal relaxation process is marginal and only the transverse ones [i.e., terms like ($\delta{\bf s}_{k}\cdot{\bf a}_{1}$) and ($\delta{\bf s}_{k}\cdot{\bf a}_{2}$)] play the important role. Particularly, there is no reason for the two transverse relaxation processes, which provide the damping of the two collective modes separately as mentioned above, to share the same rate. Most importantly, since ${\delta}s^{z}_{k}$ is related to the density fluctuation [i.e., $\delta{n}=\sum_{\bf k}{\delta}s^{z}_{k}=0$ from Eq. (7)], as a consequence of the charge conservation, the relaxation terms should not have any component along $z$ direction. All above features, unsatisfied in Ref. NL6, , are well kept in our relaxation terms in Eq. (43) , thanks to the microscopic scattering in the GIKE.

III.2.2 Damping of Higgs mode

By taking the optical response of the density matrix $\delta\rho_{\bf k}(t=0)=\rho^{\omega}_{\bf k}+\rho^{2\omega}_{\bf k}$, we first perform the numerical calculation to self-consistently solve Eq. (42) with Eqs. (15) and (16). Then, the temporal evolution of the Higgs mode $\delta|\Delta|(t)$ and NG mode $\mu_{\rm eff}(t)$ can be self-consistently obtained. We focus on the measurable Higgs mode in this part.

The temporal evolution of the Higgs mode after the optical pulse is plotted in Fig. 3 at different scattering rates. As seen from the figure, $\delta|\Delta|(t)$ exhibits an oscillatory decay behavior, in consistency with the experimental observation NL1 ; NL2 ; NL3 ; NL4 ; NL6 ; NL7 . The frequency of the oscillation is around $2\Delta_{0}$, in agreement with the energy spectrum of the Higgs mode. Moreover, it is also found that the damping of $\delta|\Delta|(t)$ shows a monotonic enhancement with the increase of the scattering rate.

To further understand the temporal evolution of $\delta|\Delta|(t)$, we analytically derive the solution of Eq. (42) by first transforming Eq. (42) into the quasiparticle space through the unitary transformation $\rho^{q}_{\bf k}=U_{k}\rho_{k}U_{k}^{\dagger}$. Then, under a weak excitation (i.e., small $\delta\rho_{\bf k}$), one has the components of the equation:

[TABLE]

with $\gamma_{k}=\Gamma_{0}{\rm sgn}(\xi_{k})\frac{2\xi_{k}}{E_{k}}$ and $a_{k}=(-\frac{\Delta_{0}}{E_{k}}\mu_{\rm eff}+\delta|\Delta|\frac{\xi_{k}}{E_{k}})\rho_{k3}^{q0}$.

An exact solution from above equations is difficult. However, at the weak scattering, similar to the Elliot-Yafet mechanism in the spin relaxation of the semiconductor spintronics spintronic , the coupling terms between $\rho^{q}_{k3}$ and $\rho^{q}_{k\pm}$ in Eq. (46) can be effectively removed through the unitary transformation as the Löwdin partition method showed dia . Then, $\rho^{q}_{\bf k}(t)$ and hence $\rho_{\bf k}(t)$ can be solved (refer to Appendix F). Consequently, from the gap equation [Eq. (15)], the temporal-evolution equation of the excited Higgs mode is given by:

[TABLE]

with $a_{c}=\sqrt{c^{2}_{1}+c^{2}_{2}}$ and $\tan\theta_{c}=c_{2}/c_{1}$ and $\delta{t}=t-t^{\prime}$. The coefficients $c_{ki}$ are determined by the initial optical excitation:

[TABLE]

As seen from the right-hand side of Eq. (47), the first and second terms are related with the initial excitation; By only considering the third term, one has $\partial_{t}\delta|\Delta|=-\left[g\sum_{\bf k}\gamma_{k}g(E_{k})\frac{\Delta_{0}^{2}}{E_{k}^{2}}\right]\delta|\Delta|$. Thus, the third term on the right-hand side of Eq. (47) causes the damping of $\delta|\Delta|$ with the relaxation rate proportional to $\gamma_{k}$. The last term show the oscillatory decay with the time evolution, and hence, directly lead to the oscillating damping of $\delta|\Delta|$ with the relaxation rate proportional to $\gamma_{k}$. The relaxation rate of the Higgs mode therefore increases by increasing the impurity density, similar to the Elliot-Yafet mechanism in the spin relaxation of the semiconductor spintronics spintronic . Comparisons between the analytic solution [Eq. (47)] and full numerical results are plotted in the insets of Fig. 3, where the results from the two sets of calculations agree well with each other.

Finally, from Eq. (47), it is found that the long-time dynamic of the Higgs mode behaves as (refer to Appendix G)

[TABLE]

exhibiting an oscillatory decay behavior with oscillating frequency at the Higgs-mode energy $2\Delta_{0}$. Here, ${\bar{\gamma}}$ is the average of ${\gamma_{\bf k}}$ in the momentum space. In the absence of disorder (${\bar{\gamma}}=0$), Eq. (51) reduces to the previous coherent BCS oscillatory decay B2 ; B3 ; OD1 ; OD2 ; OD3 ; OD4 as it should be, since our kinetic equation [Eq. (43) or Eq. (42)] without the scattering exactly recovers the linearized Bloch (i.e., Anderson-pseudospin) equations around the equilibrium state B2 ; B3 ; OD3 ; OD4 . Whereas the presence of the impurity leads to exponential decay.

IV SUMMARY AND DISCUSSION

Within the GIKE approach, we analytically investigate the influence of the scattering on the optical response of superconductors in the normal-skin-effect region ($l<\delta$). Two extreme situations: during and after a multi-cycle THz pulse pulse, are considered with a careful implementation of the Markovian approximation for the microscopic scattering. During the pulse, the multi-cycle optical field with the stable phase and narrow frequency bandwidth as applied in recent experiments NL7 , exhibits the continuous-wave-like behavior. Then, response of the density of matrix, as the solution of the free GIKE in superconductors, is forced to oscillate with the multiples of the optical frequency. Consequently, due to this forced oscillation, after the Markovian approximation, the energy conservation of the scattering is influenced by the optical frequency. Whereas after the optical pulse, the system is free from the optical field, and the density of matrix in this situation exhibits the free coherent oscillation in the clean limit. Then, after the Markovian approximation, the energy conservation of the scattering becomes free from the influence from the optical frequency. Rich physics in both extreme cases is revealed.

Specifically, during the pulse, responses of the superconductivity in linear and second-order regimes are studied. In the linear regime, we analytically derive the optical conductivity from the GIKE at the weak scattering ($l>\xi$). We show that by taking $T>T_{c}$ the optical conductivity from our theory obtained at $T<T_{c}$ exactly recovers the one in normal metals as the Drude model or conventional Boltzmann equation revealed. To the best of our knowledge, so far there is no theory in the literature that can rigorously make this recovery. Whereas in the superconducting states, we find that the optical absorption $\sigma_{1s}(\omega)$, due to the contribution of superfluid density, always exhibits a finite value when $\omega<2\Delta_{0}$ even at low temperature, and shows an upturn with the decrease of frequency below $2\Delta_{0}$, in contrast to the vanishing $\sigma_{1s}(\omega)$ in the anomalous-skin-effect region as MB theory revealed MB . Moreover, $\sigma_{1s}(\omega)$ shows a significant crossover at $\omega=2\Delta_{0}(T)$, which comes from opened interband-scattering channel for $\omega>2\Delta_{0}$. Through the full numerical calculation, we further show that both the upturn of the finite $\sigma_{1s}(\omega)$ below $2\Delta_{0}$ and the crossover point at $\omega=2\Delta_{0}(T)$ in $\sigma_{1s}(\omega)$ also appear in the dirty-limit regime ($\xi<l$), in qualitative agreement with the experimental observations in disordered type-II superconductors like Nb NSL1 ; NSL2 , NbN NSL10 , MgB2 NSL3 ; NSL4 ; NSL0 ; NSL11 , NbTiN NSL5 ; NSL6 ; NSL8 and Al NSL9 ; NSL12 .

As for the second-order regime, we study the response of the Higgs mode. We show that the scattering causes a phase shift in this second-order optical response. Particularly, we find that this phase shift exhibits a significant $\pi$-jump at $\omega=\Delta_{0}$, which provides a very clear feature for the experimental detection. Recently, thanks to the advanced pump-probe technique, a $\pi$-jump of the phase shift has been experimentally observed at $\omega=\Delta_{0}(T)$ in the second-order optical response of the disordered high-$T_{c}$ cuprates-based superconductors NL7 . The origin of this jump is still controversial. Whereas our present work suggests that the $\pi$-jump of the phase shift in the second-order optical response can also be realized in the conventional superconductors through the scattering effect.

Finally, we study the relaxation mechanism of the excited collective modes after the pulse. In this situation, based on the complete GIKE, a simplified model with the damping terms in the Anderson pseudospin picture is proposed. The damping terms in this model exactly come from the microscopic scattering, differing from and going beyond the phenomenological relaxation mechanism in the previous works NL6 ; NL7 . Particularly, both the charge conservation and the unique feature of the dominant elastic scattering in superconductors: vanishing longitudinal relaxation process, are kept in our relaxation terms, in sharp contrast to Ref. NL6, . Then, by studying the damping of the Higgs-mode excitation, we reveal an exponential relaxation mechanism due to the elastic scattering, which shows a monotonic enhancement with the increase of the impurity density. In addition, we also investigate the damping of the NG mode (refer to Appendix H). It is found that in the conventional BCS superconductors, the damping of the phase fluctuation (NG mode) is much faster than that of the amplitude fluctuation (Higgs mode) of the order parameter.

Note added: After the completion of our manuscript, we became aware of a very recent paper by SilaevCO . In that paper, by separately using Eilenberger equation and diagram formalism, the author studied the Higgs mode excitation in the presence of the scattering. This is indeed the very first paper that rigorously calculates the scattering influence on optical properties within the Eilenberger equation in the literature, even though it is too complex to obtain final analytic solution. Nevertheless, based on the following reasons, the results in that paper are not correct. Firstly, in Ref. CO, by Silaev, the claimed conclusion that the Higgs-mode generation is zero without impurity is based on the incomplete electromagnetic effect in his approach. Specifically, both the Hamiltonian used in his diagram formalism and the Eilenberger equation are not gauge invariant with vector potential ${\bf A}$ alone GEG . It is well known that the gauge invariance is the basic character of the electromagnetic field. The absence of the gauge invariance indicates that the incomplete electromagnetic effect. Secondly, another conclusion in Ref. CO, that the Higgs mode is not sensitive to disorder, is also incorrect. This can be easily seen by the following simple analysis through the general physics. In the Nambu space, the BdG Hamiltonian in the presence of the Higgs mode excitation is written as $H_{\rm BdG}=\xi_{\hat{p}}\tau_{3}+\Delta_{0}({\bf r})\tau_{1}+\delta|\Delta({\bf r})|\tau_{1}$ in the real space, and the electron-impurity interaction is given by $V({\bf r})\tau_{3}$. Then, due to the non-commutation relation

[TABLE]

the Higgs mode must be sensitive to the disorder. In fact, the scattering influence on the Higgs mode in the present work exactly comes from this non-commutation relation. Specifically, our scattering term of the isotropic part [Eq. (41) with $w_{kk^{\prime}}|_{\xi_{k}=-\xi_{k^{\prime}}}=0$] is given by

[TABLE]

in which the projection operator $\Gamma^{\eta}_{k}$ picks up the energy-conserved scattering channel. Then, it is immediately observed that the Higgs-mode part ($\tau_{1}$ component of $\delta\rho_{k}$) $[\sum_{\eta}\Gamma^{\eta}\tau_{3}\Gamma^{\eta},\delta\rho_{k1}\tau_{1}]$ has the form of Eq. (52) limited by the energy conservation.

Acknowledgements.

This work was supported by the National Natural Science Foundation of China under Grants No. 11334014 and No. 61411136001.

Appendix A Derivation of Eq. (25)

In this part, we derive Eq. (25). From Eq. (18), one has

[TABLE]

in which $e^{itH_{\bf k}}=\sum_{\eta}\Gamma^{\eta}_{k}e^{itE^{\eta}_{\bf k}}$ is used.

The $n$-th order of above equation during the optical response is written as

[TABLE]

Similarly, one also finds

[TABLE]

Consequently, Eq. (25) is derived. For completeness, The explicit expressions of $Y^{i}_{\bf kk^{\prime}}(n\omega)$ [Eq. (26)] are given by

[TABLE]

One also has $N^{1}_{\bf kk^{\prime}}(n\omega)=2Y^{1}_{\bf kk^{\prime}}(n\omega)$, $N^{2}_{\bf kk^{\prime}}(n\omega)=2Y^{2}_{\bf kk^{\prime}}(n\omega)$ and $N^{0}_{\bf kk^{\prime}}(n\omega)=N^{3}_{\bf kk^{\prime}}(n\omega)=0$.

Appendix B Derivation of Eqs. (31) and (32)

We derive Eqs. (31) and (32) in this part. At the weak scattering, substituting the solved $\rho^{\omega}_{{\bf k}0}$ [Eq. (30)] into Eq. (8), one has

[TABLE]

with

[TABLE]

Then, with the explicit expression of $Y^{0}_{\bf kk^{\prime}}$ in Eq. (57), the above equation becomes

[TABLE]

in which we have taken care of the particle-hole symmetry to remove terms with the odd orders of $\xi_{k}$ and $\xi_{k^{\prime}}$ in the summation of ${\bf k}$ and ${\bf k^{\prime}}$. After the mathematical integral, above equation becomes

[TABLE]

Consequently, the optically excited current ${\bf j}$ in the linear regime and hence the optical conductivity $\sigma_{s}(\omega)=\sigma_{1s}(\omega)+i\sigma_{2s}(\omega)$ are derived.

Appendix C Optical conductivity at $T>T_{c}$

We give the optical conductivity at $T>T_{c}$. In the normal state at $T>T_{c}$, with $\Delta_{0}=0$, one finds that $l(E)=-\partial_{E}f(E)$ and $\frac{E(E+\omega)-\Delta^{2}_{0}}{\sqrt{(E+\omega)^{2}-\Delta^{2}_{0}}\sqrt{E^{2}-\Delta^{2}_{0}}}=E(E+\omega)/(|E||E+\omega|)$. Then, thanks to the constant density of states in normal states, Eqs. (31) and (32) become

[TABLE]

which are exactly the optical conductivity in normal metals as the Drude model or conventional Boltzmann equation revealed.

Appendix D Optical absorption in the dirty limit

By full numerical calculation of Eq. (29), the frequency dependence of the optical absorption towards the dirty limit are plotted in Fig. 4. As seen from Fig. 4(a), both the upturn of the finite $\sigma_{1s}(\omega)$ below $2\Delta_{0}$ and the crossover point at $\omega=2\Delta_{0}(T)$ in $\sigma_{1s}(\omega)$ appear in the dirty-limit regime ($\xi<l$), justifying our analysis in Sec. III.1.1 and in qualitative agreement with the experimental findings NSL1 ; NSL2 ; NSL10 ; NSL3 ; NSL4 ; NSL0 ; NSL11 ; NSL5 ; NSL6 ; NSL8 ; NSL9 ; NSL12 . Thus, the GIKE provides an efficient approach to capture the optical conductivity in the normal-skin-effect region. To quantitatively fit the experimental data in the dirty limit, the specific parameters of the density, effective mass and momentum-relaxation rate are necessary, and this goes beyond the scope of the present work.

Appendix E Derivation of Eq. (34)

We derive Eq. (34) in this part. Following the approach in our previous work in the clean limit GOBE5 , the solution of $\rho^{2\omega}_{\bf k}$ from Eq. (33) in the presence of the scattering is written as

[TABLE]

where $\partial_{t}\rho_{\bf k}|^{2\omega,\tau_{i}}_{\rm sc}$ denotes the $\tau_{i}$ component of the scattering term $\partial_{t}\rho_{\bf k}|^{2\omega}_{\rm sc}$; $A_{\bf k}$, $B_{\bf k}$ and $S^{c}_{\bf k}$ are given by

[TABLE]

The scattering term $\partial_{t}\rho_{\bf k}|^{2\omega}_{\rm sc}$ [Eq. (25)] reads:

[TABLE]

At the weak scattering, by substituting the clean-limit solution of $\rho_{\bf k}^{2\omega}$ into the scattering terms as the first-order iteration, Eq. (74) becomes

[TABLE]

Then, $\rho^{2\omega}_{\bf k}$ is solved.

Consequently, with $Y^{i}_{\bf kk^{\prime}}(2\omega)$ given by Eqs. (57)-(60), substituting $\rho^{2\omega}_{{\bf k}1}$ into Eq. (15), one has

[TABLE]

where $C_{k}$ is given by

[TABLE]

Here, we have taken care of the particle-hole symmetry to remove terms with the odd orders of $\xi_{k}$ and $\xi_{k^{\prime}}$ in the summation of ${\bf k}$ and ${\bf k^{\prime}}$; we also take $\eta_{\bf k}$ in $\rho^{\omega}_{\bf k}$ [Eq. (30)] as its average value ${\bar{\eta}}_{\bf k}$ in the momentum space. Then, after the mathematical integral, Eq. (34) is obtained.

Appendix F Solution of Eqs. (44)-(46)

In this part, we analytically solve Eqs. (44)-(46). Considering the weak scattering, we only keep zeroth and first orders of the scattering strength $\Gamma_{0}$ in the following derivation. Similar to the Elliot-Yafet relaxation mechanism in the semiconductor spintronics spintronic , following the Löwdin partition method dia , through a unitary transformation $(\delta\rho^{s}_{k+},\delta\rho^{s}_{k-},\delta\rho^{s}_{k3})^{T}=(1-S)(\delta\rho^{q}_{k+},\delta\rho^{q}_{k-},\delta\rho^{q}_{k3})^{T}$ with

[TABLE]

Eqs. (44)-(46) become

[TABLE]

from which $\delta\rho_{k}^{s}$ can be directly solved:

[TABLE]

Through the inverse transformations $\delta\rho_{k}=U^{\dagger}_{k}\delta\rho^{q}_{k}U_{k}$ and $(\delta\rho^{q}_{k+},\delta\rho^{q}_{k-},\delta\rho^{q}_{k3})^{T}=(1+S)(\delta\rho^{s}_{k+},\delta\rho^{s}_{k-},\delta\rho^{s}_{k3})^{T}$, one has

[TABLE]

Finally, substituting Eq. (84) into Eq. (15), by taking care of the particle-hole symmetry to remove terms with the odd order of $\xi_{k}$ in the summation of ${\bf k}$, one obtains

[TABLE]

where the phase shift $\phi_{k}=\arctan[{\rm sgn}(\xi_{k})/\xi_{k}\Delta^{2}_{0}\Gamma_{0}/E_{k}^{2}]$ can be neglected at the weak scattering. Then, Eq. (47) is derived.

Appendix G Derivation of Eq. (51)

In this part, we derive Eq. (51). To consider the long-time dynamic behavior of the Higgs mode, by approximately taking the starting point of time as $-\infty$ in Eq. (47), one has

[TABLE]

with $A_{i}=\sum_{\bf k}\frac{\Delta_{0}}{E_{k}}c_{k3}$ and \Gamma_{H}=\sum_{\bf k}\big{[}\gamma_{k}g(E_{k})\frac{\Delta_{0}^{2}}{E_{k}^{2}}\big{]}.

In the frequency space $\delta|\Delta|(t)=\int\frac{d\Omega}{2\pi}\delta|\Delta|_{\Omega}e^{-i\Omega{t}+0^{+}t}$, the above equation becomes

[TABLE]

By using Eq. (24) to replace $g$, one has

[TABLE]

Consequently, the temporal evolution of the Higgs mode is given by

[TABLE]

Considering the weak scattering, the second term on the right-hand side of the above equation can be neglected. By keeping the zeroth and first orders of the scattering, one obtains

[TABLE]

It is noted that for the integrand in Eq. (92), in the complex plane of $\Omega$, there exist two branching points at $\Omega=\pm 2\Delta_{0}$. Then, similar to the previous workOD3 , after the standard construction of the closed contour, one obtains

[TABLE]

Appendix H Response of NG mode

As mentioned in the introduction, in our latest work for the clean limit GOBE5 , a finite second-order response of the NG mode, free from the influence of the Anderson-Higgs mechanism, is predicted as a consequence of charge conservation. An experimental scheme for this response is further proposed based on Josephson junction. In this part, for completeness, we study the influence of the scattering on this response during and after the THz pulse.

H.0.1 Excitation of NG mode in second-order response

During the pulse, substituting $\rho^{2\omega}_{{\bf k}2}$ into Eq. (16), the NG mode can be self-consistently derived:

[TABLE]

Here, $g_{\omega}={{\int}d\xi_{k}l(E_{k})z(E_{k})}/[{\int{d}\xi_{k}g(E_{k})z(E_{k})}]$ and $l_{\omega}={{\int}d\xi_{k}m(E_{k})z(E_{k})}/[{\int{d}\xi_{k}g(E_{k})z(E_{k})}]$ with $z(E_{k})=(E^{2}_{k}-\omega^{2})^{-1}$ and $m(E_{k})=(2\xi_{k}\partial_{\xi_{k}}+1)l(E_{k})$. The scattering contributions are given by

[TABLE]

and $\gamma_{1}$ and $\gamma_{2}$ are determined via replacing function $g(x)$ on the right-hand side of Eq. (96) by $l(x)$ and $m(x)$, respectively. It is noted that in the absence of the scattering, Eq. (94) exactly reduces to the clean-limit one revealed in our previous work GOBE5 .

Consequently, similar to the investigation of the Higgs mode in Sec. III.1.2, the scattering also causes a phase-shift in the second-order response of the NG mode. Nevertheless, this phase shift is hard to detect, differing from the measurable optical response of the Higgs mode in Sec. III.1.2.

H.0.2 Damping of NG mode

After the pulse, by numerically solving our simplified model in Sec. III.2.1 [Eq. (42) with Eqs. (15) and (16)], the temporal evolution of the optically excited NG mode is plotted in Fig. 5 at different scattering rates. As seen from the figure, the NG mode $\mu_{\rm eff}(t)=\partial_{t}\theta(t)$, i.e., the phase fluctuation, after the optical excitation exhibits an oscillatory decay behavior. The oscillating frequency is around $2\Delta_{0}$, and the damping shows a monotonic enhancement with the increase of the scattering rate. Particularly, by further comparing Figs. 3 and 5, it is interesting to find that the damping of the phase fluctuation (NG mode) is much faster than that of the amplitude fluctuation (Higgs mode) of the order parameter.

Substituting the analytic solution of $\delta\rho_{k2}$ [Eq. (85)] into Eq. (16), by taking care of the particle-hole symmetry to remove terms with the odd order of $\xi_{k}$ in the summation of ${\bf k}$, the analytic solution of $\mu_{\rm eff}(t)$ is derived:

[TABLE]

As seen from Eq. (97), terms on both the left- and right-hand sides show the oscillatory decay with the time evolution, and hence, directly lead to the oscillating damping of $\mu_{\rm eff}(t)$ with the relaxation rate proportional to $\Gamma_{0}$. Comparisons between the solution from Eq. (97) and the full numerical results from Eq. (42) are plotted in the insets of Fig. 5, and the results from the two sets of calculations agree with each other again.