Nonlinear electromagnetic response and Higgs mode excitation in BCS superconductors with impurities

TL;DR

This paper demonstrates that impurities in BCS superconductors enable effective excitation of the Higgs mode by electromagnetic radiation, explaining recent experimental observations of third-harmonic signals.

Contribution

It introduces a quasiclassical formalism to describe nonlinear electromagnetic response and Higgs mode excitation in disordered BCS superconductors, including analytical solutions matching impurity diagram summations.

Findings

Impurities facilitate Higgs mode excitation by electromagnetic radiation.

Analytical solutions match impurity diagram summations.

Resonant third-harmonic signals are explained by Higgs mode excitation.

Abstract

We reveal that due to the presence of disorder oscillations of the order parameter amplitude called the Higgs mode can be effectively excited by the external electromagnetic radiation in usual BCS superconductors. This mechanism works for superconductors with both isotropic s-wave and anisotropic, such as d-wave, pairings. The non-linear response in the presence of impurities is captured by the quasiclassical formalism. We demonstrate that analytical solutions of the Eilenberger equation with impurity collision integral and external field drive coincide with the exact summation of ladder impurity diagrams. Using the developed formalism we show that resonant third-harmonic signal observed in recent experiments is naturally explained by the excitation of Higgs mode mediated by impurity scattering.

Figures

b. As before we are interested in the momentum-integrated function $\hat{g}_{A\Delta}=i\pi^{-1}\int d\xi_{p}\hat{\tau}_{3}\hat{G}_{A\Delta}$ because it determines the correction to the current. Treating the diagram Fig.[6](#S4.F6)d using exactly the same procedure as above we arrive to the equation for $\hat{g}_{A\Delta}$ obtained directly from the Eilenbeger formalism. Its solution is given by the Eq.([38](#S3.E38)).

Equations

$${\hat{H}}_p={\hat{V}}_1+{\hat{V}}_2$$

$${\hat{V}}_1=-\frac{e}{c}(\mathbf{v}\mathbf{A})$$

$${\hat{V}}_2={\hat{\tau }}_3\frac{e^2}{2m{c}^2}{A}^2,$$

$$\mathbf{j}(\tau )=e\int \frac{d^3\mathbf{p}}{(2\pi {)}^3}\left(\mathbf{v}-{\tau }_3\frac{e}{c}{\hat{m}}^{-1}\mathbf{A}(\tau )\right)\hat{G}(\mathbf{p},{\tau }_{1,2}=\tau )$$

$$\hat{\Delta }(\tau )=\tilde{\lambda }\int \frac{d^3\mathbf{p}}{(2\pi {)}^3}\hat{P}[\hat{G}](\mathbf{p},{\tau }_{1,2}=\tau )$$

$${\hat{\mathcal{D}}}_0[\hat{G}]=\tilde{\lambda }T\sum _{\omega }\int \frac{d^3\mathbf{p}}{(2\pi {)}^3}\hat{P}[\hat{G}](\omega ,\Omega ,\mathbf{p}).$$

$${\hat{G}}_0(\omega ,\mathbf{p})=\frac{\tilde{\Delta }{\hat{\tau }}_1-i\tilde{\omega }{\hat{\tau }}_0-{\xi }_p{\hat{\tau }}_3}{{\tilde{\Delta }}^2+{\tilde{\omega }}^2+{\xi }_p^2}$$

$$\tilde{\omega }=\omega \frac{\tilde{s}(\omega )}{s(\omega )};\tilde{\Delta }=\Delta \frac{\tilde{s}(\omega )}{s(\omega )},$$

$$\hat{g}=\frac{i}{\pi }\int d{\xi }_p{\hat{\tau }}_3\hat{G}$$

$$\frac{ie}{c}{\mathbf{v}}_F[{\hat{\tau }}_3\mathbf{A},\hat{g}{]}_{\tau }=$$

$$-i\{{\hat{\tau }}_3{\partial }_{\tau },\hat{g}{\}}_{\tau }+i[{\hat{\tau }}_3\hat{\Delta },\hat{g}{]}_{\tau }+\frac{1}{2{\tau }_{imp}}[⟨\hat{g}⟩\circ ,\hat{g}].$$

$$\mathbf{j}(\tau )=-i\pi e\nu \mathrm{Tr}[{\hat{\tau }}_3⟨{\mathbf{v}}_F\hat{g}(\tau ,\tau )⟩]$$

$$\hat{\Delta }(\tau )=-i\lambda \hat{P}[{\hat{\tau }}_3⟨\hat{g}(\tau ,\tau )⟩],$$

$$-i\{{\hat{\tau }}_3{\partial }_{\tau },\hat{g}{\}}_{\tau }+i[{\hat{\tau }}_3\hat{\Delta },\hat{g}{]}_{\tau }=\hat{\Sigma }\circ \hat{g}-\hat{g}\circ \hat{\Sigma }$$

$${\hat{\Sigma }}_{em}({\tau }_1,{\tau }_2)=\frac{D{e}^2}{c^2}\mathbf{A}({\tau }_1){\hat{\tau }}_3\hat{g}({\tau }_1,{\tau }_2){\hat{\tau }}_3\mathbf{A}({\tau }_2)$$

$$\mathbf{j}(\tau )=\pi \frac{\sigma }{c}\times$$

$$\mathrm{Tr}[{\hat{\tau }}_3\hat{g}(\tau ,{\tau }_1)\circ \left(\hat{g}({\tau }_1,\tau )\mathbf{A}(\tau ){\hat{\tau }}_3-{\hat{\tau }}_3\mathbf{A}({\tau }_1)\hat{g}({\tau }_1,\tau )\right)]$$

$$\hat{g}({\tau }_1,{\tau }_2)=$$

$$T\sum _{\omega }{e}^{i{\omega }_1{\tau }_1}[{\hat{g}}_0(1)e^{-i{\omega }_1{\tau }_2}+{\hat{g}}_A(12)e^{-i{\omega }_2{\tau }_2}+$$

$${\hat{g}}_{AA}(13)e^{-i{\omega }_3{\tau }_2}+{\hat{g}}_{AAA}(14)e^{-i{\omega }_4{\tau }_2}]$$

$${\hat{g}}_0(\omega )=\frac{{\hat{\tau }}_3\omega -{\hat{\tau }}_2\Delta }{s(\omega )}$$

$${\hat{g}}_A=\alpha \cos \chi {\hat{g}}_{1a}(\omega )$$

$${\hat{g}}_{AA}={\alpha }^2({\cos }^2\chi -1/3){\hat{g}}_{2a}+{\alpha }^2{\hat{g}}_{2s}/3$$

$${\hat{g}}_{AAA}={\alpha }^3\cos \chi \left[({\cos }^2\chi -1/3){\hat{g}}_{3a}+{\hat{g}}_{3s}/3\right]$$

$$i\left[{\hat{\tau }}_3{\hat{g}}_0(3)-{\hat{g}}_0(2){\hat{\tau }}_3\right]=$$

$${\tilde{s}}_2{\hat{g}}_0(2){\hat{g}}_{1a}(23)-{\tilde{s}}_3{\hat{g}}_{1a}(23){\hat{g}}_0(3),$$

$$i\left[{\hat{\tau }}_3{\hat{g}}_{1a}(34)-{\hat{g}}_{1a}(23){\hat{\tau }}_3\right]=$$

$$s_2{\hat{g}}_0(2){\hat{g}}_{2s}(234)-s_4{\hat{g}}_{2s}(234){\hat{g}}_0(4).$$

$$\left[i{\hat{\tau }}_3+\frac{{\hat{g}}_{1a}(12)}{2{\tau }_{imp}}\right]{\hat{g}}_{2s}(234)-{\hat{g}}_{2s}(123)\left[i{\hat{\tau }}_3+\frac{{\hat{g}}_{1a}(34)}{2{\tau }_{imp}}\right]=$$

$${\tilde{s}}_1{\hat{g}}_0(1){\hat{g}}_{3s}-{\tilde{s}}_4{\hat{g}}_{3s}{\hat{g}}_0(4),$$

$$i\left[{\hat{\tau }}_3{\hat{g}}_{2a}(234)-{\hat{g}}_{2a}(123){\hat{\tau }}_3\right]={\tilde{s}}_1{\hat{g}}_0(1){\hat{g}}_{3a}-{\tilde{s}}_4{\hat{g}}_{3a}{\hat{g}}_0(4).$$

$${\hat{g}}_{1a}(12)=i\frac{{\hat{g}}_0(1){\hat{\tau }}_3{\hat{g}}_0(2)-{\hat{\tau }}_3}{s_1+s_2+{\tau }_{imp}^{-1}}.$$

$${\hat{g}}_{2s}(234)=i\frac{s_2{\hat{g}}_0(2)\hat{X}(234)+s_4\hat{X}(234){\hat{g}}_0(4)}{s_2^2-s_4^2},$$

$$\hat{X}(234)={\hat{\tau }}_3{\hat{g}}_{1a}(34)-{\hat{g}}_{1a}(23){\tau }_3$$

$${\hat{g}}_{3s,3a}=i\frac{{\tilde{s}}_1{\hat{g}}_0(1){\hat{Y}}_{s,a}+{\tilde{s}}_4{\hat{Y}}_{s,a}{\hat{g}}_0(4)}{(s_1-s_4)(s_1+s_4+{\tau }_{imp}^{-1})},$$

$${\hat{Y}}_s=$$

$$\left[{\hat{\tau }}_3-\frac{i{\hat{g}}_{1a}(12)}{2{\tau }_{imp}}\right]{\hat{g}}_{2s}(234)-{\hat{g}}_{2s}(123)\left[{\hat{\tau }}_3-\frac{i{\hat{g}}_{1a}(34)}{2{\tau }_{imp}}\right],$$

$${\hat{Y}}_a={\hat{\tau }}_3{\hat{g}}_{2a}(234)-{\hat{g}}_{2a}(123){\hat{\tau }}_3.$$

$${\hat{g}}_0(1){\hat{g}}_{1a}(13)+{\hat{g}}_{1a}(13){\hat{g}}_0(3)=0$$

$${\hat{g}}_0(1){\hat{g}}_{2s}(123)+{\hat{g}}_{2s}(123){\hat{g}}_0(3)+{\hat{g}}_{1a}(12){\hat{g}}_{1a}(23)=0$$

$${\hat{g}}_0(1){\hat{g}}_{3s}+{\hat{g}}_{3s}{\hat{g}}_0(4)+$$

$${\hat{g}}_{1a}(12){\hat{g}}_{2s}(234)+{\hat{g}}_{2s}(123){\hat{g}}_{1a}(34)=0$$

Full text

I Introduction

Nonlinear electromagnetic responses are ubiquitous in superconducting systems and have attracted interest for many years Gor’kov and Eliashberg (1968); Gorkov and Eliashberg (1968); Amato and McLean (1976); Gorkov and Eliashberg (1969). External field with the frequency $\Omega$ produces several important second-order corrections to the superconducting order parameter $\Delta$. The zero-frequency change of $\Delta$ leads to the critical temperature and critical current enhancements Ivlev and Eliashberg (1971); Eliashberg (1970); Tikhonov et al. (2018) known as the microwave stimulation of the superconducting state. The time-dependent corrections to $\Delta$ at the frequency $2\Omega$ produce the current oscillating with frequency Gor’kov and Eliashberg (1968) $3\Omega$. This effect called the third harmonic generation (THG) has been observed experimentally in microwaves and explained with the help of the time-dependent Ginzburg-Landau theory Amato and McLean (1976).

Recently the terahertz (THz) spectroscopy of superconducting state has become experimentally available Beck et al. (2011, 2013); Matsunaga and Shimano (2012); Matsunaga et al. (2013, 2014); Giorgianni et al. (2019). This range of frequencies is especially interesting since it overlaps with the typical gap sizes in low-temperature superconductors like NbN. Thus measuring nonlinear responses in THz domain allows for probing dynamics of the order parameter amplitude Matsunaga et al. (2013, 2014) predicted to feature oscillations with an eigen frequency $2\Delta(T)$ where $\Delta(T)$ is the the gap at a given temperature Volkov and Kogan (1974); Kulik et al. (1981); Barankov et al. (2004); Barankov and Levitov (2006). By analogy with the Higgs bosonHiggs (1964) in particle physics this type of collective excitation in condensed matter systems is called the Higgs modeVarma (2002); Podolsky et al. (2011); Pashkin and Leitenstorfer (2014); Volovik and Zubkov (2014); Pekker and Varma (2015). The order parameter amplitude oscillations excited by the short optical pulses has been observed in several pump-probe experiments Matsunaga et al. (2013, 2014). Recent measurements report the evidence of resonant Higgs mode excitation in the THG component of the THz signal transmitted through the superconducting plateMatsunaga et al. (2014, 2017). Earlier the amplitude modes have been observed by Raman scattering in superconductors with charge density wave order Sooryakumar and Klein (1980); Littlewood and Varma (1982); Grasset et al. (2018a, b) and by the nuclear magnetic resonance in superfluid 3He Paulson et al. (1973); Lawson et al. (1973); Zavjalov et al. (2016).

Despite the significant experimental advances, theoretical understanding of high-frequency non-linear properties in superconductors is still lacking. Numerical simulations presented in several works consider strongly non-equilibrium regimes without any disorder Krull et al. (2014); Papenkort et al. (2008); Fauseweh et al. (2017) and have to attribute rather large wavevector to the radiation field in order to obtain the sizable coupling with the order parameter. At the same time perturbative calculations of nonlinear responses reveal several important limitations imposed by the absence of impurity scattering. The order parameter modulation by the external radiation was studied in the pioneering work of Gor’kov and Eliashberg Gorkov and Eliashberg (1968) who considered the electron-photon coupling linear by the vector potential $\bm{A}_{\Omega}$. The contribution of such terms to the order parameter modulation $\Delta_{2\Omega}$ at the frequency $2\Omega$ is described by the second-order diagram in Fig.1a. Here current vertices $\bullet$ describe the linear coupling to the external field. This contribution disappears the absence of impurity scatteringGorkov and Eliashberg (1968). One can think of this as a consequence of the Galilean invariance featured by the superconducting condensate. Indeed switching to the moving frame can eliminate the condensate velocity induced by external field. In the moving frame the order parameter amplitude remains unaffected thus coinciding with that of the stationary condensate. Therefore, in order to perturb the amplitude through the linear electron-photon coupling terms the Galilean invariance should be broken either by the spatially-inhomogeneous field or the inhomogeneous potential which can naturally related to the presence of impurities.

Based on these arguments at zero or negligibly small wave vectors the Higgs mode excitation Tsuji and Aoki (2015) and non-linear responseCea et al. (2016) in the absence of impurities are possible only through the electron density modulation generated by the term $\propto A_{\Omega}^{2}$ in the Hamiltonian. The corresponding density vertices $\bigcirc$ are similar to those which determine Raman scattering in superconductors Abrikosov and Falkovsky (1961); Abrikosov and Genkin (1974); Sooryakumar and Klein (1980); Abrikosov and Falkovsky (1988); Falkovsky (1993). The perturbation of the order parameter due to such coupling is shown by the diagram in Fig.1b. However, this contribution to $\Delta_{2\Omega}$ vanishes within the Bardeen-Cooper-Schrieffer (BCS) model of superconductivity and becomes non-zero only due to the various extensions Tsuji and Aoki (2015); Cea et al. (2016); Matsunaga et al. (2017).

The prediction of negligible Higgs mode generation in BCS superconductors Cea et al. (2016) has been in contrast both to the recent THz probesBeck et al. (2013, 2011); Matsunaga et al. (2014, 2013); Matsunaga and Shimano (2012); Katsumi et al. (2018); Matsunaga et al. (2017) and earlier microwave experimentsAmato and McLean (1976) where significant radiation-induced oscillation $\Delta_{2\Omega}$ have been observed. Here we resolve this controversy and show that previous theoretical conclusions about the non-linear response of superconductors are drastically altered in the presence of disorder. We consider the arbitrary amount of impurity scattering treated within the self-consistent Born approximation. The calculations are implemented using both the diagrammatic technique and quasiclassical Eilenberger theory formalism Eilenberger (1968) and the agreement between these two approaches is demonstrated. We show that in the presence of disorder the non-linear process shown in Fig. 1a provides quite an effective generation of the Higgs mode which shows up through the resonant THG in agreement with several recent experiments.

In the diffusive limit typical for the superconducting thin films our findings confirm that the resonant third-harmonic generation at the frequency $\Omega=\Delta(T)$ observed in Matsunaga et al. (2014) is determined solely by the Higgs mode generation. This result is in sharp contrast with the system without impurities where the Higgs mode generation is negligible and the resonant contribution comes from the other source Cea et al. (2016, 2018). It agrees qualitatively with the studies of linear Higgs mode generation in diffusive current-carrying superconductorMoor et al. (2017); Nakamura et al. (2018).

II Formalism

II.1 General approach

We describe the interaction of electrons with electromagnetic field using the following Hamiltonian which contains two qualitatively different terms Abrikosov and Genkin (1974); Abrikosov and Falkovsky (1988); Falkovsky (1993)

[TABLE]

where $\bm{v}=\partial E_{p}/\partial\bm{p}$ is the band velocity, $m$ and $e$ are the electron mass and charge.

Here we introduce the notation $\hat{\tau}_{0,1,2,3}$ for the Pauli matrices in Nambu particle-hole space. In the diagrammatic representation the perturbation term $\hat{V}_{1}$ linear in the external field amplitude is described by the current vertex with attached single external field line. Such current vertices determine diagrams of the type shown in Fig.1a. The term $\hat{V}_{2}$ quadratic by the external field produces radiation-induced electronic density modulation. Thus light-matter coupling is described by diagrams with density vertices $\propto\hat{\tau}_{3}$ such as shown by the open circles in Fig.1b.

The charge current and order parameter as functions of the at imaginary time $\tau[0,\beta]$ are given by

[TABLE]

where $\hat{G}(\bm{p},\tau_{1,2})$ is the imaginary time Green’s function (GF), $\tilde{\lambda}$ is the pairing constant, $\hat{P}[\hat{G}]=\hat{\tau}_{1}{\rm Tr}(\hat{\tau}_{1}\hat{G})+\hat{\tau}_{2}{\rm Tr}(\hat{\tau}_{2}\hat{G})$ is the projection operator. It will be convenient to use also the frequency representation which can be defined as $\hat{G}(\omega,\Omega)=\int d\tau_{1}d\tau_{2}\hat{G}(\tau_{1},\tau_{2})e^{i\omega(\tau_{1}-\tau_{2})+i\Omega(\tau_{1}+\tau_{2})}$. In this representation the self-consistency relation can be written as $\Delta_{2\Omega}=\hat{\cal D}_{0}[\hat{G}]$, where $\hat{\cal D}_{0}$ is the order parameter vertex which appears first in the diagrams in Fig.(1). The algebraical expression for $\hat{\cal D}_{0}$ reads

[TABLE]

The stationary propagators depend only on the frequency corresponding to the relative time $\hat{G}=\hat{G}_{0}(\omega,\bm{p})$. In disordered superconductors with arbitrary amount of point-like impurity scatterers treated within the self-consistent Born approximation propagators are given by (Abrikosov et al., 1963)

[TABLE]

where $\xi_{p}=p^{2}/2m-\mu$ is the deviation of the kinetic energy from the chemical potential $\mu$ and $\hat{\tau}_{1,2,3}$ are the Pauli matrices in Nambu space. We denote $s=\sqrt{\omega^{2}+\Delta^{2}}$ and $\tilde{s}=s+1/2\tau_{imp}$.

The propagator (7) is averaged over the randomly disordered point scatterers configurations. It is parametrised by the scattering time $\tau_{imp}^{-1}=2\pi\nu nu^{2}$, where $\nu$ is the normal metal density of states at the Fermi level, $u$ is the impurity potential strength and $n$ is the density of impurities.

We are interested in the non-linear contribution to the current determined by the diagrams in Fig.2 where the coupling to the electromagnetic field is determined by the vector current-type vertices. Such diagrams are generated by the perturbation potential $\hat{V}_{1}$ given by Eq.(2). As we see below the contribution of such diagrams to the measurable quantities is captured by the quasiclassical Eilenberger equations Eilenberger (1968). We assume here that all external field lines have the same time dependence $\bm{A}_{\Omega}e^{i\Omega t}$, so that diagrams in Fig.2 yield the third harmonic response of the current $(j_{AAA}+j_{H})e^{3i\Omega t}$.

The diagrams shown in Fig.(2a) determine the current $j_{AAA}$ generated by the direct coupling to the external field through the current vertices. Dashed lines correspond to impurity scattering averaged over the random disorder configuration. Analytically one should integrate over the input/output momentum the content between the dashed line start and end points as well as multiply the result by the factor $1/(2\pi\tau_{imp})$. The shaded regions show impurity ladders discussed in detail in Sec.IV.2.

Besides the direct coupling to external field equally important is the third-harmonic generation by the current $j_{H}$ determined by the order parameter modulation $\Delta_{2\Omega}e^{i2\Omega\tau}$. Taking into account that $\Delta_{2\Omega}$ is generated by external source according to the diagram in Fig.1a the current $j_{H}$ can be represented by the third-order response diagram in Fig.2b. This contribution is of the special interest since as we will demonstrate its frequency dependence contains the information about the Higgs mode, that is the resonant enhanceable of the order parameter oscillations amplitude at $\Omega=\Delta$. Technically the resonant Higgs mode contribution is determined by the polarization operator which modifies the gap function vertex shown by the filled triangle in Fig.(2b). This vertex yields the order parameter coupling to the external source such as the second-order correction of the GF by the electromagnetic field discussed above. The diagrammatic representation of the order parameter vertex with polarization bubble insertions is shown in Fig.(3). Note that it contains impurity scattering corrections, both as the self-energies modifying individual propagator lines and the ladder dressing the $\rhd=\hat{\tau}_{1}$ vertex.

Previously, it has been shown that in the the absence of impurities $\tau_{imp}=\infty$ the contribution of diagrams with current vertices ${\Large{\bullet}}$ to the order parameter modulation $\Delta_{2\Omega}$ disappears Gorkov and Eliashberg (1968); Tsuji and Aoki (2015); Matsunaga et al. (2014); Cea et al. (2016). In this case the only non-zero contribution is given by the diagrams with density vertices $\bigcirc$ shown in Fig.4. Below we compare the contributions of these diagrams to obtain the threshold impurity concentration when the paramagnetic diagrams start to dominate.

II.2 Quasiclassical theory

In general it is believed that the quasiclassical approximation introduced by Eilenberger Eilenberger (1968) takes into account all diagrams with current vertices shown in Fig. (2) but neglects the ones with density vertices. Technically this happens because the $\hat{\tau}_{3}A^{2}$ terms in the Hamiltonian generically drop out from the quasiclassical equations. However, previously there has been done no direct comparison of the results given by quasiclassical theory with those obtained from diagrams describing coupling to the external field in the presence of impurity scattering. We will implement this check and demonstrate the summation of diagrams with impurity ladders give the same results as the quasiclassical calculation implemented according to the formalism described below.

The quasiclassical propagator is defined as

[TABLE]

In the imaginary time domain $\hat{g}=\hat{g}(\tau_{1},\tau_{2},\bm{r},\bm{v}_{F})$ is determined by the Eilenberger equation Eilenberger (1968)

[TABLE]

Here we denote the commutators $[X,g]_{\tau}=X(\tau_{1})g(\tau_{1},\tau_{2})-g(\tau_{1},\tau_{2})X(\tau_{2})$ and the convolution is given by $\langle\hat{g}\rangle\circ\hat{g}=\int_{0}^{\beta}d\tau\langle\hat{g}\rangle(\tau_{1},\tau)\hat{g}(\tau,\tau_{2})$. The angle-averaging over the Fermi surface is given by $\langle g\rangle$. The current and order parameter are given by

[TABLE]

where $\lambda=\pi\nu\tilde{\lambda}$ and $\bm{v}_{F}=\bm{v}(\bm{p}=\bm{p}_{F})$ is the Fermi velocity which is the band velocity at the momentum equal to the Fermi momentum $\bm{p}_{F}$. The quasiclassical equations are supplemented by the normalization condition $\hat{g}\circ\hat{g}=1$.

Dirty limit: Usadel theory. In the dirty limit $\tau_{imp}T_{c}\ll 1$ the calculations can be significantly simplified using the Usadel equation formalismUsadel (1970). The key idea of this approximation is to represent the GF as the superposition of isotropic $\langle\hat{g}\rangle$ and anisotropic $\hat{g}_{an}\propto(\bm{v}_{F}\bm{A})$ parts. For the validity of Usadel theory it is required that the anisotropic part is small $\hat{g}_{an}\ll\langle\hat{g}\rangle$ which is satisfied in the diffusive limit. Then using the normalization condition $\langle\hat{g}\rangle\circ\hat{g}_{an}+\hat{g}_{an}\circ\langle\hat{g}\rangle=0$ one obtain get the Usadel equation for the isotropic component $\langle\hat{g}\rangle$ which reads (we omit angular brackets)

[TABLE]

where the self-energy $\hat{\Sigma}_{em}$ describes coupling to the electromagnetic field, the the diffusion coefficient is $D=\tau_{imp}v_{F}^{2}/3$ and conductivity is $\sigma=e^{2}\nu D$.

III Quasiclassical calculations

First, we obtain time-dependent perturbations of the order parameter and the third-order current response using the Eilenberger theory for quasiclassical propagators. In the Sec.IV these results will be confirmed by the direct summation of impurity ladder diagrams with current-type vertices describing the light-matter interaction.

III.1 Nonlinear response: direct light-matter interaction

We start with calculating corrections generated directly by the time-dependent vector potential. Corrections due to the order parameter oscillations which contain the Higgs mode contribution are considered below in Sec.III.2. In the presence of the oscillating external field $\bm{A}_{\Omega}e^{i\Omega\tau}$ we can find the solution of Eilenberger equation (10) in the form of expansion by the orders of $\bm{A}_{\Omega}$ as follows

[TABLE]

where we introduce the shortened notation to define the frequency dependence such as $\hat{g}_{A}(ij)=\hat{g}_{A}(\omega_{i},\omega_{j})$ and the shifted Matsubara frequencies $\omega_{1}=\omega+2\Omega$, $\omega_{2}=\omega+\Omega$, $\omega=\omega$, $\omega_{4}=\omega-\Omega$. The zeroth-order solution given by the unperturbed propagator given from Eqs.(7,9)

[TABLE]

where $s(\omega)=\sqrt{\omega^{2}+\Delta^{2}}$. The expansion terms in Eq.(16) can be found in the form where momentum direction dependence is explicitly defined

[TABLE]

where we denote $\cos\chi=\bm{v}_{F}\bm{A}_{\Omega}/v_{F}A_{\Omega}$ and $\alpha=ev_{F}A/c$.

The impurity collision integral $[\langle\hat{g}\rangle,\hat{g}]$ disappears for the isotropic parts of the propagator, so that we have in the first-order $[\langle\hat{g}\rangle,\hat{g}]^{(1)}=[\hat{g}_{0},\hat{g}_{A}]$, in the second-order $[\langle\hat{g}\rangle,\hat{g}]^{(2)}=\alpha^{2}(\cos^{2}\chi-1/3)[\hat{g}_{0},\hat{g}_{2a}]$ and for the third-order correction $[\langle\hat{g}\rangle,\hat{g}]^{(3)}=[\hat{g}_{0},\hat{g}_{3}]+[\hat{g}_{2s},\hat{g}_{A}]/3$. Then we obtain a chain of equations to determine corrections driven by the direct coupling to the vector potential. The equation for the first-order correction reads

[TABLE]

where we introduce the notation $s_{i}=s(\omega_{i})$ and $\tilde{s}_{i}=s_{i}+1/2\tau_{imp}$. The second order correction $\hat{g}_{2s}$ is

[TABLE]

The equation for $\hat{g}_{2a}$ has the similar form with the replacements $s_{2,4}\to\tilde{s}_{2,4}$ in the right hand side. And finally the equation for the third-order correction is

[TABLE]

The solution of above equations reads as follows. The first-order correction is given by

[TABLE]

The isotropic second-order correction $\hat{g}_{2s}$ is given by

[TABLE]

and the anisotropic part $\hat{g}_{2a}$ is obtained by the replacements $s_{2,4}\to\tilde{s}_{2,4}$ in the right hand side of Eq.(26). Finally, the third-order correction is given by

[TABLE]

As discussed below in Sec.III.4 in order to use solutions (28) for the numerical calculation it is necessary to convert them into the form which does not have the factor $s_{1}-s_{4}$ in the denominators. This can be done with the help of normalization condition $\hat{g}\circ\hat{g}=1$ which yields the following relations for the corrections

[TABLE]

The commutation relations for $\hat{g}_{2a}$ and $\hat{g}_{3a}$ are obtained by substituting these functions instead of $\hat{g}_{2s}$ and $\hat{g}_{3s}$ to the relation (32,33), respectively. As the consistency check, one can prove by the direct calculation that the solutions (26,28) satisfy Eqs.(32,33). With the help of these relations one can rewrite Eq.(28) in the form suitable for numerics as described in the Appendix A. Using this form it is also possible to derive the diffusive limit consistent with the results given directly by the Usadel equations as discussed in Sec.V.4 and in the Appendix B.

III.2 Higgs mode contribution to the nonlinear response

Besides corrections to the GF induced directly by the coupling to vector potential we need to take into account the nonlinear current induced through order parameter oscillations. This process is depicted by the diagrams shown in Fig. 2a. First, let us calculate the current generated by the vector potential $\bm{A}_{\Omega}e^{i\Omega t}$ combined with the order parameter oscillating with the frequency $2\Omega$ which we denote as $\hat{\tau}_{1}\Delta_{2\Omega}e^{i\Omega t}$. We can find the corresponding solutions of Eilenberger equation in the form

[TABLE]

We use a similar approach to solve the chain of equations for corrections as in the previous section. Then the obtained first- and second-order solutions read

[TABLE]

where the first-order correction due to the vector potential $\hat{g}_{1a}$ is given by the Eq.(25). Note that the correction $\hat{g}_{\Delta}$ induced by the order parameter oscillations is completely isotropic and therefore is not affected by the impurity scattering collision integral. From this one can immediately conclude that the polarization operator which determines the Higgs mode is not sensitive to the disorder.

As discussed below in Sec.III.4 in order to use the solution (38) for the numerical calculation of the analytical continuation it is necessary to convert it into the form which does not have the factor $s_{1}-s_{4}$ in the denominator. This can be done with the help of normalization condition of the quasiclassical GF. For the corrections it yields

[TABLE]

One can check by the direct calculation with the certain analytical effort that the solutions (37,38) satisfy the conditions (41,42). Then, using these relations it is possible to convert (42) to the required form as described in the Appendix A.

The correction $\hat{g}_{A\Delta}$ is the basic building block to calculate the Higgs mode-related current $j_{H}$ according to the diagram in Fig.2b. In order to obtain $j_{H}$ as the third-order response to the external field we need to calculate the order parameter amplitude $\Delta_{2\Omega}$ excited in the second order by $\bm{A}_{\Omega}$ as shown by the diagram in Fig.1a. This can be done in two steps described below.

External perturbation of the order parameter. First, from the diagram in Fig.1a we find the source, that is the time-dependent order parameter induced directly by the external field. The amplitude $F_{\Delta}$ can be found using the isotropic second-order correction for the GF $\hat{g}_{2s}$ calculated above (26)

[TABLE]

Polarization operator Second, to find the order parameter amplitude $\Delta_{2\Omega}$ driven by the external source $F_{\Delta}$ we need to take into account the polarization corrections given by the diagrammatic series shown in Fig. (3). Thus the renormalized order parameter vertex $\hat{{\cal D}}$ denoted in Fig.(3) by the filled triangle is related to the bare one $\hat{S}_{0}$ as $\hat{{\cal D}}=\hat{{\cal D}}_{0}/[1-\Pi(2\Omega)]$. Here $\Pi$ is the polarization operator given by the single bubble in Fig. 3. Thus the total order parameter perturbation is related to the external source $F_{\Delta}$ as

[TABLE]

Previously, the polarization operator has been calculated in the clean limit Tsuji and Aoki (2015); Cea et al. (2016). In the presence of impurities the diagram summation becomes more complicated because besides the modification of propagator lines we need also to take into account the impurity ladders. However, as shown below in result the expression for $\Pi$ remains the same as in the clean limit.

Instead of the direct diagram summation it is much faster to calculate the polarization operator in the presence of disorder within the quasiclassical formalism. Let us assume that there is an external driving term in the gap function given by $e^{2i\Omega\tau}F_{\Delta}\hat{\tau}_{1}$. Then in the Eilenberger equation we have the source $i\hat{\tau}_{3}\hat{\tau}_{1}e^{2i\Omega\tau}F_{\Delta}=-\hat{\tau}_{2}ie^{2i\Omega\tau}F_{\Delta}$. Besides that there is a response of the order parameter which we denote $ie^{2i\Omega\tau}\Delta_{2\Omega}\hat{\tau}_{1}$. Then the amplitude of the total off-diagonal driving term in the Eilenberger equation is $\Delta_{2\Omega}+F_{\Delta}$. Then we can use this amplitude instead of $\Delta_{2\Omega}$ in the expression (37) for the corresponding correction to the propagator. The self-consistency equation $\Delta_{2\Omega}=-\lambda T\sum_{\omega}{\rm Tr}\hat{\tau}_{2}\hat{g}_{\Delta}(13)$ yields

[TABLE]

where we denote as before $s_{1}=\sqrt{\omega_{1}^{2}+\Delta^{2}}$ and $s_{3}=\sqrt{\omega_{3}^{2}+\Delta^{2}}$. Taking into account that $\sum_{\omega}(s_{3}^{-1}-s_{3}^{-1})=0$ the expression for polarization operator (45) can be transformed to

[TABLE]

The analytical continuation of this expression to real frequencies as explained in Sec.(III.4) yields the result coinciding with the one obtained previously in the clean limit Volkov and Kogan (1974); Kulik et al. (1981); Littlewood and Varma (1982); Tsuji and Aoki (2015); Cea et al. (2016). Note however, that here we have demonstrated that this expression is valid for arbitrary impurity scattering time.

III.3 Total current

Finally we collect expressions for the current $j=j_{AAA}+j_{H}$ given by the diagrams shown in Fig.2

[TABLE]

In order to find the current at real frequency we need to make analytical continuation of Eqs.(47,48) as described in Sec.III.4. In general, this procedure leads to the expressions which can be handled only numerically and in Sec.VI we show the characteristic dependencies of THG current on various parameters. However, in a number of limiting cases it is possible to treat these expressions analytically which we discuss in Sec.V. In the next Sec. we prove that the results of quasiclassical calculations of the non-linear response coincide with those obtained by the summation of diagrams with current vertices in the presence of impurity self-energy and ladders.

III.4 Analytical continuation

In order to find the real-frequency response we need to implement the analytical continuation of Eq.(47,48). These third-order responses are obtained by the summation of expressions which depend on the four shifted fermionic frequencies such as $g=g(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$. The analytical continuation of the sum by Matsubara frequencies is determined according to the general rule Kopnin (2001)

[TABLE]

where $n_{0}(\varepsilon)=\tanh(\varepsilon/2T)$ is the equilibrium distribution function. In the r.h.s. of (49) we substitute in each term $\omega_{k<l}=-i\varepsilon^{R}_{k}$ and $\omega_{k>l}=-i\varepsilon^{A}_{k}$ for $k=1,2,3,4$, denote $\varepsilon_{k}=\varepsilon+(3-k)\Omega$ and $\varepsilon^{R}=\varepsilon+i\Gamma$, $\varepsilon^{A}=\varepsilon-i\Gamma$. Here the term with $\Gamma>0$ is added to shift of the integration contour into the corresponding half-plane. At the same time, $\Gamma$ can be used as the Dynes parameter Dynes et al. (1984) to describe the effect of different depairing mechanisms on spectral functions in the superconductor. We implement the analytical continuation in such a way that $s(-i\varepsilon^{R,A})=-i\sqrt{(\varepsilon^{R,A})^{2}-\Delta^{2}}$ assuming that the branch cuts run from $(\Delta,\infty)$ and $(-\infty,-\Delta)$.

A special care should given to the differences $s_{i}-s_{j}$ in the denominators of Eq.(26,28,38) for the second- and third-order responses. When analytically continued to the real energies and frequencies these combinations become zero for certain energies. Indeed, e.g. for $i=1$, $j=4$ we have $\varepsilon-2\Omega=-(\varepsilon+\Omega)$ is $\varepsilon=\Omega/2$ so that $s_{1}^{A}=s_{4}^{R}$ for such energy. Thus, the numerical integration of the expressions that contain combinations like $s_{1}^{A}-s_{4}^{R}$ in the denominators is not possible. Fortunately, the expressions (26,28,38) with certain analytical effort can be written in the form which does not contain $s_{i}-s_{j}$ combinations in the denominators. The procedure of how this can be done using the commutation relations (31,32, 33) and (41, 42) for the second-order corrections $\hat{g}_{2s,2a}$ and third-order correction $\hat{g}_{3s,3a}$, $\hat{g}_{A\Delta}$ is given in the Appendix A.

IV Diagram summation

It is the goal of this section to demonstrate that corrections given by the diagrams with current vertices can be calculated using quasiclassical approximation introduced in the previous section. For this purpose we consider impurity ladder diagrams for corrections up to the third order in external field and use them to derive equations for the corresponding momentum-integrated propagators introduced according to Eq.(9).

IV.1 First-order corrections

Let us start with the simplest diagrams for the first-order corrections induced by the interaction with external field $\bm{A}_{\Omega}e^{i\Omega t}$ through the current-type vertex and by the order parameter modulation $\Delta_{2\Omega}e^{2i\Omega t}$ as shown in Fig.(5).

We aim to demonstrate the general approach for deriving equation for the momentum-integrated propagators using the simplest diagram shown in Fig.5a. Let us introduce the notation

[TABLE]

The key idea of the derivation of the simplified equation for $\hat{g}_{A}(12)$ is to use the following trick. Let us multiply the function $\hat{G}_{A}(12)$ by $\hat{G}_{0}^{-1}(1)$ from the left and by $\hat{G}_{0}^{-1}(2)$ from the right, subtract the results and integrate by $\xi_{p}$. We use that Eq.(7) yields the relations $\hat{G}_{0}^{-1}(j)=\tilde{\Delta}_{j}\hat{\tau}_{1}+i\tilde{\omega}_{j}\hat{\tau}_{0}+\xi_{p}\hat{\tau}_{3}$ and $\tilde{\Delta}_{j}\hat{\tau}_{1}+i\tilde{\omega}_{j}\hat{\tau}_{0}=i(s_{j}+1/2\tau_{imp})\hat{g}_{0}(\omega_{j})\hat{\tau}_{3}$, where $s_{j}=\sqrt{\omega_{j}^{2}+\Delta^{2}}$. Then we eliminate off-shell contributions in the momentum integrals to express the result through quasiclassical propagators

[TABLE]

where we introduce the notation $\tilde{s}_{j}=s_{j}+1/2\tau_{imp}$. The obtained expression coincides with the r.h.s. of the Eilenberger Eq.(21) for the correction $\hat{g}_{A}$. Next let us derive the l.h.s. of the equation for $\hat{g}_{A}$. Using the diagram Fig.5a we get that $\hat{G}_{0}^{-1}(1)\hat{G}_{A}(12)-\hat{\tau}_{3}\hat{G}_{A}(12)\hat{G}_{0}^{-1}(2)\hat{\tau}_{3}=\alpha\cos\chi[\hat{G}_{0}(2)-\hat{\tau}_{3}\hat{G}_{0}(1)\hat{\tau}_{3}]$. Then, integrating by $\xi_{p}$ we obtain the equation to determine $\hat{g}_{A}(12)$ coinciding with Eqs.(18, 21).

Let us now consider the momentum integrated correction

[TABLE]

where $\hat{G}_{\Delta}$ is given by the more complicated diagram Fig.5b with ladder insertion. Following the procedure described above we get the r.h.s. of the equation for $\hat{g}_{\Delta}(13)$ in the form similar to (51). To obtain the l.h.s one needs to use the diagrammatic equation Fig.(5c). Using the fact that $\hat{G}_{\Delta}$ does not depend on the momentum direction we write this equation in the algebraic form as follows

[TABLE]

From these equations we obtain $\hat{G}_{0}^{-1}(1)\hat{G}_{\Delta}(13)-\hat{\tau}_{3}\hat{G}_{A}(13)\hat{G}_{0}^{-1}(3)\hat{\tau}_{3}=(\hat{\Delta}_{2\Omega}-i\hat{g}_{\Delta}\hat{\tau}_{3}/2\tau_{imp})\hat{G}_{0}(3)-\hat{\tau}_{3}\hat{G}_{0}(1)(\hat{\Delta}_{2\Omega}\hat{\tau}_{3}-i\hat{g}_{\Delta}/2\tau_{imp})$. Then terms with impurity scattering in l.h.s. and r.h.s. appear to be exactly the same. In result we recover the quasiclassical equation for $\hat{g}_{\Delta}$ with the solution Eq.(37) which is not sensitive to impurity scattering.

IV.2 Second-order correction

In this section we calculate the basic element of the nonlinear response diagrams Fig.2, that is the second-order correction $\hat{G}_{AA}$ modified by the impurity ladder as shown in Fig.6a. We denote the frequencies $\omega_{1}=\omega+2\Omega$, $\omega_{2}=\omega+\Omega$, $\omega_{3}=\omega$, $\omega_{4}=\omega-\Omega$.

The correction $\hat{G}_{AA}$ can be though of as a result of the second-order electron-photon interaction process. It is determined by the equation shown diagrammatically in Fig.6b. Due to the presence of two current vertices the angle average associated with the dashed impurity line produces non-zero result. The corresponding integral equation for $\hat{G}_{AA}$ reads

[TABLE]

where we denote the angular average $\langle f\rangle_{\chi}=\int_{0}^{\pi}d\chi\sin\chi f(\chi)$ and the momentum-integrated GF $\hat{g}_{AA}=i\pi^{-1}\int d\xi_{p}\hat{G}_{AA}$. From the Eq.(54) we can obtain equation for the momentum-integrated GF $\hat{g}_{AA}$. Using the same method as described above in Sec.IV.1 we get the equation for $\hat{g}_{AA}=\hat{g}_{AA}(123)$ which reads as

[TABLE]

Using the ansatz (19) we get from here Eq.(26) for the components $\hat{g}_{2s}$ and $\hat{g}_{2a\\ }$.

IV.3 Third-order correction

The third-order correction to the GF is determined by the diagram shown in Fig. 6c . The corresponding analytical expression reads

[TABLE]

The straightforward calculation of $\hat{G}_{AAA}$ is rather lengthy. However, it is possible to show that the momentum-integrated function $\hat{g}_{AAA}=i\pi^{-1}\int d\xi_{p}\hat{G}_{AAA}$ satisfies Eilenberger Eqs. (23,24). The key idea of this derivation is to use the same trick as introduced above in Sec.IV.1 to rewrite the diagrammatic equation in terms of the quasiclassical propagators. First, we multiply Eq.(56) from the left by $\hat{G}_{0}^{-1}(1)$, from the right by $\hat{G}_{0}^{-1}(4)$ and then subtract the results in the same way as given by the Eq.(51) to obtain the r.h.s. of the Eilenberger equation in the form

[TABLE]

The l.h.s. of the resulting equation can be expressed in terms of the quasiclassical propagators directly from Eq.(56). In this way we obtain

[TABLE]

where $\hat{Y}_{s,a}$ are given by the Eqs.(29,30). Using the ansatz (20) we get from Eq.(58) the solutions for components (28).

IV.4 Corrections due to the Higgs mode

Correction to the GF $\hat{G}_{A\Delta}$ generated by the combined action of the vector potential and the time-dependent order parameter is given by the diagram in Fig.6d. This perturbation determine the non-linear current $j_{H}$ which is sensitive to the excitation of the Higgs mode as shown by the diagram in Fig.2b. As before we are interested in the momentum-integrated function $\hat{g}_{A\Delta}=i\pi^{-1}\int d\xi_{p}\hat{\tau}_{3}\hat{G}_{A\Delta}$ because it determines the correction to the current. Treating the diagram Fig.6d using exactly the same procedure as above we arrive to the equation for $\hat{g}_{A\Delta}$ obtained directly from the Eilenbeger formalism. Its solution is given by the Eq.(38).

V Limiting cases and estimations

V.1 Normal state or very high frequencies $\Omega\gg\Delta$

First, let us check that in the normal state the third-harmonic response disappears. This can be seen from the analytical continuations of the Eq. (47) for $j_{AAA}$ because the Higgs mode-related part $j_{H}$ does not exist in the normal state. Let us introduce the condensed notation e.g. $g(-i\varepsilon_{1},-i\varepsilon_{2},-i\varepsilon_{3},-i\varepsilon_{4})=g(1^{R},2^{R},3^{R},4^{R})$. In the normal state $\hat{g}_{0}(-i\varepsilon^{R,A})=\pm\hat{\tau}_{3}$. Besides that we have $s(-i\varepsilon_{k}^{R})=-i[\varepsilon+(k-3)\Omega]$ and $s(-i\varepsilon_{k}^{A})=i[\varepsilon+(k-3)\Omega]$ so that the sum $s(-i\varepsilon_{k}^{R})+s(-i\varepsilon_{n}^{A})=i(n-k)\Omega$, that is does not depend on energy.

Then for the third-order corrections given by Eqs.(28) one can see that $\hat{g}_{3a}(1^{R},2^{R},3^{R},4^{R})=\hat{g}_{3s}(1^{A},2^{A},3^{A},4^{A})=0$. Moreover, $\hat{g}_{3a}(1^{A},2^{A},3^{R},4^{R})=-2\hat{g}_{3a}(1^{A},2^{R},3^{R},4^{R})=-2\hat{g}_{3a}(1^{A},2^{A},3^{A},4^{R})$. The same relations are true for $\hat{g}_{3s}$. Due to this relations, the contributions from different branch cuts (49) in the expression for the normal state current $j_{AAA}$ cancel identically.

Note that for the frequencies $\Omega\gg\Delta$ the effect of superconducting correlations disappears, so one can consider the system as normal metal. Hence the non-linear response vanishes in this high-frequency limit. This is the reason why diagrams with current vertices can be neglected e.g. when studying Raman response where the frequencies associated with single external lines are assumed to be rather large.

V.2 Absence of impurities $\tau_{imp}=\infty$

Previously it has been noted that in the clean limit, in the spatially homogeneous case and finite frequency the order parameter amplitude is not affected by the irradiation Gorkov and Eliashberg (1968). In that work by Gor’kov and Eliashber only the contribution of diagrams with current vertices has been taken into account. Under the same assumptions non-linear current response also disappears. Here we check our general results for consistency against this limiting case of $\tau_{imp}=\infty$.

Order parameter perturbation. First, let us look at the second-order corrections which determine perturbation of the order parameter $F_{\Delta}$ according to the diagram Fig.1a. In the limit $\tau_{imp}=\infty$ the solutions isotropic second-order correction to GF (26) yields (see Appendix C for details)

[TABLE]

Therefore in the Eq.(43) the sum over frequencies disappears so that

[TABLE]

Now let us consider the contribution of the diagram with density modulation Fig.1b. In this case we obtain

[TABLE]

where ${\rm tr}[\hat{X}]=T\sum_{\omega}\int d^{3}\bm{p}/(2\pi)^{3}{\rm Tr}[\hat{X}]$. The momentum integral in (61) disappears due to the particle-hole symmetryTsuji and Aoki (2015); Cea et al. (2016) which holds up to the corrections of the order $T_{c}/E_{F}\ll 1$. The same conclusion holds in the presence of impurity scattering which modifies propagators and the density vertex as shown in Fig.1b by the shaded region.

Thus, the contribution of density modulation to the Higgs mode excitation always vanishes due to the particle-hole symmetry. Therefore it has been claimed that the Higgs mode does not contribute to the third harmonic generation Cea et al. (2016, 2018). Below we show that in the presence of impurities the non-zero $F_{\Delta}$ is obtained within the quasiclassical approximation and does not contain any such small prefactors .

THG current response. Let us now consider expression for the current $j_{AAA}$ which in the absence of impurities is determined by the third-order correction $\hat{g}_{AAA}$ given by Eqs.(20). The expression which determined the current can be written as

[TABLE]

where $f(123)=f(\omega_{1},\omega_{2},\omega_{3})$ is the function which exact form is rather lengthy and not particularly important. Then due to the Eq.(62) the sum over Matsubara frequencies in the Eq.(47) for current disappears since $\sum_{\omega}f(123)=\sum_{\omega}f(234)$.

V.3 Transition to the clean limit $\tau_{imp}T_{c}\gg 1$.

As shown above in Sec.V.2 the finite-frequency contribution of diagrams with current vertices Fig.1a disappears in the absence of impurity scattering and other relaxation mechanisms. The contribution of diagram with density vertices Fig.1b is zero with the accuracy of the particle-hole symmetry near the Fermi level. So that the contribution of Higgs mode excitation to the THG signal is expected to have the negligible amplitudeCea et al. (2016, 2018) as compared to the direct coupling described by the diagram in Fig.4. Thus the polarization dependence of THG signal is expected to be sensitive to the lattice anisotropy. Here we show that the above conclusions can drastically altered in the presence of rather small impurity concentrations, typical for all realistic superconducting samples, especially in thin film geometries. Our goal is to find the threshold value of impurity scattering when the diagrams with current vertices become dominant and so that the quasiclassical theory described in Sec.II.2 is applicable for calculating non-linear properties.

To understand the magnitude of this threshold in this subsection let us analyse the amplitudes in the regime of low frequencies, impurity scattering rates and small temperatures as compared to the critical temperature $T,\Omega,\tau^{-1}_{imp}\ll T_{c}$. In Sec.VI the numerical results with broad range frequency and scattering rate dependencies are presented.

Order parameter perturbation. First, let us estimate the magnitude of the external order parameter perturbation in the presence of weak disorder. The contribution of diagrams with current vertices Fig.1a to the order parameter perturbation $F_{\Delta}$ can be obtained using solution (26) expanded in the regime $\Omega,\tau^{-1}_{imp}\ll T_{c}$

[TABLE]

Then at small temperatures $T\ll T_{c}$ the first term here vanishes upon the integration over $\omega$ while the second term yields the leading order expansion by the scattering rate

[TABLE]

Thus the Higgs mode amplitude perturbation driven by the linear electron-photon coupling terms in the Hamiltonian is non-zero for the finite concentration of impurities.

THG current response. At non-zero frequencies the current $j_{AAA}$ given by diagrams with current vertices Fig.(2) disappears without impurity scattering and other relaxation mechanisms. At the same time the current $j_{AAA}^{(d)}$ determined by diagrams with density vertices in Fig.4 remains non-zero. Let us find the threshold value of impurity scattering when the latter contribution $j_{AAA}^{(d)}$ can be neglected.

It is convenient to introduce dimensionless amplitudes $I_{AAA}$, $I_{H}$ of the currents (47,48) determined by the contribution of diagrams with current vertices

[TABLE]

Here $j_{0}$ is the normalization current density. In the same way, the contribution of the density modulation diagram in Fig.4 can be written as in terms of the dimensionless amplitude

[TABLE]

Then the dimensionless amplitude is given by

[TABLE]

Hence from the Eqs.(69,65) the ratio of different contributions to the current can be expressed through the ratio of dimensionless amplitudes as follows

[TABLE]

where the prefactor is determined by $j_{0}/j_{0}^{(d)}=T_{c}^{2}/2E_{F}^{2}$ and $E_{F}=mv_{F}^{2}/2$ is the Fermi energy. As we mentioned above, in the limit $\tau_{imp}\to\infty$ the amplitude $I_{AAA}$ disappears while $I^{(d)}_{AAA}$ remains finite. However, due to the very small prefactor in Eq.(70) $(T_{c}/E_{F})^{2}\gg 1$ the threshold level of $\tau_{imp}$ occurs to be quite large.

To understand the magnitude of this threshold let us analyse the amplitudes in the regime of low frequencies, impurity scattering rates and small temperatures as compared to the critical temperature $T,\Omega,\tau^{-1}_{imp}\ll T_{c}$.

Under the above assumption one can obtain the analytical expression for the density modulation-induced amplitude $I^{(d)}_{AAA}$. For small $\Omega$ we can replace $\Omega^{-1}(s_{3}^{-1}-s_{1}^{-1})=-2ds^{-1}_{3}/d\omega=\omega(\omega^{2}+\Delta^{2})^{-3/2}$. So that the current amplitude becomes just $I^{(d)}_{AAA}=1/2$.

The calculation of quasiclassical contribution is more involved. Let us implement the Taylor expansion in terms of the frequency $\Omega$ and the scattering rate $\tau^{-1}_{imp}$ of the Eqs.(26,28). In this way we obtain the leading-order terms

[TABLE]

Since we consider low temperatures, when calculating contributions to the current the frequency summation has to be replaced by the integral, e.g. $2\pi T\sum_{\omega}\hat{g}_{3a}=\int\hat{g}_{3a}d\omega$. Then the first terms in the r.h.s. of expansions (71,72) disappears so we end up with the result

[TABLE]

Thus we get the leading-order contribution to the quasiclassical current amplitude

[TABLE]

where we took into account the relation $\Delta=1.76T_{c}$. Thus from (70) we get the ratio of the density- modulation and quasiclassical currents given by

[TABLE]

Note that Eqs.(74,75) are valid at small frequencies $\Omega\ll\Delta$ and they also agree with the static limit when $\Omega=0$. Indeed, in the static limit $j_{AAA}+j_{H}$ determines non-linear correction to the Meisser current. This correction can be shown to vanish at $T=0$ in the absence of disorder $\tau_{imp}=\infty$ in agreement with Eq..

Based on the estimation (75) the contribution of diagram with density vertex becomes dominating in the limit determined by the condition

[TABLE]

Taking into account that in usual superconductors like NbN with $T_{c}\approx 10K$ and $E_{F}\approx 10^{4}K$ the above condition yields $\tau_{imp}T_{c}>10^{3}$. This criterion means that the superconductor should be in the super-clean regime Kopnin and Volovik (1997); Kopnin and Lopatin (1995) defined as $\tau_{imp}T_{c}>E_{F}/T_{c}\approx 10^{3}$. Up to now the only known system where this regime is realizedKopnin and Lopatin (1995); Kopnin and Salomaa (1991) is the superfluid He3 which generically does not contain any impurities. In solid state systems the certain amount of disorder is always present. Besides that in thin films the scattering time is bounded from above by the time of flight of electrons between interfaces. Taking into account that in the clean limit the coherence length is $\xi=T_{c}/v_{F}$ the above criterion means that the film should be thicker than $100\xi$ which is much larger than what has been used in experiments Matsunaga et al. (2014); Katsumi et al. (2018); Nakamura et al. (2018). Besides that, typical materials used in THz spectroscopy experiments like NbN superconductors usually have strong intrinsic disorder so the dirty limit $\tau_{imp}T_{c}<1$ is realized there even without taking into account the interface scattering.

Another interesting tendency characteristic for the transition to the clean limit $\tau_{imp}T_{c}\gg 1$ is that the Higgs mode-related current is suppressed much strongly than the other component so that $j_{H}\ll j_{AAA}$. This can be understood from the expansion b small parameter $(\tau_{imp}T_{c})^{-1}\ll 1$ of the GF correction $\hat{g}_{A\Delta}$ which determine $j_{H}$ according to the Eq.(48). In this case the amplitude of the order parameter external perturbation is given by (64), that is it already contains the small parameter $(\tau_{imp}T_{c})^{-1}\ll 1$. Besides that expanding $\hat{g}_{A\Delta}$ given by Eqs.(37,38,40) we get

[TABLE]

The first term here vanishes as usual upon the integration by $\omega$ while the second term together with Eqs.(64,44) yields the amplitude of the current $j_{H}$. Here we need to take into account the low-frequency asymptotic of the polarization operator which follows from Eq.(45) $\Pi(\Omega\ll\Delta)=1-\lambda$. The we get the non-linear current generated due to the Higgs mode excitation with the dimensionless amplitude given by

[TABLE]

where at low temperatures $(T_{c}/\Delta)^{4}\approx 0.1$. Thus one can see that the suppression of $I_{H}(\tau_{imp})$ in the transition to clean case is determined by the second order of the small parameter $(\tau_{imp}T_{c})^{-1}\ll 1$. Therefore in this limit $I_{H}\ll I_{AAA}$ except of the vicinity if the Higgs mode resonance at $\Omega=\Delta$ where the amplitude $I_{H}$ is enhanced by the factor $\sqrt{\Delta/\Gamma}$, where $\Gamma$ is the Dynes parameter. However, if the impurity scattering is sufficiently weak $(\tau_{imp}T_{c})^{-1}<\sqrt{\Gamma/\Delta}$ the direct contribution to nonlinear current $I_{AAA}$ dominates for all frequencies. These different regimes are illustrated by the numerical results below in the sectionVI.

Estimations that we provided above rule out the necessity to consider the contribution of density-vertex diagrams to describe the non-linear properties of known superconducting materials. Besides that the thin-film samples used for the non-linear response studies in the THz regime are generically in the dirty regime, because the mean free path is bounded from above by thickness because the surface scattering of electrons mimics impurity scattering. Therefore in realistic superconducting samples where the impurity scattering rate is always above the superclean limit the non-linear response is determined only by the diagrams with current vertices Fig.2. As shown above this means that one can use the quasiclassical approximation with impurity collision integrals which yields major simplification as compared to the direct summation of diagrammatic series. Besides that usually the low-temperature superconductors are in the dirty regime $\tau_{imp}T_{c}\ll 1$ which can be treated within even simpler Usadel theory as discussed in the next section. However, the general solutions we obtain can be applied with some modifications to calculate nonlinear responses in clean materials like the d-wave high temperature superconductors Katsumi et al. (2018) or iron-based superconductors which are in the regime $\tau_{imp}T_{c}>1$.

V.4 Dirty limit $\tau_{imp}T_{c}\ll 1$

Previously, non-linear electromagnetic properties of superconductors have been studied mostly in the diffusive system using the Usadel formalism Gorkov and Eliashberg (1968, 1969); Gor’kov and Eliashberg (1968); Eliashberg (1970); Ivlev and Eliashberg (1971); Moor et al. (2017); Tikhonov et al. (2018). The general case with arbitrary impurity concentration has been discussed in the linear response regimeArtemenko and Volkov (1979). Here we present for the first time calculations of the non-linear responses with arbitrary impurity scattering time. Therefore it is important to establish connection with the Usadel theory results which should be obtained from our general expressions in the limit $\tau_{imp}T_{c}\ll 1$. That is, the isotropic second-order correction $\hat{g}_{2s}$ which determines the order parameter perturbation and the expression for non-linear current response can be obtained directly from Eqs.(13,14,15) as $(\alpha^{2}D/v_{F}^{2})g_{2}(123)$ where

[TABLE]

where we normalize the current by the amplitude $j_{3\Omega}=\pi e^{2}(A_{\Omega}/c)^{3}D\sigma$. The correction $\hat{g}_{\Delta}$ induced by the order parameter oscillation is given by Eq.(37) and does not depend on the impurity scattering rate. In the Appendix (B) we demonstrate that the same expressions follow from the general Eqs.(26,47) as the leading term expansions by the small parameter $\tau_{imp}T_{c}\ll 1$.

VI Numerical results and discussion

Analytical estimations in previous section are obtained at small frequencies $\Omega\ll\Delta$. To understand the full frequency dependencies we implement the numerical calculation using the analytical continuation procedure explained in Sec. III.4. At first, our aim is to compare the dimensionless amplitudes of the three contributions to the current $I_{AAA}$, $I_{H}$ and $I_{AAA}^{(d)}$ for different values of the parameters. The results are presented in Fig.7 for different values of the scattering rate varying between the clean and dirty regimes. Here we consider the regime of small temperatures $T=0.1T_{c}$ and the Dynes parameter is $\Gamma=0.01$. These plots confirm qualitative conclusions made above. In the clean case $(\tau_{imp}T_{c})^{-1}=0.5$ (Fig.7a) the Higgs contribution is much smaller than the direct coupling one $I_{H}\ll I_{AAA}$. This is despite the fact that Higgs mode contribution is resonant and its maximal value scales like $\sqrt{\Delta/\Gamma}$. The $I_{AAA}$ contribution also has peaks both as $\Omega=\Delta$ and $\Omega=2\Delta/3$ although their amplitude does not diverge with $\Gamma\to 0$. The origin of these peaks is the BCS density of states singularity at the energy $\varepsilon=\pm\Delta$. The external radiation at frequencies $\Omega=\Delta,2\Delta/3$ causes transition between these states with the enhanced probability due to the large density of states.

With increased scattering rate the general amplitude of the Higgs mode contribution rises so that its maximal value becomes much larger than the direct contribution. For the considered value of Dynes parameter $\Gamma=0.01$ this happens at rather large scattering $(\tau_{imp}T_{c})^{-1}\approx 10$ as shown in Figs.7d,e. At such parameters $I_{H}\approx 10^{-3}I_{AAA}^{(d)}$ which means that $j_{H}\gg j_{AAA}^{(d)}$ if one recalls the overal factor $(T_{c}/E_{F})^{2}\approx 10^{-6}$ in front of the density modulation current (70). The non-resonant peak in $I_{AAA}$ at $\Omega=\Delta$ remains although the one at $\Omega=2\Delta/3$ is eliminated by the impurity scattering. The overall dependencies of $I_{AAAA}$ and $I_{H}$ on frequency and scattering rates as shown in 7c,f. One can see that first they increase with $(\tau_{imp}T_{c})^{-1}$ but then start to decrease. We discuss the regime with large scattering rate below using the diffusive limit approach .

In thin-film NbN samples used for the non-linear response measurement the amount of scattering is typically rather large. Therefore we consider the diffusive limit below using the Usadel theory results from Sec.V.4 to calculate non-linear currents. Besides that since NbN has strong electron-phonon coupling Matsunaga et al. (2017) resulting in the enhanced inelastic relaxtion we use larger value of Dynes parameter $\Gamma=0.1$ as compared to the previous example. The results for temperature and frequency dependence of currents $j_{AAA}$ and $j_{H}$ are shown in Fig.8. As one can see in Fig.8a,b both currents $j_{AAA}$ and $j_{H}$ have peaks at $\Omega=\Delta(T)$. However, the peak value of Higgs contribution is several times larger. It occurs when the denominator in Eq.(44) reaches it minimal value of $1-\Pi(2\Delta)\approx\sqrt{\Gamma/\Delta}$. So the Higgs-mode related part of non-linear current has the same maximal amplitude $j_{H}(\Omega=\Delta)\propto j_{3\Omega}\sqrt{\Delta/\Gamma}$ as compared the $j_{AAA}(\Omega=\Delta)\propto j_{3\Omega}$ part. These estimation agrees with the numerical result in Fig.8c,d.

For higher values of the Dynes parameter which can mimic the enhanced inelastic relaxation in superconductors with strong electron-phonon interaction like NbN Matsunaga et al. (2017) the resonant Higgs peak broadens and decreases. This tendency is illustrated in Fig.8. The broadened peaks featured by the temperature dependencies of $j_{H}(T)$ at fixed frequencies (Fig.8f) are similar to those obtained in the experiment Matsunaga et al. (2014). At the same time the dependencies of $j_{AAA}(T)$in Fig.8e show no peaks at all.

As we discussed above, the presence of impurities triggers the non-linear response and Higgs mode generation. However, as one can see from the sequence of the plots in Fig.7c,f above certain threshold scattering rate the amplitudes $I_{H}$ and $I_{AAA}$ start to decrease with decreasing $\tau_{imp}$. This agrees with the diffusive limit results Eqs.(80, 81) where the currents are determined by the amplitude has an additional small parameter $(\tau_{imp}T_{c})^{2}$ as compared to the prefactor in (68) so that $j_{3\Omega}\propto e\nu(eA/c)^{3}(\tau_{imp}E_{F})^{2}$. At the same time the density modulation-related current $j_{AAA}^{(d)}$ is not sensitive for disorder and therefore should dominate in the very dirty system as well as in the very clean one. Thus comparing $j_{3\Omega}$ with the prefactor in Eq.(68) we obtain that in the diffusive limit $j_{AAA}^{(d)}$ is negligible as long as $\tau_{imp}E_{F}>>1$ and start to dominate in the opposite case, that is close to the localization threshold.

The technique developed in the present work for the arbitrary scattering rates can be applied as well to study Higgs mode generation in non-trivial superconductors like those with the d-wave symmetry of the order parameterKatsumi et al. (2018). Within quasiclassical formalism such states are described with the help of the anisotropic pairing constant $\lambda_{d}(\theta,\theta^{\prime})=\lambda\sin(2\theta)\sin(2\theta^{\prime})$. Here $\theta$ and $\theta^{\prime}$ are the angles corresponding to the momenta of interacting electrons. Correspondingly the order parameter acquires momentum dependence $\Delta(\theta)\propto\sin(2\theta)$ which should be taking into account when solving the Eilenbeger equation in d-wave superconductor. With that one can see that in the absence of impurities the same limitations on the non-linear response pertain as for the s-wave superconductor. That is, in the completely pure system the second-order electron-photon coupling through the linear terms $V_{1}$ in the Hamiltonian (2) does not excite Higgs mode and does not produce any THG signal. The presence of impurities certainly helps the situation although their effect is a bit more tricky than in the isotropic s-wave considered here. However since d-wave pairing can be found only in the clean regime $(T_{c}\tau_{imp})^{-1}\ll 1$ one can expect that the amplitude of Higgs mode contribution should be strongly suppressed according to the Fig.7a. However, to figure out the resulting amplitude it is necessary to figure out the value of Dynes parameter which can be much smaller than in low-temperature superconductors thus allowing for large $I_{H}$ peaks even in the clean system. At the same time we don’t expect the direct coupling current $j_{AAA}$ to feature pronounced peaks because of the lack of the density of states singularities in the d-wave superconductor. Thus the influence of impurities on the collective modes in superconductors with non-trivial pairing is potentially very interesting although its detailed study is beyond the scope of the present paper.

Another interesting direction which can be addressed using the formalism developed by us is the nonlinear response and generation of collective modes in multiband superconductors like MgB2 and iron-pnictide compounds. Here in addition to the impurity scattering important effects can be related to the interband tunnelling of quasiparticles and Cooper pairs which should significantly affect non-linear response. With that we can address recent experimental results on the THz pump-probe experiments with MgB2 Giorgianni et al. (2019).

VII Conclusions

We have studied non-linear electromagnet response of superconductors with the amount of disorder varying between the completely pure limit and the dirty regime. The impact of our study is threefold. First, it it demonstrated that the quasiclassical approximation allows for the correct description of the non-linear effects in superconductors coupled to the external electromagnetic field. Propagators obtained by solving quasiclassical Eilenberger equation with the impurity collision integral coincide with those obtained by the direct summation of diagrams with current vertices taking into account the impurity self-energies and ladders.

Second, we demonstrated that effective Higgs mode excitation is possible in usual BCS superconductors without any extensions of the model suggested in previous works. We show that the contribution of diagrams with current vertices start to dominate over the density-modulation related processes for the level of disorder above the extremely weak superclean threshold. Since the superclean regime is hardly realizable in experiments our results provide the basis for the analysis of non-linear responses in realistic superconducting samples to describe the pump-probe or the THG generation experiments in the broad range of frequencies. The same conclusion holds for compounds with unconventional pairing such as the d-wave cuprates Katsumi et al. (2018) or multiband superconductorsAkbari et al. (2013). Theory suggested in the present paper can be applied to analyse the recent data on the Higgs mode in a d-wave superconductorKatsumi et al. (2018) and the collective mode in MgB2 Giorgianni et al. (2019). In general the impurity scattering determines collective mode excitation in superconductors with both the s-wave and non-trivial pairings and should modify the Higgs mode spectroscopy approach which has been suggested recentlyFauseweh et al. (2017)

Third, we have demonstrated that in the diffusive regime which is typical for thin-film superconducting samples used for recent THz measurements the resonant contribution to THG signal is determined by the Higgs mode excitation thus providing the natural explanation of recent experimentsMatsunaga et al. (2014, 2017). The amplitude of this peak is bounded from above by the Dynes parameter which describes the quasiparticle recombination rate due to the electron-phonon interaction.

VIII Acknowledgements

I am grateful to Lara Benfatto for useful discussions. This work was supported by the Academy of Finland (Project No. 297439).

Appendix A Removing singularities in Eqs.(26,28,38)

Eqs. (26) for the second- and for the third-order (28,38) corrections contain the differences $s_{i}-s_{j}$ in the denominators. When analytically continued this difference becomes zero for certain energy. Therefore these equations cannot be used directly for the numerical integration along the real energy axis. In order to make them suitable for numerics they Eqs. should be written in such a way to eliminate singular $(s_{i}-s_{j})$ combinations in the denominators.

First, we note that one can simplify the expression for the second-order corrections as follows. Substituting Eq. for the first-order correction $g_{1a}$ into the Eqs. (26,28,29) we can rewrite them as follows

[TABLE]

[TABLE]

where we denote

[TABLE]

With the third-order it requires more effort to get rid of the differences $s_{1}-s_{4}$ in denominators. First, let us demonstrate how this can be done for the correction $\hat{g}_{\Delta A}$ that determine the Higgs mode contribution to the current (38). We rewrite it as follows $\hat{g}_{\Delta A}=\hat{g}^{(1)}_{\Delta A}+\hat{g}^{(2)}_{\Delta A}$, where

[TABLE]

where $\hat{Y}_{A\Delta}=\hat{Y}_{\Delta}-I\hat{Y}_{A}$. Now we need to treat only the term $\hat{g}_{A\Delta 2}$. In order to do that we note the relation

[TABLE]

Using the commutation relationw (41, 42) we obtain then the expression for $\hat{g}^{(2)}_{A\Delta}$ without singularity

[TABLE]

Now let us apply the same trick to the corrections $\hat{g}_{3s,3a}$. We write them as the superposition of two parts, e.g. $\hat{g}_{3s}=\hat{g}^{(1)}_{3s}+\hat{g}^{(2)}_{3s}$

[TABLE]

Using commutation relations (31,32,33) we rewrite expression for $\hat{g}^{(2)}_{3s}$ in the following form

[TABLE]

that does not have singularities. Expressions for the anisotropic part $\hat{g}_{3a}=\hat{g}^{(1)}_{3a}+\hat{g}^{(2)}_{3a}$ are similar to Eqs.(90, 92) with the change of $\hat{g}_{2s}$ by $\hat{g}_{2s}$ and $\hat{Y}_{s}$ by $\hat{Y}_{a}$ .

Appendix B Derivation of the response in diffusive limit

First, in the dirty limit we can find corrections to the propagators directly from the Usadel equation which is a simplified version of the Eilenberger equation. In the frequency domain

[TABLE]

we get from Eq.(13,14,15)

[TABLE]

where we denote again $\omega_{1}=\omega+2\Omega$, $\omega_{2}=\omega+\Omega$, $\omega_{3}=\omega$, $\omega_{4}=\omega-\Omega$. The solution can be written as (79)

[TABLE]

Taking into account the diffusion coefficient $D=v_{F}^{2}\tau_{imp}/3$ the solution for $\hat{g}_{2}$ coincides with $\hat{g}_{2s}$ obtained from the general expression (26) up to the leading term in $\tau^{-1}_{imp}$.

Within the Usadel theory the current can be calculated using general expression (15) to have the form (80,81)

[TABLE]

Taking the dirty limit for general Eq.(47) that determines the current is more tricky. The leading-order correction in the limit $\tau_{imp}\to 0$ is given by (28). First, from the Eq.(25) we find for the first-order corrects in the limit $\tau_{imp}\to 0$ given by

[TABLE]

Using this relation and the commutation relations $\hat{g}_{2s}(234)\hat{g}_{0}(4)=-\hat{g}_{0}(2)\hat{g}_{2s}(234)$ we obtain

[TABLE]

Thus from Eq.(28) and Eqs.(90,92) we get

[TABLE]

[TABLE]

We substitute the result (98,99) into the expression for current (47) and use the summation $\sum_{\omega}\hat{g}_{0}(2)\hat{g}_{2s}(234)=-\sum_{\omega}\hat{g}_{2s}(123)\hat{g}_{0}(3)$ to get the expression for the current which is equal to the Eq.(80).

The dirty limit for Higgs mode-related part of the current can be obtained from Eq.(48,37,38,40) in the similar way as above. Using the relation $\hat{g}_{\Delta}(24)\hat{g}_{0}(4)=-\hat{g}_{0}(2)\hat{g}_{\Delta}(24)$ we obtain

[TABLE]

Thus from Eq.(38) and Eqs.(86, 89) we get

[TABLE]

[TABLE]

Using the summation $\sum_{\omega}\hat{g}_{0}(2)\hat{g}_{\Delta}(24)=-\sum_{\omega}\hat{g}_{\Delta}(13)\hat{g}_{0}(3)$ and the expression for current (48) we get the dirty limit expression for the current (96).

Appendix C Absence of the Higgs mode generation without impurities

Without disorder we get from Eq.(82) the expression which according to Eq.(43) determines the order parameter (82)

[TABLE]

Using the relations

[TABLE]

we get

[TABLE]

Substituting this result into Eq.(103) we get as required by the Eq.(59)

[TABLE]

As explained in the main text since the summation by Matsubara frequencies of this expression yields zero, it yields no perturbation of the order parameter amplitude.