Higgs-mode resonance in third harmonic generation in NbN superconductors: multiband electron-phonon coupling, impurity scattering, and polarization-angle dependence

We theoretically investigate the resonance of third harmonic generation (THG) that has been observed at frequency being half of the superconducting gap in a multiband disordered superconductor NbN. The central question is whether the dominant contribution to the THG resonance comes from the Higgs mode (the collective amplitude mode of the superconducting order parameter) or quasiparticle excitations. To resolve this issue, we analyze a realistic three-band model with effective intraband and interband phonon-mediated interactions estimated from first-principles together with nonmagnetic impurity scatterings. Using the estimated pairing interactions with multiband impurity scattering rates being varied from clean to dirty regimes, we calculate the THG susceptibility for NbN in a channel-resolved manner by means of the BCS and self-consistent Born approximations. In the dirty regime, which is close to the experimental situation, the leading contribution is given by the paramagnetic channel of the Higgs mode having almost no polarization-angle dependence, while the second leading contribution comes from the paramagnetic channel of quasiparticles generally showing significant polarization-angle dependence. The result is consistent with the recent experimental observation of no polarization-angle dependence of THG, giving firm evidence that the Higgs mode dominantly contributes to the THG resonance in NbN superconductors.

We theoretically investigate the resonance of third harmonic generation (THG) that has been observed at frequency being half of the superconducting gap in a multiband disordered superconductor NbN. The central question is whether the dominant contribution to the THG resonance comes from the Higgs mode (the collective amplitude mode of the superconducting order parameter) or quasiparticle excitations. To resolve this issue, we analyze a realistic three-band model with effective intraband and interband phonon-mediated interactions estimated from first-principles together with nonmagnetic impurity scatterings. Using the estimated pairing interactions with multiband impurity scattering rates being varied from clean to dirty regimes, we calculate the THG susceptibility for NbN in a channel-resolved manner by means of the BCS and self-consistent Born approximations. In the dirty regime, which is close to the experimental situation, the leading contribution is given by the paramagnetic channel of the Higgs mode having almost no polarization-angle dependence, while the second leading contribution comes from the paramagnetic channel of quasiparticles generally showing significant polarization-angle dependence. The result is consistent with the recent experimental observation of no polarization-angle dependence of THG, giving firm evidence that the Higgs mode dominantly contributes to the THG resonance in NbN superconductors.

I. INTRODUCTION
The standard microscopic theory of superconductivity, i.e., the BCS theory, predicts the presence of the collective amplitude mode of the superconducting order parameter [1][2][3][4][5][6], which is recently referred to as the Higgs mode due to the close analogy with the Higgs boson in particle physics (for recent reviews, see [7,8]). Despite the fundamental and universal aspects of the Higgs mode, its observation in ordinary superconductors had been elusive until recently. One exception was a superconductor 2H-NbSe 2 , which is special in the sense that superconductivity and charge density wave (CDW) coexist in a single material. In this particular situation, the Higgs mode becomes Raman active, and has been observed in the early stage by Raman experiments [9,10] (see also [11][12][13] for recent studies). However, the Higgs mode itself should exist irrespective of CDW, so that its observation in superconductors without any other orders has been long awaited.
The difficulty in observing the Higgs mode in superconductors without other coexisting orders is that the Higgs mode does not linearly couple to external electromagnetic fields, and that the energy of the Higgs mode, which lies around the superconducting gap energy 2∆, is in the terahertz (THz) frequency range, for which an intense light source had been lacking for a long time. The recent development of THz laser techniques, however, has made it possible to excite the Higgs mode directly through the nonlinear light-Higgs coupling [14]. In fact, coherent oscillation of the superconducting order parameter with frequency 2∆ after irradiation with a monocycle THz pulse has been observed in a superconducting NbN [15]. Subsequently, resonant enhancement of third harmonic gen-eration (THG) at the condition of 2Ω = 2∆ with Ω being the incident light frequency has been reported for NbN using multicycle THz pulses [16].
While all these measurements are consistent with the interpretation that the Higgs mode is excited by THz laser excitations, it is not sufficient to confirm that the mode energy is 2∆, since the pair-breaking energy of quasiparticles is also equal to 2∆. This forces one to distinguish the collective Higgs mode from individual excitations of quasiparticles by properties other than the mode energy.
One way to discriminate them is to measure the polarization-angle dependence of the resonant THG [17]. According to the BCS mean-field calculation in the clean limit for a single-band model, the quasiparticle contribution has strong angle dependence in THG, whereas the Higgs-mode contribution does not. Followed by the theoretical proposal, the polarization-resolved measurement of THG has been performed for a single-crystal NbN, showing that the THG intensity at the resonance has almost no polarization-angle dependence [18]. Does this mean that the origin of the resonant THG observed in NbN is the Higgs mode?
The story is not so simple, because the BCS clean limit calculation also suggests that the absolute magnitude of the quasiparticle contribution to the THG resonance is generally much larger than that of the Higgs mode in the BCS clean limit [8,17]. Considering both the polarization-angle dependence and absolute magnitude of the Higgs and quasiparticle contributions to the THG, we come to the conclusion that at least the BCS mean-field treatment in the clean limit fails to describe the THG experiments for NbN superconductors.
There are several possibilities to circumvent this controversial situation: One is to go beyond the BCS approx-imation and include, e.g., phonon retardation effects. In fact, NbN is known to have a moderately strong electronphonon coupling (with a dimensionless coupling constant λ ∼ 1) [19][20][21]. Based on the nonequlibrium dynamical mean-field theory [22], it has been shown that the Higgs mode can contribute to THG with an order of magnitude comparable to quasiparticles [23].
Another possibility is to depart from the clean limit and consider the effect of impurity scattering. Since the optical conductivity of NbN used in the THz laser experiments agrees well [15,16] with the Mattis-Bardeen form [24], the NbN samples are close to the dirty limit. The effect of impurity scattering on THG in BCS superconductors has been studied in Refs. [25][26][27]. Strikingly, in the dirty regime the magnitude of the Higgs-mode contribution to THG can exceed by far that of quasiparticles.
These studies suggest that impurity scattering has more substantial effects on the magnitude of the THG resonance than phonon retardation. It is natural to expect so, because the impurity scattering rate generally approaches a nonzero constant value in the low-energy limit, while the electron-phonon scattering rate decays to zero in the Fermi-liquid regime. Hence, in the present study we focus on the effect of impurity scattering. A crucial key to the open issue of which of the Higgs mode or quasiparticles are dominant in the THG resonance in NbN is the polarization-angle dependence of THG in the dirty regime of superconductors, which has not been addressed so far.
In this paper, we study the polarization-angle dependence of THG in NbN superconductors with impurities. For this purpose, we use an effective three-band model including the phonon-mediated multiband pairing interactions for NbN. In particular, we take special care of the relative magnitude between the intraband and interband pairing interactions, since the polarization-angle dependence might be strongly affected by it. In the previous study [18], the calculations of the polarization angle dependence for the three-band model have been performed at the BCS clean limit without first-principles estimate of the pairing interactions. The calculations assuming the same amplitude of the intraband and interband pairing interaction have shown that the Higgs-mode contribution in THG is isotropic, while the quasiparticle contribution has significant angle dependence [18]. For other choices of the relative magnitude between the intraband and interband interaction parameters, the Higgs mode can also exhibit the polarization-angle dependence [28]. In the present study, we go beyond these previous studies by taking into account the effect of impurities with a realistic estimate of the ratio between the intraband and interband interactions.
We first estimate the pairing interaction parameters for NbN from first principles calculations of the phonon band structure and the electron-phonon couplings. Using the estimated ratio between the intraband and interband interaction parameters, we calculate the THG susceptibility for the multiband superconductor NbN within the BCS mean-field theory. The effect of nonmagnetic impurity scattering is treated by means of the self-consistent Born approximation. We consider both intraband and interband impurity scatterings [29] for multiband NbN superconductors. The calculated THG susceptibility is classified according to the physical origin (quasiparticle or Higgs mode), the coupling channel to light (diamagnetic or paramagnetic), and the diagrammatic representation in the presence of impurities.
The results show that the THG resonance is dominated by the paramagnetic channel in the dirty regime in NbN, in which the Higgs-mode contribution generally becomes larger than the quasiparticle contribution. This behavior is similar to the previous results for single-band superconductors [25,27]. With the estimated ratio between the intraband and interband pairing interactions, the quasiparticles always show clear polarization-angle dependence of THG, while the Higgs mode does not in general, except in the vicinity of the parameter region where the interband impurity scattering rate vanishes. By comparing with the polarization-resolved THG measurements for NbN [18], we conclude that the dominant contribution to the THG resonance is coming from the Higgs mode rather than quasiparticles.
The paper is organized as follows. In Sec. II, we evaluate the electron-phonon couplings and the effective pairing interactions in NbN from first principles calculations. In Sec. III, we describe the method to calculate the THG susceptibility using an effective three-band model for NbN with multiband pairing interactions and impurity scatterings. In Sec. IV, we show the numerical results for THG in NbN superconductors, focusing on its magnitude and polarization-angle dependence of the Higgs-mode and quasiparticle contributions for various impurity scattering rates. The paper is summarized in Sec. V.

II. FIRST PRINCIPLES ESTIMATION OF THE ELECTRON-PHONON COUPLING IN NbN
In this section, we evaluate the intraband and interband effective pairing interactions of NbN from first principles, which are important to determine the polarization-angle dependence of THG in multiband NbN superconductors as discussed in the introduction. The electronic band structure of NbN has been calculated from first principles in the previous literatures [18,[30][31][32][33]. The ab initio estimate of the phonon band structure and the electron-phonon coupling constant of NbN has been reported in [34][35][36][37]. While the total effective pairing interaction (summed over the band indices) has been derived in the previous calculations, here we need the band-resolved matrix elements of the pairing interaction.
Our approach is based on the ab initio construction of a low-energy effective model of NbN including the electron-phonon coupling. To this end, we perform the density functional calculation for NbN using Quantum ESPRESSO package [38,39]. We use the Troullier-Martins norm-conserving pseudopotentials [40] in the Kleinman-Bylander representation [41] with the Perdew-Burke-Ernzerhof [42] exchange-correlation functional. We set the cutoff energy for the wave functions and charge density to be 100 eV and 400 eV, respectively, and take N k = 8×8×8 k points for electron's momentum mesh.
The phonon band structure and the electron-phonon coupling constants are evaluated by the density functional perturbation theory [43], for which we use N q = 8×8×8 q points for phonon's momentum mesh. The previous phonon band calculation [35,36] shows that NbN in the NaCl-type structure has a structural instability as indicated by imaginary phonon frequencies, which is, however, not observed in experiments. To avoid such an instability, we employ a virtual crystal approximation, where we create a pseudopotential for Nb with the nuclear charge Z = +40.5. With this, we fully optimize the lattice structure, obtaining the lattice constant a = 4.497 A, which agrees well with the experimental data [44].
Our calculation of the electronic band structure of NbN (red curves in Fig. 1) well reproduces the previous results [18,[30][31][32][33]. Near the Fermi energy, there are three bands consisting of Nb's 4d t 2g orbitals (xy, yz and zx), which are occupied by two electrons in one unit cell in average. Therefore, NbN can be effectively regarded as a threeband system at one third filling at low energy. We first construct the effective three-band tight-binding Hamiltonian on the basis of the maximally localized Wannier orbitals [45,46], for which we use the open-source package RESPACK [47]. We can simplify the effective three-band model by taking the leading hopping processes [18], where c † knσ is a creation operator of electrons with momentum k, orbital n, and spin σ. The simplified energy dispersion for the d xy orbital is given by The remaining band dispersions k,yz and k,zx are given by permuting x, y, and z in k,xy . The three hopping parameters t, t , and t are fitted with the three bands constructed from the maximally localized Wannier orbitals. The results are t = −0.79 eV, t = −0.22 eV, and t = 0.19 eV, which slightly deviate from the previous result in Ref. [18]. The difference arises because the previous study fit the expression (2) directly with the original band structures (corresponding to red curves in Fig. 1), while here we fit Eq. (2) using the band dispersion of the maximally localized Wannier orbitals. In Fig. 1, we show the simplified band dispersions kn by blue dots. One can see that both of the band dispersions constructed from the density functional calculation and from the simplified model (2) agree fairly well with each other. We employ the simplified dispersion kn for the model-based calculation of THG in Sec. IV, where we need higher-order derivatives of the dispersion such as ∂ 2 kn ∂ki∂kj that can be analytically evaluated with the expression (2). I. The magnitude of the static part (ω = 0) of the Fermi-surface(FS)-averaged phonon-mediated attractions estimated from first principles for NbN. The column "Broadening" shows the broadening width η of the delta functions used to evaluate the FS average. Vintra and Vinter denote the intraband and interband attractions, respectively (see the text for details).

Broadening [Ry]
Vintra where ω qν is the phonon frequency, and b † qν is the creation operator of phonons with momentum q at νth branch. There are six phonon modes (ν = 1, . . . , 6) in total, corresponding to two atoms (Nb and N) in the unit cell each of which can oscillate along three orthogonal directions (x, y, and z). In Fig. 2, we plot the phonon band dispersion ω qν of NbN obtained from the first principles calculation. Three of them are acoustic phonons with linear dispersions around Γ point, while the rest are optical phonons with energy gaps. Here we do not see imaginary phonon frequencies, implying that the present lattice structure is dynamically stable.
The electron-phonon coupling term is written as where g ν mn (k, q) represents the matrix elements of the multiband electron-phonon coupling constant estimated from the density function perturbation theory. The calculations of the matrix elements are done in the maximally localized Wannier orbital basis [48].
The electron-phonon coupling mediates effective attraction between the electrons, whose magnitude is given by ν |g ν mn (k, q)| 2 D ph qν (ω) with the phonon propagator D ph qν (ω). We simplify the form of attraction and model it by the BCS-type attraction form as Here, V nn ≡ V intra and V mn ≡ V inter (m = n) denote the intraband and interband effective attractions, respectively. To estimate the realistic ratio between V intra and V inter , we compute the static part (ω = 0) of the Fermisurface(FS)-averaged phonon-mediated attraction as where D( F ) is the density of states for each t 2g orbital at the Fermi energy (by symmetry, the density of states is the same among t 2g orbitals), w kn is the weight of n maximally localized Wannier orbital at the Fermi energy at momentum k. In practical numerical calculations, w kn is calculated using the Gaussian smearing e −x 2 /2η 2 / √ 2πη with a broadening width η.
In Table I, we list the magnitude of the static part of the FS-averaged phonon-mediated attraction V intra FS and V inter FS for NbN. We note that the density of states D( F ) is about 0.12 states/eV and that the coupling constant λ = D( F ) ( V intra FS + 2 V inter FS ) amounts to ∼ 1, in accord with the previous estimates [34][35][36][37]. The important quantity in the following THG calculations is the ratio between the intraband and interband interactions V inter /V intra in Eq. (5). We find that the first-principles estimate of the ratio V inter FS / V intra FS is about 0.17-0.18. Referring to the ab initio value, in the following sections, we set the ratio V inter /V intra to be 0.18.

III. METHOD FOR THE CALCULATION OF THIRD HARMONIC GENERATION
Having evaluated the effective intraband and interband pairing interactions for NbN evaluated in the previous section, we now move on to the calculation and classification of the THG susceptibility for NbN with impurities. Our method is based on the BCS mean-field approximation. For the treatment of impurities, we employ the self-consistent Born approximation, which is valid in the weak disorder case (i.e., the impurity scattering rate is much smaller than the hopping but can be larger than the superconducting gap). In the calculation of the THG susceptibility, we need to take into account the vertex corrections represented by impurity ladder diagrams [25,27,49].

A. Formalism
Let us consider a multiband system described by the BCS Hamiltonian with the intraband and interband pairing interactions and nonmagnetic impurity scatterings, where v imn represents the impurity potential that hybridizes m and n bands at a lattice site i. We assume that v imn is a Gaussain random variable with the disorder average given by v imn v i m n disorder = γ 2 mn δ ii δ mm δ nn . Here the impurity scattering rate is parametrized by the intraband and interband ones, γ intra = γ nn and γ inter = γ mn (m = n), respectively.
To calculate real-frequency spectra, we introduce nonequilibrium (retarded, advanced, and lesser) Green's functions for multiband superconductors, where Ψ † kn = (c † kn↑ c −kn↓ ) is the two-component Nambu spinor, a, b = 1, 2 are the Nambu space indices, and θ(t) is the step function [θ(t) = 1 (t ≥ 0) and θ(t) = 0 (t < 0)]. The retarded Green's function satisfies the Dyson equation, where δ is a positive infinitesimal constant, ξ kn = kn −µ, µ is the chemical potential,τ 0 is the unit matrix,τ i (i = 1, 2, 3) are Pauli matrices in the Nambu space, and Σ R n (ω) is the retarded self-energy. We put a hat on a matrix which has Nambu spinor indices. The advanced and lesser Green's functions are given byĜ A f (ω) = 1/(e βω +1) is the Fermi distribution function and β = (k B T ) −1 is the inverse temperature.
In the BCS and self-consistent Born approximations, the retarded self-energy is determined bŷ where ∆ n is the superconducting gap for the band n defined by The superconducting gap, self-energy, and Green's functions are self-consistently determined by Eqs. (11), (12), and (13). In Fig. 3, we show the diagrammatic representation for the Dyson equation in the BCS mean-field and self-consistent Born approximations. In order to evaluate the third harmonic generation, we employ the field-derivative approach developed in Ref. [23], which allows one to systematically derive nonlinear optical susceptibilities. The idea is to analytically differentiate the current, with respect to the amplitude of the external field A(t) = eAe −iΩt , where v k,n = ∂ kn ∂k is the group velocity, e is the unit polarization vector ( e = 1), A and Ω are the amplitude and frequency of the field. In this process, we repeatedly differentiate the nonequilibrium extension of the self-consistent equations (11), (12), and (13) in the presence of the external field A(t).
For the THG susceptibility, we need to calculate the second derivative of the self-energy and superconducting gap,Σ n (ω) = ∂ 2 ∂A 2Σn (ω) and∆ n = ∂ 2 ∂A 2 ∆ n [23] (the derivative with respect to A is denoted by dots), which are self-consistently determined by the doubly differentiated equations. By taking the second derivative of the nonequilibrium extension of Eq. (11) one obtains the relation with α = R, A, <. The superscript should be understood according to the Langreth rule [50], i.e., (XY ) and so on. In Eq. (15), one can see that there are two types of couplings to the light field: one is the diamagnetic coupling (∝ ρA 2 ; ρ is the electron density) given through¨ k , and the other is the paramagnetic coupling (∝ j · A) given through two˙ k 's. The role of the latter paramagnetic coupling has been emphasized as a dominant interaction between the Higgs mode and electromagnetic fields [23].
The second derivative of the self-energyΣ n (ω) is determined by doubly differentiating the nonequilibrium extension of Eq. (12), where τ α 1 = τ 1 for α = R, A and τ α 1 = 0 for α =<. The first term on the right hand side of Eq. (16) represents the effect of amplitude fluctuation of the superconducting gap (i.e., the Higgs mode), while the second term corresponds to the impurity-ladder vertex corrections. Finally, the second derivative of the superconducting gap∆ n is given by doubly differentiating Eq. (13), The second derivatives,G kn (ω),Σ n (ω), and∆ n , are calculated by solving the self-consistent equations (15), (16), and (17). The results are plugged into the third derivative of the current ... j to obtain the THG susceptibility. More details are described in Appendix A.

B. Classification of THG susceptibilities
In Table II, we list all the THG diagrams including both the quasiparticle and Higgs-mode contributions. The distinction between the two contributions is defined by whether or not the diagram includes the Higgs-mode propagator depicted by the double wavy lines. The Higgs-mode propagator contains the fluctuation of the superconducting gap amplitude, which is diagrammatically shown in Fig. 4. Here the dotted lines represent the bare attractive interaction V mn , while the bold lines with arrows represent the electron propagatorĜ kn . The shaded square represents the impurity-ladder correction, as shown in Fig. 5.
There are five topologically inequivalent diagrams for quasiparticles and three diagrams for the Higgs mode. They have different couplings to external laser fields (single wavy lines) classified into the paramagnetic, diamagnetic, and mixed channels in Table II. Each outer vertex attached to photon lines in the THG diagram is assigned to the th derivative d k dA , which has the same parity as the density if is even and has the same parity as the current if is odd. Hence we call the channel of the coupling to light diamagnetic when is even and paramagnetic when is odd. There are THG diagrams in which the paramagnetic and diamagnetic couplings coexist, which we refer to as the mixed channel. We remark that the impurity correction is absent for vertices with odd number of photon lines since odd parity terms vanish after momentum summation.
In Table II, we also show which THG susceptibility has a resonance at frequency 2∆. As we will see in Sec. IV, χ H originates from the collective Higgs mode whose energy gap corresponds to 2∆. On the other hand, the quasiparticle contributions χ (i) qp also exhibit the resonance at 2∆, which is equal to the lowest pair-breaking energy. The degeneracy of the resonance energy between the Higgs mode and quasiparticles forces us to distinguish them by properties other than the resonance frequency, as discussed in Sec. I.
We also indicate in Table II which THG susceptibility is robust (insensitive) against nonmagnetic impurity scattering. Generally the THG susceptibility in the diamagnetic coupling channels do not depend on the impurity scattering rate, while the paramagnetic and mixed channels exhibit strong impurity dependence, which plays a key role in enhancing the Higgs-mode contribution in THG in dirty regimes. This fact is related to Anderson's theorem [51] (which states robustness of the supercon-TABLE II. Classification of the THG susceptibility according to the physical origin (quasiparticle or Higgs mode), the coupling channel to light (diamagnetic or paramagnetic), and the diagrammatic representation in the presence of impurities. In the diagrams, the lines with arrows, single wavy lines, double wavy lines, and shaded squares represent the electron propagators, photon propagators, Higgs-mode propagators, and impurity ladder corrections, respectively. In each diagram, three of the four photon propagators carry incoming frequencies Ω, while the rest carries outgoing frequency 3Ω. The fifth column shows whether each THG susceptibility has a resonance at frequency 2Ω = 2∆. The sixth column shows whether each THG susceptibility is robust against nonmagnetic impurities.

Higgs mode paramagnetic
ducting gap against nonmagnetic impurity scattering in equilibrium s-wave superconductors), which will be discussed in a separated publication [52]. We confirm the robustness of each THG channel against impurity scattering by numerical simulations in Sec. IV.

IV. THIRD HARMONIC GENERATION IN NbN SUPERCONDUCTOR
Based on the method described in Sec. III, we numerically evaluate the THG susceptibilities for NbN superconductors. We use the simplified band dispersion of NbN (2) derived from the first principles calculation in Sec. II. Throughout this section, we use eV as the unit of energy, and fix the ratio between the intraband and interband phonon-mediated interactions to be V inter /V intra = 0.18 (Sec. II). The absolute value of the interaction is chosen such that the superconducting gap is fixed. For the numerical feasibility, we take a relatively large superconducting gap 2∆ = 0.8 to maintain sufficiently high frequency resolution. We have checked that the results do not change qualitatively as we vary the value of 2∆. The intraband and interband impurity scattering rates γ intra and γ inter are free parameters. In order for the selfconsistent Born approximation to be valid (which is the case in the experimental situation [16,18]), the impurity scattering rates should be sufficiently smaller than the Fermi energy (γ intra , γ inter F ). Here we restrict ourselves to γ intra /2∆, γ inter /2∆ ≤ 2.5 (cf. F ∼ 3 − 4). We first focus on the case of γ intra = 2γ inter = γ with 0 ≤ γ/2∆ ≤ 2.5. Then we scan the parameter space of (γ intra , γ inter ) with 0 ≤ γ intra /2∆, γ inter /2∆ ≤ 2.5. We set the filling to be one third for NbN superconductor (two electrons in the three bands) and the inverse temperature β = 50, which is sufficiently lower than the superconducting critical temperature. In the numerical simulation, we take a finite value of the constant δ = 0.01 [which has been introduced as a positive infinitesimal in Eq. (11)].

A. Channel-resolved THG intensity
We first present the results for the THG intensity in NbN superconductors in the channel resolved manner as classified in Table II in Sec. III. In this subsection, the polarization direction e of light is set to be parallel to x crystal axis.
In Fig. 6, we plot the frequency dependence of the THG intensity for NbN superconductors with γ intra /2∆ = 2.0 and γ inter /2∆ = 1.0 (the dirty regime) in the linear (top panel) and log scale (middle). The leading contribution comes from the Higgs mode in the paramagnetic channel (χ (3) H ), showing a clear resonance peak at 2Ω = 2∆. The second dominant contribution comes from quasiparticles in the paramagnetic channel (χ (5) qp ), which has a relatively broadened resonance peak at 2Ω = 2∆. One can see that the resonance at 2Ω = 2∆ occurs in the channels χ  Table II. In the bottom panel of Fig. 6, we plot the total quasiparticle (|χ qp contributions to the THG intensity in the dirty regime of NbN superconductors. Clearly, the Higgs-mode contribution is larger than the quasiparticles. This result agrees with the previous observations that the Higgs-mode contribution is drastically enhanced due to impurity scattering [25][26][27]. Similar enhancement has been found due to phonon retardation effects [23].
In Fig. 7, we plot the impurity dependence of the THG intensity for NbN superconductors at frequency 2Ω = 2∆, where we set γ intra = 2γ inter = γ. The THG intensity in the paramagnetic and mixed channels [|χ   Table II in Sec. III. The γ dependence of the quasiparticle and Higgs-mode contributions in the diamagnetic and paramagnetic channels is qualitatively consistent with the previous order estimate in [26]. In the clean limit (γ → 0), the most dominant contribution comes from quasiparticles in the diamagnetic channel (χ The paramagnetic channel is also present at γ = 0, since we broaden the THG spectrum by taking the finite value of δ (so that the results at γ = 0 slightly deviate from the ideal clean limit). As we increase γ, the quasiparticle contribution in the paramagnetic channel (χ (5) qp ) quickly grows, and exceeds over the other components. At the same time, the Higgs-mode contribution in the paramagnetic channel (χ H ) also grows rapidly. Up to γ/2∆ 1, χ qp remains to be most dominant. When the system enters the dirty regime (γ/2∆ 1), the Higgs mode takes over the dominant part of the THG resonance, and χ becomes the largest contribution. This tendency seems to continue toward the dirty limit. The maximum magnitude of the THG intensity in the paramagnetic channel is attained around γ ∼ ∆, in agreement with the previous results [26].
In Fig. 8, we plot the ratio between the total Higgsmode (|χ H | 2 ) and quasiparticle (|χ qp | 2 ) contributions to the THG intensity for NbN superconductors as a function of γ intra = 2γ inter = γ. In the clean regime (γ 2∆) the quasiparticle contribution is dominant, whereas in the dirty regime (γ 2∆) the Higgs-mode contribution exceeds the quasiparticle one. We expect that the ratio |χ H | 2 /|χ qp | 2 continues to increase towards the dirty limit (as observed in the single-band case [27]), while our calculation is limited to γ/2∆ ≤ 2.5 in order to maintain the validity of the self-consistent Born approximation. One can see a little increase of |χ H | 2 /|χ qp | 2 around γ/2∆ = 0, which we attribute to the effect of the relatively large superconducting gap (2∆ = 0.8) and the finite broadening factor (δ = 0.01) in our simulation.
We separately change the intraband and interband impurity scattering rates to plot |χ H | 2 /|χ qp | 2 in Fig. 9. The quasiparticle dominant region is shown by blue in the color plot, while the Higgs dominant regions is shown by red. The boundary between the two regions is roughly given by γ intra /2∆ ∼ 1.75 and γ inter /2∆ ∼ 1. The interband scattering is more effective to enhance the Higgsmode contribution than the intraband one, since the interband scattering takes place more frequently in threeband systems. Each quasiparticle and Higgs-mode contribution is normalized by the value at θ = 0 • , respectively.

B. Polarization-angle dependence
Next, we study the polarization-angle dependence of the THG intensity for NbN superconductors, which may allow one to distinguish the quasiparticle and Higgs-mode contributions in experiments. The polarization angle is measured from the x crystal axis, and the polarization vector e is rotated in the xy plane.
In Fig. 10, we plot the normalized |χ qp (θ)| 2 /|χ qp (θ = 0 • )| 2 and |χ H (θ)| 2 /|χ H (θ = 0 • )| 2 for several values of γ intra /2∆ and γ inter /2∆. In the clean limit (top panel), we find that the quasiparticle contribution grows monotonically by ∼ 6% as the angle varies from θ = 0 • to 45 • , whereas the Higgs-mode contribution decreases by ∼ 7%. The angle dependence of quasiparticles arises due to the anisotropic band structure of NbN. The change of the quasiparticle contribution from θ = 0 • to 45 • is smaller than that of the previous result [18]. This is mainly due to the difference of the hopping parameters that we used in the effective model of NbN. The angle dependence of the Higgs mode in the clean limit with small V inter /V intra is consistent with the previous study [28].
As we increase γ intra = 2γ inter = γ (the middle and bottom panels in Fig. 10), we observe qualitatively different polarization-angle dependence for quasiparticles. Namely, |χ qp | 2 decreases by ∼ 20 − 30% as θ changes from 0 • to 45 • . This behavior is mostly coming from the paramagnetic channel, which becomes dominant in the dirty regime. Namely, the quasiparticle contribution in the paramagnetic channel always tends to decrease from θ = 0 • to 45 • for arbitrary impurity scattering rates. The transition from the increasing to decreasing dependence on the polarization angle is very rapid, taking place around γ/2∆ ∼ 0.1 where the paramagnetic channel starts to exceed the diamagnetic one. Contrary to the significant angle dependence for quasiparticles, the Higgs-mode contribution quickly becomes isotropic as one deviates from the clean limit. The angle dependence of the Higgs-mode contribution is no larger than 1.5% at γ/2∆ ≥ 1.
In Fig. 11, we plot |χ qp (θ = 45 • )| 2 /|χ qp (θ = 0 • )| 2 (left panel) and |χ H (θ = 45 • )| 2 /|χ H (θ = 0 • )| 2 (right) in the space of γ intra /2∆ and γ inter /2∆. The quasiparticle contribution generally shows clear angle dependence for arbitrary impurity scattering rates. The increasing behavior of |χ qp | 2 as a function of θ is seen only in the vicinity of the clean limit, apart from which |χ qp | 2 decreases by 5 − 50%. It seems that the angle dependence of the quasiparticle contribution does not vanish in the large intraband and/or interband impurity scattering. On the other hand, the angle dependence of the Higgs mode is suppressed for general impurity scattering rates as compared to quasiparticles. At the vanishing of the interband impurity scattering, the Higgs-mode contribution shows slight angle dependence of few %. One can also see that the angle dependence of both the quasiparticle and Higgs-mode contributions is sensitive to the interband impurity scattering rather than the intraband one in the dirty regime.
We expect that there are generally nonvanishing interband impurity scatterings in NbN. Then, one can use the polarization-angle dependence of THG to discrimi- nate the Higgs-mode and quasiparticle contributions in the dirty regime of multiband superconductors. The optical conductivity measurement [16,18] suggests that the NbN samples used in the THG experiment is close to the dirty limit (γ/2∆ 1). The experimental observation of no angle dependence of THG in NbN superconductors [18] together with our results on the channel-resolved THG intensity in the dirty regime imply that the dominant contribution to the THG resonance originates from the Higgs mode.
Finally, let us comment on the behavior of the polarization-angle dependence on the ratio V inter /V intra . While we used the realistic value of V inter /V intra = 0.18 for NbN estimated from first principles calculations throughout the paper, we have checked the angle dependence for several other values of V inter /V intra (not shown). In general, the angle dependence of the Higgs mode tends to be strongly suppressed as one increases V inter /V intra (for the case of V inter /V intra = 1 in the clean limit, see [18]), whereas the angle dependence of quasiparticles remains almost unchanged. Although the realistic value of V inter /V intra = 0.18 that we obtained in the present paper is not so large, we find that the Higgs-mode contribution is almost polarization-angle independent in the dirty regime for NbN superconductors. This suggests that the effect of impurity scattering plays an important role in understanding the behavior of the THG resonance in NbN superconductors.

V. SUMMARY
To summarize, we study the resonance of third harmonic generation and its polarization-angle dependence in disordered NbN superconductors based on the effec-tive three-band model constructed from first principles calculations on the electron and phonon band structures of NbN. Using the density functional perturbation theory, we evaluate the band-resolved matrix elements of the electron-phonon coupling constants for NbN, and the ratio between the intraband and interband pairing interactions, V inter /V intra , is found to be about 0.17-0.18.
We input the evaluated pairing interaction parameters in the effective model, whose THG susceptibility is calculated in the channel-resolved manner with the BCS mean-field and self-consistent Born approximations. The results show that in the dirty regime the dominant contribution to the THG resonance is given by the Higgs mode in the paramagnetic channel, which does not have polarization-angle dependence with nonvanishing interband impurity scattering. The second dominant one is given by quasiparticles in the paramagnetic channel, which exhibit clear polarization-angle dependence in the dirty regime. Our results are quite consistent with the polarization-resolved THG experiment on NbN superconductors [18], which have found no polarization-angle dependence in the THG resonance.
Our scheme of classifying and calculating THG susceptibilities for disordered superconductors from first principles can be applied not only to NbN but also to other superconductors. Interesting future applications include the THG resonance in MgB 2 , a multigap superconductor with multiple Higgs modes as well as the Leggett mode [53][54][55][56], and NbSe 2 , where superconductivity and charge density wave coexist. For unconventional superconductors such as cuprates [57][58][59][60][61] and iron-based superconductors, we need to extend the present formalism to take into account strong correlation effects beyond the BCS approximation in THG, which we leave as a future problem.
which are again solved numerically by matrix inversion.
Finally, the THG susceptibility is determined from the vertex functions. The classification of the THG susceptibility in terms of the diagrammatic topology has been given in Table II (Ω)τ 1Ĝkn (ω + Ω)˙ knĜkn (ω)] < .