Quasiclassical Theory on Third-Harmonic Generation in Conventional Superconductors with Paramagnetic Impurities
Takanobu Jujo
cond-mat.supr-con



-23.5pt
Physical Properties and Evolutionary States of EA-type Eclipsing Binaries Observed by LAMOST
    Qian S.-B.1,2,3,4, Zhang J.1,2,3, He J.-J.1,2,3,
Zhu L.-Y.1,2,3,4, Zhao E.-G.1,2,3,4, Shi X.-D.1,2,3, Zhou X.1,2,3 and Han Z.-T.1,2,3,4
    December 30, 2023
===========================================================================================================================================





§ INTRODUCTION


In recent years, studies on the nonlinear optical 
response in superconductors have advanced 
with the development of microwave spectroscopy. 
Many of these studies were on the dynamics 
of transient responses using pump-probe spectroscopy, 
and the main objective was to elucidate the interaction 
effect in superconductors including 
strongly correlated electron systems. 
(For a review, see Ref. 1

for example.) 

Another important aspect of the nonlinear optical response 
in the research of superconductors is that we can 
investigate phenomena that do not appear in the linear response. 
The amplitude fluctuation of the superconducting 
order parameter <cit.> 
is a typical example because the amplitude mode 
is not reflected in the linear absorption spectrum. 
(This is different from the phenomenon that 
the phase mode is pushed up to a high energy 
by the long-range Coulomb interaction. <cit.>) 

So far, the amplitude mode has been observed in the phonon spectrum 
through the coupling of the superconductivity and 
the charge density wave in a superconductor 
coexisting with the charge density wave state. <cit.> 
Recently, it has been reported that 
the amplitude mode can be observed in ordinary superconductors 
by pump-probe spectroscopy <cit.> and 
third-harmonic generation (THG). <cit.> 
The former experimental study asserted that the observed 
transient oscillation 
of optical conductivity is an expression of the amplitude mode. 
The latter showed that the peak appears at 
a temperature at which the superconducting gap is 
equal to the frequency of the incident wave. 
This may indicate the existence of the amplitude mode. <cit.> 

Theoretically, the excitation of the amplitude mode 
overlaps the quasiparticle excitation 
by two-photon absorption in the spectrum. <cit.> 
Thus, it is required to distinguish these two phenomena 
in order to confirm the amplitude mode. 
As a proposal for identifying the amplitude mode, 
in this paper we investigate the THG intensity 
in superconductors with paramagnetic impurities. 
Time-reversal symmetry breaking due to paramagnetic impurities 
leads to instability of the amplitude mode. <cit.> 
Therefore, the peak of the THG intensity will decrease 
when the amplitude mode is dominant in the THG intensity. 
In addition, in order to clarify the existence of 
the amplitude mode, 
we conduct a semiquantitative evaluation, 
which has not been done so far 
(for example, supplementary materials of 
Refs. 8 and 12 

and references therein), 
and compare it with the experimental result. 


The following is the structure of this paper. 
Section 2 gives the formulation by the quasi-classical 
approximation for calculating the response function. 
Section 3 gives an expression for the THG intensity 
under pulse irradiation and shows 
the results of a numerical calculation such as 
its temperature dependence. 
A semiquantitative evaluation is also given. 
Section 4 gives an expression for the diamagnetic term 
and shows that this term is negligible in 
the THG intensity. 






§ FORMULATION


A nonlinear current is written using the 
quasiclassical Green function (g^K) as follows: <cit.> 

    J^μ_ω=e∫_ FSv^μ_k
    mk_F/2π∫dϵ/4π i∑_ν=x,y,z Tr[v_k^νĝ^K(3)_ν(ϵ+ω,ϵ)].

Here, v_k and m are 
the velocity and mass of electrons, respectively, and 
k_F is the Fermi wave number. (We put ħ=c=1 in this paper with 
c the velocity of light.) 
∫_ FS means integration on the Fermi surface. 
The superscripts K and (3) in ĝ^K(3) indicate 
the Keldysh Green function <cit.> 
and the third order of the external fields, 
respectively. The modifier “  ”
means a matrix in Nambu representation. 
μ,ν=x,y,z indicate spatial directions. 
We assume the isotropic case and omit this index 
hereafter. 

ĝ^K(3)(ϵ+ω,ϵ) 
is derived from kinetic equations in the dirty limit. <cit.> 
In this limit 
the quasiclassical Green function can be divided into 
an odd order and an even order with respect to the external field. 
In the kinetic equation for the former, 
the effect of scattering by nonmagnetic 
impurities does not vanish in the collision integral and 
is predominant over other terms such 
as the superconducting gap. 
Thus, the kinetic equation in this case 
can be solved, and its solution 
is written as 

    ĝ^K(3)(ϵ+ω,ϵ)
    =e∫dω'/2π
    A_ω'-1/2α(
    ĝ^+_ϵ+ωĝ^K_ϵ+ω-ω',ϵ
    +
    ĝ^K_ϵ+ωĝ^-_ϵ+ω-ω',ϵ
    +
    ĝ^+_ϵ,ϵ-ω+ω'ĝ^K_ϵ-ω
    +
    ĝ^K_ϵ,ϵ-ω+ω'ĝ^-_ϵ-ω
       -
    τ̂_3ĝ^+_ϵτ̂_3ĝ^K_ϵ,ϵ-ω+ω'
    -
    ĝ^K_ϵτ̂_3
    ĝ^-_ϵ,ϵ-ω+ω'τ̂_3
    -
    τ̂_3ĝ^+_ϵ,ϵ-ω+ω'τ̂_3ĝ^K_ϵ-ω+ω' 
    -
    ĝ^K_ϵ,ϵ-ω+ω'τ̂_3
    ĝ^-_ϵ-ω+ω'τ̂_3
    ).

Here, α=(mk_F/2π)n_i u_i^2 with 
n_i the concentration of nonmagnetic impurities and 
u_i the magnitude of the potential. 
The dirty limit means that α≫Δ. 
A_ω is the external vector potential. 
τ̂_3=(
[  1  0;  0 -1 ]). 
ĝ_ϵ is the quasiclassical Green function 
in the equilibrium state and 
ĝ_ϵ,ϵ' is a second-order function on external 
fields (+ and - in the superscript indicate 
the retarded and advanced Green function, respectively). 
The kinetic equations for the latter are written as 

    τ̂_3ϵĝ^+_ϵ,ϵ'
    -ĝ^+_ϵ,ϵ'ϵ'τ̂_3
    -
    [
    τ̂_3
    Σ̂^+_ϵĝ^+_ϵ,ϵ'
    -ĝ^+_ϵ,ϵ'Σ̂^+_ϵ'τ̂_3
    +
    τ̂_3
    Σ̂^+_ϵ,ϵ'ĝ^+_ϵ'
    -ĝ^+_ϵΣ̂^+_ϵ,ϵ'τ̂_3
    ] 
    
       -e^2D_α∫dω_1dω_2/(2π)^2
    A_ω_1A_ω_2[
    τ̂_3ĝ^+_ϵ-ω_1ĝ^+_ϵ'
    -ĝ^+_ϵĝ^+_ϵ'+ω_1τ̂_3
    ]δ(ϵ-ϵ'-ω_1-ω_2)
    =0

and 

    τ̂_3ϵĝ^(a)_ϵ,ϵ'
    -ĝ^(a)_ϵ,ϵ'ϵ'τ̂_3
    -
    [
    τ̂_3
    Σ̂^+_ϵĝ^(a)_ϵ,ϵ'
    -ĝ^(a)_ϵ,ϵ'Σ̂^-_ϵ'τ̂_3
    +
    τ̂_3
    Σ̂^(a)_ϵ,ϵ'ĝ^-_ϵ'
    -ĝ^+_ϵΣ̂^(a)_ϵ,ϵ'τ̂_3
    ] 
       -e^2D_α∫dω_1dω_2/(2π)^2
    A_ω_1A_ω_2{τ̂_3
    [
    (t^h_ϵ-ω_1-t^h_ϵ-ω_1-ω_2)
    ĝ^+_ϵ-ω_1
    +
    (t^h_ϵ-t^h_ϵ-ω_1)
    ĝ^-_ϵ-ω_1]
    ĝ^-_ϵ'
       -
    ĝ^+_ϵ[
    (t^h_ϵ'+ω_1-t^h_ϵ')
    ĝ^+_ϵ'+ω_1
    +
    (t^h_ϵ'+ω_1-t^h_ϵ)
    ĝ^-_ϵ'+ω_1]
    τ̂_3
    }δ(ϵ-ϵ'-ω_1-ω_2)
    =0.

Here, t^h_ϵ:= tanh(ϵ/2T) (T is temperature), 
δ(·) is the delta function, and 
ĝ^(a)_ϵ,ϵ':=
ĝ^K_ϵ,ϵ'
-t^h_ϵ'ĝ^+_ϵ,ϵ'
+t^h_ϵĝ^-_ϵ,ϵ' 
(ĝ^(a)=Σ̂^(a)=0̂ 
in the equilibrium state). <cit.> 
The effect of impurity scatterings is calculated 
with the Born approximation, <cit.> 
and the interaction between electrons and phonons 
is treated with the weak-coupling approximation. 
The self-energy is written as 
Σ̂^s_ϵ,ϵ'=
Σ̂^(ep)s_ϵ,ϵ'
+Σ̂^(ni)s_ϵ,ϵ'
+Σ̂^(pi)s_ϵ,ϵ' 
with 
Σ̂^(ni)s_ϵ,ϵ'=
ατ̂_3ĝ^s_ϵ,ϵ'τ̂_3 
(the effect of nonmagnetic impurity scattering) and 
Σ̂^(pi)s_ϵ,ϵ'=
α_p ĝ^s_ϵ,ϵ'
(the effect of paramagnetic impurity scattering <cit.>) 
[s=+, - or (a)]. 
Σ̂^(ep)+_ϵ,ϵ'
=Σ̂^(ep)-_ϵ,ϵ'
=g_0∫dϵ_1/2π iτ̂_3ĝ^K_ϵ_1+(ϵ-ϵ')/2,
ϵ_1-(ϵ-ϵ')/2τ̂_3 
(the electron-phonon interaction) 
and 
Σ̂^(ep)(a)_ϵ,ϵ'=
(t^h_ϵ-t^h_ϵ')
Σ̂^(ep)+_ϵ,ϵ'. 
α_p=(mk_F/2π)n'_i u'^2 
and g_0=(mk_F/2π)(g_ph^2/ω_D) 
with ω_D the Debye frequency and 
g_ph the coupling constant between 
electrons and phonons. 
D_α=v_F^2/6α=v_F^2τ/3 
(τ=1/2α is the relaxation time) is 
the diffusion constant. 

Equations (<ref>) and (<ref>) 
are solved by introducing 

    ĝ^s_ϵ,ϵ'=
    g^s_ϵ,ϵ'τ̂_0
    +f^s_ϵ,ϵ'τ̂_1

and 

    Σ̂^s_ϵ,ϵ'=
    Σ^n,s_ϵ,ϵ'τ̂_0
    +Σ^a,s_ϵ,ϵ'τ̂_1.

Here, 
n and a indicate the normal and anomalous 
self-energy, respectively. 
τ̂_0=(
[ 1 0; 0 1 ]) and 
τ̂_1=(
[ 0 1; 1 0 ]). 

By solving the above equations using 
Σ^(ep)+_ϵ=Σ^(ep)-_ϵ
=Δτ̂_1 
(Δ is the superconducting gap 
with the effect of impurity scattering included), 
the solution is written as 

    [ g^s_ϵ,ϵ'; f^s_ϵ,ϵ' ]
    :=
    e^2D_α∫dω_1dω_2/(2π)^2
    A_ω_1A_ω_2δ(ϵ-ϵ'-ω_1-ω_2)
    [ g̅^s_ϵ,ϵ'(ω_1,ω_2); f̅^s_ϵ,ϵ'(ω_1,ω_2) ].

Here, s=+, -, or (a), and 

    [ g̅^±_ϵ,ϵ'(ω_1,ω_2); f̅^±_ϵ,ϵ'(ω_1,ω_2) ]
    =
    M̂^±±_ϵ,ϵ'1/2[
    [ g^±_ϵ-ω_1; f^±_ϵ-ω_1 ]
    +
    [ g^±_ϵ-ω_2; f^±_ϵ-ω_2 ]]
    +
    M̂^±±_ϵ,ϵ'[                    0; Σ̅^a,+_ϵ,ϵ'(ω_1,ω_2) ]

(the double signs correspond), 
and 

    [ g̅^(a)_ϵ,ϵ'(ω_1,ω_2); f̅^(a)_ϵ,ϵ'(ω_1,ω_2) ]
    =   M̂^+-_ϵ,ϵ'1/2∑_ω'=ω_1,ω_2[
    (t^h_ϵ-ω'-t^h_ϵ')
    [ g^+_ϵ-ω'; f^+_ϵ-ω' ]
    +
    (t^h_ϵ-t^h_ϵ-ω')
    [ g^-_ϵ-ω'; f^-_ϵ-ω' ]]  
       +
    M̂^+-_ϵ,ϵ'
    (t^h_ϵ-t^h_ϵ')
    [                    0; Σ̅^a,+_ϵ,ϵ'(ω_1,ω_2) ].

Here, Σ̅ indicates the vertex correction, 
which is calculated in the next subsection, 
and the diagonal element (Σ̅^n,+) 
is shown to vanish. 

    M̂^ab_ϵ,ϵ'
    :=
    i[
    τ̂_3
    -X^ab_ϵ,ϵ'τ̂_0
    -Y^ab_ϵ,ϵ'τ̂_1
    ] 
    /z_ϵ^a +z_ϵ'^b
    +2iα_p X^ab_ϵ,ϵ'.

Here, 
X^ab_ϵ,ϵ':= 
(ϵ^a_pϵ'^b_p
+Δ^a_ϵΔ^b_ϵ')/
z^a_ϵz^b_ϵ', 
Y^ab_ϵ,ϵ':= 
(ϵ^a_p Δ_ϵ'^b
+Δ^a_ϵϵ^'b_p)/
z^a_ϵz^b_ϵ', 
ϵ^±_p:=ϵ-α_p g^±_ϵ, 
Δ^±_ϵ:=Δ+α_p f^±_ϵ, 
and 
z^±_ϵ:=√((ϵ_p^±)^2
-(Δ^±_ϵ)^2). 
Quasiclassical Green functions in the equilibrium state
are written as 
g_ϵ^±=-iϵ_p^±/z^±_ϵ
and
f_ϵ^±=-iΔ_ϵ^±/z^±_ϵ. 




By using the above results, the nonlinear current 
is written as 

    J_ω=   -(e^2D_α)^2
    mk_F/2π∫dw dw'/(2π)^3
    A_ω-2wA_w+w'/2A_w-w'/2
    I^(3)_ω,w,w'.

Here, 

    I^(3)_ω,w,w':=   ∫dϵ/4π i Tr[ 
    (
    ĝ^+_ϵ+ω-w
    -τ̂_3ĝ^+_ϵ+wτ̂_3
    +ĝ^-_ϵ-ω+w
    -τ̂_3ĝ^-_ϵ-wτ̂_3)
    ĝ̅̂^K_ϵ+w,ϵ-w(w+w'/2,w-w'/2)
       +
    (
    ĝ^K_ϵ+ω-w
    -τ̂_3ĝ^K_ϵ+wτ̂_3
    )
    ĝ̅̂^-_ϵ+w,ϵ-w(w+w'/2,w-w'/2)
    +
    (
    ĝ^K_ϵ-ω+w
    -τ̂_3ĝ^K_ϵ-wτ̂_3
    )
    ĝ̅̂^+_ϵ+w,ϵ-w(w+w'/2,w-w'/2)].





 §.§ Vertex correction


We obtain Σ^(ep)+_ϵ,ϵ' 
by solving 

    [ Σ^(ep)n,+_ϵ,ϵ'; Σ^(ep)a,+_ϵ,ϵ' ]
    =
    g_0∫dϵ_1/2π iτ̂_3
    [
    [ g^(a)_ϵ_1+w,ϵ_1-w; f^(a)_ϵ_1+w,ϵ_1-w ]
    +
    t^h_ϵ_1-w[ g^+_ϵ_1+w,ϵ_1-w; f^+_ϵ_1+w,ϵ_1-w ]
    -
    t^h_ϵ_1+w[ g^-_ϵ_1+w,ϵ_1-w; f^-_ϵ_1+w,ϵ_1-w ]]

[w:=(ϵ-ϵ')/2]. 
Using Eqs. (<ref>)-(<ref>) and 

    [ Σ^(ep)n,+_ϵ,ϵ'; Σ^(ep)a,+_ϵ,ϵ' ]
    :=e^2D_α∫dω_1dω_2/(2π)^2
    A_ω_1A_ω_2δ(ϵ-ϵ'-ω_1-ω_2)
    [ Σ̅^n,+_ϵ,ϵ'(ω_1,ω_2); Σ̅^a,+_ϵ,ϵ'(ω_1,ω_2) ],

the solution is written as 

    [ Σ̅^n,+_ϵ,ϵ'(ω_1,ω_2); Σ̅^a,+_ϵ,ϵ'(ω_1,ω_2) ]
    =
    g_0/2D_2w∫dϵ_1/2π i[  0  0;  0 -1 ]∑_s,s'=±
    s[    
    (t^h_ϵ_1-s'w'-t^h_ϵ_1-sw)
    M̂^+-_ϵ_1+w,ϵ_1-w
       +
    t^h_ϵ_1-swM̂^ss_ϵ_1+w,ϵ_1-w]
    [ g^s_ϵ_1-s'w'; f^s_ϵ_1-s'w' ]

with w=ϵ-ϵ'=(ω_1+ω_2)/2 
and w'=(ω_1-ω_2)/2. 
Here, 

    D_2w=
    1-g_0∫dϵ_1/2π i[
    (t^h_ϵ_1+w-t^h_ϵ_1-w)
    m^+-_ϵ_1+w,ϵ_1-w
    +
    t^h_ϵ_1-w
    m^++_ϵ_1+w,ϵ_1-w
    -
    t^h_ϵ_1+w
    m^–_ϵ_1+w,ϵ_1-w]

with 

    m^ab_ϵ,ϵ':=
    i(1+X^ab_ϵ,ϵ')
    /z_ϵ^a +z_ϵ'^b
    +2iα_pX^ab_ϵ,ϵ'.





§ THIRD-HARMONIC GENERATION




 §.§ Nonlinear response under incident pulse


The electric field is written as 
Ẽ_̃t̃=E̅_̅t̅e^-iΩ t
+E̅_̅t̅^*e^iΩ t. 
Using of this field, 
the vector potential is written as 
A_ω=E_ω/iω 
with 
E_ω=∫dt/2πẼ_̃t̃e^iω t. 
We assume a Gaussian pulse for E̅_̅t̅: 
|E̅_̅t̅|=E̅_̅0̅e^-(t/t_0)^2. 
Then

    A_ω=√(π)t_0/iω[
    E̅_̅0̅e^-(ω-Ω)^2/4c+
    E̅_̅0̅^*e^-(ω+Ω)^2/4c],

and we introduce a dimensionless external field, 

    A̅_ω:=
    Δ_0/√(π)t_0|E̅_̅0̅|A_ω.

(Δ_0 is the superconducting gap 
at T=0 without impurity scatterings.) 
Using of A̅_ω, 
the current [Eq. (<ref>)] is 
rewritten as 

    J_ω=σ_0 |E_Ω|
    (e|E̅_̅0̅|ξ_0/Δ_0)^2
    l/ξ_0j_ω

with 

    j_ω:=-1/48∫ dw dw't_0^2
    A̅_ω-2wA̅_w+w'/2A̅_w-w'/2
    I^(3)_ω,w,w'.

Here,
|E_Ω|=√(π)t_0|E̅_̅0̅+E̅_̅0̅^*
e^-(Ω t_0)^2|≃√(π)t_0|E̅_̅0̅| and 
σ_0=n_e e^2 τ/m=e^2 D_α mk_F/π^2 
are used. 
ξ_0=v_F/πΔ_0 and 
l=v_Fτ 
are the coherence length and 
the mean free path, respectively. 

Equation (<ref>) does not 
depend on α. 
This is because 
ĝ^s_ϵ 
and ĝ̅̂^s_ϵ,ϵ' 
do not include α owing to 
Anderson's theorem <cit.> and 
the absence of a collision term derived from nonmagnetic impurity 
scattering <cit.>, respectively. 
The dependence of the nonlinear current 
[Eq. (<ref>)] on nonmagnetic impurity scattering 
is included only in l. 
Thus, the following numerical calculations are performed 
without specifying the value of α. 






 §.§ Numerical calculations


We calculate the THG intensity using 
the dimensionless j_ω from Eq. (<ref>). 
We take Δ_0 as the unit of energy (Δ_0=1), 
and in the variables in the subsequent figures, 
the notation of Δ_0 
is omitted for the sake of simplicity. 

The dependence of the THG intensity (|j_3Ω|^2 
with Ω the frequency of the incident pulse) 
on temperature is shown in Fig. <ref>. 

There is a sharp peak for small values of α_p. 
This peak becomes blurred and shifts to low temperatures 
with increasing α_p. 
As shown below, this behavior is caused by the fact that 
this peak originates from the amplitude mode, 
and the blurring occurs owing to 
the vanishing of the well-defined mode for 
finite values of α_p. 

The amplitude mode described by D_2ω 
[Eq. (<ref>)] 
is included in the vertex correction 
Σ̅^a [Eq. (<ref>)]. 
The THG intensities with and 
without the vertex correction are shown 
in Figs. <ref>(a) and 2(b), respectively. 
 
Figure 2(a) [2(b)] is calculated using g̅_ω, 
in which only the second [first] terms 
in the left-hand side of Eqs. (<ref>) and (<ref>) 
are included. 
The THG intensity is proportional to the square of 
the absolute value of j_3Ω, and thus, 
the result in Fig. 1 is not a simple summation of 
the results in Figs. 2(a) and 2(b). 
Figure 2(a) indicates that 
the vertex correction causes a sharp peak 
for small values of α_p. 

The dependence of the amplitude mode (D_2ω) 
on the frequency ω at T/T_C=0.5 is shown in Fig. <ref>.

Figure 3(a) shows that a well-defined amplitude mode 
disappears with increasing α_p. 
- ImD_2ω is the damping rate of 
the amplitude mode and takes finite values 
in the range ω<Δ for finite values of α_p. 
This behavior originates from the fact that 
the gap edge (E_g) of a one-particle spectrum 
is smaller than Δ for finite values of α_p: 
E_g=Δ[1-(2α_p/Δ)^2/3]^3/2 
(see Ref. 19). 


We introduce ω_am as the frequency 
at which ReD_2ω takes a local minimum. 
The dependence of ω_am on temperature 
is shown in Fig. <ref>. 

The values of ω_am are numerically calculated 
from the dependences of ReD_2ω on frequency. 
Figure 4 shows that 
ω_am≃Δ regardless of the value of α_p. 
In contrast, 
E_g deviates from Δ with increasing α_p as noted 
above. 
The reason for including horizontal lines in Fig. 4 is 
to investigate how the peak of the THG intensity is related 
to the temperature at which the frequency (Ω) crosses Δ 
or E_g. 

The temperature (normalized by T_C) 
at which Δ (or E_g) and Ω (=0.46Δ_0, 0.69Δ_0, 
0.92Δ_0) intersect is shown in Fig. <ref>. 

(T_Δ and T_E_g are dimensionless quantities.) 
T_Δ and T_E_g are set to 0 
when Δ and E_g do not 
cross Ω, respectively. 
Figures 1 and 5 show that 
the THG intensity does not show a peak 
at T/T_C=T_Δ or T/T_C=T_E_g. 
The temperature at which the THG intensity shows a peak 
is between T_Δ and T_E_g. 

This is verified by calculating the 
THG intensity for other values of Ω. 
The THG intensities for Ω=0.69Δ_0 
and Ω=0.92Δ_0 are shown 
in Figs. <ref>(a) and 6(b), respectively. 

Figures 6(a) and 6(b) show that the peak of the THG intensity 
almost vanishes in the case of T_E_g=0. 
The deviation of the peak from T_Δ 
is similar to the case of the two-photon absorption (TPA) spectrum. 
In the latter case, 
the frequency at which 
the TPA spectrum shows a peak 
is between E_g and Δ. <cit.> 
This behavior is caused by a shift of the spectrum 
to a lower energy with increasing α_p. 

The THG intensities for various values of 1/t_0^2 are shown 
in Fig. <ref>. 

A small value of 1/t_0^2 corresponds to 
a small width of the incident pulse in the frequency space. 
The temperature at which the THG intensity shows a peak 
does not vary with 1/t_0^2, 
but the peak becomes sharp with decreasing 1/t_0^2. 
This behavior supports an expectation about the peak 
in the THG intensity suggested in Ref. 8. 

 



 §.§ Quantitative evaluation


The THG intensity divided by the incident pulse 
is approximately written as 

    |
    4π d/2n_3Ω c|^2
    |J_3Ω|^2/|E_Ω|^2
    =
    |4π dσ_0/2n_3Ω c|^2
    (e|E̅_̅0̅|ξ_0/Δ_0)^4
    (l/ξ_0)^2
    |j_3Ω|^2

(see, for example, Ref. 24). 

Here, n_3Ω is the refractive index, 
|n_3Ω|^2
≃ 4π |σ_3Ω|/3Ω, and 
d is the thickness of the film. 
When we set ξ_0=5 nm, <cit.> Δ_0=2.7 meV, 
d=24 nm, and σ_0=1.5× 10^4 Ω^-1 cm^-1, 
and assume |σ_3Ω|≃σ_0 because we consider 
the case of 3Ω>Δ, we obtain 

    |4π dσ_0/2n_3Ω c|^2
    (e|E̅_̅0̅|ξ_0/Δ_0)^4
    (l/ξ_0)^2
    ≃
    1.37× 10^-5

for 
l=0.8 nm (l≃ 0.58 ∼ 0.83 nm <cit.>), 
|E̅_̅0̅|=3.5 kV/cm, and Ω=0.92Δ_0. 
Then 
|4π dσ_0/2n_3Ω c|^2
(e|E̅_̅0̅|ξ_0/Δ_0)^4
(l/ξ_0)^2
|j_3Ω|^2 ≃ 6.85× 10^-5 
using the calculated results of |j_3Ω|^2 
in Fig. 6(b) for Ω=0.92Δ_0. 
This is about the same value as an experimental value in 
Ref. 8, 

where 8× 10^-5 was reported as 
a value of the THG intensity normalized by that of 
the pump pulse. 
 










§ DIAMAGNETIC TERM


When we take account of the diamagnetic coupling 
in the interaction between electrons and 
external fields <cit.>, there is an additional term 
in a nonlinear current as follows: 

    J^d_ω=-e^2/m∫dω'/2π
    A_ω'∫_ FSmk_F/2π∫dϵ/4π i Tr[
    τ̂_3ĝ^K_ϵ+ω-ω',ϵ].

In this case, 
there are also additional terms in the kinetic equations. 
These are written as 

    ∫dω_1dω_2/(2π)^2e^2/2
    A_ω_1A_ω_21/m[
    ĝ^+_ϵ-ĝ^+_ϵ']δ(ϵ-ϵ'-ω_1-ω_2)

and 

    ∫dω_1dω_2/(2π)^2e^2/2
    A_ω_1A_ω_21/m(t^h_ϵ-t^h_ϵ')
    [
    ĝ^+_ϵ-ĝ^-_ϵ']δ(ϵ-ϵ'-ω_1-ω_2)

for the left-hand sides of 
Eqs. (<ref>) and (<ref>), 
respectively. 
The additional term arises in the solution 
[Eq. (<ref>)] 
and is written as <cit.>  

    g'_ϵ,ϵ'τ̂_3
    +if'_ϵ,ϵ'τ̂_2
 
with iτ̂_2=(
[  0  1; -1  0 ]). 

Since (g',f') are different matrix elements from 
(g,f) in Sects. 2 and 3, 
they can be obtained separately. 
The solution is written as 

    [ g^'s_ϵ,ϵ'; f^'s_ϵ,ϵ' ]
    =
    ∫dω_1dω_2/(2π)^2e^2/2
    A_ω_1A_ω_2δ(ϵ-ϵ'-ω_1-ω_2)
    1/m[ g̅^'s_ϵ,ϵ'; f̅^'s_ϵ,ϵ' ]

with 

    [ g̅^'±_ϵ,ϵ'; f̅^'±_ϵ,ϵ' ]
    :=
    M̂^'±±_ϵ,ϵ'[ 1; 0 ]

and 

    [ g̅^'(a)_ϵ,ϵ'; f̅^'(a)_ϵ,ϵ' ]
    :=
    (t^h_ϵ-t^h_ϵ')
    M̂^'+-_ϵ,ϵ'[ 1; 0 ].

Here, 

    M̂^'ab_ϵ,ϵ':=
    i[
    τ̂_3
    -X^'ab_ϵ,ϵ'τ̂_0
    -Y^'ab_ϵ,ϵ'τ̂_1
    ] 
    /z_ϵ^a +z_ϵ'^b
    +2iα_p X^'ab_ϵ,ϵ',

where 
X^'ab_ϵ,ϵ':= 
(ϵ^a_pϵ'^b_p
-Δ^a_ϵΔ^b_ϵ')
/z^a_ϵz^b_ϵ' 
and 
Y^'ab_ϵ,ϵ':= 
(ϵ^a_p Δ_ϵ'^b
-Δ^a_ϵϵ^'b_p)
/z^a_ϵz^b_ϵ'. 
Using the above quantities, 
Eq. (<ref>) is rewritten as 

    J^d_ω=
    -1/2(e^2/m)^2
    ∫dw dw'/(2π)^3
    A_ω-2wA_w+w'/2A_w-w'/2mk_F/2π∫dϵ/4π i Tr[τ̂_3
    ĝ̅̂^'K_ϵ+w,ϵ-w].


By using the electric field introduced in Sect. 3, 
the nonlinear current is written as 

    J^d_ω=
    σ_0|E_Ω|
    (e|E̅_̅0̅|ξ_0/Δ_0)^2
    (Δ_0/E_F)^2
    9π^2ξ_0/8l
    j^d_ω.

Here, E_F=k_F^2/2m is the Fermi energy and 

    j^d_ω=
    -1/48∫ dw dw' t_0^2
    A̅_ω-2wA̅_w+w'/2A̅_w-w'/2∫dϵ/4π i Tr[τ̂_3
    ĝ̅̂^'K_ϵ+w,ϵ-w].





We consider the ratio |J_3Ω^d|^2/|J_3Ω|^2
in order to evaluate this term quantitatively. 
From Eqs. (<ref>) and (<ref>), 

    |J_3Ω^d|^2/|J_3Ω|^2
    =(Δ_0/E_F)^4
    (3πξ_0/2√(2)l)^4
    |j_3Ω^d|^2/|j_3Ω|^2.

Here, 
3πξ_0/2√(2)l≃ 20.8 
for ξ_0=5 nm and l=0.8 nm, and 
then 
|J_3Ω^d|^2/|J_3Ω|^2≃ 2× 10^-3
when Δ_0/E_F=0.01 
and |j_3Ω^d|^2/|j_3Ω|^2=1 
(this last value is an overestimated value as mentioned below). 
Thus, the diamagnetic term is negligible 
as suggested in the annotation of Ref. 10 

(Ref. 14 therein). 
This result is consistent with the experimental result 
that the diamagnetic term is not observed. <cit.> 

The result of a numerical calculation shows that 
|j_3Ω^d|^2/|j_3Ω|^2 is smaller than 1. 
The reason for this is that 
there is no increase due to the amplitude mode in 
|j^d_3Ω|^2. 
The absence of the amplitude mode in |j^d_3Ω|^2 
is the result of the properties of 
the diamagnetic coupling term, which is represented by 
τ̂_3 and iτ̂_2 
in Eq. (<ref>). 
(τ̂_3 and iτ̂_2 indicate 
the density fluctuation and the phase mode, respectively. 
The amplitude mode in Sect. 2 comes from τ̂_1.) 










§ SUMMARY AND DISCUSSION


In this paper a theoretical study of the third-harmonic generation 
in dirty BCS superconductors was carried out. 
We calculated the temperature dependence and frequency dependence 
of the THG intensity, and showed that 
the vertex correction term including the amplitude mode 
is dominant. 
We introduced the effect of paramagnetic impurities, 
and showed that time-reversal symmetry breaking 
destabilizes the amplitude mode and that this effect is 
reflected in the THG intensity. 

We showed that the dependences of the THG intensity on temperature 
and pulse width reproduce experimental results 
at the limit where the effect of paramagnetic impurities 
is small. 
Quantitatively, almost the same result as the experimental 
result was obtained. 
In addition, it is known that the diamagnetic term is 
small in the experiment. <cit.> 
The calculated result of Sect. 4 gives an explanation 
of this. 
This result is based on the fact that the reduction 
effect due to the superconducting gap being smaller 
than the Fermi energy is more dominant than the increase caused 
by the mean free path being shorter than the coherence length. 


It is known that the effect of time-reversal symmetry breaking 
(due to paramagnetic impurities in this paper) 
can also be obtained by applying an in-plane magnetic field. <cit.> 
Referring to Figs. 2(a) and 2(b), 
the dependence of the THG intensity on α_p differs between 
the term including the amplitude mode and that not including 
this mode. 
Therefore, the results of our calculation show that 
the presence of the amplitude mode can be 
confirmed experimentally by applying 
a magnetic field in the plane. 











§ ACKNOWLEDGEMENT


The numerical computation in this work was carried out 
at the Yukawa Institute Computer Facility. 







9

giannetti C. Giannetti, M. Capone, D. Fausti, 
M. Fabrizio, F. Parmigiani, and D. Mihailovic, 
Adv. Phys. 65, 58 (2016). 

schmid A. Schmid, 
Phys. Kondens. Mater. 8, 129 (1968). 

volkov A. F. Volkov and S. M. Kogan, 
Sov. Phys. JETP 38, 1018 (1974). 

anderson58 P. W. Anderson, 
Phys. Rev. 112, 1900 (1958).

sk R. Sooryakumar and M. V. Klein, 
Phys. Rev. Lett. 45, 660 (1980). 

littlewood P. B. Littlewood and C. M. Varma, 
Phys. Rev. B 26, 4883 (1982). 

matsunaga R. Matsunaga, Y. I. Hamada, K. Makise, Y. Uzawa, 
H. Terai, Z. Wang, and R. Shimano, 
Phys. Rev. Lett. 111, 057002 (2013).

matsunaga14 R. Matsunaga, N. Tsuji, H. Fujita, A. Sugioka, 
K. Makise, Y. Uzawa, H. Terai, Z. Wang, H. Aoki, and R. Shimano, 
Science 345, 1145 (2014).

note0 
In Ref. 8, 

a theoretical work <cit.> is cited as an explanation 
of the amplitude mode in the THG intensity. 
In this theory, however, the paramagnetic coupling 
term <cit.> is not included because 
this calculation takes the clean and local limit. 
The linear absorption spectrum in the experiment <cit.> 
indicates that this coupling term is effective 
because there is a clear gap structure in the spectrum. 
Thus, the paramagnetic coupling term cannot be neglected 
in calculations of response functions, <cit.> and 
it makes a predominant contribution to the THG intensity 
as shown in Sect. 3. 



jujo15 T. Jujo, 
J. Phys. Soc. Jpn. 84, 114711 (2015). 

ivlev B. I. Ivlev, 
JETP Lett. 15, 313 (1972). 

matsunaga17 R. Matsunaga, N. Tsuji, 
K. Makise, H. Terai, H. Aoki, and R. Shimano, 
Phys. Rev. B 96, 020505(R) (2017). 


eliashberg G. M. Éliashberg, 
Zh. Eksp. Teor. Fiz. 61, 1254 (1971) 
[Sov. Phys. JETP 34, 668 (1972)]. 

keldysh L. V. Keldysh, 
Sov. Phys. JETP 20, 1018 (1965). 

usadel K. D. Usadel, 
Phys. Rev. Lett. 25, 507 (1970). 

kopnin N. B. Kopnin, 
Theory of Nonequilibrium Superconductivity 
(Oxford University Press, New York, 2001), Chaps. 5 and 9. 

abrikosov59 A. A. Abrikosov and L. P. Gor'kov, 
Sov. Phys. JETP 8, 1090 (1959). 

abrikosov61 A. A. Abrikosov and L. P. Gor'kov, 
Sov. Phys. JETP 12, 1243 (1961). 

skalski S. Skalski, O. Betbeder-Matibet, and P. R. Weiss, 
Phys. Rev. 136, A1500 (1964). 

anderson59 P. W. Anderson, 
J. Phys. Chem. Solids 11, 26 (1959). 

vollhardt D. Vollhardt and P. Wölfle, 
Phys. Rev. B 22, 4666 (1980). 

jujo17 T. Jujo, 
J. Phys. Soc. Jpn. 86, 024709 (2017). 

jujolt28 T. Jujo, 
presented at LT28, 28th Int. Conf. on Low Temperature Physics, 
2017. 

shen 
Y. R. Shen, The Principles of Nonlinear Optics 
(Wiley, New Jersey, 2003) Chap. 7. 

ikebe Y. Ikebe, R. Shimano, M. Ikeda, T. Fukumura, 
and M. Kawasaki, 
Phys. Rev. B 79, 174525 (2009). 

semenov A. Semenov, B. Günther, U. Böttger, 
H.-W. Hübers, H. Bartolf, A. Engel, A. Schilling, 
K. Ilin, M. Siegel, R. Schneider, D. Gerthsen, 
and N. A. Gippius, 
Phys. Rev. B 80, 054510 (2009). 

note1 Terms such as 
the paramagnetic coupling and the diamagnetic coupling 
between electrons and electromagnetic fields 
are the same as those in Ref. 31. 


note2 It is known that 
only the anisotropic term remains 
in the diamagnetic term. <cit.> 
The value of the anisotropic term 
with vertex corrections included 
is at most several times the value of the isotropic term without vertex 
corrections. 
Thus, it is sufficient to evaluate the latter case 
in order to find out whether the diamagnetic term 
or the paramagnetic term dominates.

maki
K. Maki, in Superconductivity, 
ed. R. D. Parks (Dekker, New York, 1969) Chap. 18.

tsuji 
N. Tsuji and H. Aoki, 
Phys. Rev. B 92, 064508 (2015).

schrieffer J. R. Schrieffer, 
Theory of Superconductivity 
(Addison-Wesley, Redwood City, CA, 1983) revised ed., Chap. 8. 

abrikosov74 A. A. Abrikosov and V. M. Genkin, 
Sov. Phys. JETP 38, 417 (1974). 










