# Optical response of Higgs mode in superconductors at clean limit:
formulation through Eilenberger equation and Ginzburg-Landau Lagrangian

F. Yang<EMAIL_ADDRESS>Hefei National Laboratory for Physical
Sciences at Microscale, Department of Physics, and CAS Key Laboratory of
Strongly-Coupled Quantum Matter Physics, University of Science and Technology
of China, Hefei, Anhui, 230026, China  M. W. Wu<EMAIL_ADDRESS>Hefei
National Laboratory for Physical Sciences at Microscale, Department of
Physics, and CAS Key Laboratory of Strongly-Coupled Quantum Matter Physics,
University of Science and Technology of China, Hefei, Anhui, 230026, China

###### Abstract

Both macroscopic Ginzburg-Landau Lagrangian and microscopic gauge-invariant
kinetic equation suggest a finite Higgs-mode generation in the second-order
optical response of superconductors at clean limit, whereas the previous
derivations through the path-integral approach and Eilenberger equation within
the Matsubara formalism failed to give such generation. The crucial treatment
leading to this controversy lies at an artificial scheme that whether the
external optical frequency is taken as continuous variable or bosonic
Matsubara frequency to handle the gap dynamics within the Matsubara formalism.
To resolve this issue, we derive the effective action of the superconducting
gap near $T_{c}$ in the presence of the vector potential through the path-
integral approach, to fill the long missing blank of the microscopic
derivation of the Ginzburg-Landau Lagrangian in superconductors. It is shown
that only by taking optical frequency as continuous variable within the
Matsubara formalism, can one achieve the fundamental Ginzburg-Landau
Lagrangian, and in particular, the finite Ginzburg-Landau kinetic term leads
to a finite Higgs-mode generation at clean limit. To further eliminate the
confusion of the Matsubara frequency through a separate framework, we apply
the Eilenberger equation within the Keldysh formalism, which is totally
irrelevant to the Matsubara space. By calculating the gap dynamics in the
second-order response, it is analytically proved that the involved optical
frequency is a continuous variable rather than bosonic Matsubara frequency,
causing a finite Higgs-mode generation at clean limit.

###### pacs:

74.40.Gh, 74.25.Gz, 74.25.N-

##  I Introduction

In the past few decades, the Higgs mode in the field of superconductivity,
which describes the amplitude fluctuation $\delta|\Delta|$ of the
superconducting order parameter $\Delta$, has attracted much attention. This
collective excitation corresponds to the radial excitation in the Mexican-hat
potential of free energyAm0 , and hence, exhibits a gapful energy spectrum at
a long wavelengthOD1 ; OD2 ; OD3 ; pm5 ; Am0 ; Am5 ; Am6 ; Am12 ; AK2 . Owing
to the advanced ultrafast terahertz technique in nonlinear optics, the Higgs
mode has been experimentally observed and identified as the origin of the
excited superfluid-density oscillation $\delta\rho_{s}$ in the second-order
harmonic generationNL7 ; NL8 ; NL9 ; NL10 ; NL11 ; DHM2 ; DHM3 . The nonlinear
optics in superconductors has since stimulated a lot of experimental interest
and inspired a great deal of theoretical studies.

Whereas the existing and growing experimental observations exhibit a very
convincing evidence, the theoretical descriptions concerning the Higgs-mode
excitation in the literature are filled with controversies, which are
detrimental to the understanding of the related experimental findings. The
central issue lies at a question that with the conventional light-matter
interactions $H_{d}={{\hat{\bf p}}\cdot{e{\bf A}}}\tau_{0}/{m}$ (current-
vertex-related drive effect) and $H_{p}=e^{2}A^{2}\tau_{3}/(2m)$ (density-
vertex-related pump effect) by the vector potential ${\bf A}$Ba0 ; G1 ; GIKE2
, whether the Higgs mode can be optically excited at clean limit. Here,
$\tau_{i=0,1,2,3}$ denote the Pauli matrices in Nambu space.

Concerning this issue, although the early stage of works through the BlochAm1
; Am2 ; Am7 ; Am9 ; Am11 ; Am14 ; Am15 ; NL7 ; NL8 ; NL9 ; NL10 ; NL11  or
LiouvilleAm3 ; Am4 ; Am8 ; Am10 ; Am16  equation within the Anderson
pseudospin pictureAs  revealed an excited fluctuation of the order parameter
by pump effect $H_{p}$, a later symmetry analysissymmetry  implies that the
pumped order-parameter fluctuation by density-vertex-related $H_{p}$ is a
phase fluctuation rather than the claimed amplitude oneGIKE2 , as the revealed
correlation between amplitude (phase) mode and $H_{p}$ in Ref. symmetry,  is
zero (nonzero). The path-integral approach is naturally capable of
distinguishing the excitations of phase and Higgs modes by deriving the
corresponding effective actions, and includes both pump and drive effects.
Using this approach, Cea et al. derived a vanishing Higgs-mode generation in
second-order response at clean limitCea1 ; Cea2 ; Cea3 , in inconsistency with
the previous experimental understandingNL7 ; NL8 ; NL9 . To explain the
experimental findings, Cea et al. pointed outCea1 ; Cea2 ; Cea3  that the
density-vertex-related pump effect $H_{p}$ can excite a finite fluctuation
$\delta{n}$ of charge density $n$, so they speculated that the experimentally
observed superfluid-density oscillation $\delta\rho_{s}$ is attributed to
charge-density fluctuation $\delta{n}$ rather than the Higgs mode
$\delta|\Delta|$, since the superfluid density
$\rho_{s}\propto{n}|\Delta|^{2}$.

In the several polarization-resolved measurements afterwardsNL10 ; NL11 ; DHM2
; DHM3 , an isotropic second-order harmonic signal was timely reported and
provides a firm evidence to rule out the possible charge-density fluctuation,
as the theoretically predicted response of the Higgs mode (charge-density
fluctuation) is isotropic (anisotropic)Cea1 . Since then, it is believed that
the Higgs-mode generation is zero at clean limit and one has to reply on
impurity scattering to mediate the Higgs-mode generationFHM ; Am16 ; ImR1 ;
ImR2 ; ImR3 ; Silaev . In this situation, to take account of the microscopic
scattering, Silaev applied the Eilenberger equationEilen  within the Matsubara
formalismSilaev ; Silaev0 , which solely includes the current-vertex-related
drive effect $H_{d}$. He also derived a vanishing Higgs-mode generation at
clean limit, but with impurities, a finite one to dominate over the charge-
density fluctuation is obtained. Particularly, Silaev showed that the impurity
scattering only mediate the Higgs-mode excitation and is incapable of causing
the damping of this collective excitationSilaev , so the increase of the
impurity density can enhance the optical signal of Higgs mode.

Meanwhile, using a gauge-invariant kinetic equation approachGIKE1 ; GIKE2 ;
GIKE3  with complete electromagnetic effect and microscopic scattering, Yang
and Wu derived totally opposite results. They obtain a finite Higgs-mode
generation contributed by the drive effect at clean limitGIKE2 , and show that
the charge-density fluctuation in fact vanishes in the second-order response,
as a consequence of the charge conservation and forbidden second-order
harmonic current in systems with the spacial inversion symmetry. The revealed
Higgs-mode generation at clean limit can capture the experimental observation
wellGIKE2 , and in particular, a finite damping/lifetime of the Higgs-mode
excitation by impurities is also derivedGIKE3 , providing a possible origin
for the experimentally observed broadening of the Higgs-mode resonance signal
as well as the fast Higgs-mode damping after optical excitation. This damping
agrees with the analysis of Heisenberg equation of motion since the Higgs-mode
excitation and electron-impurity interaction are non-commutative in Nambu
space.

This finite Higgs-mode generation has therefore been in sharp contrast to the
aforementioned vanishing one from Eilenberger equation and path-integral
approach. Actually, at clean limit, the finite Higgs-mode generation in
second-order response of superconductors is a direct consequence of the
Ginzburg-Landau LagrangianEPM . This is because that from the time-dependent
Ginzburg-Landau superconducting Lagrangian at clean limit, which was proposed
by Pekker and Varma through the symmetry analysis and Lorentz invariance from
the Landau phase-transition theoryAm6 , one can directly reveal the equation
of motion of the Higgs mode by considering the amplitude fluctuation of the
Landau order parameter. Then, a finite Higgs-mode generation in the second-
order response is immediately obtainedEPM . This directly leads to a
particularly bizarre question that why both path-integral approachPI2GL  and
Eilenberger equationBa20  can recover the Ginzburg-Landau equation but reach a
zero Higgs-mode generation.

To resolve this controversy, by re-examining the previous derivations within
the path-integral approach in Refs. Cea1, ; Cea2, ; Cea3,  and Eilenberger
equation in Ref. Silaev, , it is pointed out by Yang and WuEPM  that both
previous derivations contain flaws. Specifically, the previous worksCea1 ;
Cea2 ; Cea3  within the path-integral approach only kept the perturbation
expansion of the action up to second order, whereas the essential coupling of
the Higgs mode to the second order of the drive effect $H_{d}$ emerges in the
third-order perturbation expansion of the actionEPM . So this coupling and
hence finite second-order harmonic generation of Higgs mode by the drive
effect are excessively overlooked in Refs. Cea1, ; Cea2, ; Cea3, . Within the
Matsubara formalism, picking up this coupling in the path-integral approach,
one can find the exactly same amplitude-response coefficient as the one
derived from Eilenberger equationSilaev :

$\displaystyle\lambda_{E}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\frac{T}{2\Omega^{2}}\sum_{p_{n}}\Big{[}\frac{2}{\sqrt{(p_{n}\\!-\\!i\Omega)^{2}\\!+\\!\Delta_{0}^{2}}}\\!-\\!\frac{1}{\sqrt{(p_{n}\\!-\\!2i\Omega)^{2}\\!+\\!\Delta_{0}^{2}}}$
(1)
$\displaystyle\mbox{}-\\!\frac{1}{\sqrt{p_{n}^{2}\\!+\\!\Delta_{0}^{2}}}\Big{]},$

where $p_{n}=(2n+1)\pi{T}$ represents the fermionic Matsubara frequencies.

Nevertheless, in Ref. Silaev, , the involved optical frequency $\Omega$ is
taken as bosonic Matsubara frequency $i\Omega_{m}$, leading to a vanishing
response coefficient $\lambda_{E}$ in Eq. (1) strongly against the finite one
from gauge-invariant kinetic equationEPM ; GIKE2  and Ginzburg-Landau
LagrangianEPM . Moreover, because of this treatment, the prefactor
$1/\Omega^{2}$ in Eq. (1) causes an undefined singularity at zero frequency,
and an unphysical discontinuity between cases at $\Omega=0$ and
$\Omega\rightarrow 0$ emergesEPM . In contrast, the previous work in Ref. EPM,
takes $\Omega$ as continuous variable. A finite Higgs-mode generation at clean
limit in agreement with the Ginzburg-Landau Lagrangian is then derived, and
the obtained $\lambda_{E}$ from both Eilenberger equation and path-integral
approach becomes exactly same as the one from gauge-invariant kinetic
equation.

Actually, the Matsubara formalism is developed as an auxiliary-function
technique in the finite-temperature Green function approach. In this
framework, whether taking the external optical frequency $\Omega$ as
continuous variable or bosonic Matsubara frequency $i\Omega_{m}$ can not be
self-justified by method itself. The treatment of $\Omega$ as $i\Omega_{m}$ in
the calculations of conductivity and dielectric function in normal metalsmt
can be cross-justified by many other methods irrelevant to Matsubara space,
such as equation of motion, zero-temperature Green function, Boltzmann
equation as well as Keldysh Green function approaches, whereas such
justification in the gap dynamics of superconductors has not been performed in
the literature so far. Physically, any treatment leading to result strongly
against the Ginzburg-Landau superconducting Lagrangian can not be correct.
Nevertheless, this justification in superconductors has been challenged,
arguing that the Ginzburg-Landau superconducting Lagrangian is a
phenomenological model near $T_{c}$ and is unimportant in microscopic studies,
as the microscopic derivation of this Lagrangian is absent in the literature.
Theoretically, this Lagrangian is a fundamental model by symmetry analysis and
Lorentz invariance as well as Landau phase-transition theory, whereas to fill
the long missing blank of the microscopic derivation, in the present work, we
derive the Lagrangian in superconductors through the basic path-integral
approach. In addition to this physical justification, a natural and rigorous
framework developed in the literature to eliminate the confusion of the
Matsubara frequency in superconductors is to perform the formulation within
the Keldysh formalismQA1 , which is a well-established systematic approach for
studying non-equilibrium properties and is totally irrelevant to Matsubara
space. We therefore apply the Eilenberger equation within the Keldysh
formalism to calculate the gap dynamics at clean limit to cross-justify the
treatment of the external optical frequency.

Specifically, through the path-integral approach within the Matsubara
formalism, we derive the effective action of superconducting gap near $T_{c}$
in the presence of the vector potential. We show that to recover the Ginzburg-
Landau kinetic term, one needs to keep the perturbation expansion of the
action up to the fourth order and formulate the fourth-order correlation
coefficient. Particularly, during our calculation of the correlation
coefficient, it is shown that only by taking optical frequency as continuous
variable, one can recover the finite coefficient in the Ginzburg-Landau
kinetic term. Then, if one considers the gap fluctuation to obtain the
equation of motion of the Higgs mode from the Ginzburg-Landau Lagrangian, it
is clearly seen that the finite coefficient in the Ginzburg-Landau kinetic
term directly leads to the finite response coefficient in the equation of
motion of the Higgs mode.

Furthermore, with the vector potential alone, we apply the Eilenberger
equation within the Keldysh formalism to derive the equation of motion of the
Higgs mode at clean limit. It is established in the literature that the
Keldysh Green function can be written as a function of the retarded and
advanced ones through a general relation via the distribution function. We
prove that this relationship makes the Keldysh Green function directly
satisfying the normalization condition of the Eilenberger equation. Then, we
solve the retarded Green function and distribution function to obtain the
Keldysh Green function. With the derived Keldysh Green function, we obtain the
equation of motion of the Higgs mode, which is exactly same as the one derived
through Eilenberger equation within the Matsubara formalismSilaev . In
contrast, as our derivation is performed in the Keldysh formalism and is
irrelevant to the Matsubara space, it is clearly shown that the involved
optical frequency in the response coefficient in the equation of motion of the
Higgs mode is a continuous variable rather than the Matsubara frequency.
Particularly, with the continuous optical frequency, the response coefficient
in the equation of motion of the Higgs mode does not vanish, leading to a
finite Higgs mode generation at clean limit.

##  II Hamiltonian

We first present the general Bogoliubov-de Gennes Hamiltonian of the
conventional $s$-wave superconductors in the presence of the electromagnetic
potentialG1 ; Ba0 ; AK ; GIKE2 :

$H\\!=\\!{\int}{d{\bf
x}}~{}\psi^{\dagger}(x)\left(\begin{array}[]{cc}\xi_{{\hat{\bf p}}-e{\bf
A}}+e\phi&\Delta(x)\\\ \Delta^{*}(x)&-\xi_{{\hat{\bf p}}+e{\bf
A}}-e\phi\end{array}\right)\psi(x).$ (2)

Here, $\psi(x)=[\psi_{\uparrow}(x),\psi^{\dagger}_{\downarrow}(x)]^{T}$ is the
field operator in the Nambu space and $x=(t,{\bf x})$ represents the space-
time vector; $\phi$ and ${\bf A}$ denote the scalar and vector potentials,
respectively; the momentum operator ${\hat{\bf p}}=-i\hbar{\bm{\nabla}}$ and
$\xi_{\hat{\bf p}}={{\bf{\hat{p}}}^{2}}/({2m})-\mu$ with $m$ being the
effective mass and $\mu$ denoting the chemical potential. It is noted that the
scalar potential can be generally written as $\phi={\bar{\phi}}_{0}+{\bf
E}_{\phi}\cdot{\bf x}+\delta\phi$, where ${\bar{\phi}}_{0}$ denotes the effect
of the electric voltage; ${\bf E}_{\phi}\cdot{\bf x}$ concerns the drive
effect by electric field; $\delta\phi$ is the induced scalar potential related
to the long-range Coulomb interaction (i.e., Hartree field caused by charge
density fluctuation)Ba0 ; AK ; GIKE2 .

In consideration of the phase $\delta\theta(x)$ and amplitude
$\delta|\Delta(x)|$ fluctuations around the equilibrium gap $\Delta_{0}$, the
superconducting order parameter reads:

$\Delta(x)=|\Delta(x)|e^{i\delta\theta(x)}=[\Delta_{0}+\delta|\Delta(x)|]e^{i\delta\theta(x)}.$
(3)

It is established that the phase mode $\delta\theta$ in Hamiltonian above can
be effectively removed by a unitary transformationgi0 ; gi1 ; AK

$\psi(x){\rightarrow}e^{i\tau_{3}\delta\theta(x)/2}\psi(x),$ (4)

and then, one has

$H=\\!\\!\\!{\int}{d{\bf x}}~{}\psi^{\dagger}(x)(H_{0}+H_{\rm LM})\psi(x),$
(5)

where the free BCS Hamiltonian $H_{0}$ is written as

$H_{0}=\xi_{\bf{\hat{p}}}\tau_{3}+|\Delta(x)|\tau_{1},$ (6)

and the light-matter interaction reads

$H_{\rm LM}=\frac{{\bf p}_{s}\cdot{\hat{\bf
p}}}{m}+\frac{p_{s}^{2}}{2m}\tau_{3}+\mu_{\rm eff}\tau_{3}$ (7)

with the gauge-invariant superconducting momentum ${\bf
p}_{s}={\bm{\nabla}}\delta\theta/2-e{\bf A}$ and effective field $\mu_{\rm
eff}=e\phi+\partial_{t}\delta\theta/2$.

It has been revealed in the literatureAK ; GIKE2  that thanks to the Coulomb
screening, in the linear regime, at long-wavelength limit, the induced scalar
potential $e\delta\phi$ cancels the original longitudinal part
$e{\bar{\phi}}_{0}+\partial_{t}\delta\theta/2$ in the effective field
$\mu_{\rm eff}=e{\bar{\phi}}_{0}+e{\bf E}_{\phi}\cdot{\bf
x}+e\delta\phi+\partial_{t}\delta\theta/2$, leaving only the drive effect of
scalar-potential-induced electric field $E_{\phi}$. In general, the inclusion
of this retained effect from scalar potential is essential to capture the
optical-electric-field effect in a gauge-invariant manner. Whereas considering
the fact that the vector potential characterizes the Meissner effect/Ginzburg-
Landau kinetic term in addition to the optical-electric-field effectGIKE1 ;
GIKE2 ; EPM ; PYW , in the present work, we only focus on the electromagnetic
effect from vector potential, similar to the previous works by Cea et al.Cea1
; Cea2 ; Cea3  and SilaevSilaev . The light-matter interaction then becomes

$H_{\rm LM}=\frac{{\bf p}_{s}\cdot{\hat{\bf
p}}}{m}+\frac{p_{s}^{2}}{2m}\tau_{3}.$ (8)

Moreover, it has been established in the literatureAm0 ; Ba0 ; AK ; AK2 ; pm0
; pm5 ; pi1 ; pi4 ; GIKE2  that the linear response of the phase mode, as a
scalar quantity, responds to the longitudinal electromagnetic field and hence
experiences the Coulomb screening. The resonance pole of this response is then
effectively lifted from the original gapless spectrum up to the high-energy
plasma frequency $\omega_{p}$ as a consequence of the Anderson-Higgs
mechanismAHM , and hence, no effective linear response of the phase mode
occurs at frequency $\Omega\ll\omega_{p}$. Because of this effect, the linear
response of the phase mode cancels the unphysical longitudinal vector
potential in ${\bf p}_{s}$, and the superconducting momentum that appears in
the previous theoretical descriptions in the literature only involves the
physical transverse vector potential. The light-matter interaction in Eq. (8)
then consists of the drive effect $H_{d}$ and pump one $H_{p}$

##  III Derivation of Ginzburg-Landau Lagrangian

In this section, we present the derivations of the time-dependent Ginzburg-
Landau Lagrangian and its non-equilibrium variation (equation of motion of the
Higgs mode). Specifically, the action of superconductors in the presence of
vector potential after the Hubbard-Stratonovich transformation is written
asBa0 ; pi1 ; pi4 ; PI2GL

$\displaystyle\\!\\!\\!\\!\\!S[\psi,\psi^{*}]=\\!\\!\int{dx}\bigg{[}\sum_{s=\uparrow,\downarrow}\\!\\!\psi^{*}_{s}(x)(i\partial_{t}\\!-\\!\xi_{\hat{\bf
p}-e{\bf A}})\psi_{s}(x)$
$\displaystyle\mbox{}+\\!\psi^{*}_{\uparrow}(x)\psi^{*}_{\downarrow}(x)\Delta(x)\\!+\\!\psi_{\downarrow}(x)\psi_{\uparrow}(x)\Delta^{*}(x)\\!-\\!\frac{|\Delta(x)|^{2}}{U}\bigg{]}.~{}~{}~{}~{}~{}~{}$
(9)

Considering the spatial dependence of the gap, one has $|\Delta(x)|=\sum_{\bf
q}|\Delta_{\bf q}|e^{i{\bf q}\cdot{\bf x}}$ in Fourier space. Then, applying
the unitary transformation in Eq. (4), the above action in momentum space
becomes

$\displaystyle\\!\\!\\!\\!S[\psi,\psi^{*}]=\int{dt}\Big{\\{}\sum_{\bf
kq}\big{[}\psi^{*}_{{\bf k+\frac{q}{2}}\uparrow}(i\partial_{t}\\!-\\!\xi_{\bf
k+\frac{q}{2}+p_{s}})\psi_{{\bf k+\frac{q}{2}}\uparrow}$
$\displaystyle\mbox{}+\psi^{*}_{{\bf-k+\frac{q}{2}}\downarrow}(i\partial_{t}\\!-\\!\xi_{\bf-k+\frac{q}{2}+p_{s}})\psi_{{\bf-k+\frac{q}{2}}\downarrow}+(\psi_{{\bf-k+\frac{q}{2}}\downarrow}\psi_{{\bf
k+\frac{q}{2}}\uparrow}$ $\displaystyle\mbox{}+\psi_{{\bf
k+\frac{q}{2}}\uparrow}^{*}\psi_{{\bf-k+\frac{q}{2}}\downarrow}^{*})|\Delta_{\bf
q}|\big{]}-\sum_{\bf q}\frac{|\Delta_{\bf q}|^{2}}{U}\Big{\\}},$ (10)

and in Nambu space, one finds the action related to gap:

$\displaystyle S[\psi,\psi^{*}]=\\!\\!\int{dt}\sum_{\bf q}\\!\Big{[}\sum_{\bf
k}\psi^{\dagger}_{\bf kq}(G_{0}^{-1}\\!-\\!\Sigma)\psi_{\bf
kq}\\!-\\!\frac{|\Delta_{\bf q}|^{2}}{U}\Big{]}.$ (11)

Here, the field operator $\psi_{\bf kq}^{\dagger}=(\psi_{{\bf
k+\frac{q}{2}}\uparrow}^{*},\psi_{{\bf-k+\frac{q}{2}}\downarrow})$;
$G_{0}^{-1}=i\partial_{t}-\xi_{\bf k}\tau_{3}$ which gives the Green function
$G_{0}(p)=(p_{0}-\xi_{\bf k}\tau_{3})^{-1}$ in frequency space with the four-
vector momentum $p=(p_{0},{\bf k})$; the self-energy reads

$\Sigma=|\Delta_{\bf q}|\tau_{1}+\frac{{\bf k}\cdot({\bf
q/2+p_{s}})}{m}+\frac{({\bf q/2+p_{s}})^{2}}{2m}\tau_{3}.$ (12)

Considering the small gap near critical temperature $T_{c}$ as well as the
weak strength and spatial variation of the vector potential, the self-energy
can be treated as small quantity. Then, after the standard integration over
the Fermi field, one has

$S=S_{0}\\!-\sum_{n=1}^{\infty}\frac{1}{n}{\rm{\bar{T}r}}[(G_{0}\Sigma)^{n}]\\!-\\!\int{dt}\sum_{\bf
q}\frac{|\Delta_{\bf q}|^{2}}{U}.$ (13)

To derive the Lagrangian related to the superconductivity, one needs to
formulate the expansions of the action with respect to the fourth order of the
self-energy (i.e., keep the expansions up to $n=4$). Then, with expansions up
to $n=4$, by only keeping the terms related to the gap, one can obtain the
effective action $S_{s}=\int{dq}\mathscr{L}$ with the frequency-momentum four
vector $q=(\Omega,{\bf q})$ and the Lagrangian of the superconductivity:

$\displaystyle\mathscr{L}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!-\chi_{1}|\Delta_{\bf
q}|\\!-\\!\big{(}\frac{1}{2}\chi_{11}+\frac{1}{U}\big{)}|\Delta_{\bf
q}|^{2}\\!-\\!\chi_{13}\frac{({\bf q/2\\!+\\!p_{s}})^{2}|\Delta_{\bf q}|}{2m}$
(14) $\displaystyle\\!\\!\mbox{}-\\!\chi_{111}\frac{|\Delta_{\bf
q}|^{3}}{3}\\!-\\!\chi_{100}\frac{k_{F}^{2}({\bf
q/2\\!+\\!p_{s}})^{2}|\Delta_{\bf
q}|}{3m^{2}}\\!-\\!\frac{\chi_{1111}|\Delta_{\bf q}|^{4}}{4}$
$\displaystyle\\!\\!\mbox{}-\\!\chi_{113}\frac{({\bf
q/2\\!+\\!p_{s}})^{2}|\Delta_{\bf q}|^{2}}{2m}\\!-\\!\chi_{1113}\frac{({\bf
q/2\\!+\\!p_{s}})^{2}|\Delta_{\bf q}|^{3}}{2m}$
$\displaystyle\\!\\!\mbox{}-\\!(\chi_{1100}+\chi_{0110}+\chi_{1010})\frac{k_{F}^{2}({\bf
q/2\\!+\\!p_{s}})^{2}|\Delta_{\bf q}|^{2}}{6m^{2}}.$

Here, we have neglected the terms proportional to odd orders of ${{\bf
k}\cdot({\bf q/2+p_{s}})}$ as these anisotropic terms vanish after the
summation of the momentum ${\bf k}$. The correlation coefficients are
determined by

$\displaystyle\chi_{i}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{p}{\rm Tr}[G_{0}(p)\tau_{j}],$ (15)
$\displaystyle\chi_{ij}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{p}{\rm
Tr}[G_{0}(p\\!+\\!q)\tau_{i}G_{0}(p)\tau_{j}],$ (16)
$\displaystyle\chi_{ijk}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{p}{\rm
Tr}[G_{0}(p\\!+\\!2q)\tau_{i}G_{0}(p\\!+\\!q)\tau_{j}G_{0}(p)\tau_{k}],$ (17)
$\displaystyle\chi_{ijkl}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{p}{\rm
Tr}[G_{0}(p\\!+\\!q)\tau_{i}G_{0}(p)\tau_{j}G_{0}(p\\!+\\!q)\tau_{k}G(p)\tau_{l}].~{}~{}~{}~{}~{}$
(18)

Since the Green function only consists of the $\tau_{0}$ and $\tau_{3}$
components, one immediately finds
$\chi_{1}=\chi_{13}=\chi_{111}=\chi_{100}=\chi_{1113}=0$. Moreover, the
coefficient $\chi_{113}$ vanishes due to the particle-hole symmetry, and the
forth-order correlation coefficient $\chi_{1010}=\chi_{0110}+\chi_{1100}$
(refer to Appendix). Then, one finds an embryonic form of the Ginzburg-Landau
Lagrangian of superconductors:

$\mathscr{L}=-\chi_{p}|\Delta_{\bf q}|^{2}-\frac{\beta_{p}|\Delta_{\bf
q}|^{4}}{2}-\frac{\lambda_{p}({\bf q/2\\!+\\!p_{s}})^{2}|\Delta_{\bf
q}|^{2}}{m}.$ (19)

with the parameters $\chi_{p}={\chi_{11}}/{2}+{1}/{U}$ and
$\beta_{p}=\chi_{1111}/2$ as well as $\lambda_{p}=k_{F}^{2}\chi_{1010}/(3m)$ .

Within the Matsubara formalism $[p=(ip_{n},{\bf k})]$, the retained
correlation coefficients in Eq. (19) are given by (refer to Appendix A)

$\displaystyle\chi_{p}$ $\displaystyle=$
$\displaystyle\alpha_{p}-\Omega^{2}\gamma_{p}/2,$ (20)
$\displaystyle\beta_{p}$ $\displaystyle=$
$\displaystyle-\\!\\!\\!\\!\sum_{ip_{n>0},\eta=\pm}\frac{2\pi{iD}}{\beta}\frac{2}{(2ip_{n}\\!+\\!\eta\Omega)^{3}},$
(21) $\displaystyle\lambda_{p}$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\sum_{ip_{n>0},\eta=\pm}\\!\\!\frac{k_{F}^{2}}{3m}\frac{2\pi{i}D}{\beta\Omega^{2}}\Big{[}\frac{4}{2ip_{n}\\!+\\!\eta\Omega}\\!-\\!\frac{1}{ip_{n}}\\!-\\!\frac{1}{ip_{n}\\!+\\!\eta\Omega}\Big{]},~{}~{}~{}~{}~{}$
(22)

with parameters:

$\displaystyle\alpha_{p}$ $\displaystyle=$
$\displaystyle\\!\\!D\int^{\omega_{D}}_{-\omega_{D}}d\xi_{\bf
k}\frac{\tanh(\frac{\beta_{c}\xi_{\bf k}}{2})\\!-\\!\tanh(\frac{\beta\xi_{\bf
k}}{2})}{2\xi_{\bf k}}\\!=\\!D\ln\frac{T}{T_{c}},~{}~{}~{}~{}$ (23)
$\displaystyle\gamma_{p}$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\sum_{ip_{n>0},\eta=\pm}\\!\\!\frac{\pi{i}D}{\beta\Omega^{2}}\Big{[}\frac{4}{2ip_{n}\\!+\\!\eta\Omega}\\!-\\!\frac{1}{ip_{n}}\\!-\\!\frac{1}{ip_{n}\\!+\\!\eta\Omega}\Big{]}.$
(24)

Here, $D$ and $\omega_{D}$ denote the density of states and Debye frequency,
respectively; $f(x)$ stands for the Fermi distribution; $\beta=k_{B}T$ and
$\beta_{c}=k_{B}T_{c}$ with $k_{B}$ being the Boltzmann constant.

At low frequency ($\Omega<{k_{B}T_{c}}$), one finds the specific parameters:

$\displaystyle\gamma_{p}$ $\displaystyle=$
$\displaystyle\beta_{p}\approx\sum_{n>0}\frac{\pi{D}}{\beta{p_{n}^{3}}}=\frac{7D\zeta(3)}{8({\pi}T)^{2}},~{}~{}~{}~{}$
(25) $\displaystyle\lambda_{p}$ $\displaystyle\approx$
$\displaystyle\sum_{n>0}\frac{k_{F}^{2}}{3m}\frac{2\pi{D}}{4\beta}\frac{4}{p_{n}^{3}}=\frac{k_{F}^{2}}{3m}\frac{7D\zeta(3)}{(2{\pi}T)^{2}},$
(26)

with $\zeta(x)$ being the Riemann zeta function. Then, through the
transformation into time-space coordinate, the time-dependent Ginzburg-Landau
Lagrangian of superconductors is obtained:

$\mathscr{L}\\!=\\!\frac{\gamma_{p}|i\partial_{t}\Delta|^{2}}{2}\\!-\\!\Big{(}\alpha_{p}|\Delta|^{2}\\!+\\!\frac{\beta_{p}|\Delta|^{4}}{2}\\!+\\!\frac{\lambda_{p}|({\bm{\nabla}}\\!-\\!2ie{\bf
A})\Delta|^{2}}{4m}\Big{)},$ (27)

which is exactly same as the one obtained by Pekker and Varma through symmetry
analysis and Lorentz invariance from the general Ginzburg-Landau free
energyAm6  as it should be. Moreover, the parameters $\alpha_{p}$ [Eq. (23)]
and $\beta_{p}$ [Eq. (25)] as well as $\lambda_{p}$ [Eq. (26)] here are
exactly same as the obtained Landau parameters in the previous work by
GorkovG1  in which the Ginzburg-Landau equation is derived from Gorkov
equation.

It is noted that in the derivation of the gap dynamics within the Matsubara
formalism in the present work, we take the optical frequency $\Omega$ as
continuous variable rather than bosonic Matsubara frequency $i\Omega_{m}$, and
only with continuous $\Omega$ in this circumstance, can one recover the
Ginzburg-Landau Lagrangian as demonstrated above. Actually, as seen from Eq.
(22), if $\Omega$ is taken as $i\Omega_{m}=2mi\pi{T}$, different irrational
response coefficients are obtained at cases with odd and even $m$. For odd
$m$, a divergent pole ($2n+1=m$ and $\eta=-1$) which is unable to circumvent
emerges in the formulation of the first term on the right-hand side of Eq.
(22). As for even $m$, through the frequency displacement in the Matsubara
frequency summation, the response coefficient $\lambda_{p}$ directly vanishes
for nonzero $m$, but due to the prefactor $1/\Omega^{2}$ in Eq. (22), there
exists an undefined singularity at zero frequency. Then, an unphysical
discontinuity between $\Omega=0$ and $\Omega\rightarrow 0$ emerges. The
difference between cases with odd and even $m$ in bosonic Matsubara frequency
$i\Omega_{m}=2mi\pi{T}$ is totally unreasonable within the Matsubara
formalism, and neither of them can recover the finite Landau parameter at low
frequency. Clearly, any treatment that leads to consequence strongly against
the fundamental model of symmetry analysis and Lorentz invariance as well as
Landau phase-transition theory can not be correct. All of the irrationalities
here simply suggest that the treatment of taking $\Omega$ as $i\Omega_{m}$ can
not be correct in the derivation of gap dynamics in superconductors.

Furthermore, from the Lagrangian above, with the equilibrium gap
$\Delta_{0}=\sqrt{-\alpha_{p}/\beta_{p}}$, by considering the gap fluctuation
$\delta|\Delta|$ through $|\Delta|=\Delta_{0}+\delta|\Delta|$, its equation of
motion at long-wavelength limit is directly obtained:

$\big{[}(2\Delta_{0})^{2}-\partial_{t}^{2}\big{]}\delta|\Delta|=-\frac{\lambda_{p}}{\gamma_{p}}\frac{e^{2}A^{2}}{m}2\Delta_{0}.$
(28)

Then, one immediately finds the Higgs-mode energy spectrum (i.e., resonance
pole) $\omega_{H}=2\Delta_{0}$ on the left-hand side of the equation above,
and in particular, a finite second-order response of Higgs mode at clean limit
on the right-hand side of the equation above. It is pointed out that although
the derivation of the equation of motion of the Higgs mode here is based on
the small gap and only holds near $T_{c}$, the serious derivation in regime
extending to $T=0$ through three different microscopic approaches including
the gauge-invariant kinetic equation, Eilenberger equation as well as path-
integral approach also obtains the similar equation of motion as a consequence
of the renormalization group (scaling) theory or basic local Abelian $U(1)$
model (complex scalar field coupled to an electromagnetic potential) in the
field theory.

##  IV Derivation of Higgs-mode response through Eilenberger equation

To eliminate the confusion of the Matsubara frequency from a separate
framework in addition to the physical justification above, in this section, we
perform the formulation of the gap dynamics within the Keldysh formalism,
which is totally irrelevant to Matsubara space. Specifically, we apply the
Eilenberger equation within the Keldysh formalism to derive the second-order
optical response of the Higgs mode at clean limit. The Eilenberger
equationEilen ; Eilen1  is derived from the basic Gorkov equation of
$\tau_{3}$-Green function through the quasiclassical approximationQA1 :

$g^{R/K/A}_{{\bf R,k_{F}}}(t,t^{\prime})=\frac{i}{\pi}\\!\\!\int{d\xi_{\bf
k}}\\!\\!\int{d{\bf r}}\tau_{3}G^{R/K/A}(x,x^{\prime})e^{-i{\bf k}\cdot({\bf
x}-{\bf x^{\prime}})}.$ (29)

Here, ${\bf R}=({\bf x+x^{\prime}})/{2}$ represents the center-of-mass spatial
coordinate; the retarded (R), advanced (A) and Keldysh (K) Green functions are
defined byQA1 ; Eilen1

$\displaystyle G^{R}(x,x^{\prime})$ $\displaystyle=$
$\displaystyle-i\langle\\{\psi(x),\psi^{\dagger}(x^{\prime})\\}\rangle\theta(t-t^{\prime}),$
(30) $\displaystyle G^{A}(x,x^{\prime})$ $\displaystyle=$ $\displaystyle
i\langle\\{\psi(x),\psi^{\dagger}(x^{\prime})\\}\rangle\theta(t^{\prime}-t),$
(31) $\displaystyle G^{K}(x,x^{\prime})$ $\displaystyle=$
$\displaystyle-i\langle[\psi(x),\psi^{\dagger}(x^{\prime})].$ (32)

The Eilenberger equation within the Keldysh formalism at clean limit
readsEilen ; Eilen1 :

$i\\{\tau_{3}\partial_{t},{\hat{g}}\\}_{t}-[(\Delta_{0}+\delta|\Delta|)\tau_{1}\tau_{3},{\hat{g}}]_{t}+[e{\bf
A}\cdot{\bf v}_{F}\tau_{3},{\hat{g}}]_{t}\\!=\\!0,$ (33)

where the green function matrices ${\hat{g}}$ is defined as

${\hat{g}}=\left(\begin{array}[]{cc}g^{R}&g^{K}\\\ 0&g^{A}\end{array}\right).$
(34)

Here, the operators
$[X,{\hat{g}}]_{t}=X(t_{1}){\hat{g}}(t_{1},t_{2})-{\hat{g}}(t_{1},t_{2})X(t_{2})$
and
$\\{X,{\hat{g}}\\}_{t}=X(t_{1}){\hat{g}}(t_{1},t_{2})+{\hat{g}}(t_{1},t_{2})X(t_{2})$.

Moreover, to guarantee the unique solution, the Eilenberger equation must be
supplemented by the normalization conditionEilen ; Eilen1 :

${\hat{g}}\circ{\hat{g}}=1,$ (35)

where the operator $\circ$ is defined by relation
$A\circ{B}=\int{dt}A(t_{1},t)B(t,t_{2})$.

The corresponding gap equation is written as

$\Delta_{0}+\delta|\Delta|=-iU{\rm Tr}[\langle{g^{K}_{\bf
R,k_{F}}(t,t)}\rangle_{F}\tau_{2}/2],$ (36)

with $\langle...\rangle_{F}$ denoting the angular average over the Fermi
surface.

Considering an external optical field with ${\bf A}(t)={\bf
A}_{0}e^{-i\Omega{t}}$, by self-consistently solving Eqs. (33)-(36), one can
formulate the Higgs-mode generation at clean limit. Specifically, in this
circumstance, one can expand the quasiclassical Green function matrices as

${\hat{g}}={\hat{g}}^{(0)}+\delta{\hat{g}}^{(1)}+\delta{\hat{g}}^{(2)},$ (37)

with the $m$-th order response $\delta{\hat{g}}^{(m)}$ on the initial state
${\hat{g}}^{(0)}$. Correspondingly, the Higgs-mode generation
$\delta|\Delta|=\delta|\Delta|^{(1)}e^{-i\Omega{t}}+\delta|\Delta|^{(2)}e^{-2i\Omega{t}}$
from Eq. (36) with $\delta|\Delta|^{(1)}$ and $\delta|\Delta|^{(2)}$ being the
excitations in first- and second-order optical responses, respectively.
Particularly, the first-order response of Keldysh Green function must be
anisotropic in momentum space, leading to a vanishing $\delta|\Delta|^{(1)}$
after the angular average over the Fermi surface in the gap equation. We then
directly take $\delta|\Delta|^{(1)}=0$ for convenience.

Consequently, the Eilenberger equation in Eq. (33) becomes a chain of
equations:

$\\{\tau_{3}\partial_{t},\delta{\hat{g}}^{(1)}\\}_{t}\\!+\\![\Delta_{0}\tau_{2},\delta{\hat{g}}^{(1)}]_{t}\\!-\\!i[e{\bf
A}\cdot{\bf v}_{F}\tau_{3},{\hat{g}}^{(0)}]_{t}\\!=\\!0,$ (38)
$\displaystyle\\{\tau_{3}\partial_{t},\delta{\hat{g}}^{(2)}\\}_{t}+[\Delta_{0}\tau_{2},\delta{\hat{g}}^{(2)}]_{t}-i[e{\bf
A}\cdot{\bf v}_{F}\tau_{3},\delta{\hat{g}}^{(1)}]_{t}$
$\displaystyle\mbox{}+[\delta|\Delta|^{(2)}e^{-2i\Omega{t}}\tau_{2},{\hat{g}}^{(0)}]_{t}\\!=\\!0,$
(39)

and one can solve $\delta{\hat{g}}^{(1)}$ and $\delta{\hat{g}}^{(2)}$ in
sequence with the given initial state ${\hat{g}}^{(0)}$. Then, with the
obtained solution of the Keldysh Green function, one can further derive the
response of the Higgs mode from the gap equation.

###  IV.1 Solution of retarded Green function

In this part, we first solve the retarded Green functions from Eq. (33). From
the expansion of Green function matrices in Eq. (37), the retarded Green
function is written as

$\displaystyle g^{R}(t,t^{\prime})$ $\displaystyle=$ $\displaystyle
g^{R(0)}(t,t^{\prime})+{\delta}g^{R(1)}(t,t^{\prime})+{\delta}g^{R(2)}(t,t^{\prime})$
(40) $\displaystyle=$
$\displaystyle\\!\\!\\!\int\\!\frac{dE}{2\pi}e^{iE(t^{\prime}-t)}[g^{R(0)}(E)\\!+\\!{\delta}g^{R(1)}(E)e^{-i{\Omega}t}$
$\displaystyle\mbox{}+{\delta}g^{R(2)}(E){e^{-2i{\Omega}t}}].~{}~{}~{}~{}~{}$

The initial state of the retarded Green function has been established in the
literature from Gorkov equationEilen1 ; G1 ; Silaev , and is written as

$g^{R(0)}(E)=\int\frac{d\xi_{\bf k}}{\pi}\frac{i\tau_{3}(E\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1})}{(E\\!+\\!i0^{+})^{2}\\!-\\!\xi_{\bf
k}^{2}\\!-\\!\Delta_{0}^{2}}=\frac{E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}}{S^{R}(E)},$
(41)

with $S^{R}(E)=\sqrt{(E+i0^{+})^{2}-\Delta_{0}^{2}}$.

By defining $E_{1}=E+\Omega$ and $E_{2}=E+2\Omega$, from Eq. (38), the
equation of the first order of retarded Green function is written as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{R(1)}(E)\\!-\\!{\delta}g^{R(1)}(E)(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}=e{\bf A}_{0}\cdot{\bf
v}_{F}\Pi^{R(0)}_{3},$ (42)

from which one finds the first-order solution (refer to Appendix B):

${\delta}g^{R(1)}(E)=(e{\bf A}_{0}\cdot{\bf
v}_{F})\frac{\tau_{3}-g^{R(0)}(E_{1})\tau_{3}g^{R(0)}(E)}{S^{R}(E_{1})+S^{R}(E)}.$
(43)

Here, $\Pi^{R(i)}_{3}=g^{R(i)}(E_{1})\tau_{3}-\tau_{3}g^{R(i)}(E)$. Similarly,
the equation of the second order of retarded Green function from Eq. (39)
reads

$\displaystyle\\!\\!\\!\\!\\!\\!\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{R(2)}(E)\\!-\\!{\delta}g^{R(2)}(E)(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}=e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}\Pi^{R(1)}_{3}\\!+\\!i\delta|\Delta|^{(2)}[g^{R(0)}(E_{2})\tau_{2}\\!-\\!\tau_{2}g^{R(0)}(E)],~{}~{}~{}~{}~{}$
(44)

and gives the second-order solution (refer to Appendix B):

$\displaystyle\\!\\!\\!\\!\\!\\!\\!{\delta}g^{R(2)}(E)=i\delta|\Delta|^{(2)}\frac{\tau_{2}\\!-\\!g^{R(0)}(E_{2})\tau_{3}g^{R(0)}(E)}{S^{R}(E_{2})+S^{R}(E)}\\!+\\!\frac{e{\bf
A}_{0}\\!\cdot\\!{\bf v}_{F}}{E_{2}^{2}\\!-\\!E^{2}}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}{\times}\big{[}S^{R}(E_{2})g^{R(0)}(E_{2})\Pi^{R(1)}_{3}\\!-\\!\Pi^{R(1)}_{3}S^{R}(E)g^{R(0)}(E)\big{]}.$
(45)

Considering the response expansions, the normalization condition [Eq. (35)]
for the retarded Green function is written as

$\displaystyle[g^{R(0)}(E)]^{2}=1,$ (46) $\displaystyle
g^{R(0)}(E_{1}){\delta}g^{R(1)}(E)\\!+\\!{\delta}g^{R(1)}(E)g^{R(0)}(E)=0,$
(47) $\displaystyle
g^{R(0)}(E_{2}){\delta}g^{R(2)}(E)\\!+\\!{\delta}g^{R(2)}(E)g^{R(0)}(E)$
$\displaystyle\mbox{}\\!+\\!{\delta}g^{R(1)}(E_{1}){\delta}g^{R(1)}(E)=0.$
(48)

The initial-state $g^{R(0)}(E)$ in Eq. (41) naturally satisfies Eq. (46).
Facilitated with Eq. (46), correspondingly substituting the solutions in Eqs.
(43) and (45), one can easily demonstrate the normalization conditions in Eqs.
(47) and (48). Therefore, as the self-consistent crosscheck, the obtained
solutions of the retarded Green function satisfy the normalization condition.

Further substituting Eq. (41) into Eqs. (42) and (44), the specific solution
of the $\tau_{2}$ component of $\delta{g^{R(2)}}(E)$ is given by (refer to
Appendix B)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!{\delta}g^{R(2)}_{2}(E)\\!=\\!i\delta|\Delta|^{(2)}\Big{[}\frac{4\Delta_{0}^{2}-(2\Omega)^{2}}{{\Gamma}^{(2)}(E)}\\!+\\!\frac{1}{S^{R}(E)}\Big{]}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}+\frac{i(e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F})^{2}\Delta_{0}}{2\Omega^{2}}\Big{[}\frac{1}{S^{R}(E_{2})}\\!+\\!\frac{1}{S^{R}(E)}\\!-\\!\frac{2}{S^{R}(E_{1})}\Big{]},$
(49)

with $\Gamma^{(2)}(E)=2S^{R}(E_{2})S^{R}(E)[S^{R}(E_{2})\\!+\\!S^{R}(E)]$.

It is noted that the involved external optical frequency in the $\tau_{2}$
component of the second order of the retarded Green function in Eq. (IV.1) is
a continuous variable. To further consider the gap dynamics at nonzero
temperature and eliminate the confusion of the auxiliary Matsubara frequency,
we next derive the Keldysh Green function.

###  IV.2 Solution of Keldysh Green function

In this part, we derive the Keldysh Green function. We start with the
normalization condition [Eq. (35)] for the Keldysh Green function:

$g^{R}{\circ}g^{K}+g^{K}{\circ}g^{A}=0.$ (50)

It is established that the Keldysh Green function can be written as a function
of the retarded and advanced ones through a general relationQA1 ; Eilen1 :

$g^{K}=g^{R}\circ{h}-h\circ{g^{A}},$ (51)

where $h(t,t^{\prime})$ denotes the distribution function. Substituting Eq.
(51) to Eq. (50), one finds that the normalization condition for the Keldysh
Green function is immediately satisfied. Consequently, with the obtained
retarded and hence advanced Green function in Sec. IV.1, to solve the Keldysh
Green function, one only needs to solve the distribution function.

In the previous work to derive the Ginzburg-Landau equation from Eilenberger
equation within the Keldysh formalismBa20 , the distribution function
$h(t,t^{\prime})$ is directly taken as the equilibrium one
$\int\frac{dE}{2\pi}h(E)e^{-iE(t-t^{\prime})}$ with the Fourier component
written as

$h(E)=\tanh\Big{(}\frac{\beta{E}}{2}\Big{)}.$ (52)

This treatment usually concerns the case near equilibrium or in strongly-
interacting systems as applied in the transport theory of normal and
superconducting metalsQA1 , and has also been widely used in previous studies
through the Eilenberger equationEilen1 ; Ba8  and diffusive Usadel oneUsadel ;
Usadel1 . In the present work, at clean limit, with a weak external
excitation, we demonstrate Eq. (52) by seriously taking account of the
distribution function.

Specifically, from the expansion of Green function matrices in Eq. (37), the
Keldysh Green function reads

$\displaystyle g^{K}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\int\frac{dE}{2\pi}e^{iE(t^{\prime}-t)}[g^{K(0)}(E)\\!+\\!{\delta}g^{K(1)}(E)e^{-i{\Omega}t}$
(53)
$\displaystyle\mbox{}+{\delta}g^{K(2)}(E){e^{-2i{\Omega}t}}],~{}~{}~{}~{}~{}$

and with the general relation in Eq. (51), one has

$\displaystyle g^{K(0)}(E)\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!g^{R(0)}(E)h^{(0)}(E)-h^{(0)}(E)g^{A(0)}(E),$ (54)
$\displaystyle{\delta}g^{K(1)}(E)\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!g^{R(0)}(E_{1}){\delta}h^{(1)}(E)\\!+\\!{\delta}g^{R(1)}(E)h^{(0)}(E)$
(55)
$\displaystyle\mbox{}\\!\\!\\!-\\!h^{(0)}(E_{1}){\delta}g^{A(1)}(E)\\!-\\!{\delta}h^{(1)}(E)g^{A(0)}(E),~{}~{}~{}~{}~{}$
$\displaystyle{\delta}g^{K(2)}(E)\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!{\delta}g^{R(1)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E_{1}){\delta}g^{A(1)}(E)$
(56)
$\displaystyle\mbox{}\\!\\!\\!+\\!g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!+\\!{\delta}g^{R(2)}(E){h^{(0)}(E)}$
$\displaystyle\mbox{}\\!\\!\\!-\\!h^{(0)}(E_{2}){\delta}g^{A(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E).$

Here, ${\delta}h^{(1)}(E)$ and ${\delta}h^{(2)}(E)$ stand for the first- and
second-order responses on the initial-state $h^{(0)}(E)$, respectively.

According to Eq. (33), the Keldysh Green function satisfies the same equation
as the retarded/advanced one. Consequently, by correspondingly replacing
$g^{R(0)}$ and $\delta{g^{R(1)}}$ by $g^{K(0)}$ and $\delta{g^{K(1)}}$ in Eq.
(42) and then substituting Eq. (55), one finds the equation of the first-order
distribution function (refer to Appendix C):

$\displaystyle\\!\\!\\!\\!\\!\\!\\!(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})[g^{R(0)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E)g^{A(0)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}-[g^{R(0)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E)g^{A(0)}(E)](E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}=e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}[h^{(0)}(E_{1})\\!-\\!h^{(0)}(E)][g^{R(0)}(E_{1})\tau_{3}\\!-\\!\tau_{3}g^{A(0)}(E)].$
(57)

From above equation, the solution of the first-order distribution function
reads (refer to Appendix C)

$\displaystyle\delta{h}^{(1)}(E)$ $\displaystyle=$ $\displaystyle e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F}\frac{h^{(0)}(E_{1})\\!-\\!h^{(0)}(E)}{E_{1}^{2}\\!-\\!E^{2}}[(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{3}$
(58) $\displaystyle\mbox{}+\tau_{3}(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})]$
$\displaystyle=$ $\displaystyle e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}\frac{h^{(0)}(E+\Omega)\\!-\\!h^{(0)}(E)}{\Omega}.$

Similarly, by correspondingly replacing $g^{R(0)}$ and $\delta{g^{R(i=1,2)}}$
by $g^{K(0)}$ and $\delta{g^{K(i=1,2)}}$ in Eq. (44), with Eqs. (55) and (56)
as well as the help of Eq. (58), the equation of the second-order distribution
function reads (refer to Appendix C)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})[g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}-[g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E)](E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}=\\!e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}[g^{R(0)}(E_{2})\tau_{3}\\!-\\!\tau_{3}g^{A(0)}(E)][{\delta}h^{(1)}(E_{1})\\!-\\!{\delta}h^{(1)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\mbox{}+i\delta|\Delta|^{(2)}[h^{(0)}(E_{2})\\!-\\!h^{(0)}(E)][g^{R(0)}(E_{2})\tau_{2}\\!-\\!\tau_{2}g^{A(0)}(E)],$
(59)

from which the solution of the second-order distribution function is obtained
as (refer to Appendix C)

$\displaystyle{\delta}h^{(2)}(E)$ $\displaystyle=$ $\displaystyle(e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F})\frac{{\delta}h^{(1)}(E_{1})\\!-\\!{\delta}h^{(1)}(E)}{2\Omega}$ (60)
$\displaystyle\mbox{}-\delta|\Delta|^{(2)}\frac{\Delta_{0}}{E\\!+\\!\Omega}\frac{h^{(0)}(E_{2})\\!-\\!h^{(0)}(E)}{2\Omega}.$

Then, both the first- and second-order distribution functions are diagonal as
they should be.

The initial-state distribution function can be obtained from the Hamiltonian
in Eq. (2) in self-consistent consideration of the Higgs mode and vector
potential, and is written as

$\displaystyle h^{(0)}(E_{\bf k})$ $\displaystyle=$
$\displaystyle\tanh\Big{\\{}\frac{\beta}{2}\big{[}\sqrt{\xi_{\bf
k}^{2}\\!+\\!(\Delta_{0}\\!+\\!\delta|\Delta|)^{2}}\\!-\\!e{\bf
A}\\!\cdot\\!{\bf v}_{F}\big{]}\Big{\\}}$ (61) $\displaystyle\approx$
$\displaystyle\tanh\Big{\\{}\frac{\beta}{2}\Big{(}E_{\bf
k}+\frac{\Delta_{0}\delta|\Delta|}{E_{\bf k}}-e{\bf A}\\!\cdot\\!{\bf
v}_{F}\Big{)}\Big{\\}},$

where $E_{\bf k}=\sqrt{\xi_{\bf k}^{2}+\Delta_{0}^{2}}$ denotes the Bogoliubov
quasiparticle energy. Following the standard treatment of energy $E=E_{\bf k}$
as the previous work in superconducting stateBo , with the weak excitation, at
low frequency ($\Omega<E=E_{\bf k}$), with Eqs. (58) and (60)-(61), the total
distribution function in relative-frequency space reads

$\displaystyle h(E)$ $\displaystyle=$ $\displaystyle
h^{(0)}(E)+e^{-i\Omega{t}}\delta{h}^{(1)}(E)+e^{-2i\Omega{t}}\delta{h}^{(2)}(E)$
(62) $\displaystyle=$ $\displaystyle\Big{[}1\\!+\\!(e{\bf A}\\!\cdot\\!{\bf
v}_{F})\partial_{E}\\!+\\!\frac{(e{\bf A}\\!\cdot\\!{\bf
v}_{F})^{2}\partial_{E}^{2}}{2}\Big{]}h^{(0)}(E)$
$\displaystyle\mbox{}-\delta|\Delta|\frac{\Delta_{0}}{E}\partial_{E}h^{(0)}(E)\approx\tanh\Big{(}\frac{\beta{E}}{2}\Big{)}.$

Then, the drive effect of vector potential and Higgs-mode part in the initial-
state distribution are exactly canceled by the first- and second-order
distribution functions, leading to the widely applied distribution function
[Eq. (52)] in the literature.

###  IV.3 Higgs-mode generation

In this part, with the obtained distribution function in Eq. (52) and the
second-order retarded Green function in Eq. (IV.1), from the gap equation [Eq.
(36)], the second-order optical response of the Higgs mode at clean limit is
determined by

$\displaystyle\delta|\Delta|^{(2)}$ $\displaystyle=$ $\displaystyle-
iU\int\frac{dE}{2\pi}\langle[h(E)\delta{g}^{R(2)}_{2}(E)-h(E_{2})\delta{g}^{A(2)}_{2}(E)]\rangle_{F}=-iU\int\frac{dE}{2\pi}2h(E)\langle\delta{g}^{R(2)}_{2}(E)\rangle_{F}$
(63) $\displaystyle=$ $\displaystyle
U\int{dE}2h(E)\Big{\\{}\delta|\Delta|^{(2)}\Big{[}\frac{4\Delta_{0}^{2}-(2\Omega)^{2}}{{\Gamma}^{(2)}(E)}\\!+\\!\frac{1}{S^{R}(E)}\Big{]}+\frac{(e{A}_{0}v_{F})^{2}\Delta_{0}}{6\Omega^{2}}\Big{[}\frac{1}{S^{R}(E_{2})}\\!+\\!\frac{1}{S^{R}(E)}\\!-\\!\frac{2}{S^{R}(E_{1})}\Big{]}\Big{\\}}.$

Consequently, one arrives at the equation of motion of the Higgs mode at clean
limit:

$[4\Delta_{0}^{2}-(2\Omega)^{2}]\delta|\Delta|^{(2)}=-\frac{(e{A}_{0}v_{F})^{2}2\Delta_{0}}{3}\frac{\lambda_{E}}{\beta_{E}},$
(64)

similar to the one [Eq. (28)] obtained from the Ginzburg-Landau Lagrangian.
Here, through the standard contour integral, the amplitude-correlation
coefficient reads

$\displaystyle\beta_{E}$ $\displaystyle=$
$\displaystyle-\int\frac{dE}{2\pi}\frac{h(E)}{S^{R}(E_{2})S^{R}(E)[S^{R}(E_{2})\\!+\\!S^{R}(E)]}$
$\displaystyle=$
$\displaystyle\sum_{n>0}\frac{1/(4\beta{i\Omega})}{p_{n}\\!-\\!i\Omega}\Big{[}\frac{1}{\sqrt{(p_{n}\\!-\\!2i\Omega)^{2}\\!+\\!\Delta_{0}^{2}}}\\!-\\!\frac{1}{\sqrt{p_{n}^{2}\\!+\\!\Delta_{0}^{2}}}\Big{]},$

and the essential response coefficient is given by

$\displaystyle\lambda_{E}$ $\displaystyle=$
$\displaystyle\\!\\!\\!-\int\frac{dE}{2\pi}\frac{h(E)}{2\Omega^{2}}\Big{[}\frac{1}{S^{R}(E_{2})}\\!+\\!\frac{1}{S^{R}(E)}\\!-\\!\frac{2}{S^{R}(E_{1})}\Big{]}$
(66) $\displaystyle=$
$\displaystyle\\!\\!\\!\frac{1}{2\beta\Omega^{2}}\sum_{n>0}\Big{[}\frac{2}{\sqrt{(p_{n}\\!-\\!i\Omega)^{2}\\!+\\!\Delta_{0}^{2}}}\\!-\\!\frac{1}{\sqrt{(p_{n}\\!-\\!2i\Omega)^{2}\\!+\\!\Delta_{0}^{2}}}$
$\displaystyle\mbox{}\\!-\\!\frac{1}{\sqrt{p_{n}^{2}\\!+\\!\Delta_{0}^{2}}}\Big{]}.$

It is noted that both amplitude-correlation coefficient $\beta_{E}$ [Eq.
(LABEL:betae)] and response one $\lambda_{E}$ [Eq. (66)] derived here are
exactly same as the ones obtained in the previous workSilaev  by Silaev
through Eilenberger equation within Matsubara formalism. However, in Eqs.
(LABEL:betae) and (66), the fermionic Matsubara frequencies $ip_{n}$ arise
from the singularities in the distribution function $h(E)$ during the standard
contour integral in the complex plane, whereas the involved external optical
frequency $\Omega$ within the Keldysh formalism is always a continuous
variable, in contrast to the treatment of taking optical frequency as bosonic
Matsubara frequency in Ref. Silaev, . As mentioned in the introduction, the
treatment of taking optical frequency as bosonic Matsubara frequency leads to
the vanishing response coefficient $\lambda_{E}$ (i.e., zero Higgs-mode
generation) at all $\Omega\neq 0$, strongly against the Ginzburg-Landau
Lagrangian, whereas the prefactor $1/\Omega^{2}$ in $\lambda_{E}$ causes an
undefined singularity at zero frequency, and hence, an unphysical
discontinuity between cases at $\Omega=0$ and $\Omega\rightarrow 0$. Actually,
even for $\Omega\neq 0$, from Eqs. (LABEL:betae)-(66), near $T_{c}$ with a
weak gap, one has

$\beta_{E}\approx\frac{1}{4\beta\Omega^{2}}\sum_{n>0}\Big{[}\frac{2}{p_{n}\\!-\\!i\Omega}\\!-\\!\frac{1}{p_{n}\\!-\\!2i\Omega}\\!-\\!\frac{1}{{p_{n}}}\Big{]},$
(67)

and

$\lambda_{E}\approx\frac{1}{2\beta\Omega^{2}}\sum_{n>0}\Big{[}\frac{2}{p_{n}\\!-\\!i\Omega}\\!-\\!\frac{1}{p_{n}\\!-\\!2i\Omega}\\!-\\!\frac{1}{{p_{n}}}\Big{]}.$
(68)

In this circumstance, as $\beta_{E}=\lambda_{E}/2$, the Higgs-mode generation
$\delta|\Delta|^{(2)}$, proportional to $\lambda_{E}/\beta_{E}$ from the
equation of motion in Eq. (64), becomes undefined at Matsubara frequency
$i\Omega_{m}$ which leads to $\lambda_{E}=\beta_{E}=0$. This directly poses a
sharp challenge to the study in Ref. Silaev, .

The derivation in the present study, which is performed in the Keldysh
formalism and totally irrelevant to Matsubara space, naturally and
analytically proves the continuous variable of the optical frequency in this
situation. With the continuous optical frequency, near $T_{c}$, from Eqs.
(67)-(68), one finds a finite Higgs-mode generation at all $\Omega$:

$\delta|\Delta|^{(2)}=-\frac{(e{A}_{0}v_{F})^{2}}{3}\frac{4\Delta_{0}}{[4\Delta_{0}^{2}-(2\Omega)^{2}]},$
(69)

which exactly recovers the one [Eq. (28)] derived from the Ginzburg-Landau
Lagrangian. As for the regime with temperature far below $T_{c}$, with the
continuous optical frequency, at low frequency ($\Omega<\Delta_{0}$), one
finds the coefficient
$\beta_{E}=\frac{1}{\beta}\sum_{n>0}({p_{n}^{2}+\Delta_{0}^{2}})^{-3/2}$ and
in particular, a finite response coefficient:

$\lambda_{E}\approx\frac{1}{2\beta}\sum_{n>0}\partial_{p_{n}}^{2}\Big{[}\frac{1}{\sqrt{p_{n}^{2}\\!+\\!\Delta_{0}^{2}}}\Big{]},$
(70)

implying a finite Higgs-mode generation at clean case.

It is also noted that although the Eilenberger equation with the continuous
optical frequency can recover the finite Higgs-mode generation at clean limit
revealed by Ginzburg-Landau Lagrangian and gauge-invariant kinetic
equationGIKE2 ; EPM , this approach fails to derive the Higgs-mode damping by
impurity scattering due to the generically incomplete scattering
integralSilaev . As proved in Ref. ESYW, , because of the quasiclassical
approximation on $\tau_{3}$-Green function, the scattering integral in
Eilenberger equation only involves the anisotropic part of the Green function
that is related to the transport property, but generically drops out the
isotropic one which determines the Higgs-mode lifetime. In this circumstance,
the path-integral approachPYW  and gauge-invariant kinetic equationGIKE3
provide efficient and separate approaches to derive the induced damping of the
Higgs mode by impurities, which agrees with the analysis through Heisenberg
equation of motion as mentioned in the introduction and provides a possible
origin for the experimentally observed broadening of the Higgs-mode resonance
signal as well as the fast Higgs-mode damping after excitationGIKE3 ; PYW .

##  V summary

In summary, we have resolved the current controversy in the literature that
why the previous derivations at clean limit through the path-integral
approachCea1 ; Cea2 ; Cea3  and Eilenberger equationSilaev  within the
Matsubara formalism failed to reach the Higgs-mode generation revealed by
Ginzburg-Landau LagrangianEPM  and gauge-invariant kinetic equationGIKE2 ; EPM
. The crucial treatment leading to this controversy lies at an artificial
scheme within the Matsubara formalism that whether the involved external
optical frequency $\Omega$ in the gap dynamics is taken as continuous variable
or bosonic Matsubara frequency $i\Omega_{m}$. To resolve this confusion, we
derive the effective action of superconducting gap near $T_{c}$ in the
presence of the vector potential through the path-integral approach, and show
that only by taking $\Omega$ as continuous variable within Matsubara
formalism, one can achieve the fundamental Ginzburg-Landau superconducting
Lagrangian in agreement with Landau phase-transition theory and symmetry
analysis. In addition to this physical justification, we also perform the
formulation of the gap dynamics within a separate and rigorous
framework—Keldysh formalism, which is totally irrelevant to Matsubara space.
By applying the Eilenberger equation in Keldysh space to calculate the second-
order response of the Higgs mode, it is analytically proved that the involved
optical frequency is always a continuous variable, leading to finite response
coefficient at clean limit.

Consequently, the present study confirms the unified conclusion, i.e., a
finite Higgs-mode generation at clean limit in the second-order response of
superconductors from three different microscopic approaches (including the
gauge-invariant kinetic equation, Eilenberger equation and path-integral
approach) as well as from Ginzburg-Landau Lagrangian, and can therefore help
understanding the experimental findings of the observed Higgs-mode
excitationNL7 ; NL8 ; NL9 ; NL10 ; NL11 ; DHM2 ; DHM3 .

###### Acknowledgements.

The authors acknowledge financial support from the National Natural Science
Foundation of China under Grants No. 11334014 and No. 61411136001.

##  Appendix A Derivation of correlation coefficients

In this part, we present the specific expressions of the related correlation
coefficients in the superconducting Lagrangian in Eq. (14). Firstly, as the
Green function $G_{0}(p)=\frac{p_{0}+\xi_{\bf k}\tau_{3}}{p_{0}^{2}-\xi_{\bf
k}^{2}}$ only consists of the $\tau_{0}$ and $\tau_{3}$ components, from Eqs.
(15)-(18), one immediately finds
$\chi_{1}=\chi_{13}=\chi_{111}=\chi_{100}=\chi_{1113}=0$. Moreover, within the
Matsubara formalism $[p=(ip_{n},{\bf k})]$, one has

$\displaystyle\chi_{113}$ $\displaystyle=$ $\displaystyle\sum_{p}{\rm
Tr}[G_{0}(p\\!+\\!2q)\tau_{1}G_{0}(p\\!+\\!q)\tau_{1}G_{0}(p)\tau_{3}]=\sum_{p}\frac{2\xi_{\bf
k}[(ip_{n}\\!+\\!2\Omega)(ip_{n}\\!+\\!\Omega)\\!-\\!\xi_{\bf
k}^{2}\\!-\\!ip_{n}\Omega]}{[(ip_{n}\\!+\\!2\Omega)^{2}\\!-\\!\xi_{\bf
k}^{2}][(ip_{n}\\!+\\!\Omega)^{2}\\!-\\!\xi_{\bf
k}^{2}][(ip_{n})^{2}\\!-\\!\xi_{\bf k}^{2}]}=0,~{}~{}~{}~{}~{}$ (71)
$\displaystyle\chi_{p}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{p}{\rm
Tr}[G_{0}(p\\!+\\!q)\tau_{1}G_{0}(p)\tau_{1}]\\!+\\!\frac{1}{U}=\sum_{p}\frac{\\!(ip_{n}\\!+\\!\Omega)^{2}\\!+\\!(ip_{n})^{2}\\!-\\!(ip_{n}\\!+\\!\Omega\\!-\\!ip_{n})^{2}\\!-\\!2\xi_{\bf
k}^{2}}{2[(ip_{n}\\!+\\!\Omega)^{2}\\!-\\!\xi_{\bf
k}^{2}][(ip_{n})^{2}\\!-\\!\xi_{\bf k}^{2}]}\\!+\\!\frac{1}{U}$ (72)
$\displaystyle\approx$
$\displaystyle-\frac{\Omega^{2}}{2}\sum_{p}\frac{1}{[(ip_{n}\\!+\\!\Omega)^{2}\\!-\\!\xi_{\bf
k}^{2}][(ip_{n})^{2}\\!-\\!\xi_{\bf k}^{2}]}\\!+\\!\sum_{\bf
k}\frac{f(\xi_{\bf k})-f(-\xi_{\bf k})}{2\xi_{\bf k}}\\!+\\!\frac{1}{U}$
$\displaystyle=$
$\displaystyle-\frac{\Omega^{2}}{2}\sum_{p}\frac{1}{[(ip_{n}\\!+\\!\Omega)^{2}\\!-\\!\xi_{\bf
k}^{2}][(ip_{n})^{2}\\!-\\!\xi_{\bf
k}^{2}]}\\!+\\!D\int^{\omega_{D}}_{-\omega_{D}}d\xi_{\bf
k}\frac{\tanh({\beta_{c}\xi_{\bf k}}/{2})\\!-\\!\tanh({\beta\xi_{\bf
k}}/{2})}{2\xi_{\bf k}},$ $\displaystyle\chi_{1111}$ $\displaystyle=$
$\displaystyle\sum_{p}{\rm
Tr}[G_{0}(p\\!+\\!q)\tau_{1}G_{0}(p)\tau_{1}G_{0}(p\\!+\\!q)\tau_{1}G(p)\tau_{1}]=\sum_{p}\frac{2}{(ip_{n}\\!+\\!\Omega\\!-\\!\xi_{\bf
k})^{2}(ip_{n}\\!+\\!\xi_{\bf k})^{2}},$ $\displaystyle\chi_{1010}$
$\displaystyle=$ $\displaystyle\sum_{p}{\rm
Tr}[G_{0}(p\\!+\\!q)\tau_{1}G_{0}(p)\tau_{0}G_{0}(p\\!+\\!q)\tau_{1}G(p)\tau_{0}]=\sum_{p}\frac{2}{[(ip_{n}\\!+\\!\Omega)^{2}\\!-\\!\xi^{2}_{\bf
k}][(ip_{n})^{2}\\!-\\!\xi^{2}_{\bf k}]},$ (73)
$\displaystyle\chi_{1100}\\!+\\!\chi_{0110}$ $\displaystyle=$
$\displaystyle\sum_{p}\Big{[}\frac{2}{(ip_{n}\\!+\\!\Omega\\!-\\!\xi_{\bf
k})^{2}[(ip_{n})^{2}\\!-\\!\xi^{2}_{\bf
k}]}+\frac{2}{[(ip_{n}\\!+\\!\Omega)^{2}\\!-\\!\xi^{2}_{\bf
k}](ip_{n}\\!-\\!\xi_{\bf k})^{2}}\Big{]}.$ (74)

Here, we have used the gap equation
$\frac{1}{U}=D\int^{\omega_{D}}_{-\omega_{D}}d\xi_{\bf
k}\frac{\tanh({\beta_{c}\xi_{\bf k}}/{2})}{2\xi_{\bf k}}$ at the critical
temperature in the BCS theoryG1 . It is noted that $\chi_{113}$ vanishes as
the consequence of the particle-hole symmetry, which eliminates the terms with
the odd order of $\xi_{\bf k}$ in the summation of ${\bf k}$.

Then, further using the facts:

$\displaystyle\sum_{p}\frac{1}{[(ip_{n}\\!+\\!\Omega)^{2}\\!-\\!\xi^{2}_{\bf
k}][(ip_{n})^{2}\\!-\\!\xi^{2}_{\bf
k}]}=\frac{D}{\beta}\sum_{ip_{n>0},\eta=\pm}\int{d\xi_{\bf
k}}\frac{1}{[(ip_{n}\\!+\\!\eta\Omega)^{2}\\!-\\!\xi_{\bf
k}^{2}][(ip_{n})^{2}\\!-\\!\xi_{\bf k}^{2}]}$
$\displaystyle\mbox{}=\frac{2\pi{i}D}{\beta\Omega}\sum_{ip_{n>0},\eta=\pm}\Big{[}\frac{1}{(2ip_{n}\\!+\\!2\Omega)(2ip_{n}\\!+\\!\Omega)}\\!-\\!\frac{1}{2ip_{n}(2ip_{n}\\!+\\!\Omega)}\Big{]}=\frac{\pi{i}D}{\beta\Omega^{2}}\sum_{ip_{n>0},\eta=\pm}\Big{[}\frac{4}{2ip_{n}\\!+\\!\Omega}\\!-\\!\frac{1}{ip_{n}}\\!-\\!\frac{1}{ip_{n}\\!+\\!\Omega}\Big{]},$
(75) $\displaystyle\sum_{p}\frac{2}{(ip_{n}\\!+\\!\Omega\\!-\\!\xi_{\bf
k})^{2}(ip_{n}\\!+\\!\xi_{\bf
k})^{2}}=\frac{D}{\beta}\sum_{ip_{n>0},\eta=\pm}\int{d\xi_{\bf
k}}\frac{2}{(ip_{n}\\!+\\!\eta\Omega\\!-\\!\xi_{\bf
k})^{2}(ip_{n}\\!+\\!\xi_{\bf
k})^{2}}=-\sum_{ip_{n>0},\eta=\pm}\frac{8\pi{i}D/\beta}{(2ip_{n}\\!+\\!\eta\Omega)^{3}},$
(76) $\displaystyle\sum_{p}\Big{[}\frac{2}{(ip_{n}\\!+\\!\Omega\\!-\\!\xi_{\bf
k})^{2}[(ip_{n})^{2}\\!-\\!\xi^{2}_{\bf
k}]}+\frac{2}{[(ip_{n}\\!+\\!\Omega)^{2}\\!-\\!\xi^{2}_{\bf
k}](ip_{n}\\!-\\!\xi_{\bf k})^{2}}\Big{]}$
$\displaystyle\mbox{}=\frac{D}{\beta}\sum_{ip_{n>0},\eta=\pm}\int{d\xi_{\bf
k}}\Big{[}\frac{2}{(ip_{n}\\!+\\!\eta\Omega\\!-\\!\xi_{\bf
k})^{2}[(ip_{n})^{2}\\!-\\!\xi^{2}_{\bf
k}]}+\frac{2}{[(ip_{n}\\!+\\!\eta\Omega)^{2}\\!-\\!\xi^{2}_{\bf
k}](ip_{n}\\!-\\!\xi_{\bf k})^{2}}\Big{]}$
$\displaystyle\mbox{}=\frac{2\pi{i}D}{\beta\Omega^{2}}\sum_{ip_{n>0},\eta=\pm}\Big{[}\frac{4}{2ip_{n}\\!+\\!\eta\Omega}\\!-\\!\frac{1}{ip_{n}}\\!-\\!\frac{1}{ip_{n}\\!+\\!\eta\Omega}\Big{]},$
(77)

the Landau parameters $\beta_{p}$ [Eq. (21)], $\lambda_{p}$ [Eq. (22)],
$\alpha_{p}$ [Eq. (23)] and $\gamma_{p}$ [Eq. (24)] are derived.

##  Appendix B Derivation of retarded Green function from Eilenberger equation

In this part, we present the derivation of the retarded Green function from
the Eilenberger equation. With the initial state of the retarded Green
function in Eq. (41), the equation of the first order of retarded Green
function in Eq. (42) is re-written as

$S^{R}(E_{1})g^{R(0)}(E_{1}){\delta}g^{R(1)}(E)\\!-\\!{\delta}g^{R(1)}(E)S^{R}(E)g^{R(0)}(E)=e{\bf
A}_{0}\cdot{\bf v}_{F}\Pi^{R(0)}_{3}.$ (78)

From above equation, considering the normalization condition of $g^{R(0)}(E)$
in Eq. (46), one has

$\displaystyle[S^{R}(E_{1})]^{2}{\delta}g^{R(1)}(E)\\!-\\!S^{R}(E_{1})S^{R}(E)g^{R(0)}(E_{1}){\delta}g^{R(1)}(E)g^{R(0)}(E)=S^{R}(E_{1})e{\bf
A}_{0}\cdot{\bf v}_{F}[\tau_{3}\\!-\\!g^{R(0)}(E_{1})\tau_{3}g^{R(0)}(E)],$
(79) $\displaystyle
S^{R}(E)S^{R}(E_{1})g^{R(0)}(E_{1}){\delta}g^{R(1)}(E)g^{R(0)}(E)\\!-\\![S^{R}(E)]^{2}{\delta}g^{R(1)}(E)=S^{R}(E)e{\bf
A}_{0}\cdot{\bf v}_{F}[g^{R(0)}(E_{1})\tau_{3}g^{R(0)}(E)\\!-\\!\tau_{3}].$
(80)

Then, the solution of ${\delta}g^{R(1)}(E)$ in Eq. (43) can be easily obtained
by adding Eqs. (79) and (80). Moreover, substituting Eq. (41) into Eq. (43),
the specific expression of ${\delta}g^{R(1)}(E)$ is given by

${\delta}g^{R(1)}(E)=(e{\bf A}_{0}\cdot{\bf
v}_{F})\frac{\tau_{3}-[S^{R}(E_{1})S^{R}(E)]^{-1}[E_{1}E\tau_{3}+i\Delta_{0}(E_{0}+E_{1})\tau_{2}+\Delta_{0}^{2}\tau_{3}]}{S^{R}(E_{1})+S^{R}(E)}.$
(81)

Similarly, the equation of the second order of retarded Green function in Eq.
(44) is re-written as

$S^{R}(E_{2})g^{R(0)}(E_{2}){\delta}g^{R(2)}(E)\\!-\\!{\delta}g^{R(2)}(E)S^{R}(E)g^{R(0)}(E)=e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F}\Pi^{R(1)}_{3}\\!+\\!i\delta|\Delta|^{(2)}[g^{R(0)}(E_{2})\tau_{2}\\!-\\!\tau_{2}g^{R(0)}(E)],$
(82)

and using the normalization condition of $g^{R(0)}(E)$ in Eq. (46), one easily
gets the solution of ${\delta}g^{R(1)}(E)$ in Eq. (45). Substituting Eqs. (41)
and (81) into Eq. (45), the specific expression of the $\tau_{2}$ component of
${\delta}g^{R(2)}(E)$ reads

$\displaystyle{\delta}g^{R(2)}_{2}(E)$ $\displaystyle=$ $\displaystyle{i(e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F})^{2}\Delta_{0}}\Big{\\{}\frac{(EE_{2}\\!+\\!E_{1}E\\!+E_{1}E_{2}\\!+\\!\Delta_{0}^{2})[S^{R}(E_{1})\\!+\\!S^{R}(E_{2})\\!+\\!S^{R}(E)]}{S^{R}(E_{1})S^{R}(E_{2})S^{R}(E)}-1\Big{\\}}\frac{1}{[S^{R}(E)\\!+\\!S^{R}(E_{2})]}$
(83)
$\displaystyle\mbox{}\times\frac{1}{[S^{R}(E_{1})\\!+\\!S^{R}(E_{2})][S^{R}(E_{1})\\!+\\!S^{R}(E)]}\\!+\\!i\delta|\Delta|^{(2)}\frac{2\Delta_{0}^{2}\\!+\\!2EE_{2}\\!+\\!2S^{R}(E)S^{R}(E_{2})}{\Gamma^{(2)}(E)}.$

Further considering

$\frac{[S^{R}(E_{1})\\!+\\!S^{R}(E_{2})\\!+\\!S^{R}(E)][S^{R}(E)\\!-\\!S^{R}(E_{1})][S^{R}(E_{1})\\!-\\!S^{R}(E_{2})][S^{R}(E)\\!-\\!S^{R}(E_{2})]}{S^{R}(E_{1})S^{R}(E_{2})S^{R}(E)}=\frac{E_{1}^{2}\\!-\\!E_{2}^{2}}{S^{R}(E)}\\!+\\!\frac{E_{2}^{2}\\!-\\!E^{2}}{S^{R}(E_{1})}\\!+\\!\frac{E^{2}\\!-\\!E_{1}^{2}}{S^{R}(E_{2})},~{}$
(84)

one has

$\displaystyle{\delta}g^{R(2)}_{2}(E)$ $\displaystyle=$
$\displaystyle\frac{i(e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F})^{2}\Delta_{0}}{(E^{2}_{1}-E^{2})(E_{2}^{2}-E^{2})(E_{1}^{2}-E_{2}^{2})}\Big{\\{}{(EE_{2}\\!+\\!E_{1}E\\!+E_{1}E_{2}\\!+\\!\Delta_{0}^{2})}{}\Big{[}\frac{E_{1}^{2}\\!-\\!E_{2}^{2}}{S^{R}(E)}\\!+\\!\frac{E_{2}^{2}\\!-\\!E^{2}}{S^{R}(E_{1})}\\!+\\!\frac{E^{2}\\!-\\!E_{1}^{2}}{S^{R}(E_{2})}\Big{]}$
(85)
$\displaystyle\mbox{}-{[S^{R}(E)\\!-\\!S^{R}(E_{2})][S^{R}(E_{1})\\!-\\!S^{R}(E_{2})][S^{R}(E)\\!-\\!S^{R}(E_{1})]}\Big{\\}}+i\delta|\Delta|^{(2)}\frac{2\Delta_{0}^{2}\\!+\\!2EE_{2}\\!+\\!2S^{R}(E)S^{R}(E_{2})}{\Gamma^{(2)}(E)}$
$\displaystyle=$ $\displaystyle\frac{i(e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F})^{2}\Delta_{0}}{(E^{2}_{1}\\!-\\!E^{2})(E_{2}^{2}\\!-\\!E^{2})(E_{1}^{2}\\!-\\!E_{2}^{2})}\Big{\\{}(E_{1}^{2}\\!-\\!E_{2}^{2})\frac{(E\\!+\\!E_{2})(E_{1}\\!+\\!E)\\!-\\![S^{R}(E)]^{2}}{S^{R}(E)}\\!+\\!\frac{(E\\!+\\!E_{1})(E_{1}\\!+\\!E_{2})\\!-\\![S^{R}(E_{1})]^{2}}{S^{R}(E_{1})}$
$\displaystyle\mbox{}\times(E_{2}^{2}\\!-\\!E^{2})\\!+\\!(E^{2}\\!-\\!E_{1}^{2})\frac{(E\\!+\\!E_{2})(E_{2}\\!+\\!E_{1})\\!-\\![S^{R}(E_{2})]^{2}}{S^{R}(E_{2})}\\!+\\!S^{R}(E)\big{\\{}[S^{R}(E_{1})]^{2}\\!-\\![S^{R}(E_{2})]^{2}\big{\\}}\\!+\\!\big{\\{}[S^{R}(E_{2})]^{2}$
$\displaystyle\mbox{}-[S^{R}(E)]^{2}\big{\\}}S^{R}(E_{1})\\!+\\!\\!S^{R}(E_{2})\big{\\{}[S^{R}(E)]^{2}\\!-\\![S^{R}(E_{1})]^{2}\big{\\}}\Big{\\}}\\!+\\!i\delta|\Delta|^{(2)}\frac{4\Delta_{0}^{2}\\!-\\!(E\\!-\\!E_{2})^{2}\\!+\\![S^{R}(E)\\!+\\!S^{R}(E_{2})]^{2}}{\Gamma^{(2)}(E)}$
$\displaystyle=$ $\displaystyle\frac{i(e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F})^{2}\Delta_{0}}{2\Omega^{2}}\Big{[}\frac{1}{S^{R}(E_{2})}\\!+\\!\frac{1}{S^{R}(E)}\\!-\\!\frac{2}{S^{R}(E_{1})}\Big{]}+i\delta|\Delta|^{(2)}\Big{[}\frac{4\Delta_{0}^{2}-(2\Omega)^{2}}{{\Gamma}^{(2)}(E)}\\!+\\!\frac{1}{S^{R}(E)}\Big{]}.$

Then, Eq. (IV.1) is derived.

##  Appendix C Derivation of distribution function from Eilenberger equation

In this part, we present the derivation of the distribution function from the
Eilenberger equation. As mentioned in Sec. IV.2, according to Eq. (33), the
Keldysh Green function satisfies the same equation as the retarded/advanced
one. Therefore, the equation of the first order of the Keldysh Green function
reads

$(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{K(1)}(E)\\!-\\!{\delta}g^{K(1)}(E)(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})=e{\bf
A}_{0}\cdot{\bf v}_{F}\Pi^{K(0)}_{3}.$ (86)

Substituting $g^{K(0)}$ [Eq. (54)] and ${\delta}g^{K(1)}(E)$ [Eq. (55)], the
above equation becomes

$\displaystyle\\!\\!\\!\\!\\!\\!\\!(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})[g^{R(0)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E)g^{A(0)}(E)]\\!-\\![g^{R(0)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E)g^{A(0)}(E)](E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\\!=\\!h^{(0)}(E)$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}\times[{\delta}g^{R(1)}(E)(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\\!-\\!(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{R(1)}(E)]\\!-\\!h^{(0)}(E_{1})[{\delta}g^{A(1)}(E)(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\\!-\\!(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{A(1)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}-e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}[h^{(0)}(E)\tau_{3}g^{R(0)}(E)\\!-\\!g^{R(0)}(E_{1})h^{(0)}(E_{1})\tau_{3}\\!-\\!h^{(0)}(E)\tau_{3}g^{A(0)}(E)\\!+\\!g^{A(0)}(E_{1})h^{(0)}(E_{1})\tau_{3}].~{}~{}~{}~{}$
(87)

Then, facilitated with the equations of the first order of the retarded [Eq.
(42)] and advanced Green functions, one arrives at the equation of the first
order of the distribution function in Eq. (57). By first multiplying Eq. (57)
by $(E_{1}\tau_{3}+i\Delta_{0}\tau_{2})$ from the left side and
$(E\tau_{3}+i\Delta_{0}\tau_{2})$ from the right side respectively and adding
the obtained two equations afterwards, one has

$\displaystyle[(E_{1}^{2}\\!+\\!\Delta_{0}^{2})\\!-\\!(E^{2}\\!+\\!\Delta_{0}^{2})][g^{R(0)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E)g^{A(0)}(E)]=e{\bf
A}_{0}\\!\cdot\\!{\bf v}_{F}[h^{(0)}(E_{1})\\!-\\!h^{(0)}(E)]$
$\displaystyle\mbox{}\times\big{\\{}g^{R(0)}(E_{1})[(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{3}\\!+\\!\tau_{3}(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})]\\!-\\![(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{3}\\!+\\!\tau_{3}(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})]g^{A(0)}(E)\big{\\}},$
(88)

from which the solution of the first order of the distribution function reads

${\delta}h^{(1)}(E)=e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}[h^{(0)}(E_{1})\\!-\\!h^{(0)}(E)]\frac{[(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{3}\\!+\\!\tau_{3}(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})]}{(E_{1}^{2}\\!+\\!\Delta_{0}^{2})\\!-\\!(E^{2}\\!+\\!\Delta_{0}^{2})},$
(89)

and then, Eq. (58) is derived.

Similarly, the equation of the second order of the Keldysh Green function
reads

$(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{K(2)}(E)\\!-\\!{\delta}g^{K(2)}(E)(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})=e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F}\Pi^{K(1)}_{3}\\!+\\!i\delta|\Delta|^{(2)}[g^{K(0)}(E_{2})\tau_{2}\\!-\\!\tau_{2}g^{K(0)}(E)].$
(90)

Substituting $g^{K(0)}$ [Eq. (54)] and ${\delta}g^{K(1)}(E)$ [Eq. (55)] as
well as ${\delta}g^{K(2)}(E)$ [Eq. (56)], the above equation becomes

$\displaystyle\\!\\!\\!\\!\\!\\!\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})[g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E)]\\!-\\![g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E)](E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!=[{\delta}g^{R(1)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E_{1}){\delta}g^{A(1)}(E)](E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\\!-\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})[{\delta}g^{R(1)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E_{1}){\delta}g^{A(1)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}+h^{(0)}(E)[{\delta}g^{R(2)}(E)(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\\!-\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{R(2)}(E)]\\!-\\![{\delta}g^{A(2)}(E)(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\\!-\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{A(2)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}{\times}h^{(0)}(E_{2})\\!+\\!e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F}[{\delta}g^{R(1)}(E_{1})h^{(0)}(E_{1})\tau_{3}\\!-\\!\tau_{3}{\delta}g^{R(1)}(E)h^{(0)}(E)\\!+\\!g^{R(0)}(E_{2}){\delta}h^{(1)}(E_{1})\tau_{3}\\!-\\!\tau_{3}g^{R(0)}(E_{1}){\delta}h^{(1)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}-e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}[h^{(0)}(E_{2}){\delta}g^{A(1)}(E_{1})\tau_{3}\\!-\\!\tau_{3}h^{(0)}(E_{1}){\delta}g^{A(1)}(E)\\!+\\!{\delta}h^{(1)}(E_{1})g^{A(0)}(E_{1})\tau_{3}\\!-\\!\tau_{3}{\delta}h^{(1)}(E)g^{A(0)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}-i\delta|\Delta|^{(2)}[\tau_{2}g^{R(0)}(E)h^{(0)}(E)\\!-\\!g^{R(0)}(E_{2})h^{(0)}(E_{2})\tau_{2}\\!-\\!\tau_{2}g^{A(0)}(E)h^{(0)}(E)\\!+\\!g^{A(0)}(E_{2})h^{(0)}(E_{2})\tau_{2}].$
(91)

Facilitated with the equations of the second order of the retarded [Eq. (44)]
and advanced Green functions, one has

$\displaystyle\\!\\!\\!\\!\\!\\!\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})[g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E)]\\!-\\![g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E)](E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!=[{\delta}g^{R(1)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E_{1}){\delta}g^{A(1)}(E)](E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\\!-\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})[{\delta}g^{R(1)}(E_{1}){\delta}h^{(1)}(E)\\!-\\!{\delta}h^{(1)}(E_{1}){\delta}g^{A(1)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}+e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}[h^{(0)}(E_{1})-h^{(0)}(E)]{\delta}g^{R(1)}(E_{1})\tau_{3}\\!+\\!e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F}[g^{R(0)}(E_{2}){\delta}h^{(1)}(E_{1})\tau_{3}\\!-\\!\tau_{3}g^{R(0)}(E_{1}){\delta}h^{(1)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}-e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}[h^{(0)}(E_{2})\\!-\\!h^{(0)}(E_{1})]\tau_{3}{\delta}g^{A(1)}(E)\\!-\\!e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F}[{\delta}h^{(1)}(E_{1})g^{A(0)}(E_{1})\tau_{3}\\!-\\!\tau_{3}{\delta}h^{(1)}(E)g^{A(0)}(E)]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}+i\delta|\Delta|^{(2)}[h^{(0)}(E_{2})\\!-\\!h^{(0)}(E)][g^{R(0)}(E_{2})\tau_{2}\\!-\\!\tau_{2}g^{A(0)}(E)].$
(92)

Further using Eq. (58) to replace $e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}[h^{(0)}(E_{1})-h^{(0)}(E)]$ with $\Omega\delta{h^{1}}(E)$, the above
equation is simplified as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})[g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E)]\\!-\\![g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E)](E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!={\delta}h^{(1)}(E)[{\delta}g^{R(1)}(E_{1})(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\\!-\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{R(1)}(E_{1})\\!-\\!e{\bf
A}_{0}\\!\cdot\\!{\bf v}_{F}\tau_{3}g^{R(0)}(E_{1})]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}-{\delta}h^{(1)}(E_{1})[{\delta}g^{A(1)}(E_{1})(E_{1}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\\!-\\!(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2}){\delta}g^{A(1)}(E_{1})\\!+\\!e{\bf
A}_{0}\\!\cdot\\!{\bf v}_{F}g^{A(0)}(E_{1})\tau_{3}]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mbox{}+e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}[g^{R(0)}(E_{2}){\delta}h^{(1)}(E_{1})\tau_{3}\\!+\\!\tau_{3}{\delta}h^{(1)}(E)g^{A(0)}(E)]\\!+\\!i\delta|\Delta|^{(2)}[h^{(0)}(E_{2})\\!-\\!h^{(0)}(E)][g^{R(0)}(E_{2})\tau_{2}\\!-\\!\tau_{2}g^{A(0)}(E)].~{}~{}~{}~{}$
(93)

Consequently, substituting equations of the first order of the retarded [Eq.
(42)] and advanced Green functions to above equation, one arrives at the
equation of the second order of the distribution function in Eq. (59). By
first multiplying Eq. (57) by $(E_{2}\tau_{3}+i\Delta_{0}\tau_{2})$ from the
left side and $(E\tau_{3}+i\Delta_{0}\tau_{2})$ from the right side
respectively and adding the obtained two equations afterwards, one has

$\displaystyle\\!\\!\\!\\!\\!\\![g^{R(0)}(E_{2}){\delta}h^{(2)}(E)\\!-\\!{\delta}h^{(2)}(E)g^{A(0)}(E)]\\!=\\!\frac{e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F}[{\delta}h^{(1)}(E_{1})\\!-\\!{\delta}h^{(1)}(E)]}{E_{2}^{2}\\!-\\!E^{2}}\big{\\{}g^{R(0)}(E_{2})[(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{3}\\!+\\!\tau_{3}(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})]$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}-[(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{3}\\!+\\!\tau_{3}(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})]g^{A(0)}(E)\big{\\}}\\!+\\!i\delta|\Delta|^{(2)}\frac{h^{(0)}(E_{2})\\!-\\!h^{(0)}(E)}{E_{2}^{2}\\!-\\!E^{2}}\big{\\{}g^{R(0)}(E_{2})[(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{2}\\!+\\!\tau_{2}(E\tau_{3}$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}+i\Delta_{0}\tau_{2})]\\!-\\![(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{2}\\!+\\!\tau_{2}(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})]g^{A(0)}(E)\big{\\}},$
(94)

from which the solution of the first order of the distribution function reads

$\displaystyle{\delta}h^{(2)}(E)$ $\displaystyle=$ $\displaystyle e{\bf
A}_{0}\\!\cdot\\!{\bf
v}_{F}\frac{[{\delta}h^{(1)}(E_{1})\\!-\\!{\delta}h^{(1)}(E)][(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{3}\\!+\\!\tau_{3}(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})]}{E_{2}^{2}\\!-\\!E^{2}}$
(95)
$\displaystyle+i\delta|\Delta|^{(2)}\frac{[h^{(0)}(E_{2})\\!-\\!h^{(0)}(E)][(E_{2}\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})\tau_{2}\\!+\\!\tau_{2}(E\tau_{3}\\!+\\!i\Delta_{0}\tau_{2})]}{E_{2}^{2}\\!-\\!E^{2}}$
$\displaystyle=$ $\displaystyle e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}\frac{{\delta}h^{(1)}(E_{1})\\!-\\!{\delta}h^{(1)}(E)}{E_{2}-E}\\!+\\!i\delta|\Delta|^{(2)}\frac{[h^{(0)}(E_{2})\\!-\\!h^{(0)}(E)](2i\Delta_{0}\\!-\\!2i\Omega{\tau_{1}})}{E_{2}^{2}\\!-\\!E^{2}}$
$\displaystyle\approx$ $\displaystyle e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F}\frac{{\delta}h^{(1)}(E_{1})\\!-\\!{\delta}h^{(1)}(E)}{E_{2}-E}\\!+\\!i\delta|\Delta|^{(2)}2i\Delta_{0}\frac{h^{(0)}(E_{2})\\!-\\!h^{(0)}(E)}{E_{2}^{2}\\!-\\!E^{2}},$

and then, Eq. (60) is derived.

## References

  * (1) P. B. Littlewood and C. M. Varma, Phys. Rev. Lett. 47, 811 (1981); Phys. Rev. B 26, 4883 (1982).
  * (2) A. F. Volkov and S. M. Kogan, Zh. Eksp. Teor. Fiz 65, 2038 (1974) [Sov. Phys. JETP 38, 1018 (1974)].
  * (3) E. A. Yuzbashyan and M. Dzero, Phys. Rev. Lett 96, 230404 (2006).
  * (4) V. Gurarie, Phys. Rev. Lett. 103, 075301 (2009).
  * (5) A. Moor, P. A. Volkov, A. F. Volkov, and K. B. Efetov, Phys. Rev. B 90, 024511 (2014).
  * (6) D. Pekker and C. Varma, Annu. Rev. Condens. Matter Phys. 6, 269 (2015).
  * (7) N. Tsuji, Y. Murakami, and H. Aoki, Phys. Rev. B 94, 224519 (2016).
  * (8) T. Yanagisawa, Commun. Comput. Phys. 23, 459 (2017).
  * (9) Z. Y. Sun, M. M. Fogler, D. N. Basov, and Andrew J. Millis, Phys. Rev. Research 2, 023413 (2020).
  * (10) R. Matsunaga and R. Shimano, Phys. Rev. Lett. 109, 187002 (2012).
  * (11) R. Matsunaga, Y. I. Hamada, K. Makise, Y. Uzawa, H. Terai, Z. Wang, and R. Shimano, Phys. Rev. Lett. 111, 057002 (2013).
  * (12) 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).
  * (13) R. Matsunaga, N. Tsuji, K. Makise, H. Terai, H. Aoki, and R. Shimano, Phys. Rev. B 96, 020505 (2017).
  * (14) K. Katsumi, N. Tsuji, Y. I. Hamada, R. Matsunaga, J. Schneeloch, R. D. Zhong, G. D. Gu, H. Aoki, Y. Gallais, and R. Shimano, Phys. Rev. Lett. 120, 117001 (2018).
  * (15) H. Chu, M. J. Kim, K. Katsumi, S. Kovalev, R. D. Dawson, L. Schwarz, N. Yoshikawa, G. Kim, D. Putzky, Z. Z. Li, H. Raffy, S. Germanskiy, J. C. Deinert, N. Awari, I. Ilyakov, B. Green, M. Chen, M. Bawatna, G. Christiani, G. Logvenov, Y. Gallais, A. V. Boris, B. Keimer, A. Schnyder, D. Manske, M. Gensch, Z. Wang, R. Shimano, and S. Kaiser, Nat. Commun. 11, 1793 (2020).
  * (16) K. Katsumi, Z. Z. Li, H. Raffy, Y. Gallais, R. Shimano, Phys. Rev. B 102, 054510 (2020).
  * (17) J. R. Schrieffer, Theory of Superconductivity (W. A. Benjamin, New York, 1964).
  * (18) A. A. Abrikosov, L. P. Gor’kov, and I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics (Prentice Hall, Englewood Cliffs, 1963).
  * (19) F. Yang and M. W. Wu, Phys. Rev. B 100, 104513 (2019).
  * (20) R. A. Barankov, L. S. Levitov, and B. Z. Spivak, Phys. Rev. Lett. 93, 160401 (2004).
  * (21) R. A. Barankov and L. S. Levitov, Phys. Rev. Lett. 96, 230403 (2006).
  * (22) N. Tsuji and H. Aoki, Phys. Rev. B 92, 064508 (2015).
  * (23) M. Dzero, M. Khodas, and A. Levchenko, Phys. Rev. B 91, 214505 (2015).
  * (24) M. Lu, H. W. Liu, P. Wang, and X. C. Xie, Phys. Rev. B 93, 064516 (2016).
  * (25) Y. Murotani, N. Tsuji, and H. Aoki, Phys. Rev. B 95, 104503 (2017).
  * (26) L. Schwarz and D. Manske, Phys. Rev. B 101, 184519 (2020).
  * (27) T. Papenkort, V. M. Axt, and T. Kuhn, Phys. Rev. B 76, 224522 (2007).
  * (28) T. Papenkort, T. Kuhn, and V. M. Axt, Phys. Rev. B 78, 132505 (2008).
  * (29) A. F. Kemper, M. A. Sentef, B. Moritz, J. K. Freericks, and T. P. Devereaux, Phys. Rev. B 92, 224517 (2015).
  * (30) H. Krull, N. Bittner, G. S. Uhrig, D. Manske, and A. P. Schnyder, Nat. Commun. 7, 11921 (2016).
  * (31) G. Seibold, M. Udina, C. Castellani, and L. Benfatto, Phys. Rev. B 103, 014512 (2017).
  * (32) P. W. Anderson, Phys. Rev. 112, 1900 (1958).
  * (33) S. Tsuchiya, D. Yamamoto, R. Yoshii, and M. Nitta, Phys. Rev. B 98, 094503 (2018).
  * (34) T. Cea, C. Castellani, and L. Benfatto Phys. Rev. B 93, 180507(R) (2016).
  * (35) T. Cea and L. Benfatto, Phys. Rev. B 94, 064512 (2016).
  * (36) T. Cea, P. Barone, C. Castellani, and L. Benfatto, Phys. Rev. B 97, 094516 (2018).
  * (37) K. Isoyama, N. Yoshikawa, K. Katsumi, J. Wong, N. Shikama, Y. Sakishita, F. Nabeshima, A. Maeda, and R. Shimano, Commun. Phys. 4, 160 (2021).
  * (38) Y. Murotani and R. Shimano, Phys. Rev. B 99, 224510 (2019).
  * (39) N. Tsuji and Y. Nomura, Phys. Rev. Research 2, 043029 (2020).
  * (40) R. Haenel, P. Froese, D. Manske, and L. Schwarz, Phys. Rev. B 104, 134504 (2021).
  * (41) M. Silaev, Phys. Rev. B 99, 224511 (2019).
  * (42) G. Eilenberger, Z. Phys. 214, 195 (1968).
  * (43) To justify the result from the Eilenberger equation, Ref. Silaev,  also performed the diagrammatic formulation. Nevertheless, as referred to Sec. IV A in Ref. Silaev, , a special procedure that follows the exact step to derive the Eilenberger equation by first taking the commutation between the self-energy and $\tau_{3}$-Green function and applying the quasiclassical approximation afterwards, was applied to handle the calculation within the diagrammatic formalism. Consequently, the applied diagrammatic formalism approach with this procedure becomes exactly same as the Eilenberger equation and is not a separate nontrivial techinique.
  * (44) F. Yang and M. W. Wu, Phys. Rev. B 98, 094507 (2018).
  * (45) F. Yang and M. W. Wu, Phys. Rev. B 102, 144508 (2020).
  * (46) F. Yang and M. W. Wu, arXiv:2111.06128.
  * (47) A. Tagliacozzo and F. Ventriglia, Il Nuovo Cimento D 11, 141 (1989).
  * (48) T. Kita, Statistical Mechanics of Superconductivity (Springer, Berlin, 2015).
  * (49) J. Rammer and H. Smith, Rev. Mod. Phys. 58, 323 (1986).
  * (50) G. D. Mahan, Many Particle Physics (Plenum, New York, 1990).
  * (51) V. Ambegaokar and L. P. Kadanoff, Nuovo Cimento 22, 914 (1961).
  * (52) Y. Nambu, Phys. Rev. 117, 648 (1960).
  * (53) Y. Nambu, Rev. Mod. Phys. 81, 1015 (2009).
  * (54) H. A. Fertig and S. D. Sarma, Phys. Rev. Lett. 65, 1482 (1990).
  * (55) I. J. R. Aitchison, P. Ao, D. J. Thouless, and X. M. Zhu, Phys. Rev. B 51, 6531 (1995).
  * (56) I. J. R. Aitchison, G. Metikas, and D. J. Lee, Phys. Rev. B 62, 6638 (2000).
  * (57) P. W. Anderson, Phys. Rev. 130, 439 (1963).
  * (58) H. G. Hugdal, J. Linder, and S. H. Jacobsen, Phys. Rev. B 95, 235403 (2017).
  * (59) A. I. Buzdin, Rev. Mod. Phys. 77, 935 (2005).
  * (60) K. D. Usadel, Phys. Rev. Lett. 25, 507 (1970).
  * (61) C. Espeday, T. Yokoyama, and J. Linder, Phys. Rev. Lett. 116, 127002 (2016).
  * (62) S. Li, A. V. Andreev, and B. Z. Spivak, Phys. Rev. B 92, 100506 (2015).
  * (63) F. Yang and M. W. Wu, arXiv:1912.09172.
  * (64) F. Yang and M. W. Wu, Phys. Rev. B 106, 144509 (2022).