# Optical response of Higgs mode in superconductors at clean limit

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

The phenomenological Ginzburg-Landau theory and the charge conservation
directly lead to the finite Higgs-mode generation and vanishing charge-density
fluctuation in the second-order optical response of superconductors at clean
limit. Nevertheless, recent microscopic theoretical studies of the second-
order optical response, apart from the one through the gauge-invariant kinetic
equation [Yang and Wu, Phys. Rev. B 100, 104513 (2019)], have derived a
vanishing Higgs-mode generation but finite charge-density fluctuation at clean
limit. We resolve this controversy by re-examining the previous derivations
with the vector potential alone within the path-integral and Eilenberger-
equation approaches, and show that both previous derivations contain flaws.
After fixing these flaws, a finite Higgs-mode generation through the drive
effect of vector potential is derived at clean limit, exactly recovering the
previous result from the gauge-invariant kinetic equation as well as Ginzburg-
Landau theory. By further extending the path-integral approach to include
electromagnetic effects from the scalar potential and phase mode, in the
second-order response, a finite contribution from the drive effect of scalar
potential to the Higgs-mode generation at clean limit as well as the vanishing
charge-density fluctuation are derived, also recovering the results from the
gauge-invariant kinetic equation. Particularly, we show that the phase mode is
excited in the second-order response, and exactly cancels the previously
reported unphysical excitation of the charge-density fluctuation, guaranteeing
the charge conservation.

###### pacs:

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

##  I Introduction

Recently, the Higgs mode in superconductors has attracted extensive
experimental and theoretical interest. Specifically, the angular and radial
excitations that emerge in the Mexican-hat free energy of superconductors,
describe the phase and amplitude fluctuations of the superconducting order
parameter $\Delta$Am0 , respectively. The gapless phase mode $\delta\theta$,
referred to as Nambu-Goldstone modegi0 ; AK ; Gm1 ; Gm2 ; Ba0 ; pm0 ; pi1 ;
pm1 ; pm2 ; pi4 ; gi1 ; Ba9 ; Ba10 ; pm5 ; AK2 , corresponds to gapless
Goldstone boson due to the spontaneous breaking of continuous $U(1)$
symmetryGm1 ; Gm2 . Whereas the amplitude mode, exhibiting a gapful energy
spectrum $\omega_{H}$, is referred to as the Higgs modeOD1 ; OD2 ; OD3 ; pm5 ;
Am0 ; Am5 ; Am6 ; Am12 ; AK2 , because of the similarity to Higgs boson in the
field theoryHiggs1 ; Higgs2 ; Higgs3 . Early theoretical works in conventional
$s$-wave superconductors have reported $\omega_{H}=2|\Delta_{0}|$ at long-wave
limitpm5 ; Am0 ; OD1 ; OD2 ; OD3 ; Am6 ; AK2 , with $|\Delta_{0}|$ being the
superconducting gap. Nevertheless, being charge neutral and spinless, the
Higgs mode has long been experimentally elusive. Until recently, thanks to the
advanced technique in ultrafast nonlinear optics, it is experimentally
realizedNL7 ; NL8 ; NL9 ; NL10 ; NL11 ; DHM2 ; DHM3  that an intense terahertz
optical field can excite the fluctuation $\delta\rho_{s}$ of the superfluid
density $\rho_{s}$ in the second-order optical response, which manifests
itself in the third-harmonic current. This fluctuation is attributed to the
excitation of the Higgs mode, due to an observed resonance when twice of the
optical frequency is tuned at $\omega_{H}$NL8 ; NL9 ; NL10 . Inspired by the
experimental finding, a great deal of theoretical efforts have been devoted to
the nonlinear optical response of superconductors. However, rather than
straightening out the situation, these theoretical descriptions make the
understanding of the existing and growing experimental findings muddled.

Specifically, early theoretical studies of the nonlinear optical response in
superconductors have used the BlochAm1 ; Am2 ; Am7 ; Am9 ; Am11 ; Am14 ; Am15
; NL7 ; NL8 ; NL9 ; NL10 ; NL11  or LiouvilleAm3 ; Am4 ; Am8 ; Am10 ; Am16
equation derived in Anderson pseudospin pictureAs , with the vector potential
${\bf A}$ alone. The second-order light-matter interaction
$H_{p}=e^{2}A^{2}\tau_{3}/(2m)$ naturally emerges in these descriptionsAm1 ;
Am2 ; Am3 ; Am4 ; Am7 ; Am8 ; Am9 ; Am10 ; Am11 ; Am14 ; NL7 ; NL8 ; NL9 ;
NL10 ; NL11  as a pseudo field along $z$ direction, with $\tau_{i}$ being the
Pauli matrices in Nambu space. This interaction can pump the quasiparticle
correlation (i.e., pseudospin precession) and then causes the fluctuation of
the order parameter, which was directly considered as the Higgs mode to
explain the experimental findingsAm1 ; Am2 ; Am3 ; Am4 ; Am7 ; Am8 ; Am9 ;
Am10 ; Am11 ; Am14 ; NL7 ; NL8 ; NL9 ; NL10 ; NL11 . Nevertheless, a latter
symmetry analysis in Anderson pseudospin picture finds a vanishing (finite)
correlation between amplitude (phase) mode and external pseudo fieldsymmetry .
This implies that the theoretically obtained order-parameter fluctuation in
Bloch or Liouville equation is a phase fluctuation rather than the claimed
amplitude oneGIKE2 . Moreover, with the isotropic pump effect alone, the Bloch
or Liouville equation in the literatureAm1 ; Am2 ; Am3 ; Am4 ; Am7 ; Am8 ; Am9
; Am10 ; Am11 ; Am14 ; NL7 ; NL8 ; NL9 ; NL10 ; NL11  fails to drive the
optical current and hence derive the superfluid density, since no drive effect
(i.e., linear light-matter interaction) is included.

By applying the standard path-integral approach with the vector potential
alone, Cea et al.Cea1 ; Cea2 ; Cea3  further considered the linear-order
light-matter interaction $H_{d}={\hat{\bf p}}\cdot{e}{\bf A}/m$ (i.e., drive
effect from vector potentialG1 ; GIKE1 , with ${\hat{\bf p}}$ being the
momentum operator), in addition to second-order one $H_{p}$ (pump effect). In
the second-order response, they foundCea1 ; Cea2 ; Cea3  that neither drive
nor pump effects can excite the Higgs mode $\delta|\Delta|$. Starting from
this theoretical investigation, it is believed afterwards that the Higgs-mode
generation is zero at clean limit. Meanwhile, Cea et al.Cea1 ; Cea2 ; Cea3
found that the pump effect $H_{p}$ can cause the fluctuation $\delta{n}$ of
the charge density $n$. As the superfluid density $\rho_{s}$ is proportional
to $n|\Delta|^{2}$G1 , it is therefore speculatedCea1 ; Cea2 ; Cea3  that the
experimentally observed $\delta\rho_{s}$ is attributed to the charge-density
fluctuation $\delta{n}$ rather than the Higgs mode $\delta|\Delta|$. Several
polarization-resolved measurements were performed afterwardsNL10 ; NL11 ; DHM2
; DHM3 ; FHM , as the theoretically predicted signal of the Higgs mode
(charge-density fluctuation) is isotropic (anisotropic)Cea1 . However, an
isotropic optical response is experimentally observedNL10 ; NL11 ; DHM2 ; DHM3
, giving firm evidence to support the previous observation of the Higgs mode.
Recent theoretical attention then tends to focus on and emphasize the
important role of the impurity scattering to mediate the Higgs-mode
generationAm16 ; ImR1 ; ImR2 ; ImR3 ; Silaev . To handle the microscopic
scattering seriously, SilaevSilaev ; Silaev0  used the Eilenberger
equationEilen ; Ba20 ; Eilen1  that only involves the drive effect $H_{d}$ by
vector potential. In the second-order response, he also found a vanishing
Higgs-mode excitation at clean limit, but derived a finite one in dirty case
to dominate over the charge-density fluctuationSilaev .

Through a gauge-invariant kinetic equation (GIKE) approach with complete
electromagnetic effectGIKE1 ; GIKE2 ; GIKE3 ; GIKE4 , our recent studyGIKE2
that calculates the amplitude and phase modes on an equal footing obtained
totally different results in the second-order optical response at clean limit:
a vanishing charge-density fluctuation, and a finite Higgs-mode generation,
contributed by the drive effect $H_{d}$ of vector potential. Physically, both
results can be understood as follows. Firstly, it is well known from the
symmetry analysis that there is no second-order current $j^{(2)}$ in systems
with the inversion symmetry. Then, from the charge conservation
$\partial_{t}e\delta{n}+{\bm{\nabla}_{R}}\cdot{\bf j}=0$ in the second-order
regime ($2{\Omega}e{\delta}n^{(2)}+2{\bf q}\cdot{\bf j}^{(2)}=0$, with
$\Omega$ and ${\bf q}$ being the optical frequency and momentum,
respectively), the second-order charge-density fluctuation $e{\delta}n^{(2)}$
is forbidden. Secondly, according to the Ginzburg-Landau theory, the general
superconducting Lagrangian at clean limit readsAm6 ; PI2GL ; PI2GL2 ; G1

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

with $\Psi$ denoting the order parameter and $\gamma_{L}$, $\alpha_{L}$,
$\beta_{L}$ as well as $\lambda_{L}$ representing the Landau parameters. From
the Lagrangian above, the equilibrium order parameter
$\Psi_{0}=\sqrt{-\alpha_{L}/\beta_{L}}$. Whereas by only considering the
amplitude fluctuation $\delta|\Psi|$ with $\Psi=\Psi_{0}+\delta|\Psi|$, one
can directly obtain its equation of motion:

$\big{[}(2|\Psi_{0}|)^{2}+\partial_{t}^{2}\big{]}\delta|\Psi|=-\frac{\lambda_{L}}{\beta_{L}}\frac{e^{2}A^{2}}{m}2\Psi_{0},$
(2)

where $\beta_{L}=\gamma_{L}=\frac{7R(3)}{8(\pi{T})^{2}}$ and
$\lambda_{L}=\frac{k_{F}^{2}}{3m}\frac{7R(3)}{4(\pi{T})^{2}}$ are used for
conventional $s$-wave superconductorsG1 ; PI2GL , with $T$ denoting
temperature and $R(x)$ being the Riemann zeta function. Then, from above
equation, one can immediately find the Higgs-mode energy spectrum
$\omega_{H}=2|\Psi_{0}|$ on the left-hand side of the equation, and in
particular, a finite second-order response of Higgs mode at clean limit on the
right-hand side of the equation. Moreover, it is established in the derivation
of the Ginzburg-Landau equation from the basic Gorkov equationG1  that the
kinetic term [i.e., the last term in Eq. (1)] is solely attributed to the
drive effect $H_{d}$, implying a finite (zero) contribution from the drive
(pump) effect to Higgs-mode excitation.

Consequently, in the second-order response at clean limit, the finite Higgs-
mode excitation and vanishing charge-density fluctuation derived from
GIKEGIKE2  and justified by the phenomenological Ginzburg-Landau theory and
the charge conservation mentioned above, pose a sharp contrast to the previous
derivations from the path-integral approach by Cea et al.Cea1 ; Cea2 ; Cea3
and Eilenberger equation by SilaevSilaev  where zero Higgs-mode generation and
finite charge-density fluctuation are obtained. Moreover, it also becomes
particularly bizarre that both path-integral approachPI2GL  and Eilenberger
equationBa20  can recover the Ginzburg-Landau equation, but obtained a zero
Higgs-mode generationCea1 ; Cea2 ; Cea3 ; Silaev  which holds against the
Ginzburg-Landau theory.

In the present work, we resolve this controversy by re-examining the previous
derivations with the vector potential alone within the path-integral approach
by Cea et al.Cea1 ; Cea2 ; Cea3  and Eilenberger equation by SilaevSilaev .
While we have successfully recovered their results, it is found that both
derivations contain flaws. In the path-integral approach, after integration
over the Fermi field, the coupling of the Higgs mode to the pump effect
$H_{p}$ and the second order of the drive effect $H_{d}$ emerge in the second-
and third-order perturbation expansions of the action, respectively.
Nevertheless, in Refs. Cea1, ; Cea2, ; Cea3, , only the second-order expansion
is kept, leading to zero Higgs-mode generation because of the vanishing
correlation between amplitude mode and pump effectsymmetry . The essential
third-order expansion, which is related to the Ginzburg-Landau kinetic
termPI2GL  and hence finite Higgs-mode generation, is excessively overlooked.
In Ref. Silaev,  within the Eilenberger equation, in the summation over the
Fermion Matsubara frequency, the involved continuous optical frequency
$\Omega$ over $\pi{T}$ is considered as discrete even integer, leading to
vanishing amplitude-response coefficient. We prove that after fixing these
flaws, in the second-order optical response at clean limit, one can find a
finite Higgs-mode generation contributed by the drive effect $H_{d}$, exactly
recovering the results from GIKEGIKE2  as well as the Ginzburg-Landau theory.

Generally, according to the gauge structure in superconductors first revealed
by Nambugi0 , among the scalar potential $\phi$, vector potential ${\bf A}$ as
well as the superconducting-phase-related effective electromagnetic potential
$\partial_{\mu}\delta\theta$, one can not choose two quantities simultaneously
to be zero in superconductors, e.g., considering the vector potential alone.
Therefore, we further extend the previous path-integral approach to include
the electromagnetic effects from the scalar potential and phase mode, which
have been overlooked in the previous theoretical descriptionsAm1 ; Am2 ; Am3 ;
Am4 ; Am7 ; Am8 ; Am9 ; Am10 ; Am11 ; Am14 ; NL7 ; NL8 ; NL9 ; NL10 ; NL11 ;
Cea1 ; Cea2 ; Cea3 ; Silaev . Then, in the second-order response at clean
limit, the finite contribution in the Higgs-mode generation from the scalar
potential as well as the vanishing charge-density fluctuation, both of which
have previously been obtained from GIKEGIKE2 , are recovered. On one hand,
differing from the contribution of the vector potential that emerges only at
finite temperatureDS1 , the one of the scalar potential is finite upon cooling
to zero temperature, and is essential because of the gauge structure in
superconductorsgi0 . On the other hand, we show that a spatially uniform phase
mode is generated in the second-order response, and exactly cancels the
unphysical excitation of the charge-density fluctuation reported by Cea et
al.Cea1 ; Cea2 ; Cea3 , guaranteeing the charge conservation. Consequently,
the present study arrives at unified conclusions about the finite Higgs-mode
generation and vanishing charge-density fluctuation in the second-order
optical response at clean limit, in consistency with the phenomenological
Ginzburg-Landau theory and the charge conservation, and hence, can help
understanding the experimental findings. Furthermore, a disscussion about the
application of Matsubara formalism in the derivation of superconducting gap
dynamics is presented.

##  II Hamiltonian

We begin with the Bogoliubov-de Gennes Hamiltonian of the conventional
$s$-wave superconducting states in the presence of the electromagnetic
potential $A_{\mu}=(\phi,{\bf A})$G1 :

$H={\int}{d{\bf x}}\psi^{\dagger}(x)[(\xi_{{\hat{\bf p}}-e{\bf
A}\tau_{3}}+e\phi+\mu_{H})\tau_{3}+{\hat{\Delta}}(x)]\psi(x),$ (3)

where $\psi(x)=[\psi_{\uparrow}(x),\psi^{\dagger}_{\downarrow}(x)]^{T}$ is the
field operator in the Nambu space; $x=(x_{0},{\bf x})$ denotes the space-time
vector; $\xi_{\hat{\bf p}}={{\bf{\hat{p}}}^{2}}/({2m})-\mu$ with $m$ and $\mu$
being the effective mass and chemical potential; the momentum operator
${\hat{\bf p}}=-i\hbar{\bm{\nabla}}$;
$\mu_{H}(x)=\sum_{x^{\prime}}V(x-x^{\prime})\delta{n}(x^{\prime})$ denotes the
Hartree field, which is equivalent to the Poisson equation and characterizes
the induced scalar potential by density fluctuation $\delta{n}(x^{\prime})$;
$V(x-x^{\prime})$ represents the Coulomb potential;
${\hat{\Delta}}(x)=\Delta(x)\tau_{+}+{\Delta}^{*}(x)\tau_{-}$, while
considering the phase and amplitude modes, the order parameter reads
$\Delta(x)=[\Delta_{0}+\delta|\Delta|(x)]e^{i\delta\theta(x)}$.

It is noted that the phase mode $\delta\theta$ in Eq. (3) can be effectively
removed by a unitary transformation:

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

and then, one hasgi0 ; AK

$H=\\!\\!\\!{\int}{d{\bf x}}\psi^{\dagger}(x)(H_{0}+H_{\rm
LM}+\delta|\Delta|\tau_{1})\psi(x),$ (5)

with the free BCS Hamiltonian

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

and the light-matter interaction

$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)

Here, the gauge-invariant superconducting momentum ${\bf
p}_{s}={\bm{\nabla}}_{\bf x}\delta\theta/2-e{\bf A}$ and effective field
$\mu_{\rm eff}=e\phi+\partial_{x_{0}}\delta\theta/2+\mu_{H}$gi0 ; AK . Then,
it is clearly seen that the phase mode provides an effective electromagnetic
potential $eA^{\rm
eff}_{\mu}=(\partial_{x_{0}}\delta\theta/2,-{\bm{\nabla}_{\bf
x}}\delta\theta/2)$, in consistency with the gauge structure in
superconductors first revealed by Nambugi0 :

$\displaystyle eA_{\mu}$ $\displaystyle\rightarrow$ $\displaystyle
eA_{\mu}-\partial_{\mu}\chi,$ (8) $\displaystyle\delta\theta$
$\displaystyle\rightarrow$ $\displaystyle\delta\theta+2\chi.$ (9)

Here, $\partial_{\mu}=(\partial_{x_{0}},-{\bm{\nabla}}_{\bf x})$. It is noted
that the gauge invariance is essential for the theoretical descriptions, since
it guarantees the charge conservation, as first proved by Nambu via the
generalized Ward’s identitygi0 ; Ba0 .

By assuming the electromagnetic potential ${\bf A}(x)={\bf
A}_{0}e^{i\Omega{x_{0}}-i{\bf q}\cdot{\bf x}}$ and $\phi(x)=\phi_{0}({\bf
x})e^{i\Omega{x_{0}}-i{\bf q}\cdot{\bf x}}$ with ${\phi_{0}({\bf
x})}={\bar{\phi}_{0}}+{\bf E}_{\phi}\cdot{\bf x}$ and $E_{\phi}$ being the
transverse field, one has

$\displaystyle e\delta{n}$ $\displaystyle=$ $\displaystyle
e{\delta}n^{(1)}e^{i\Omega{x_{0}}-i{\bf q}\cdot{\bf
x}}+{\delta}n^{(2)}e^{2i\Omega{x_{0}}-2i{\bf q}\cdot{\bf x}},$ (10)
$\displaystyle\delta\theta$ $\displaystyle=$
$\displaystyle\delta\theta^{(1)}e^{i\Omega{x_{0}}-i{\bf q}\cdot{\bf
x}}+\delta\theta^{(2)}e^{2i\Omega{x_{0}}-2i{\bf q}\cdot{\bf x}},$ (11)
$\displaystyle\delta|\Delta|$ $\displaystyle=$
$\displaystyle\delta|\Delta|^{(1)}e^{i\Omega{x_{0}}-i{\bf q}\cdot{\bf
x}}+\delta|\Delta|^{(2)}e^{2i\Omega{x_{0}}-2i{\bf q}\cdot{\bf x}},$ (12)

where $e{\delta}n^{(l)}$ as well as $\delta\theta^{(l)}$ and
$\delta|\Delta|^{(l)}$ denote the $l$-th order responses of the charge
density, phase and Higgs modes, respectively. Then, one correspondingly finds
the amplitudes of $\mu_{H}$ and $\mu_{\rm eff}$ as well as ${\bf p}_{s}$ in
the $l$-th-order response as

$\displaystyle\mu_{H}^{(l)}$ $\displaystyle=$ $\displaystyle
V_{q}\delta{n}^{(l)},$ (13) $\displaystyle\mu_{\rm eff}^{(l)}$
$\displaystyle=$ $\displaystyle e\phi_{0}({\bf
x})\delta_{l,1}+il\Omega\delta\theta^{(l)}/2+\partial_{x_{0}}\delta\theta^{(l)}/2+\mu_{H}^{(l)},~{}~{}~{}$
(14) $\displaystyle{\bf p}^{(l)}_{s}$ $\displaystyle=$ $\displaystyle-i{\bf
q}\delta\theta^{(l)}/2+{\bm{\nabla}}_{\bf x}\delta\theta^{(l)}/2-e{\bf
A}_{0}\delta_{l,1}.$ (15)

Particularly, it is noted that based on the gauge structure in Eqs. (8) and
(9) of superconductors, in the $l$-th-order response, one can choose the
phase-related effective electromagnetic potential
$\partial_{\mu}\delta\theta^{(l)}$ to be zero. Then, the amplitude
$\delta\theta^{l}$ of the $l$-th-order response of the phase mode is spatially
uniform and time-independent as a background. In this situation, it has been
established in the literatureAK ; AK2 ; Ba0 ; pm0 ; Am0 ; Ba9 ; Ba10 ; pm5 ;
pi1 ; pi4 ; GIKE2  that the phase mode in the linear regime
$\delta\theta^{(1)}$, as a scalar quantity, responds to the longitudinal
electromagnetic field solely, which experiences the Coulomb screening.
Consequently, the uniform linear response of the phase mode
${\delta\theta^{(1)}}/{2}=\frac{i\Omega{e{\bar{\phi}}_{0}-\omega^{2}_{p}{i{\bf
q}{\cdot}e{\bf A}_{0}}/{q^{2}}}}{\Omega^{2}-v^{2}_{p}q^{2}}$ becomes
${\delta\theta^{(1)}}/{2}=\frac{i\Omega{e{\bar{\phi}}_{0}-\omega^{2}_{p}{i{\bf
q}{\cdot}e{\bf
A}_{0}}/{q^{2}}}}{(\Omega^{2}-v^{2}_{p}q^{2})(1-\omega^{2}_{p}/\Omega^{2})}\approx\frac{i\Omega{e{\bar{\phi}}_{0}-\omega^{2}_{p}{i{\bf
q}{\cdot}e{\bf A}_{0}}/{q^{2}}}}{\Omega^{2}-\omega^{2}_{p}}$ after considering
the long-range Coulomb interactionAK ; Ba0 ; pm0 ; Am0 ; Ba9 ; Ba10 ; pm5 ;
AK2 ; GIKE2 , with $v_{p}$ being the velocity of the phase mode and
$\omega_{p}$ denoting the plasma frequency. The original gapless spectrum
(resonance pole) is then effectively lifted up to the high-energy plasma
frequency as a consequence of the Anderson-Higgs mechanismAHM . At this case,
with $\Omega\ll\omega_{p}$, one has $\delta\theta^{(1)}/2\approx{i{\bf
q}{\cdot}e{\bf A}_{0}}/{q^{2}}$, which cancels the unphysical longitudinal
vector potential in ${\bf p}_{s}^{(1)}={\bf q}({\bf q}\cdot{e}{\bf
A}_{0})/q^{2}-e{\bf A}_{0}$, and then, the superconducting momentum ${\bf
p}_{s}$ that appears in the previous theoretical descriptions such as
Ginzburg-Landau equationG1  and Meissner supercurrentG1  as well as Anderson-
pump effectAm1 ; Am2 ; Am3 ; Am4 ; Am7 ; Am8 ; Am9 ; Am10 ; Am11 ; Am14 ; NL7
; NL8 ; NL9 ; NL10 ; NL11  only involves the physical transverse vector
potential. Moreover, thanks to the Coulomb screening (i.e.,
$\mu_{H}^{(1)}=-2DV_{q}\frac{e{\bar{\phi}}_{0}+i\Omega\delta\theta^{(1)}/2}{1+2DV_{q}}$
with $D$ being the density of states)AK ; GIKE2 , at long-wave limit, one
finds $\mu_{\rm eff}^{(1)}=e{\bf E}_{\phi}\cdot{\bf x}$, in which the original
longitudinal part $e{\bar{\phi}}_{0}+i\Omega\delta\theta^{(1)}$ vanishes.
Consequently, considering the spatially uniform transverse fields for the
optical response (i.e., ${\bf q}\rightarrow 0$ and the optical electric field
${\bf E}_{0}=-{\bf E}_{\phi}-i\Omega{\bf A}_{0}$ are spatially uniform and
transverse one), the linear-order component of the light-matter interaction
$H_{\rm LM}$ in Eq. (7) is written as

$H^{(1)}_{\rm LM}=\Big{[}-\frac{{\hat{\bf p}}{\cdot}{e{\bf A}_{0}}}{m}+(e{\bf
E}_{\phi}\cdot{\bf x})\tau_{3}\Big{]}e^{i\Omega{x_{0}}},$ (16)

whereas the second-order one reads

$H^{(2)}_{\rm
LM}=\Big{(}\frac{e^{2}A^{2}_{0}}{2m}+i\Omega\delta\theta^{(2)}+\mu_{H}^{(2)}\Big{)}e^{2i\Omega{x_{0}}}\tau_{3}.$
(17)

It is noted that $H^{(1)}_{\rm LM}$ represents the drive effects of the vector
and scalar potentialsGIKE1 . $H^{(2)}_{\rm LM}$ denotes the pump effect, in
which besides the conventional contribution $H_{p}$ from the vector potential
as mentioned in the introduction, the second-order response of the phase mode
and Hartree field also play an important role.

##  III Analytic Derivation

In this section, for the convenience of the comparison and understanding, we
first briefly introduce the results of the second-order response of the
collective modes from GIKE at clean limitGIKE2 , and then, separately use the
Eilenberger equation and path-integral approach to derive the second-order
response of Higgs mode at clean limit.

###  III.1 GIKE

In this part, we briefly introduce the results of the second-order responses
of the collective modes from GIKE at clean limitGIKE2 . Particularly, we
extend our previous results in Ref. GIKE2,  at low temperature up to $T_{c}$.
Specifically, the GIKEGIKE1 ; GIKE2  is derived from the basic Gorkov equation
of $\tau_{0}$-Green function
$G_{0}(x,x^{\prime})=-i\langle{\hat{T}}\psi(x)\psi^{\dagger}(x^{\prime})\rangle$
based on equal-time scheme ($t=t^{\prime}$)GQ2 ; GQ3 , with ${\hat{T}}$ being
the chronological ordering. To retain the gauge invariance, the gauge-
invariant $\tau_{0}$-Green function is constructed through the Wilson
lineWilson . As a result of the gauge invariance, the complete electromagnetic
effects are includedGIKE1  and the charge conservation is naturally
satisfiedGIKE2  in the GIKE.

In this microscopic approach, the response of system is described by density
matrix $\rho_{\bf k}=\rho^{(0)}_{\bf k}+\delta\rho_{\bf k}(R)$ in Nambu space,
which consists of the equilibrium part $\rho^{(0)}_{\bf
k}=\frac{1}{2}+\frac{f(E_{\bf k}^{+})-f(E_{\bf
k}^{-})}{2}(\frac{\xi_{k}}{E_{\bf k}}\tau_{3}+\frac{\Delta_{0}}{E_{\bf
k}}\tau_{1})$ and nonequilibrium one $\delta\rho_{\bf k}(R)$. Here,
$R=(x+x^{\prime})/2=(t,{\bf R})$ represents the center-of-mass coordinate;
$f(x)$ denotes the Fermi-distribution function; $E_{\bf k}^{\pm}$ represents
the quasi-electron and quasi-hole energies, which in the presence of
superconducting momentum is written asFF4 ; FF5 ; FF6 ; FF8 ; FF9 ; GIKE1

$E_{\bf k}^{\pm}={\bf v}_{\bf k}\cdot{\bf p}_{s}\pm{E_{\bf k}}.$ (18)

Here, $E_{\bf k}=\sqrt{\xi_{\bf k}^{2}+\Delta_{0}^{2}}$ is the Bogoliubov
quasiparticle energy and ${\bf v}_{\bf k}\cdot{\bf p}_{s}$ denotes the Doppler
shiftFF4 ; FF5 ; FF6 ; GIKE1 , with the group velocity ${\bf v}_{\bf
k}=\partial_{\bf k}\xi_{\bf k}$. The nonequilibrium $\delta\rho_{\bf k}$ can
be solved from the GIKEGIKE2 :

$\displaystyle\\!\\!\\!\\!\\!\\!\partial_{t}\rho_{\bf
k}\\!+\\!i\Big{[}\Big{(}\xi_{k}\\!+\\!\frac{p_{s}^{2}}{2m}\\!+\\!\mu_{\rm
eff}\\!+\\!\mu_{H}\Big{)}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}\\!+\\!\delta|\Delta|\tau_{1},\rho_{\bf
k}\Big{]}$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}\\!-\\!\frac{i}{8}\left[({\bm{\nabla}}_{\bf
R}\\!+\\!2i{\bf p}_{s}\tau_{3})({\bm{\nabla}}_{\bf R}\\!\\!+\\!\\!2i{\bf
p}_{s}\tau_{3})|\Delta|\tau_{1},\partial_{\bf k}\partial_{\bf k}\rho_{\bf
k}\right]$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}\\!+\\!\frac{1}{2}\left\\{e{\bf
E}\tau_{3}\\!\\!-\\!\\!({\bm{\nabla}}_{\bf R}\\!\\!+\\!\\!2i{\bf
p}_{s}\tau_{3})|\Delta|\tau_{1},\partial_{\bf k}\rho_{\bf
k}\right\\}\\!-\\!\Big{[}\frac{i\nabla^{2}_{\bf R}}{8m}\tau_{3},\rho_{\bf
k}\Big{]}$ $\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}\\!+\\!\Big{\\{}\frac{{\bf
k}\\!\cdot\\!{\bm{\nabla}}_{\bf R}}{2m}\tau_{3},\rho_{\bf
k}\Big{\\}}\\!-\\!\Big{[}\frac{{\nabla_{\bf R}}\\!\circ\\!{\bf
p}_{s}}{4m}\tau_{3},\tau_{3}\rho_{\bf k}\Big{]}=\partial_{t}\rho_{\bf k}|_{\rm
scat},$ (19)

where we have applied the unitary transformation in Eq. (4) to effectively
remove the phase mode from the order parameter. Here, ${\nabla_{\bf
R}}\\!\circ\\!{\bf p}_{s}=(2{\bf p}_{s}\\!\cdot\\!{\bm{\nabla}}_{\bf
R}\\!+\\!{\bm{\nabla}}_{\bf R}\\!\cdot\\!{\bf p}_{s})$; the electric field
$e{\bf E}=-{\bm{\nabla}}_{\bf R}(e\phi+\mu_{H})-\partial_{t}e{\bf
A}=-{\bm{\nabla}}_{\bf R}\mu_{\rm eff}+\partial_{t}{\bf p}_{s}$.

The gauge-invariant density and current read $n=\sum_{\bf k}(1+2\rho_{{\bf
k}3})$ and ${\bf j}=\sum_{\bf k}\big{(}\frac{e{\bf k}}{m}\rho_{{\bf
k}0}\big{)}$, respectively. Moreover, after the unitary transformation, the
equation of the order parameter becomesGIKE2

$\displaystyle U{\sum_{\bf k}}^{\prime}\rho_{{\bf k}1}$ $\displaystyle=$
$\displaystyle-|\Delta|,$ (20) $\displaystyle U{\sum_{\bf
k}}^{\prime}\rho_{{\bf k}2}$ $\displaystyle=$ $\displaystyle 0,$ (21)

where $U$ denotes the pairing potential and $\rho_{{\bf k}i}$ stands for the
$\tau_{i}$ component of $\rho_{\bf k}$; $\sum_{\bf k}^{\prime}$ here and
hereafter stands for the summation restricted in the spherical shell
($|\xi_{\bf k}|\leq\omega_{D}$) with $\omega_{D}$ being the Debye frequency.
It is noted that Eq. (20) gives the gap equation and hence the Higgs mode,
whereas Eq. (21) determines the phase fluctuation as revealed in our previous
workGIKE2 .

For the weak probe, by expanding $\delta\rho_{\bf k}=\delta\rho^{(1)}_{\bf
k}+\delta\rho^{(2)}_{\bf k}$ with $\delta\rho^{(1)}_{\bf k}$ and
$\delta\rho^{(2)}_{\bf k}$ being the first and second order responses to
optical probe, the GIKE becomes a chain of equations, as its first order only
involves $\delta\rho^{(1)}_{\bf k}$ and equilibrium $\rho^{(0)}_{\bf k}$ and
its second order involves $\delta\rho^{(2)}_{\bf k}$, $\delta\rho^{(1)}_{\bf
k}$ and $\rho^{(0)}_{\bf k}$. Consequently, starting from the lowest order,
one can calculate $\delta\rho^{(1)}_{\bf k}$ and $\delta\rho^{(2)}_{\bf k}$ in
sequence, and then, obtain the linear and second-order responses of the Higgs
(phase) mode by substituting the solved $\delta\rho^{(1)}_{\bf k}$ and
$\delta\rho^{(2)}_{\bf k}$ into Eq. (20) [Eq. (21)], respectively.

As revealed in our previous workGIKE2 , the linear response of the Higgs mode
vanishes in the long-wave limit, whereas the linear response of the phase mode
recovers the previous resultsAK ; Ba0 ; pm0 ; Am0 ; Ba9 ; Ba10 ; pm5  of the
Anderson-Higgs mechanismAHM  mentioned in Sec. II. Here, we present the
second-order responses of the Higgs $\delta|\Delta|^{(2)}$ and phase
$\delta\theta^{(2)}$ modes as well as charge-density fluctuation
$\delta{n}^{(2)}$ derived from the GIKE at clean limit (the specific
derivation can be found in Ref. GIKE2, ):

$4(\Delta_{0}^{2}-\Omega^{2})\beta_{g}\delta|\Delta|^{(2)}\\!\\!=\\!-\frac{2\Delta_{0}v_{F}^{2}}{3}[\gamma_{g}{e^{2}_{\mu_{\rm
eff}}}\\!+\\!(p_{s}^{(1)})^{2}\lambda_{g}],$ (22)

and

$-i\Omega\delta\theta^{(2)}\\!-\mu_{\rm
H}^{(2)}\\!=\\!\frac{{p_{s}^{2}}}{2m}\\!-\\!\frac{1}{3mu_{g}}({w_{g}}{e^{2}_{\mu_{\rm
eff}}}\\!+t_{g}p_{s}^{(1)}e_{\mu_{\rm eff}}),$ (23)

as well as

$\delta{n^{(2)}}=0,$ (24)

with the amplitude-correlation coefficient

$\beta_{g}={\sum_{\bf k}}^{\prime}\Big{[}\frac{1-2f(E_{\bf k})}{2E_{\bf
k}}\frac{1}{E_{\bf k}^{2}-\Omega^{2}}+\frac{\partial_{E_{\bf k}}f(E_{\bf
k})}{E^{2}_{\bf k}}\Big{]},$ (25)

amplitude-response coefficients

$\displaystyle\gamma_{g}$ $\displaystyle=$ $\displaystyle{\sum_{\bf
k}}^{\prime}\frac{\xi_{\bf k}}{E_{\bf k}^{2}-\Omega^{2}}\partial^{2}_{\xi_{\bf
k}}\Big{[}{\xi_{\bf k}}\frac{2f(E_{\bf k})-1}{2E_{\bf k}}\Big{]},$ (26)
$\displaystyle\lambda_{g}$ $\displaystyle=$ $\displaystyle{\sum_{\bf
k}}^{\prime}\frac{\partial^{2}_{E_{\bf k}}f(E_{\bf k})}{E_{\bf k}},$ (27)

and the phase-response coefficients

$\displaystyle u_{g}$ $\displaystyle=$ $\displaystyle{\sum_{\bf
k}}^{\prime}\frac{\Delta^{2}_{0}}{E_{\bf k}^{2}-\Omega^{2}}\frac{2f(E_{\bf
k})-1}{E_{\bf k}},$ (28) $\displaystyle t_{g}$ $\displaystyle=$
$\displaystyle{\sum_{\bf k}}^{\prime}\frac{\Delta^{2}_{0}}{E_{\bf
k}^{2}-\Omega^{2}}\partial_{\xi_{\bf k}}\Big{[}\xi_{\bf k}\frac{2f(E_{\bf
k})-1}{E_{\bf k}}\Big{]},$ (29) $\displaystyle w_{g}$ $\displaystyle=$
$\displaystyle{\sum_{\bf k}}^{\prime}\frac{2\Delta^{2}_{0}}{E_{\bf
k}^{2}-\Omega^{2}}(\xi_{\bf k}\partial^{2}_{\xi_{\bf k}}+\partial_{\xi_{\bf
k}})\Big{[}\xi_{\bf k}\frac{2f(E_{\bf k})-1}{E_{\bf k}}\Big{]},~{}~{}~{}$ (30)

as well as $e_{\mu_{\rm eff}}=\nabla_{\bf R}{\mu^{(1)}_{\rm eff}}/(i\Omega)$.
It is noted that in our previous workGIKE2  which considers the low
temperature and weak optical probe, we neglect the Doppler shift in the
quasiparticle energy [Eq. (18)] by assuming $v_{k}p_{s}<\Delta_{0}$. In the
present work, we sublate this approximation in order to extend the calculation
to the entire temperature regime. The considered Doppler shift does not
influence the previous calculation, but causes a additional contribution
[i.e.,the second term on the right-hand side of Eq. (22)] in the second-order
response of the Higgs mode $\delta|\Delta|^{(2)}$ through $\rho_{\bf k}^{(0)}$
in Eq. (20). From Eqs. (22) and (23), one finds that the second-order
responses of the Higgs $\delta|\Delta|^{(2)}$ and phase $\delta\theta^{(2)}$
modes are decoupled, as they represent mutually orthogonal excitations in the
Mexican-hat potential of free energy. It is also noted that all source terms
on the right-hand side of both Eqs. (22) and (23) are gauge-invariant.

Higgs-mode generation.—Considering the spatially uniform transverse fields for
the optical response and choosing the phase-related effective electromagnetic
potential $\partial_{\mu}\delta\theta^{(1)}$ to be zero, Eq. (22) following
the analysis of the light-matter interaction $H_{\rm LM}$ in Sec. II becomes

$4(\Delta_{0}^{2}-\Omega^{2})\beta_{g}\delta|\Delta|^{(2)}\\!\\!=\\!-\frac{e^{2}v_{F}^{2}}{3}2\Delta_{0}\Big{[}\frac{\gamma_{g}E^{2}_{\phi}}{(i\Omega)^{2}}\\!+\\!{A_{0}^{2}}\lambda_{g}\Big{]},$
(31)

It is pointed out that the right-hand side of Eq. (31) arises from the second
order of the drive effects $H_{LM}^{(1)}$ of vector potential and scalar
potential, whereas the pump effect $H_{LM}^{(2)}$ makes no contribution, in
consistency with the vanishing correlation between amplitude mode and pump
effect as revealed by previous symmetry analysissymmetry .

Consequently, from Eq. (31), one immediately finds a finite second-order
response of the Higgs mode at clean limit, contributed solely by the drive
effects. Actually, according to the analysis of the Ginzburg-Landau theory as
mentioned in the introduction, this finite second-order response derived from
the GIKE is expected, since the GIKE near $T_{c}$ can recover the Ginzburg-
Landau equationGIKE1 . Particularly, from Eqs. (25) and (27), at low
frequency, near $T_{c}$, one has

$\displaystyle\beta_{g}\\!\\!$ $\displaystyle\approx$
$\displaystyle\\!\\!D\int{d\xi_{\bf k}}\frac{1}{|\xi_{\bf
k}|}\partial_{|\xi_{\bf k}|}\Big{[}\frac{2f(|\xi_{\bf k}|)-1}{2|\xi_{\bf
k}|}\Big{]}$ (32) $\displaystyle=$ $\displaystyle
DT\sum_{\omega_{n}}\int{d\xi_{\bf k}}\frac{1}{|\xi_{\bf
k}|}\partial_{|\xi_{\bf k}|}\left[\frac{1}{(i\omega_{n})^{2}-\xi_{\bf
k}^{2}}\right]$ $\displaystyle=$ $\displaystyle
DT\sum_{\omega_{n}}\int{d\xi_{\bf k}}\frac{2}{[(\omega_{n})^{2}+\xi_{\bf
k}^{2}]^{2}}=\frac{7DR(3)}{4(\pi{T})^{2}},~{}~{}~{}$
$\displaystyle\lambda_{g}\\!\\!$ $\displaystyle\approx$
$\displaystyle\\!\\!D\int{d\xi_{\bf k}}\frac{\partial_{|\xi_{\bf
k}|}^{2}f(|\xi_{\bf k}|)}{|\xi_{\bf k}|}\\!=\\!\sum_{\omega_{n}}\int{d\xi_{\bf
k}}\frac{2DT}{(i\omega_{n}\\!-\\!|\xi_{\bf k}|)^{3}|\xi_{\bf k}|}$ (33)
$\displaystyle=$ $\displaystyle 4DT\sum_{\omega_{n}}\int{d\xi_{\bf
k}}\frac{3\omega_{n}^{2}-\xi_{\bf k}^{2}}{(\omega^{2}_{n}+\xi^{2}_{\bf
k})^{3}}=\frac{7DR(3)}{2(\pi{T})^{2}},~{}~{}~{}$

with $\omega_{n}=(2n+1)\pi{T}$ being the Matsubara frequency. Then,
considering the vector potential alone, Eq. (31) derived from the GIKE exactly
recovers Eq. (2) derived from the Ginzburg-Landau theory.

Moreover, it is noted in Eq. (31) that the drive effect $E_{\phi}$ of the
scalar potential also contributes to the Higgs-mode generation. This
contribution, being finite from $T=0$ to $T=T_{c}$, is different from the one
of the vector potential that emerges only at finite temperatureDS1 . Actually,
this difference is natural, since the superconductors can directly respond to
vector potential ${\bf A}$ (Meissner effect/Ginzburg-Landau kinetic term) in
addition to the electric field ${\bf E}=-{\nabla}_{\bf R}\phi-\partial_{t}{\bf
A}$ (optical-electric-field effect), differing from normal metals that solely
respond to electric field. Consequently, the contribution of the scalar
potential captures the optical-electric-field effect, whereas the contribution
of the vector potential characterizes the electromagnetic effects including
the Meissner effect/Ginzburg-Landau kinetic term as well as the optical-
electric-field effect. As mentioned in Sec. II, according to the gauge
structure [Eqs. (8) and (9)] in superconductors, among the scalar potential,
vector potential as well as the phase-related effective electromagnetic
potential $\partial_{\mu}\delta\theta$, one can only choose one quantity to be
zero in superconductors. Consequently, the inclusion of the contribution from
the scalar potential here is essential, since we have chosen zero
$\partial_{\mu}\delta\theta^{(1)}$.

Phase-mode generation.—For spatially uniform transverse optical fields and
zero $\partial_{\mu}\delta\theta^{(1)}$, Eq. (23) becomes

$-i\Omega\delta\theta^{(2)}\\!-\mu_{\rm
H}^{(2)}\\!=\\!\frac{{e^{2}A_{0}^{2}}}{2m}\\!-\\!\frac{e^{2}}{3m}\Big{[}\frac{w_{g}}{u_{g}}\frac{E^{2}_{\phi}}{(i\Omega)^{2}}\\!+\\!\frac{t_{g}}{u_{g}}\frac{E_{\phi}A_{0}}{i\Omega}\Big{]}.$
(34)

The first term on the right-hand side of Eq. (34) arises from the pump effect
$H_{p}$ of vector potential and the last two terms come from the drive effect
$H_{\rm LM}^{(1)}$. The finite contribution of the pump effect here agrees
with the finite correlation between phase mode and $H_{p}$ as revealed by
previous symmetry analysissymmetry .

Consequently, from Eqs. (24) and (34), one finds in the second-order optical
response a vanishing charge-density fluctuation $\delta{n^{(2)}}$ but a finite
phase-mode generation $\delta\theta^{(2)}$, respectively. As mentioned in the
introduction, the vanishing charge-density fluctuation $\delta{n^{(2)}}$
agrees with the inversion symmetry and charge conservation. As for the phase-
mode generation, the Hartree field on the left-hand side of Eq. (34) vanishes
as $\mu_{H}^{(2)}=2V_{q}\delta{n}^{(2)}=0$, whereas the right-hand side of the
equation is determined by the transverse optical field and hence free from the
influence of the Coulomb screening. As we pointed out in Ref. GIKE2, , this
phase-mode generation, showing a spatially uniform but temporally oscillating
phase, is a unique feature of the optical properties in the second-order
response, and does not manifest itself or incur any consequence in the
thermodynamic, electric or magnetic properties. Nevertheless, we show in the
following Sec. III.3 that this phase-mode generation $\delta\theta^{(2)}$ that
has long been overlooked in the literature is essential in the theoretical
description of the second-order optical response, since $\delta\theta^{(2)}$
provides an effective field to exactly cancel the unphysical excitation of the
charge-density fluctuation reported by Cea et al.Cea1 ; Cea2 ; Cea3  and hence
guarantee the charge conservation.

It is stressed that all results from the GIKE, including the Higgs-mode
generation [Eq. (31)] from both contributions of the drive effects of scalar
and vector potentials and phase mode generation [Eq. (34)] as well as
vanishing charge-density fluctuation [Eq. (24)] can be exactly recovered from
the path-integral approach within the gauge-invariant manner, to be shown in
the following Sec. III.3.

###  III.2 Eilenberger equation

Following the previous work by SilaevSilaev , we next use the Eilenberger
equationEilen ; Ba20 ; Eilen1  to derive the second-order optical response of
the Higgs mode. The Eilenberger equationEilen ; Ba20 ; Eilen1  is derived from
the basic Gorkov equation of $\tau_{3}$-Green function
$G_{3}(x,x^{\prime})=-i\tau_{3}\langle{\hat{T}}\psi(x)\psi^{\dagger}(x^{\prime})\rangle$
through the quasiclassical approximationQA1 :

$g(x_{0},x_{0}^{\prime},{\bf R},{{\bf k}}_{F})=\frac{i}{\pi}\int{d\xi_{\bf
k}}[G_{3}(x_{0},x_{0}^{\prime},{\bf R},{{\bf k}})].$ (35)

Here, $G(x_{0},x_{0}^{\prime},{\bf R},{{\bf k}})=\int{d{\bf
r}}G(x,x^{\prime})e^{-i{\bf k}\cdot({\bf x}-{\bf x^{\prime}})}$. In the
imaginary time domain ($x_{0}\rightarrow{i\tau_{1}}$,
$x_{0}^{\prime}\rightarrow{i\tau_{2}}$), in consideration of the spatially
uniform transverse vector potential alone, the Eilenberger equation at clean
limit readsSilaev :

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

where
$[X,{\hat{g}}]_{\tau}=X(\tau_{1})g(\tau_{1},\tau_{2})-g(\tau_{1},\tau_{2})X(\tau_{2})$
and
$\\{X,{\hat{g}}\\}_{\tau}=X(\tau_{1})g(\tau_{1},\tau_{2})+g(\tau_{1},\tau_{2})X(\tau_{2})$.
Moreover, the Eilenberger equation is supplemented by the normalization
conditionSilaev ; Ba20 ; Eilen1
$\int{d\tau}g(\tau_{1},\tau)g(\tau,\tau_{2})=1$. While the corresponding gap
equation is written asEilen ; Silaev ; Ba20 ; Eilen1

$\Delta=U{\rm Tr}[\langle{g(\tau,\tau)}\rangle_{F}\tau_{2}/2],$ (37)

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

By self-consistently solving Eqs. (36)-(37), one can formulate the Higgs-mode
generation at clean limit. Specifically, in the optical response with ${\bf
A}(\tau)={\bf A}_{0}e^{i\Omega\tau}$, the quasiclassical Green function is
given by ${\hat{g}}={\hat{g}}^{(0)}+{\hat{g}}^{(1)}+{\hat{g}}^{(2)}$ with the
$m$-th order response ${\hat{g}}^{(m)}$ written as

${\hat{g}}^{(m)}\\!\\!=\\!T\\!\sum_{\omega_{n}}g^{(m)}(i\omega_{n}\\!+\\!im\Omega,i\omega_{n})e^{i(\omega_{n}+m\Omega)\tau_{1}-i\omega_{n}\tau_{2}}.$
(38)

Consequently, the Eilenberger equation in Eq. (36) becomes a chain of
equations, whose first order only involves $g^{(0)}$ and $g^{(1)}$ and second
order involves $g^{(0)}$ and $g^{(1)}$ as well as $g^{(2)}$. Then, with the
equilibrium
$g^{(0)}(i\omega_{n})=\frac{\omega_{n}\tau_{3}-\Delta_{0}\tau_{2}}{\sqrt{\omega_{n}^{2}+\Delta_{0}^{2}}}$,
one can solve $g^{(1)}$ and $g^{(2)}$ in sequence, whose specific expressions
are given by (the detailed derivation can be found in Ref. Silaev, )

$\displaystyle\\!\\!\\!\\!\\!\\!g^{(1)}(i\omega_{n}\\!+\\!i\Omega,i\omega_{n})=\\!-\delta\Delta^{(1)}\Gamma_{2}(i\omega_{n}\\!+\\!i\Omega,i\omega_{n})$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}+\\!i(e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F})\Gamma_{3}(i\omega_{n}\\!+\\!i\Omega,i\omega_{n}),$ (39)
$\displaystyle\\!\\!\\!\\!\\!\\!g^{(2)}(i\omega_{n}\\!\\!+\\!2i\Omega,i\omega_{n})=-\delta\Delta^{(2)}\Gamma_{2}(i\omega_{n}\\!+\\!2i\Omega,i\omega_{n})\\!\\!$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}+\\!(e{\bf A}_{0}\\!\cdot\\!{\bf
v}_{F})^{2}\frac{[(\omega_{n}\\!\\!+\\!2\Omega)\tau_{3}\\!+\\!\Delta_{0}\tau_{2}]{\hat{\Xi}}\\!+\\!{\hat{\Xi}}(\omega_{n}\tau_{3}\\!+\\!\Delta_{0}\tau_{2})}{(\omega_{n}\\!+\\!2\Omega)^{2}\\!-\\!\omega^{2}_{n}},~{}~{}~{}~{}~{}$
(40)

with the correlation function

$\Gamma_{j}(i\omega_{n}\\!+\\!i\Omega,i\omega_{n})=\frac{g^{(0)}(i\omega_{n}+i\Omega)\tau_{j}g^{(0)}(i\omega_{n})\\!-\\!\tau_{j}}{\sqrt{(\omega_{n}+\Omega)^{2}\\!+\\!\Delta_{0}^{2}}\\!+\\!\sqrt{\omega^{2}_{n}\\!+\\!\Delta_{0}^{2}}},$
(41)

and
${\hat{\Xi}}=\Gamma_{3}(i\omega_{n}\\!+\\!2i\Omega,i\omega_{n}\\!+\\!i\Omega)\tau_{3}\\!-\\!\tau_{3}\Gamma_{3}(i\omega_{n}\\!+\\!i\Omega,i\omega_{n})$.

Higgs-mode generation.—Substituting the solved $g^{(1)}$ and $g^{(2)}$ into
the gap equation [Eq. (37)], one can obtain the linear and second-order
response of the Higgs mode, respectively. The linear response
$\delta\Delta^{(1)}=0$, as the anisotropic source term from the vector
potential on the right-hand side of Eq. (39) vanishes after the angular
average over the Fermi surface. The obtained second-order response of the
Higgs mode is written as

$(4\Delta_{0}^{2}\\!+\\!4\Omega^{2})\delta|\Delta|^{(2)}\beta_{E}=-\frac{e^{2}A_{0}^{2}v_{F}^{2}2\Delta_{0}}{3}\lambda_{E},$
(42)

with the amplitude-correlation coefficient $\beta_{E}$ and the amplitude-
response coefficient $\lambda_{E}$ given by

$\displaystyle\beta_{E}$ $\displaystyle=$ $\displaystyle
2\Big{\\{}\frac{{U^{-1}\\!+\\!T\sum_{\omega_{n}}{\rm
Tr}[\Gamma_{2}(i\omega_{n}\\!+\\!i2\Omega,i\omega_{n})\tau_{2}/2]}}{(4\Delta_{0}^{2}\\!+\\!4\Omega^{2})}\Big{\\}}=T\sum_{\omega_{n}}\frac{1/[\sqrt{(\omega_{n}\\!+\\!2\Omega)^{2}\\!+\\!\Delta_{0}^{2}}\sqrt{\omega^{2}_{n}\\!+\\!\Delta_{0}^{2}}]}{\sqrt{(\omega_{n}\\!+\\!2\Omega)^{2}\\!+\\!\Delta_{0}^{2}}\\!+\\!\sqrt{\omega^{2}_{n}\\!+\\!\Delta_{0}^{2}}},$
(43) $\lambda_{E}\\!=\\!-T\sum_{\omega_{n}}\frac{{\rm
Tr}[i(\omega_{n}\\!\\!+\\!2\Omega)\tau_{1}{\hat{\Xi}}\\!-\\!i\omega_{n}{\hat{\Xi}}\tau_{1}\\!+\\!2\Delta_{0}{\hat{\Xi}}]}{\Delta_{0}[(\omega_{n}\\!+\\!2\Omega)^{2}\\!-\\!\omega^{2}_{n}]}\\!=\\!-\frac{T}{2\Omega^{2}}\sum_{\omega_{n}}\Big{[}\frac{2}{\sqrt{(\omega_{n}\\!+\\!\Omega)^{2}\\!+\\!\Delta_{0}^{2}}}\\!-\\!\frac{1}{\sqrt{(\omega_{n}\\!+\\!2\Omega)^{2}\\!+\\!\Delta_{0}^{2}}}\\!-\\!\frac{1}{\sqrt{\omega_{n}^{2}\\!+\\!\Delta_{0}^{2}}}\Big{]}.$
(44)

It is noted that Eq. (42) is exactly same as the one obtained in the previous
work by SilaevSilaev . Nevertheless, in Ref. Silaev, , the amplitude-response
coefficient $\lambda_{E}$ is directly considered to disappear after the
summation over the Matsubara frequency, leading to a zero Higgs-mode
generation. However, in contrast to the discrete $\omega_{n}/(\pi{T})=2n+1$,
the optical frequency $\Omega$ must be continuous in this circumstance (refer
to Sec. IV). Considering this point, the amplitude-response coefficient
$\lambda_{E}$ does not vanish. In fact, using the fact:

$\frac{1}{\sqrt{\Delta_{0}^{2}\\!+\\!(\omega_{n}\\!+\\!m\Omega)^{2}}}=\\!\\!\int\frac{d\xi_{\bf
k}}{\pi}\frac{1}{E_{\bf k}^{2}\\!-\\!(i\omega_{n}\\!+\\!im\Omega)^{2}}$ (45)

and considering low-frequency regime ($\Omega<E_{\bf k}$), after the standard
Matsubara-frequency summations, one can find a nonzero amplitude-response
coefficient

$\displaystyle\\!\\!\\!\\!\\!\\!\lambda_{E}=-\frac{1}{2\Omega^{2}}\frac{1}{\pi}\int{d\xi_{\bf
k}}\sum_{\eta=\pm}\Big{[}\frac{f({\eta}E_{\bf k})}{2{\eta}E_{\bf
k}}+\frac{f({\eta}E_{\bf k}-2i\Omega)}{2{\eta}E_{\bf k}}$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}-\frac{2f({\eta}E_{\bf
k}\\!-\\!i\Omega)}{2{\eta}E_{\bf
k}}\Big{]}\\!\approx\\!\frac{1}{\pi}\\!\\!\int\\!\\!{d\xi_{\bf
k}}\Big{[}\frac{\partial_{E_{\bf k}}^{2}f(E_{\bf k})}{2E_{\bf
k}}\\!+\\!O\Big{(}\frac{\Omega^{3}}{E^{3}_{\bf k}}\Big{)}\Big{]}.~{}~{}~{}~{}$
(46)

Consequently, a finite second-order response of the Higgs mode, contributed by
the drive effect of the vector potential, is achieved at clean limit by using
Eilenberger equation, in contrast to the previous work by SilaevSilaev .
Actually, according to the analysis of the Ginzburg-Landau theory as mentioned
in the introduction, this finite second-order response of the Higgs mode
derived from the Eilenberger equation is expected, since the Eilenberger
equation near $T_{c}$ can recover the Ginzburg-Landau equationBa20 .
Particularly, at low frequency, the amplitude-correlation coefficient

$\displaystyle\beta_{E}$ $\displaystyle\approx$
$\displaystyle\\!\frac{T}{2}\\!\sum_{\omega_{n}}\frac{1}{(\Delta_{0}^{2}+\omega^{2}_{n})^{3/2}}\\!=\\!-\frac{\partial_{\Delta_{0}}}{2\Delta_{0}}\Big{[}T\\!\sum_{\omega_{n}}\frac{1}{\sqrt{\Delta_{0}^{2}+\omega^{2}_{n}}}\Big{]}$
(47) $\displaystyle=$ $\displaystyle\\!2\int\frac{d\xi_{\bf
k}}{\pi}\frac{\partial_{E_{\bf k}}}{4E_{\bf k}}\Big{[}\frac{2f(E_{\bf
k})\\!-\\!1}{2E_{\bf k}}\Big{]}.$

Then, by comparing Eqs. (32)-(33) and Eqs. (47)-(46), one has
$\beta_{E}=\beta_{g}/(2D\pi)$ and $\lambda_{E}=\lambda_{g}/(2D\pi)$.
Therefore, for real optical frequency ($i\Omega\rightarrow\Omega$), Eq. (42)
derived from the Eilenberger equation exactly recovers Eq. (31) derived from
the GIKE, and hence, near $T_{c}$, can also recover Eq. (2) derived from the
Ginzburg-Landau theory.

Nevertheless, it is noted from Eq. (36) that the quasiclassical Eilenberger
equationEilen ; Ba20 ; Eilen1  only involves the drive effect of the vector
potential, i.e., the first term of $H^{(1)}_{\rm LM}$ in Eq. (16). Whereas the
drive effect of the scalar potential, i.e., the second term of $H^{(1)}_{\rm
LM}$ in Eq. (16), is hard to handle in the quasiclassical formalism due to its
spatial dependence (${\hat{\bf x}}\rightarrow{-i{\hat{\partial}}_{\bf k}}$).
Hence, the finite contribution from the drive effect of the scalar potential
to the Higgs-mode generation, which is nonzero at $T=0$, is generically
dropped out in this approach. Furthermore, even in consideration of the vector
potential alone, the density-related pump effect $H_{p}$ is generically
dropped out in the quasiclassical Eilenberger equation. Accordingly, the
response of the density-related phase mode, i.e., the second term of
$H^{(2)}_{\rm LM}$ in Eq. (17), as well as the Hartree field that is related
to the charge-density fluctuation and long-range Coulomb interaction, are also
dropped out. In fact, these deficiencies are because that the Eilenberger
equation in Eq. (36) is not gauge-invariant, and hence, the contained
electromagnetic effect is incomplete.

###  III.3 Path-integral approach

Following the previous work by Cea et al.Cea1 ; Cea2 ; Cea3 , we next use the
path-integral approach to derive the second-order optical response of the
Higgs mode. We start with the generalized action of superconductors in the
presence of electromagnetic potential $A_{\mu}$ Ba0 ; G1 :

$\displaystyle\\!\\!\\!\\!\\!\\!S[\psi,\psi^{*}]=\int{dx}\Big{[}\sum_{s=\uparrow,\downarrow}\psi^{*}_{s}(x)(i\partial_{t}-\xi_{\hat{\bf
p}-e{\bf A}}-e\phi)\psi_{s}(x)$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{}+\\!U\psi^{*}_{\uparrow}(x)\psi^{*}_{\downarrow}(x)\psi_{\downarrow}(x)\psi_{\uparrow}(x)\\!-\\!\frac{1}{2}\\!\\!\int\\!{dx^{\prime}}V(x\\!-\\!x^{\prime})n(x)n(x^{\prime})\Big{]},$

differing from the one used in Refs. Cea1, ; Cea2, ; Cea3,  with the vector
potential alone. Here, the density
$n(x)=\sum_{s=\uparrow,\downarrow}\psi^{*}_{s}(x)\psi_{s}(x)$. After the
Hubbard-Stratonovich transformation, one has

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

Here, $\mu_{H}$ stands for the auxiliary field, i.e., the Hartree field that
reflects the density fluctuation; $V_{Q}$ denotes the Fourier component of
$V(x-x^{\prime})$. The action in Eq. (49) satisfies the gauge structure in
Eqs. (8)-(9) revealed by Nambugi0  and hence is gauge invariant.

By further using the unitary transformation in Eq. (4) to effectively remove
the phase mode from the order parameter, the action in Eq. (49) becomes

$\displaystyle S[\psi,\psi^{*}]$ $\displaystyle=$
$\displaystyle\int{dx}\Big{\\{}\psi^{*}(x)[G_{0}^{-1}(x)\\!-\\!\Sigma(x)]\psi(x)\\!-\\!\frac{|\Delta(x)|^{2}}{U}$
(50) $\displaystyle\mbox{}-\eta_{f}{\rm
Tr}[\Sigma(x)\tau_{3}/2]\Big{\\}}+\frac{1}{2}\sum_{Q}\frac{|\mu_{H}(Q)|^{2}}{V_{Q}}.$

where $\eta_{f}=\sum_{\bf k}1$ emerges due to the anti-commutation of the
Fermi field; the Green function $G_{0}^{-1}(x)=i\partial_{t}-H_{0}$, which in
frequency-momentum space [$x\rightarrow{p=(p_{0},{\bf k})}$] reads
$G_{0}(p)=(p_{0}+\xi_{\bf k}\tau_{3}+\Delta_{0}\tau_{1})/(p_{0}^{2}-E_{\bf
k}^{2})$ and the self-energy $\Sigma(x)=H_{\rm
LM}+\delta|\Delta|\tau_{1}+\mu_{H}\tau_{3}$.

After the standard integration over the Fermi field, one has $S=S_{0}+{S}_{\rm
ne}[A_{\mu},\delta|\Delta|,\delta\theta]$, consisting of the equilibrium part
$S_{0}$ and the non-equilibrium one:

$\displaystyle\\!\\!\\!\\!\\!S_{\rm
ne}[A_{\mu},\delta|\Delta|,\delta\theta]\\!=\\!-\sum_{n=1}^{\infty}\frac{1}{n}{\rm{\bar{T}r}}[(G_{0}\Sigma)^{n}]\\!+\\!\frac{1}{2}\sum_{Q}\frac{|\mu_{H}(Q)|^{2}}{V_{Q}}$
$\displaystyle\\!\\!\\!\\!\mbox{}-\\!\\!\int{dx}\eta_{f}{\rm
Tr}[\Sigma(x)\tau_{3}/2]-\\!\int{dx}\frac{(\delta|\Delta|)^{2}+2\Delta_{0}\delta|\Delta|}{U}.$
(51)

In the non-equilibrium $S_{\rm ne}[A_{\mu},\delta|\Delta|,\delta\theta]$, the
anisotropic linear and third orders with respect to the electromagnetic
potential vanish after the angular integration of momentum, whereas the
second-order part that corresponds to the linear current excitation has been
well established in the literaturepi1 ; pi4 . To discuss the experimentally
observed third-harmonic current, one needs to formulate the expansion of the
action with respect to the fourth order of the electromagnetic potential, and
hence, keeps the expansions up to $n=4$. For the convenience of the
derivation, we consider the spatially uniform transverse optical fields and
choose zero phase-related effective electromagnetic potential
$\partial_{\mu}\delta\theta^{(l)}$. Then, following the analysis of the light-
matter interaction $H_{\rm LM}$ in Sec. II, the related action with expansions
up to $n=4$ is written as

$\displaystyle S\big{[}A^{4}_{\mu}\big{]}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}{\rm{\bar{T}r}}[G_{0}(H^{(2)}_{\rm
LM}\\!+\\!\delta|\Delta|\tau_{1}\\!+\\!\mu_{H}\tau_{3})G_{0}(H^{(2)}_{\rm
LM}\\!+\\!\delta|\Delta|\tau_{1}\\!+\\!\mu_{\rm
H}\tau_{3})]\\!-\\!{\rm{\bar{T}r}}[G_{0}H_{\rm LM}^{(1)}G_{0}H_{\rm
LM}^{(1)}G_{0}(H^{(2)}_{\rm LM}\\!+\\!\delta|\Delta|\tau_{1}\\!+\\!\mu_{\rm
H}\tau_{3})]$ (52)
$\displaystyle\mbox{}\\!-\\!\frac{1}{4}{\rm{\bar{T}r}}[G_{0}H_{\rm
LM}^{(1)}G_{0}H_{\rm LM}^{(1)}G_{0}H_{\rm LM}^{(1)}G_{0}H_{\rm
LM}^{(1)}]\\!-\\!\sum_{Q}\frac{|\delta|\Delta|_{Q}|^{2}}{U}\\!+\\!\frac{1}{2}\sum_{Q}\frac{|\mu_{H}(Q)|^{2}}{V_{Q}}$
$\displaystyle=$
$\displaystyle-\\!\sum_{Q}\Big{[}\chi_{H}|\delta|\Delta|_{2Q}|^{2}\\!+\\!\chi_{33}\Big{|}\frac{e^{2}A_{0}^{2}}{2m}\\!+\\!\mu_{H}\\!+\\!iQ_{0}\delta\theta^{(2)}\Big{|}^{2}\\!-\\!\frac{|\mu_{H}|^{2}}{2V_{Q}}\Big{]}\\!-\\!\sum_{Q}\Bigg{\\{}\\!2\chi_{13}\Big{(}\frac{e^{2}A_{0}^{2}}{2m}\\!+\\!\mu_{H}\\!+\\!iQ_{0}\delta\theta^{(2)}\Big{)}_{2Q}\delta|\Delta|_{-2Q}$
$\displaystyle\mbox{}\\!+\\!\Big{[}\chi_{001}\frac{e^{2}A^{2}_{0}v_{F}^{2}}{3}\\!+\\!(\chi_{0{\bar{3}}1}-\chi_{3{\bar{0}}1})\frac{eA_{0}eE_{\phi}v_{F}^{2}}{3i}+\chi_{3{\bar{\bar{3}}}1}\frac{e^{2}E_{\phi}^{2}v_{F}^{2}}{3}\Big{]}_{2Q}\delta\Delta_{-2Q}\\!+\\!\Big{(}\frac{e^{2}A_{0}^{2}}{2m}\\!+\\!\mu_{H}\\!+\\!iQ_{0}\delta\theta^{(2)}\Big{)}_{-2Q}$
$\displaystyle\mbox{}\\!\times\\!\Big{[}\chi_{003}\frac{e^{2}A^{2}_{0}v_{F}^{2}}{3}\\!+\\!(\chi_{0{\bar{3}}3}-\chi_{3{\bar{0}}3})\frac{eA_{0}eE_{\phi}v_{F}^{2}}{3i}+\chi_{3{\bar{\bar{3}}}3}\frac{e^{2}E_{\phi}^{2}v_{F}^{2}}{3}\\!+\\!(\chi_{0{3}3}-\chi_{303})\frac{eA_{0}eE_{\phi}}{3mi}+\chi_{3{\bar{3}}3}\frac{e^{2}E_{\phi}^{2}}{3m}\Big{]}_{-2Q}\\!+\\!h.c.\Bigg{\\}}$
$\displaystyle\mbox{}\\!-\\!\frac{e^{4}v_{F}^{4}}{20}{\rm{\bar{T}r}}\\{[G_{0}(A_{0}-E_{\phi}\tau_{3}\partial_{\xi_{\bf
k}})]^{4}\\},$

in which we have considered a large $v_{F}$ (i.e., neglected the terms
proportional to $v_{F}^{2}$ and $v_{F}^{0}$ and only kept the ones
proportional to $v_{F}^{4}$) in the expansion of $n=4$. Here, the frequency-
momentum vector $Q=(Q_{0},{\bf Q})$; $\chi_{H}=\chi_{11}+1/U$ denotes the
energy-spectrum function of the Higgs mode; the correlation coefficients are
written as

$\displaystyle\chi_{ij}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\frac{1}{2}\sum_{p}{\rm
Tr}[G_{0}(p+2Q)\tau_{i}G_{0}(p)\tau_{j}],$ (53)
$\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}],$ (54)
$\displaystyle\chi_{i{\bar{j}}k}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{p}{\rm Tr}[G_{0}(p+2Q)\tau_{i}\partial_{\xi_{\bf
k}}G_{0}(p+Q)\tau_{j}G_{0}(p)\tau_{k}],$ (55)
$\displaystyle\chi_{i{\bar{\bar{3}}}k}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\frac{1}{3}\sum_{p}{\rm
Tr}[G_{0}(p\\!+\\!2Q)\tau_{i}{\partial^{2}_{\xi_{\bf
k}}}G_{0}(p\\!+\\!Q)\tau_{3}G_{0}(p)\tau_{k}].~{}~{}~{}~{}~{}~{}$ (56)

It is noted that the action in Eq. (52) exactly recovers the one in the
previous work by Cea et al.Cea1 ; Cea2 ; Cea3 , if one only keeps the second-
order perturbation (i.e., $n=2$) expansion and neglects the third- and forth-
order perturbation (i.e., $n=3$ and $n=4$) expansions. As revealed in Refs.
Cea1, ; Cea2, ; Cea3, , the second-order correlations $\chi_{11}$ and
$\chi_{33}$ characterize the amplitude-amplitude and density-density
correlations, respectively. The density-amplitude correlation $\chi_{13}$ is
zero as a consequence of the particle-hole symmetry, and hence, the only
coupling between the Higgs mode and second-order optical field in the second-
order perturbation expansion, i.e., the coupling between the Higgs mode and
pump effect, vanishes.

As for the third-order correlations, one can prove that $\chi_{003}$,
$\chi_{0{\bar{3}}3}$, $\chi_{3{\bar{0}}3}$ and $\chi_{3{\bar{\bar{3}}}3}$
vanish as a consequence of the particle-hole symmetry, and
$\chi_{0{\bar{3}}1}-\chi_{3{\bar{0}}1}=0$ (refer to Appendix A). Then, the
action in Eq. (52) is simplified as

$\displaystyle S\big{[}A^{4}_{\mu}\big{]}$ $\displaystyle=$
$\displaystyle-\\!\sum_{Q}\Big{\\{}\chi_{H}|\delta|\Delta|_{2Q}|^{2}\\!+\\!\Big{[}\Big{(}\chi_{001}\frac{e^{2}A^{2}_{0}v_{F}^{2}}{3}\\!+\\!\chi_{3{\bar{\bar{3}}}1}\frac{e^{2}E_{\phi}^{2}v_{F}^{2}}{3}\Big{)}_{2Q}\delta\Delta_{-2Q}+h.c.\Big{]}\Big{\\}}$
(57)
$\displaystyle\mbox{}\\!-\\!\sum_{Q}\Big{(}\chi_{33}\Big{|}\frac{e^{2}A_{0}^{2}}{2m}\\!+\\!\mu_{H}\\!+\\!iQ_{0}\delta\theta^{(2)}\\!-\\!\frac{\zeta_{1}eA_{0}eE_{\phi}}{3miQ_{0}}\\!+\\!\frac{\zeta_{2}e^{2}E_{\phi}^{2}}{3mQ_{0}^{2}}\Big{|}^{2}\\!+\\!\frac{|\mu_{H}|^{2}}{2V_{Q}}\Big{)}\\!-\\!\frac{e^{4}v_{F}^{4}}{20}{\rm{\bar{T}r}}\\{[G_{0}(A_{0}-E_{\phi}\tau_{3}\partial_{\xi_{\bf
k}})]^{4}\\},~{}~{}~{}$

with $\zeta_{1}=Q_{0}\frac{\chi_{0{3}3}-\chi_{303}}{\chi_{33}}$ and
$\zeta_{2}=Q_{0}^{2}\frac{\chi_{3{\bar{3}}3}}{\chi_{33}}$. In Eq. (57), the
first term is related to the Higgs mode, and the second one is related to the
phase mode and charge-density fluctuation. Whereas the third one denotes the
fourth order of the drive effect, which is related to the thermal effect.

####  III.3.1 Higgs-mode generation

By using the action in Eq. (57), we discuss the Higgs-mode generation in the
second-order optical response and its contribution to the third-harmonic
current. Considering the optical response
$A_{\mu}(Q)\rightarrow{A_{\mu}}\delta(Q_{0}-\Omega)\delta({\bf q})$, one has
$\chi_{H}(2Q)=(4\Delta_{0}^{2}-4\Omega^{2})\beta_{p}$ and
$\chi_{001}=2\Delta_{0}\lambda_{p}$ as well as
$\chi_{3{\bar{\bar{3}}}1}=\frac{2\Delta_{0}}{(i\Omega)^{2}}\gamma_{p}$, with
the coefficients $\beta_{p}$ and $\lambda_{p}$ as well as $\gamma_{p}$ at low
frequency given by (refer to Appendix A)

$\displaystyle\beta_{p}$ $\displaystyle\approx$ $\displaystyle D\int{d\xi_{\bf
k}}\frac{\partial_{E_{\bf k}}}{2E_{\bf k}}\Big{[}\frac{2f(E_{\bf
k})\\!-\\!1}{2E_{\bf k}}\Big{]},$ (58) $\displaystyle\lambda_{p}$
$\displaystyle=$ $\displaystyle\frac{1}{2Q_{0}^{2}}\sum_{\bf
k}\Big{\\{}\Big{[}\frac{f(E_{\bf k}\\!-\\!2Q_{0})}{2E_{\bf
k}}\\!+\\!\frac{f(E_{\bf k})}{2E_{\bf k}}\\!-\\!\frac{2f(E_{\bf
k}\\!-\\!Q_{0})}{2E_{\bf k}}\Big{]}$ (59) $\displaystyle\mbox{}+\\{E_{\bf
k}\rightarrow-E_{\bf k}\\}\Big{\\}},$ $\displaystyle\gamma_{p}$
$\displaystyle\approx$ $\displaystyle D\int{d\xi_{\bf k}}\frac{\xi_{\bf
k}^{4}}{E_{\bf k}^{5}}\partial_{E_{\bf k}}\Big{[}\frac{2f(E_{\bf
k})-1}{2E_{\bf k}}\Big{]}.$ (60)

Then, from Eq. (57), the equation of motion of the Higgs mode reads

$4(\Delta_{0}^{2}\\!-\\!\Omega^{2})\beta_{p}\delta|\Delta|^{(2)}=-2\Delta_{0}\frac{e^{2}v_{F}^{2}}{3}\Big{[}\lambda_{p}{A^{2}_{0}}\\!+\\!\frac{\gamma_{p}E^{2}_{\phi}}{(i\Omega)^{2}}\Big{]}.$
(61)

Consequently, one immediately finds a finite second-order response of the
Higgs mode from the path-integral approach at clean limit, contributed by both
drive effects of the vector and scalar potentials. Particularly, at low
frequency, by comparing Eqs. (32)-(33) and Eqs. (58)-(59), one finds
$\beta_{p}=\beta_{g}/2$ and

$\lambda_{p}~{}{\approx}~{}D\int{d\xi_{\bf k}}\frac{\partial_{E_{\bf
k}}^{2}f(E_{\bf k})}{2E_{\bf k}}=\lambda_{g}/2.$ (62)

In addition, from Eq. (26) at low frequency, one has

$\displaystyle\gamma_{g}$ $\displaystyle\approx$
$\displaystyle{D}\int{d\xi_{\bf k}}\frac{\xi_{\bf k}}{E_{\bf
k}^{2}}\partial^{2}_{\xi_{\bf k}}\Big{[}{\xi_{\bf k}}\frac{2f(E_{\bf
k})-1}{2E_{\bf k}}\Big{]}$ (63) $\displaystyle=$ $\displaystyle
D\int{d\xi_{\bf k}}\Big{(}\frac{3\xi_{\bf k}^{2}}{E_{\bf
k}^{3}}\\!-\\!\frac{\xi_{\bf k}^{4}}{E_{\bf k}^{5}}\\!+\\!\frac{\xi_{\bf
k}^{3}}{E_{\bf k}^{3}}\partial_{\xi_{\bf k}}\Big{)}\partial_{E_{\bf
k}}\Big{[}\frac{2f(E_{\bf k})\\!-\\!1}{2E_{\bf k}}\Big{]}$ $\displaystyle=$
$\displaystyle D\int{d\xi_{\bf k}}\Big{[}\frac{3\xi_{\bf k}^{2}}{E_{\bf
k}^{3}}\\!-\\!\frac{\xi_{\bf k}^{4}}{E_{\bf k}^{5}}\\!-\\!\partial_{\xi_{\bf
k}}\Big{(}\frac{\xi_{\bf k}^{3}}{E_{\bf k}^{3}}\Big{)}\Big{]}\partial_{E_{\bf
k}}\Big{[}\frac{2f(E_{\bf k})\\!-\\!1}{2E_{\bf k}}\Big{]}$ $\displaystyle=$
$\displaystyle 2\gamma_{p}.$

Consequently, Eq. (61) derived from the path-integral approach exactly
recovers Eq. (31) derived from the GIKE, with both contributions of the vector
and scalar potentials incorporated.

The obtained finite second-order response of the Higgs mode through the
vector-potential drive effect within the path-integral approach is quite
natural, since this approach near $T_{c}$ can also recover the Ginzburg-Landau
equationPI2GL . Near $T_{c}$, Eq. (61) with the vector potential alone
recovers Eq. (2) derived from the Ginzburg-Landau theory. Nevertheless, in the
previous works by Cea et al.Cea1 ; Cea2 ; Cea3  with the vector potential
alone, the expansions of the action for $n>2$ are excessively overlooked, and
then, only the coupling $\chi_{13}H_{p}\delta|\Delta|$ between the pump effect
and Higgs mode in $n=2$ perturbation expansion is left, leading to a zero
Higgs-mode generation at clean limit due to $\chi_{13}=0$. But in fact, in
$n=3$ expansion, there exists the essential coupling
$\chi_{001}H_{d}^{2}\delta|\Delta|$ between the second order of the vector-
potential drive effect and the Higgs mode, which leads to the finite Higgs-
mode generation in Eq. (61). Previous calculations in Refs. Cea1, ; Cea2, ;
Cea3,  that overlooked $n=3$ expansion therefore missed the finite Higgs-mode
generation in the second-order response at clean limit.

It is also noted that the previous worksCea1 ; Cea2 ; Cea3  with the vector
potential alone overlooked the drive effect of the scalar potential. Whereas
the finite contribution from this effect to the Higgs-mode generation, which
has obtained from GIKE in Eq. (31), is exactly recovered here in Eq. (61), as
we handle the path-integral approach within the gauge-invariant manner in the
present work. As mentioned in Sec. III.1, this contribution is finite upon
cooling to $T=0$, in contrast to the one of the vector potential which that
emerges only at finite temperatureDS1 . Whereas because of the gauge structure
in superconductors, the consideration of the contribution from the scalar
potential is essential here as we have chosen zero
$\partial_{\mu}\delta\theta^{(1)}$.

From the action in Eq. (57), after the integration out the Higgs mode, one
obtains the finite contribution from the Higgs mode to the third-harmonic
current:

$S[A^{4}_{\mu}]|_{\rm
Higgs}=4\Delta^{2}_{0}\frac{e^{4}v_{F}^{4}}{9}\frac{\Big{[}\lambda_{p}{A^{2}_{0}}\\!+\\!\gamma_{p}\Big{(}\frac{E_{\phi}}{i\Omega}\Big{)}^{2}\Big{]}^{2}}{(4\Delta_{0}^{2}\\!-\\!4\Omega^{2})\beta_{p}},$
(64)

which shows an isotropic signal and exhibits a resonance when
$2\Omega=2\Delta_{0}$, in consistency with the experimental findingsNL8 ; NL9
; NL10 .

####  III.3.2 Vanishing charge-density fluctuation

By using the action in Eq. (57), we next discuss the phase mode and charge-
density fluctuation in the second-order optical response and its contribution
to the third-harmonic current. In the action, for the part that is related to
the phase mode and charge-density fluctuation [i.e., the second term in Eq.
(57)], after the integration out the phase mode, one has

$\displaystyle\\!\\!\\!\\!\\!\\!S[A^{4}_{\mu}]|_{\rm
phase}=-\sum_{Q}\Big{[}\Big{|}\frac{e^{2}A_{0}^{2}}{2m}\\!+\\!\mu_{H}\\!-\\!\frac{\zeta_{1}eA_{0}eE_{\phi}}{3mi\Omega}\\!+\\!\frac{\zeta_{2}e^{2}E_{\phi}^{2}}{3m\Omega^{2}}\Big{|}^{2}$
$\displaystyle\mbox{}\\!\\!\\!\\!\\!\times\Big{(}\chi_{33}-\frac{\chi_{33}^{2}}{\chi_{33}}\Big{)}\\!+\\!\frac{|\mu_{H}|^{2}}{2V_{Q}}\Big{]}=-\sum_{Q}\frac{|\mu_{H}|^{2}}{2V_{Q}}.$
(65)

In above action, there is no coupling term between the charge-density
fluctuation and second-order optical fields. Consequently, there is no
generation of the charge-density fluctuation in the second-order response to
contribute to the third-harmonic current, exactly same as the result [Eq.
(24)] from the GIKEGIKE2  and in consistency with the analysis based on the
inversion symmetry and charge conservation as mentioned in the introduction.

Particularly, with $\mu_{H}=0$, from the action in Eq. (57), the equation of
motion of the phase mode is written as

$-i\Omega\delta\theta^{(2)}=\frac{e^{2}A_{0}^{2}}{2m}\\!-\\!\frac{\zeta_{1}eA_{0}eE_{\phi}}{3mi\Omega}\\!+\\!\frac{\zeta_{2}e^{2}E_{\phi}^{2}}{3m\Omega^{2}}.$
(66)

On the right-hand side of above equation, the pump effect $H_{p}$ of the
vector potential (the first term) from $n=2$ perturbation expansion and the
second order of the drive effect $H_{\rm LM}^{(1)}$ (the last two terms) from
$n=3$ perturbation expansion provide the source terms. Then, one finds a
finite second-order response of the phase mode, contributed by both pump and
drive effects, in agreement with the result [Eq. (34)] from the GIKEGIKE2 .

It is noted that in the previous works by Cea et al.Cea1 ; Cea2 ; Cea3 , the
excitation of the phase mode in the second-order response is overlooked. Then,
from the second term in Eq. (57), one obtains a finite coupling between the
second-order optical fields and the charge-density fluctuation $\mu_{H}$,
leading to the finite charge-density fluctuation in the second-order response.
Whereas this result violates the analysis based on the inversion symmetry and
charge conservation mentioned in the introduction. In fact, by calculating the
amplitude and phase modes on an equal footing through the path-integral
approach, the present work shows that the phase mode $\delta\theta^{(2)}$ is
excited in the second-order response. Substituting the phase-mode generation
in Eq. (66) into the action in Eq. (57), one can immediately finds that the
charge-density-fluctuation part becomes

$S[A^{4}_{\mu}]|_{\rm
CDF}=-\sum_{Q}\Big{(}\chi_{33}|\mu_{H}|^{2}\\!+\\!\frac{|\mu_{H}|^{2}}{2V_{Q}}\Big{)}.$
(67)

Consequently, the phase-mode generation $\delta\theta^{(2)}$ provides an
effective field to exactly cancel the unphysical excitation of the charge-
density fluctuation reported by Cea et al.Cea1 ; Cea2 ; Cea3 , guaranteeing
the charge conservation.

##  IV Discussion on Matsubara formalism in derivation of superconducting gap
dynamics

In this part, we address a specail issue in the application of Matsubara
formalism in superconductors. We show that in the derivation of the
superconducting gap dynamics, treating the optical frequency $i\Omega$ as
bosonic Matsubara frequencies leads to results against Ginzburg-Landau
equation. One has to take $\Omega$ as continuous variable in order to
recover/derive Ginzburg-Landau equation.

Specifically, with the vector potential alone, either from the
phenomenological time-dependent Ginzburg-Landau superconducting Lagrangian
[Eq. (1)] or through the microscopic Eilenberger equation (Sec. III.2) and
path-integral (Sec. III.3) as well as gauge-invariant kinetic equation (Sec.
III.1) approaches, the Higgs-mode dynamics in the second-order response reads:

$\beta_{H}(4\Delta_{0}^{2}-4\Omega^{2})\delta|\Delta|=-2\Delta_{0}\frac{{\lambda}k_{F}^{2}}{3m}\frac{e^{2}A^{2}}{m},$
(68)

where the response coefficients $\lambda$ from the Ginzburg-Landau theory,
Eilenberger equation and path-integral as well as gauge-invariant kinetic
equation approaches are given by
$\lambda_{l}=2m\lambda_{L}/k_{F}^{2}={7R(3)}/{(2\pi{T})^{2}}$ as well as
$\lambda_{E}$ [Eq. (44)], $\lambda_{p}$ [Eq. (59)] and $\lambda_{g}$ [Eq.
(27)], respectively. The Landau parameter $\lambda_{L}$ was derived by Gorkov
near $T_{c}$ from basic Gorkov equationG1 . Whereas as mentioned in Secs.
III.2 and III.3, $\lambda_{E}$ in Eq. (44) has been exactly derived by Silaev
in Ref. Silaev, , but $\lambda_{p}$ was directly missed in the previous works
by Cea et al.Cea1 ; Cea2 ; Cea3 .

Particularly, using the fact in Eq. (45), one finds that $\lambda_{E}$ [Eq.
(44)] from Eilenberger equation is directly equivalent to $\lambda_{p}$ [Eq.
(59)] from path-integral approach in Matsubara formalism. However, in Ref.
Silaev, , $\lambda_{E}$ derived from Eilenberger equation is considered as
zero by taking the optical frequency $i\Omega$ as bosonic Matsubara
frequencies $i\Omega_{m}$, and for bosonic Matsubara frequencies
$i\Omega_{m}$, $\lambda_{p}$ in Eq. (59) derived from path-integral approach
also vanishes. This treatment about the optical frequency is indeed the
conventional one applied in conductivity calculation of normal metals.
Nevertheless, in the derivation of the superconducting gap dynamics here, the
vanishing $\lambda_{E}$ and $\lambda_{p}$ are strongly against the finite
$\lambda_{l}$ from Ginzburg-Landau theory. Actually, it can be easily
demonstrated that both the Eilenberger equationBa20  and path-integral
approach at the stationary case ($\Omega=0$) can exactly recover the Ginzburg-
Landau equation and derive the Ginzburg-Landau kinetic term
$\frac{\lambda_{L}({\nabla}-2ie{\bf A})^{2}\Delta}{4m}$ (the detailed
derivation is given in Appendixes B and C for the sake of completeness). In
other words, due to the treatment of taking $i\Omega$ as bosonic Matsubara
frequencies, which eliminates the response coefficient $\lambda$ and hence
Ginzburg-Landau kinetic term at finite $\Omega$, an unphysical abrupt change
between $\Omega=0$ and $\Omega\rightarrow 0$ emerges. This demonstrates that
the application of Matsubara formalism in superconductors should be carefully
examined, since the treatment that fails to recover the Ginzburg-Landau
equation in conventional superconductors can not be correct.

In fact, only by taking $\Omega$ as continuous variable in this circumstance,
can one recover/derive Ginzburg-Landau equation. Similarly, only with the
continuous $\Omega$, the coefficients $\lambda_{E}$ [Eq. (44)] from
Eilenberger equation and $\lambda_{p}$ [Eq. (59)] from path-integral approach
can exactly recover the finite $\lambda_{g}$ from gauge-invariant kinetic
equation and $\lambda_{l}$ from Ginzburg-Landau theory at low frequency, as
demonstrated in Secs. III.2 and III.3. It is noted that the gauge-invariant
kinetic equation is developed by using Keldysh Green function approach,
irrelevant to the Matsubara formalism, and the coefficient $\lambda_{g}$ [Eq.
(27)] from this approach can naturally recover the one from Ginzburg-Landau
equationGIKE1 .

The finite second-order response of Higgs mode can also be understood as
follows. Specifically, with the vector potential alone at low frequency, it is
established that the vector potential drives the Doppler shiftFF4 ; FF5 ; FF6
to influence the gap equation:

${\Delta}=g\\!\sum_{\bf k}\Delta\frac{f(-{\bf v}_{\bf k}\cdot{e}{\bf A}-E_{\bf
k})-f(-{\bf v}_{\bf k}\cdot{e}{\bf A}+E_{\bf k})}{2E_{\bf k}},$ (69)

which can be directly derived according to the Hamiltonian in Eq. (3). Then,
considering the gap variation and weak field, the above equation directly
becomes

$4\Delta_{0}^{2}\delta|\Delta|\sum_{\bf k}\frac{\partial_{E_{\bf
k}}[\frac{2f(E_{\bf k})-1}{2E_{\bf k}}]}{2E_{\bf
k}}\\!=\\!-2\Delta_{0}\frac{v_{F}^{2}e^{2}A^{2}\sum_{\bf
k}\frac{\partial^{2}_{E_{\bf k}}f(E_{\bf k})}{2E_{\bf k}}}{3},$ (70)

which is exactly same as the ones derived from gauge-invariant kinetic
equation [Eq. (31)] and path-integral approach [Eq. (61)] at low continuous
optical frequency.

##  V Summary

In summary, through three different microscopic approaches including the GIKE,
Eilenberger equation as well as path-integral approach, the present study
arrives at unified conclusion about the finite Higgs-mode generation in the
second-order optical response of superconductors at clean limit, in
consistency with the phenomenological Ginzburg-Landau theory. Moreover, while
the density-related effect is hard to handle in the quasiclassical Eilenberger
equation, the vanishing charge-density fluctuation in the second-order
response, which agrees with the charge conservation, is obtained within the
GIKE and path-integral approach. Consequently, the present work resolves the
controversies among various studies in the literature (whether the
experimentally observed third-harmonic currentNL7 ; NL8 ; NL9 ; NL10 ; NL11 ;
DHM2 ; DHM3  is attributed to the Higgs-modeNL10 ; NL11 ; DHM2 ; DHM3 ; FHM ;
GIKE2  or charge-density-fluctuationCea1 ; Cea2 ; Cea3  generation, and
whether one has to rely on impurity to explain the experimentally observed
Higgs-mode generationCea1 ; Cea2 ; Cea3 ; GIKE2 ; Silaev ; FHM ; Am16 ; ImR1 ;
ImR2 ; ImR3 ), and can therefore help understanding the experimental findings.

Specifically, with the vector potential alone, by separately using the GIKE as
well as Eilenberger equation and path-integral approach, we obtain the exactly
same finite Higgs-mode generation in the second-order optical response at
clean limit. This finite Higgs-mode generation is solely contributed by the
drive effect ${\bf p}\cdot{e{\bf A}}/m$ of the vector potential, and near
$T_{c}$, exactly recovers the one from the phenomenological Ginzburg-Landau
theory. Nevertheless, the previous works within the path-integral approach by
Cea et al.Cea1 ; Cea2 ; Cea3  and Eilenberger equation by SilaevSilaev  missed
this finite generation because of the flaws in their derivations. A
disscussion about the application of Matsubara formalism in the derivation of
superconducting gap dynamics is given, and it is demonstrated that taking the
optical frequency as continuous variable is essential to recover/derive the
Ginzburg-Landau equation.

According to the gauge structure of superconductorsgi0 , among the scalar
potential $\phi$, vector potential ${\bf A}$ as well as phase-related
effective electromagnetic potential $\partial_{\mu}\delta\theta$, one can not
choose two quantities simultaneously to be zero, e.g., considering the vector
potential alone. We therefore extend the path-integral approach to include
electromagnetic effects from the scalar potential and phase mode. Then, in the
second-order response at clean limit, a finite contribution in the Higgs-mode
generation from the drive effect of scalar potential as well as the vanishing
charge-density fluctuation, which have previously been obtained from GIKEGIKE2
, are exactly recovered.

The contribution of the scalar potential is finite upon cooling to $T=0$,
differing from the one of the vector potential that emerges only at finite
temperatureDS1 . This difference is due to the fact that the superconductors
can directly respond to vector potential ${\bf A}$ (Meissner effect/Ginzburg-
Landau kinetic term) in addition to the electric field ${\bf E}=-{\nabla}_{\bf
R}\phi-\partial_{t}{\bf A}$ (optical-electric-field effect), differing from
normal metals that solely respond to ${\bf E}$. Consequently, in contrast to
the scalar potential that only captures the optical-electric-field effect, the
vector potential also characterizes the Meissner effect/Ginzburg-Landau
kinetic term in addition to the optical-electric-field effect. Particularly,
because of the gauge structure of superconductorsgi0 , the inclusion of the
contribution from the scalar potential here is essential, since we have chosen
zero $\partial_{\mu}\delta\theta^{(l)}$ (i.e., a spatially uniform background
phase mode $\delta\theta^{(l)}$).

Although the uniform background phase mode $\delta\theta^{(l)}$ does not
manifest itself in the measurable optical properties, the inclusion of this
mode is essential for theoretical description to cancel the unphysical effect,
as pointed out by Nambu in his Nobel lecturegi1 . The linear response of the
background phase mode $\delta\theta^{(1)}$ cancels the unphysical longitudinal
vector potential in ${\bf p}_{s}$, as established in the literatureAK ; AK2 ;
Ba0 ; pm0 ; Am0 ; Ba9 ; Ba10 ; pm5 ; pi1 ; pi4 ; GIKE2 . Then, the
superconducting momentum ${\bf p}_{s}$ that appears in the previous
theoretical descriptions such as Ginzburg-Landau equationG1  and Meissner
supercurrentG1  as well as Anderson-pump effectAm1 ; Am2 ; Am3 ; Am4 ; Am7 ;
Am8 ; Am9 ; Am10 ; Am11 ; Am14 ; NL7 ; NL8 ; NL9 ; NL10 ; NL11  only involves
the physical transverse vector potential. The previous works by Cea et al.Cea1
; Cea2 ; Cea3  overlooked the phase mode, and obtained the excited charge-
density fluctuation in the second-order optical response of superconductors.
Whereas this result in systems with the inversion symmetry violates the law of
charge conservation. We show in the present work that a background phase mode
$\delta\theta^{(2)}$ is actually generated in the second-order response, and
exactly cancel the unphysical excitation of the charge-density fluctuation
reported in Refs. Cea1, ; Cea2, ; Cea3, , guaranteeing the charge
conservation.

###### 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, from Eqs. (53)-(56), we present the specific expressions of the
related correlation coefficients at low frequency ($Q_{0}<E_{\bf k}$) and
long-wave limit (${\bf Q}\rightarrow 0$):

$\displaystyle\chi_{11}\\!-\\!U^{-1}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\frac{1}{2}\sum_{p}{\rm
Tr}\Big{[}\frac{p_{0}\\!+\\!2Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}}\tau_{1}\frac{p_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf
k}}\tau_{1}\Big{]}\\!-\\!U^{-1}=\sum_{p}\frac{(p_{0}+2Q_{0})p_{0}+\Delta_{0}^{2}-\xi^{2}_{\bf
k}}{[(p_{0}+2Q_{0})^{2}-E^{2}_{\bf k}](p_{0}^{2}-E^{2}_{\bf k})}-U^{-1}$ (71)
$\displaystyle=$
$\displaystyle\\!\\!\sum_{p}\Big{\\{}\frac{4\Delta_{0}^{2}-(p_{0}+2Q_{0}-p_{0})^{2}}{2[(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}](p_{0}^{2}-E^{2}_{\bf k})}+\frac{1}{2(p_{0}^{2}-E^{2}_{\bf
k})}+\frac{1}{2[(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}]}\Big{\\}}-\sum_{p}\frac{1}{p_{0}^{2}-E^{2}_{\bf k}}$
$\displaystyle\approx$
$\displaystyle\\!\\!\sum_{p}\frac{4\Delta_{0}^{2}\\!-\\!4Q_{0}^{2}}{2(p_{0}^{2}-E^{2}_{\bf
k})^{2}}\\!=\\!4({\Delta_{0}^{2}\\!-\\!Q_{0}^{2}})\sum_{\bf k}\frac{1}{2E_{\bf
k}^{2}}\Big{[}\frac{1\\!-\\!2f(E_{\bf k})}{2E_{\bf k}}\\!+\\!\partial_{E_{\bf
k}}f(E_{\bf k})\Big{]}=4({\Delta_{0}^{2}\\!-\\!Q_{0}^{2}})\sum_{\bf
k}\frac{\partial_{E_{\bf k}}}{2E_{\bf k}}\Big{[}\frac{2f(E_{\bf
k})\\!-\\!1}{2E_{\bf k}}\Big{]},~{}~{}~{}~{}~{}$ $\displaystyle\chi_{33}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{p}{\rm
Tr}\Big{[}\frac{p_{0}\\!+\\!2Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}}\tau_{3}\frac{p_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf
k}}\tau_{3}\Big{]}=\sum_{p}\frac{(p_{0}\\!+\\!2Q_{0})p_{0}\\!+\\!\xi^{2}_{\bf
k}\\!-\\!\Delta_{0}^{2}}{[(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}](p_{0}^{2}-E^{2}_{\bf k})}$ (72) $\displaystyle=$ $\displaystyle\sum_{\bf
k}\Big{\\{}\Big{[}\frac{f(E_{\bf k})}{2E_{\bf k}}\frac{2\xi_{\bf
k}^{2}\\!+\\!2Q_{0}E_{\bf k}}{4Q_{0}(E_{\bf
k}\\!+\\!Q_{0})}\\!-\\!\frac{f(E_{\bf k}\\!-\\!2Q_{0})}{2E_{\bf
k}}\frac{2\xi_{\bf k}^{2}\\!-\\!2Q_{0}E_{\bf k}}{4Q_{0}(E_{\bf
k}\\!-\\!Q_{0})}\Big{]}\\!+\\![E_{\bf k}\\!\rightarrow\\!-E_{\bf
k}]\Big{\\}}\approx\sum_{\bf k}\frac{\Delta^{2}_{0}}{E_{\bf
k}^{2}\\!-\\!Q_{0}^{2}}\frac{2f(E_{\bf k})\\!-\\!1}{2E_{\bf k}},$
$\chi_{13}=\frac{1}{2}\sum_{p}{\rm
Tr}\Big{[}\frac{p_{0}\\!+\\!2Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}}\tau_{1}\frac{p_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf
k}}\tau_{3}\Big{]}=\sum_{p}\frac{2\xi_{\bf
k}\Delta_{0}}{[(p_{0}+2Q_{0})^{2}-E^{2}_{\bf k}](p_{0}^{2}-E^{2}_{\bf k})}=0,$
(73) $\displaystyle\chi_{003}$ $\displaystyle=$ $\displaystyle\sum_{p}{\rm
Tr}\Big{[}\frac{p_{0}\\!+\\!2Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}}\frac{p_{0}\\!+\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+Q_{0})^{2}-E^{2}_{\bf
k}}\frac{p_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf k}}\tau_{3}\Big{]}$
(74) $\displaystyle=$ $\displaystyle\sum_{p}2\xi_{\bf
k}\frac{(p_{0}+2Q_{0})(p_{0}+Q_{0})+\xi_{\bf
k}^{2}+\Delta_{0}^{2}+2p^{2}_{0}+3p_{0}Q_{0})}{[(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}][(p_{0}+Q_{0})^{2}-E^{2}_{\bf k}](p_{0}^{2}-E^{2}_{\bf k})}=0.$
$\displaystyle\chi_{0{\bar{3}}3}$ $\displaystyle=$ $\displaystyle\sum_{p}{\rm
Tr}\Big{\\{}\frac{p_{0}\\!+\\!2Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}}\partial_{\xi_{\bf k}}\Big{[}\frac{p_{0}\\!+\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+Q_{0})^{2}-E^{2}_{\bf
k}}\Big{]}\tau_{3}\frac{p_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf
k}}\tau_{3}\Big{\\}}$ (75) $\displaystyle=$
$\displaystyle\sum_{p}2\Big{\\{}\frac{2\xi_{\bf
k}(p_{0}+Q_{0})}{[(p_{0}\\!+\\!2Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}](p_{0}^{2}\\!-\\!E^{2}_{\bf
k})}\\!+\\!\frac{[(p_{0}\\!+\\!2Q_{0})(p_{0}\\!+\\!Q_{0})\\!+\\!E_{\bf
k}^{2}]p_{0}\\!+\\!(2p_{0}\\!+\\!3Q_{0})(\xi_{\bf
k}^{2}-\Delta_{0}^{2})}{[(p_{0}\\!+\\!2Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}](p_{0}^{2}\\!-\\!E^{2}_{\bf k})}\partial_{\xi_{\bf
k}}\Big{\\}}\Big{[}\frac{1}{(p_{0}\\!+\\!Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}}\Big{]}$ $\displaystyle=$ $\displaystyle 0,~{}~{}~{}~{}$
$\displaystyle\chi_{3{\bar{0}}3}$ $\displaystyle=$ $\displaystyle\sum_{p}{\rm
Tr}\Big{\\{}\frac{p_{0}\\!+\\!2Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}}\tau_{3}\partial_{\xi_{\bf k}}\Big{[}\frac{p_{0}\\!+\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+Q_{0})^{2}-E^{2}_{\bf
k}}\Big{]}\frac{p_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf
k}}\tau_{3}\Big{\\}}$ (76) $\displaystyle=$
$\displaystyle\sum_{p}2\Big{\\{}\frac{2\xi_{\bf
k}(p_{0}+Q_{0})}{[(p_{0}\\!+\\!2Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}](p_{0}^{2}\\!-\\!E^{2}_{\bf
k})}\\!+\\!\frac{(p_{0}\\!+\\!2Q_{0})[(p_{0}\\!+\\!Q_{0})p_{0}\\!+\\!E_{\bf
k}^{2}]\\!+\\!(2p_{0}\\!+\\!Q_{0})(\xi_{\bf
k}^{2}-\Delta_{0}^{2})}{[(p_{0}\\!+\\!2Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}](p_{0}^{2}\\!-\\!E^{2}_{\bf k})}\partial_{\xi_{\bf
k}}\Big{\\}}\Big{[}\frac{1}{(p_{0}\\!+\\!Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}}\Big{]}$ $\displaystyle=$ $\displaystyle 0,~{}~{}~{}~{}$
$\displaystyle\chi_{3{\bar{\bar{3}}}3}$ $\displaystyle=$
$\displaystyle\sum_{p}{\rm Tr}\Big{\\{}\frac{p_{0}\\!+\\!2Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}}\tau_{3}\partial^{2}_{\xi_{\bf
k}}\Big{[}\frac{p_{0}\\!+\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+Q_{0})^{2}-E^{2}_{\bf
k}}\Big{]}\tau_{3}\frac{p_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf
k}}\tau_{3}\Big{\\}}$ (77) $\displaystyle=$
$\displaystyle\sum_{p}2\Big{\\{}\frac{2[(p_{0}\\!+\\!2Q_{0})p_{0}\\!+\\!\xi^{2}_{\bf
k}\\!-\\!\Delta_{0}^{2}]\partial_{\xi_{\bf
k}}}{[(p_{0}\\!+\\!2Q_{0})^{2}\\!-\\!E^{2}_{\bf k}](p_{0}^{2}\\!-\\!E^{2}_{\bf
k})}\\!+\\!\frac{\xi_{\bf
k}[(p_{0}\\!+\\!2Q_{0})(p_{0}\\!+\\!Q_{0})\\!+\\!\xi_{\bf
k}^{2}\\!-\\!3\Delta_{0}^{2}\\!+\\!(2p_{0}\\!+\\!3Q_{0})p_{0}]\partial^{2}_{\xi_{\bf
k}}}{[(p_{0}+2Q_{0})^{2}-E^{2}_{\bf k}](p_{0}^{2}-E^{2}_{\bf
k})}\Big{\\}}\Big{[}\frac{1}{(p_{0}\\!+\\!Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}}\Big{]}$ $\displaystyle=$ $\displaystyle 0,$ $\displaystyle\chi_{001}$
$\displaystyle=$ $\displaystyle\sum_{p}{\rm
Tr}\Big{\\{}\frac{p_{0}\\!+\\!2Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}}\frac{p_{0}\\!+\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+Q_{0})^{2}-E^{2}_{\bf
k}}\frac{p_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf
k}}\tau_{1}\Big{\\}}$ (78) $\displaystyle=$
$\displaystyle\sum_{p}2\Delta_{0}\frac{E^{2}_{\bf
k}\\!+\\!(p_{0}\\!+\\!2Q_{0})p_{0}\\!+\\!(p_{0}\\!+\\!Q_{0})p_{0}\\!+\\!(p_{0}\\!+\\!2Q_{0})(p_{0}\\!+\\!Q_{0})}{[(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}][(p_{0}+Q_{0})^{2}-E^{2}_{\bf k}](p_{0}^{2}-E^{2}_{\bf k})}$
$\displaystyle=$ $\displaystyle\sum_{\bf k}\Big{\\{}\Big{[}\frac{f(E_{\bf
k}-2Q_{0})\delta_{ip_{0}+2Q_{0},E_{\bf
k}}}{[(ip_{0}\\!+\\!Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}][(ip_{0})^{2}\\!-\\!E^{2}_{\bf k}]}\\!+\\!\frac{f(E_{\bf
k}-Q_{0})\delta_{ip_{0}+Q_{0},E_{\bf
k}}}{[(ip_{0}\\!+\\!2Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}][(ip_{0})^{2}\\!-\\!E^{2}_{\bf k}]}\\!+\\!\frac{f(E_{\bf
k})\delta_{ip_{0},E_{\bf k}}}{[(ip_{0}\\!+\\!2Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}][(ip_{0}\\!+\\!Q_{0})^{2}\\!-\\!E^{2}_{\bf k}]}\Big{]}$
$\displaystyle\mbox{}\times\frac{2\Delta_{0}[E^{2}_{\bf
k}\\!+\\!(ip_{0}\\!+\\!2Q_{0})ip_{0}\\!+\\!(ip_{0}\\!+\\!Q_{0})ip_{0}\\!+\\!(ip_{0}\\!+\\!2Q_{0})(ip_{0}\\!+\\!Q_{0})]}{2E_{\bf
k}}+\\{E_{\bf k}\rightarrow-E_{\bf k}\\}\Big{\\}}$ $\displaystyle=$
$\displaystyle 2\Delta_{0}\sum_{\bf k}\Big{\\{}\Big{[}\frac{f(E_{\bf
k}\\!-\\!2Q_{0})(ip_{0}\\!+\\!Q_{0}\\!+\\!ip_{0}\\!+\\!2Q_{0})(ip_{0}\\!+\\!ip_{0}\\!+\\!2Q_{0})\delta_{ip_{0}+2Q_{0},E_{\bf
k}}}{2E_{\bf
k}[(ip_{0}\\!+\\!Q_{0})^{2}\\!-\\!(ip_{0}\\!+\\!2Q_{0})^{2}][(ip_{0})^{2}\\!-\\!(ip_{0}\\!+\\!2Q_{0})^{2}]}\\!+\\!\frac{f(E_{\bf
k}\\!-\\!Q_{0})(ip_{0}\\!+\\!2Q_{0}\\!+\\!ip_{0}\\!+\\!Q_{0})}{2E_{\bf
k}[(ip_{0}\\!+\\!2Q_{0})^{2}\\!-\\!(ip_{0}\\!+\\!Q_{0})^{2}]}$
$\displaystyle\mbox{}\times\frac{(ip_{0}\\!+\\!Q_{0}\\!+\\!ip_{0})\delta_{ip_{0}+Q_{0},E_{\bf
k}}}{[(ip_{0})^{2}\\!-\\!(ip_{0}\\!+\\!Q_{0})^{2}]}\\!+\\!\frac{f(E_{\bf
k})\delta_{ip_{0},E_{\bf
k}}(ip_{0}\\!+\\!2Q_{0}\\!+\\!ip_{0})(ip_{0}\\!+\\!Q_{0}\\!+\\!ip_{0})}{2E_{\bf
k}[(ip_{0}\\!+\\!2Q_{0})^{2}\\!-\\!(ip_{0})^{2}][(ip_{0}\\!+\\!Q_{0})^{2}\\!-\\!(ip_{0})^{2}]}\Big{]}+[E_{\bf
k}\rightarrow-E_{\bf k}]\Big{\\}}$ $\displaystyle=$
$\displaystyle\frac{2\Delta_{0}}{2Q_{0}^{2}}\sum_{\bf
k}\Big{\\{}\Big{[}\frac{f(E_{\bf k}\\!-\\!2Q_{0})}{2E_{\bf
k}}\\!+\\!\frac{f(E_{\bf k})}{2E_{\bf k}}\\!-\\!\frac{2f(E_{\bf
k}\\!-\\!Q_{0})}{2E_{\bf k}}\Big{]}+\\{E_{\bf k}\rightarrow-E_{\bf
k}\\}\Big{\\}}.$ $\displaystyle\chi_{3{\bar{\bar{3}}}1}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\sum_{p}{\rm
Tr}\Big{\\{}\frac{p_{0}\\!+\\!2Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+2Q_{0})^{2}-E^{2}_{\bf
k}}\tau_{3}\partial^{2}_{\xi_{\bf
k}}\Big{[}\frac{p_{0}\\!+\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+Q_{0})^{2}-E^{2}_{\bf
k}}\Big{]}\tau_{3}\frac{p_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf
k}}\tau_{1}\Big{\\}}$ (79) $\displaystyle=$
$\displaystyle\frac{1}{3}\sum_{p}\partial_{\xi_{\bf k^{\prime}}}\Big{\\{}{\rm
Tr}\Big{[}\frac{p_{0}\\!+\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+Q_{0})^{2}-E^{2}_{\bf
k}}\frac{p_{0}\\!+\\!\xi_{\bf
k^{\prime}}\tau_{3}\\!-\\!\Delta_{0}\tau_{1}}{p_{0}^{2}-E^{2}_{\bf
k^{\prime}}}\frac{p_{0}\\!-\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}-Q_{0})^{2}-E^{2}_{\bf
k}}\tau_{1}\Big{]}\Big{\\}}_{\bf k^{\prime}{\rightarrow}k}$ $\displaystyle=$
$\displaystyle\frac{2\Delta_{0}}{3}\sum_{p}\partial_{\xi_{\bf
k^{\prime}}}\Big{\\{}\frac{p_{0}^{2}+Q_{0}^{2}+\xi_{\bf k}^{2}+2\xi_{\bf
k}\xi_{\bf k^{\prime}}-\Delta_{0}^{2}}{[(p_{0}+Q_{0})^{2}-E^{2}_{\bf
k}](p_{0}^{2}-E_{\bf k^{\prime}}^{2})[(p_{0}-Q_{0})^{2}-E^{2}_{\bf
k}]}\Big{\\}}_{\bf k^{\prime}{\rightarrow}k}$ $\displaystyle=$
$\displaystyle\frac{2\Delta_{0}}{3}\sum_{p}\partial_{\xi_{\bf
k^{\prime}}}\Big{\\{}\frac{(\xi_{\bf k}+\xi_{\bf
k^{\prime}})^{2}+Q_{0}^{2}}{[(p_{0}+Q_{0})^{2}-E^{2}_{\bf k}](p_{0}^{2}-E_{\bf
k^{\prime}}^{2})[(p_{0}-Q_{0})^{2}-E^{2}_{\bf
k}]}+\frac{1}{[(p_{0}+Q_{0})^{2}-E^{2}_{\bf k}][(p_{0}-Q_{0})^{2}-E^{2}_{\bf
k}]}\Big{\\}}_{\bf k^{\prime}{\rightarrow}k}$ $\displaystyle=$
$\displaystyle\frac{2\Delta_{0}}{3(iQ_{0})^{2}}\\!\sum_{\bf
k}\big{\\{}\big{[}2\\!+\\!4\xi_{\bf k}\partial_{\xi_{\bf
k^{\prime}}}\\!+\\!(4\xi_{\bf k}^{2}\\!+\\!Q_{0}^{2})\partial_{\xi_{\bf
k^{\prime}}}^{2}\big{]}I_{\bf k}(\xi_{\bf k^{\prime}})\big{\\}}_{\bf
k^{\prime}{\rightarrow}k}\\!\approx\\!\frac{2\Delta_{0}}{(iQ_{0})^{2}}\\!\sum_{\bf
k}\frac{4\xi^{2}_{\bf k}}{3}\Big{\\{}\frac{3\xi_{\bf k}^{2}[1\\!-\\!2f(E_{\bf
k})]}{8E^{7}_{\bf k}}\\!+\\!\frac{3\xi_{\bf k}^{2}}{4E_{\bf
k}^{6}}\partial_{E_{\bf k}}f(E_{\bf k})\Big{\\}}$ $\displaystyle=$
$\displaystyle\frac{2\Delta_{0}}{(iQ_{0})^{2}}\\!\sum_{\bf k}\frac{\xi_{\bf
k}^{4}}{E_{\bf k}^{5}}\partial_{E_{\bf k}}\Big{[}\frac{2f(E_{\bf
k})-1}{2E_{\bf k}}\Big{]},$
$\displaystyle\chi_{3{\bar{0}}1}\\!-\\!\chi_{0{\bar{3}}1}\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\sum_{p}\partial_{\xi_{\bf
k^{\prime}}}\Big{\\{}{\rm Tr}\Big{[}\frac{p_{0}\\!+\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}+Q_{0})^{2}-E^{2}_{\bf
k}}\frac{[\tau_{3},p_{0}\\!+\\!\xi_{\bf
k^{\prime}}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}]}{p_{0}^{2}-E^{2}_{\bf
k^{\prime}}}\frac{p_{0}\\!-\\!Q_{0}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(p_{0}-Q_{0})^{2}-E^{2}_{\bf
k}}\tau_{1}\Big{]}\Big{\\}}_{\bf k^{\prime}{\rightarrow}k}$ (80)
$\displaystyle=$ $\displaystyle\\!\\!\sum_{p}\partial_{\xi_{\bf
k^{\prime}}}\Big{\\{}\frac{-8Q_{0}\Delta_{0}\xi_{\bf
k}}{[(p_{0}\\!+\\!Q_{0})^{2}\\!-\\!E^{2}_{\bf k}](p_{0}^{2}\\!-\\!E_{\bf
k^{\prime}}^{2})[(p_{0}\\!-\\!Q_{0})^{2}\\!-\\!E^{2}_{\bf k}]}\Big{\\}}_{\bf
k^{\prime}{\rightarrow}k}\\!\\!=\\!\frac{8i\Delta_{0}}{iQ_{0}}\sum_{\bf
k}\xi_{\bf k}\big{[}\partial_{\xi_{\bf k^{\prime}}}I_{\bf k}(\xi_{\bf
k^{\prime}})\big{]}_{\bf
k^{\prime}{\rightarrow}k}\\!\\!\approx\\!0,~{}~{}~{}~{}~{}$

with

$\displaystyle I_{\bf k}(\xi_{\bf k^{\prime}})$ $\displaystyle=$
$\displaystyle\int\frac{dp_{0}}{2\pi}\frac{(iQ_{0})^{2}}{[(p_{0}\\!+\\!Q_{0})^{2}\\!-\\!E^{2}_{\bf
k}](p_{0}^{2}\\!-\\!E_{\bf
k^{\prime}}^{2})[(p_{0}\\!-\\!Q_{0})^{2}\\!-\\!E^{2}_{\bf k}]}$ (81)
$\displaystyle=$ $\displaystyle\frac{f(E_{\bf k^{\prime}})}{2E_{\bf
k^{\prime}}}\frac{1}{(4E_{\bf
k^{\prime}}^{2}-Q^{2}_{0})}\Big{[}1\\!+\\!\frac{\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}}}{(2E_{\bf
k^{\prime}}\\!+\\!Q_{0})Q_{0}}-\frac{\xi_{\bf k}^{2}\\!-\\!\xi^{2}_{\bf
k^{\prime}}}{(2E_{\bf k^{\prime}}\\!-\\!Q_{0})Q_{0}}\\!+\\!\frac{(\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}})^{2}}{(2E_{\bf
k^{\prime}}\\!+\\!Q_{0})^{2}Q^{2}_{0}}\\!+\\!\frac{(\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}})^{2}}{(2E_{\bf
k^{\prime}}\\!-\\!Q_{0})^{2}Q^{2}_{0}}$
$\displaystyle\mbox{}\\!-\\!\frac{(\xi_{\bf k}^{2}\\!-\\!\xi^{2}_{\bf
k^{\prime}})^{2}}{(4E^{2}_{\bf
k^{\prime}}\\!-\\!Q^{2}_{0})Q^{2}_{0}}\Big{]}\\!-\\!\frac{f(E_{\bf
k}\\!-\\!Q_{0})}{4E_{\bf k}}\frac{1}{(2E_{\bf k}\\!-\\!2Q_{0})(2E_{\bf
k}\\!-\\!Q_{0})}\Big{[}1\\!+\\!\frac{\xi_{\bf k}^{2}\\!-\\!\xi^{2}_{\bf
k^{\prime}}}{(2E_{\bf k}\\!-\\!Q_{0})Q_{0}}\\!+\\!\frac{(\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}})^{2}}{(2E_{\bf
k}\\!-\\!Q_{0})^{2}Q^{2}_{0}}\Big{]}$ $\displaystyle\mbox{}-\frac{f(E_{\bf
k}\\!+\\!Q_{0})}{4E_{\bf k}}\frac{1}{(2E_{\bf k}\\!+\\!2Q_{0})(2E_{\bf
k}\\!+\\!Q_{0})}\Big{[}1\\!-\\!\frac{\xi_{\bf k}^{2}\\!-\\!\xi^{2}_{\bf
k^{\prime}}}{(2E_{\bf k}\\!-\\!Q_{0})Q_{0}}\\!+\\!\frac{(\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}})^{2}}{(2E_{\bf
k}\\!+\\!Q_{0})^{2}Q^{2}_{0}}\Big{]}+\\{E_{\bf k}\rightarrow-E_{\bf k},E_{\bf
k^{\prime}}\rightarrow-E_{\bf k^{\prime}}\\}$ $\displaystyle\approx$
$\displaystyle\frac{f(E_{\bf k^{\prime}})}{2E_{\bf k^{\prime}}(4E_{\bf
k^{\prime}}^{2}-Q^{2}_{0})}\Big{[}1\\!+\\!\frac{\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}}}{(2E_{\bf
k^{\prime}}\\!+\\!Q_{0})Q_{0}}-\frac{\xi_{\bf k}^{2}\\!-\\!\xi^{2}_{\bf
k^{\prime}}}{(2E_{\bf k^{\prime}}\\!-\\!Q_{0})Q_{0}}\\!+\\!\frac{(\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}})^{2}}{(2E_{\bf
k^{\prime}}\\!+\\!Q_{0})^{2}Q^{2}_{0}}\\!+\\!\frac{(\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}})^{2}}{(2E_{\bf
k^{\prime}}\\!-\\!Q_{0})^{2}Q^{2}_{0}}$
$\displaystyle\mbox{}\\!-\\!\frac{(\xi_{\bf k}^{2}\\!-\\!\xi^{2}_{\bf
k^{\prime}})^{2}}{(4E^{2}_{\bf
k^{\prime}}\\!-\\!Q^{2}_{0})Q^{2}_{0}}\Big{]}\\!-\\!\frac{[1\\!-\\!Q_{0}\partial_{E_{\bf
k}}+\frac{Q_{0}^{2}}{2}\partial_{E_{\bf k}}^{2}]f(E_{\bf k})}{4E_{\bf
k}(2E_{\bf k}\\!-\\!2Q_{0})(2E_{\bf
k}\\!-\\!Q_{0})}\Big{[}1\\!+\\!\frac{\xi_{\bf k}^{2}\\!-\\!\xi^{2}_{\bf
k^{\prime}}}{(2E_{\bf k}\\!-\\!Q_{0})Q_{0}}\\!+\\!\frac{(\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}})^{2}}{(2E_{\bf
k}\\!-\\!Q_{0})^{2}Q^{2}_{0}}\Big{]}$
$\displaystyle\mbox{}-\frac{[1\\!+\\!Q_{0}\partial_{E_{\bf
k}}+\frac{Q_{0}^{2}}{2}\partial_{E_{\bf k}}^{2}]f(E_{\bf k})}{4E_{\bf
k}(2E_{\bf k}\\!+\\!2Q_{0})(2E_{\bf
k}\\!+\\!Q_{0})}\Big{[}1\\!-\\!\frac{\xi_{\bf k}^{2}\\!-\\!\xi^{2}_{\bf
k^{\prime}}}{(2E_{\bf k}\\!-\\!Q_{0})Q_{0}}\\!+\\!\frac{(\xi_{\bf
k}^{2}\\!-\\!\xi^{2}_{\bf k^{\prime}})^{2}}{(2E_{\bf
k}\\!+\\!Q_{0})^{2}Q^{2}_{0}}\Big{]}+\\{E_{\bf k}\rightarrow-E_{\bf k},E_{\bf
k^{\prime}}\rightarrow-E_{\bf k^{\prime}}\\},$

and hence

$[I_{\bf k}(\xi_{\bf k^{\prime}})]_{{\bf k}\rightarrow{\bf
k^{\prime}},\Omega\ll{E_{\bf k}}}=\Big{[}\frac{2f(E_{\bf k})-1}{8E^{3}_{\bf
k}}-\frac{2f(E_{\bf k})-1}{16E^{3}_{\bf k}}-\frac{2f(E_{\bf
k})-1}{16E^{3}_{\bf k}}\Big{]}=0,$ (82) $\displaystyle[\partial_{\xi_{\bf
k^{\prime}}}I_{\bf k}(\xi_{\bf k^{\prime}})]_{{\bf k}\rightarrow{\bf
k^{\prime}},\Omega\ll{E_{\bf k}}}$ $\displaystyle=$
$\displaystyle\frac{\xi_{\bf k}}{Q_{0}(2E_{\bf
k})^{4}}\Big{[}\frac{2Q_{0}f(E_{\bf k})}{E_{\bf k}}\\!+\\!f(E_{\bf
k})\Big{(}1\\!+\\!\frac{2Q_{0}}{E_{\bf k}}\Big{)}\\!-\\!Q_{0}\partial_{E_{\bf
k}}f(E_{\bf k})-f(E_{\bf k})\Big{(}1\\!-\\!\frac{2Q_{0}}{E_{\bf
k}}\Big{)}\\!-\\!Q_{0}\partial_{E_{\bf k}}f(E_{\bf k})\Big{]}$ (83)
$\displaystyle+\frac{\xi_{\bf k}}{E_{\bf k}}\partial_{E_{\bf
k}}\Big{[}\frac{f(E_{\bf k})}{8E_{\bf k}^{3}}\Big{]}+\\{E_{\bf k}\rightarrow-
E_{\bf k}\\}=0,$ $\displaystyle[\partial_{\xi_{\bf k^{\prime}}}^{2}I_{\bf
k}(\xi_{\bf k^{\prime}})]_{{\bf k}\rightarrow{\bf k^{\prime}},\Omega\ll{E_{\bf
k}}}$ $\displaystyle=$ $\displaystyle\frac{1}{\xi_{\bf k}}\partial_{\xi_{\bf
k}}\Big{[}\frac{f(E_{\bf k})}{8E_{\bf k}^{3}}\Big{]}\\!-\\!\partial_{\xi_{\bf
k}}^{2}\Big{[}\frac{f(E_{\bf k})}{8E_{\bf k}^{3}}\Big{]}\\!-\\!\xi_{\bf
k}\partial_{\xi_{\bf k}}\Big{[}\frac{f(E_{\bf k})}{4E_{\bf
k}^{5}}\Big{]}\\!+\\!\frac{\xi_{\bf k}^{2}f(E_{\bf k})}{4E_{\bf
k}^{7}}\\!-\\!\frac{\xi_{\bf k}^{2}}{4E_{\bf k}^{6}}\partial_{E_{\bf
k}}f(E_{\bf k})\\!+\\!\frac{\xi_{\bf k}^{2}}{8E_{\bf
k}^{5}}\partial^{2}_{E_{\bf k}}f(E_{\bf k})$ (84)
$\displaystyle\mbox{}+\\{E_{\bf k}\rightarrow-E_{\bf k}\\}=\frac{3\xi_{\bf
k}^{2}[1-2f(E_{\bf k})]}{8E^{7}_{\bf k}}+\frac{3\xi_{\bf k}^{2}}{4E_{\bf
k}^{6}}\partial_{E_{\bf k}}f(E_{\bf k}).$

Here, we have applied the Wick rotation [i.e.,
$\int\frac{dp_{0}}{2\pi}F(p_{0})\rightarrow{T}\sum_{\omega_{n}}F(i\omega_{n})$]
to map the frequency integral into the Matsubara frequency summationG1 . It is
noted that $\chi_{13}$ [Eq. (73)], $\chi_{003}$ [Eq. (74)],
$\chi_{0{\bar{3}}3}$ [Eq. (75)], $\chi_{3{\bar{0}}3}$ [Eq. (76)] and
$\chi_{3{\bar{\bar{3}}}3}$ [Eq. (77)] vanish 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}$. The second-order correlation
coefficients $\chi_{11}$ [Eq. (71)] and $\chi_{33}$ [Eq. (72)] and $\chi_{13}$
[Eq. (73)] here are exactly same as the ones obtained in the previous works by
Cea et al.Cea1 ; Cea2 ; Cea3 . Moreover, one also finds that the phase-related
coefficient $\chi_{33}$ [Eq. (72)] from the path-integral approach is exactly
same as the one $u_{g}$ [Eq. (28)] from GIKE.

##  Appendix B Derivation of Ginzburg-Landau equation from Eilenberger
equation

In this part, we present the derivation of Ginzburg-Landau equation from
Eilenberger equation at clean and stationary caseBa20 . In this circumstance,
the Eilenberger equation in Matsubara formalism reads

$[i\omega_{n}\tau_{3}-{\hat{\Delta}}({\bf R})\tau_{3},g]\\!+\\!i{\bf
v}_{F}\cdot{\bf\nabla}_{\bf R}g\\!+\\!{e{\bf A}\cdot{\bf
v}_{F}}[\tau_{3},g]=0.$ (85)

At the weak field, the quasiclassical $\tau_{3}$-Green function can be
expanded as $g=g^{(0)}+\sum_{n=1}\delta{g}^{(n)}$, with $\delta{g}^{(n)}$
being the $n$-th order response.

Considering the anomalous Green function (i.e., off-diagonal part
$\delta{g}^{(n)}_{12}$) of $\delta{g}^{(n)}$, from Eq. (85), one has

$2i\omega_{n}\delta{g}^{(n)}_{12}=-i{\bf v}_{F}\cdot(\partial_{\bf R}-2ie{\bf
A})\delta{g}^{(n-1)}_{12}.$ (86)

Then, keeping the expansions up to the second-order response, one finds the
solution:

$g_{12}=g^{(0)}_{12}-\frac{{\bf v}_{F}\cdot(\partial_{\bf R}-2ie{\bf
A})}{2\omega_{n}}{g}^{(0)}_{12}+\frac{[{\bf v}_{F}\cdot(\partial_{\bf
R}-2ie{\bf A})]^{2}}{4(\omega_{n})^{2}}{g}^{(0)}_{12}.$ (87)

The equilibrium quasiclassical $\tau_{3}$-Green function can be derived by
Gorkov equation and its anomalous Green function is written asSilaev

$g^{(0)}_{12}=\frac{i\Delta}{\sqrt{(\omega_{n})^{2}+|\Delta|^{2}}}.$ (88)

Then, with the solved $g^{(0)}_{12}$ and hence $g_{12}$ in Eq. (87), from the
corresponding gap equation $\Delta=-iUN(0)\langle{g_{12}}\rangle_{F}$, near
$T_{c}$, one obtains

$\sum_{\omega_{n}>0}\Big{[}\frac{v_{F}^{2}(\partial_{\bf R}-2ie{\bf
A})^{2}}{6(\omega_{n})^{3}}-\frac{|\Delta|^{2}}{(\omega_{n})^{3}}\Big{]}\Delta+\Big{[}\sum_{\omega_{n}>0}\frac{2}{\omega_{n}}-\frac{1}{UN(0)}\Big{]}\Delta=0.$
(89)

Consequently, through the mathematical calculation, the above equation becomes

$\Big{[}\frac{7R(3)k_{F}^{2}}{24(\pi{T})^{2}}\frac{(\partial_{\bf R}-2ie{\bf
A})^{2}}{2m}-\frac{7R(3)}{8({\pi}T)^{2}}{|\Delta|^{2}}+\ln\Big{(}\frac{T_{c}}{T}\Big{)}\Big{]}\Delta=0,$
(90)

which exactly recovers the Ginzburg-Landau equation.

##  Appendix C Derivation of Ginzburg-Landau equation within the path-integral
approach

In this part, we present the derivation of Ginzburg-Landau equation within the
path-integral approach at clean and stationary case. Specifically, with vector
potential alone, after the integration over Fermi field within the path-
integral approach, the gap-variation related part from Eq. (52) is written as

$S_{H}=-\Big{[}(\chi_{11}-U^{-1})\delta|\Delta|^{2}+\Big{(}\chi_{001}\delta|\Delta|\frac{e^{2}A^{2}_{0}v_{F}^{2}}{3}+\chi_{13}\delta|\Delta|\frac{e^{2}A^{2}}{2m}+h.c.\Big{)}\Big{]}.$
(91)

with

$\displaystyle\chi_{ij}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\frac{1}{2}\sum_{p}{\rm
Tr}[G_{0}(p)\tau_{i}G_{0}(p)\tau_{j}],$ (92) $\displaystyle\chi_{ijk}\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!\sum_{p}{\rm
Tr}[G_{0}(p)\tau_{i}G_{0}(p)\tau_{j}G_{0}(p)\tau_{k}],$ (93)

In Eq. (91), the second and third terms denote the couplings of Higgs mode to
second order of light-matter interaction, i.e., second order of ${\bf
k}\cdot{e{\bf A}}/m$ and linear one of $\frac{e^{2}A^{2}}{2m}\tau_{3}$,
respectively.

Through the mathematical calculation, the coupling coefficients read

$\displaystyle\chi_{11}\\!-\\!U^{-1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{p}{\rm Tr}\Big{[}\frac{ip_{n}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(ip_{n})^{2}-E^{2}_{\bf
k}}\tau_{1}\frac{ip_{n}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(ip_{n})^{2}\\!-\\!E^{2}_{\bf
k}}\tau_{1}\Big{]}\\!-\\!U^{-1}\\!\\!=\\!\\!\sum_{p}\frac{(ip_{n})^{2}\\!+\\!\Delta_{0}^{2}\\!-\\!\xi^{2}_{\bf
k}}{[(ip_{n})^{2}\\!-\\!E^{2}_{\bf k}]^{2}}\\!-\\!U^{-1}$ (94)
$\displaystyle=$
$\displaystyle\sum_{p}\Big{\\{}\frac{2\Delta_{0}^{2}}{[(ip_{n})^{2}-E^{2}_{\bf
k}]^{2}}+\frac{1}{(ip_{n})^{2}-E^{2}_{\bf
k}}\Big{\\}}-\sum_{p}\frac{1}{(ip_{n})^{2}-E^{2}_{\bf k}}$ $\displaystyle=$
$\displaystyle 2{\Delta_{0}^{2}}\sum_{\bf k}\frac{1}{2E_{\bf
k}^{2}}\Big{[}\frac{1\\!-\\!2f(E_{\bf k})}{2E_{\bf k}}\\!+\\!\partial_{E_{\bf
k}}f(E_{\bf k})\Big{]}=2{\Delta_{0}^{2}}\sum_{\bf k}\frac{\partial_{E_{\bf
k}}}{E_{\bf k}}\Big{[}\frac{2f(E_{\bf k})\\!-\\!1}{2E_{\bf k}}\Big{]},$
$\displaystyle\chi_{13}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{p}{\rm Tr}\Big{[}\frac{ip_{n}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(ip_{n})^{2}-E^{2}_{\bf
k}}\tau_{1}\frac{ip_{n}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(ip_{n})^{2}-E^{2}_{\bf
k}}\tau_{3}\Big{]}=\sum_{p}\frac{2\xi_{\bf
k}\Delta_{0}}{[(ip_{n})^{2}-E^{2}_{\bf k}]^{2}}=0,$ $\displaystyle\chi_{001}$
$\displaystyle=$ $\displaystyle\sum_{p}{\rm
Tr}\Big{[}\frac{ip_{n}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(ip_{n})^{2}-E^{2}_{\bf
k}}\frac{ip_{n}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(ip_{n})^{2}-E^{2}_{\bf
k}}\frac{ip_{n}\\!+\\!\xi_{\bf
k}\tau_{3}\\!+\\!\Delta_{0}\tau_{1}}{(ip_{n})^{2}-E^{2}_{\bf
k}}\tau_{1}\Big{]}=\sum_{p}2\Delta_{0}\frac{E^{2}_{\bf
k}\\!+\\!3(ip_{n})^{2}}{[(ip_{n})^{2}-E^{2}_{\bf k}]^{3}}$ (95)
$\displaystyle=$ $\displaystyle\Delta_{0}\sum_{\bf
k}\Big{\\{}\Big{[}\Big{(}\frac{6}{8E_{\bf k}^{3}}\\!-\\!\frac{36E_{\bf
k}}{16E_{\bf k}^{4}}\\!+\\!\frac{48E_{\bf k}^{2}}{32E_{\bf
k}^{5}}\Big{)}f(E_{\bf k})\\!+\\!2\Big{(}\frac{6E_{\bf k}}{8E_{\bf
k}^{3}}\\!-\\!\frac{12E_{\bf k}^{2}}{16E_{\bf k}^{4}}\Big{)}\partial_{E_{\bf
k}}f(E_{\bf k})\\!+\\!\frac{4E_{\bf k}^{2}}{8E_{\bf k}^{3}}\partial_{E_{\bf
k}}^{2}f(E_{\bf k})\Big{]}$ $\displaystyle\mbox{}\\!+\\![E_{\bf k}\rightarrow-
E_{\bf k}]\Big{\\}}=2\Delta_{0}\sum_{\bf k}\frac{\partial_{E_{\bf
k}}^{2}f(E_{\bf k})}{2E_{\bf k}}.$

Then, further following the derivation of Eqs. (32) and (33) and notation, one
has
$\chi_{001}=\Delta_{0}\frac{7DR(3)}{2(\pi{T})^{2}}=2D\Delta_{0}\lambda_{l}$
and
$\chi_{11}-U^{-1}=4\Delta_{0}^{2}\frac{7DR(3)}{8(\pi{T})^{2}}=\beta_{L}4\Delta_{0}^{2}D$
near $T_{c}$. In this circumstance, using
$\Delta_{0}^{2}=-\alpha_{L}/\beta_{L}$, Eq. (91) becomes

$S_{H}=-2D\Big{(}\alpha_{L}\delta|\Delta|^{2}+\frac{\beta_{L}6\Delta_{0}^{2}\delta|\Delta|^{2}}{2}+\lambda_{L}\frac{|2ieA_{0}|^{2}2\Delta_{0}\delta|\Delta|}{4m}\Big{)},$
(96)

which exactly recovers the nonequilibrium variation of Ginzburg-Landau
Lagrangian at stationary case. In fact, within the path-integral approach, the
derivation of this Lagrangian at stationary case is exactly same as the one of
effective action for Higgs mode near $T_{c}$ when $\Omega\rightarrow 0$,
guaranteeing the physical continuity between $\Omega\rightarrow 0$ and
$\Omega=0$.

Moreover, it is noted that mathematically, the finite response coefficient
$\chi_{001}$ in Eq. (95), which is derived at stationary situation, arises
from a third-order residue, and both response coefficients $\lambda_{E}$ [Eq.
(44)] derived from Eilenberger equation and $\lambda_{p}$ [Eq. (59)] derived
from path-integral approach at finite optical frequency $\Omega$, for a
continuous $\Omega\rightarrow 0$, can exactly recover this result with
$\chi_{001}=2\Delta_{0}\lambda_{E/p}$. But if the optical frequency $i\Omega$
in Eqs. (44) and (59) is taken as bosonic Matsubara frequencies, one only
encounters a first-order residue, and hence, finds vanishing response
coefficients at all $\Omega\neq 0$, leading to an unphysical abrupt change
between results at $\Omega=0$ and $\Omega\rightarrow 0$.

## References

  * (1) P. B. Littlewood and C. M. Varma, Phys. Rev. Lett. 47, 811 (1981); Phys. Rev. B 26, 4883 (1982).
  * (2) Y. Nambu, Phys. Rev. 117, 648 (1960).
  * (3) V. Ambegaokar and L. P. Kadanoff, Nuovo Cimento 22, 914 (1961).
  * (4) J. Goldstone, Nuovo Cimento 19, 154 (1961).
  * (5) J. Goldstone, A. Salam, and S. Weinberg, Phys. Rev. 127, 965 (1962).
  * (6) J. R. Schrieffer, Theory of Superconductivity (W. A. Benjamin, New York, 1964).
  * (7) H. A. Fertig and S. D. Sarma, Phys. Rev. Lett. 65, 1482 (1990).
  * (8) I. J. R. Aitchison, P. Ao, D. J. Thouless, and X. M. Zhu, Phys. Rev. B 51, 6531 (1995).
  * (9) K. Kadowaki, I. Kakeya, M. B. Gaifullin, T. Mochiku, S. Takahashi, T. Koyama, and M. Tachiki, Phys. Rev. B 56, 5617 (1997).
  * (10) K. Kadowaki, I. Kakeya, and K. Kindo, Europhys. Lett. 42, 203 (1998).
  * (11) I. J. R. Aitchison, G. Metikas, and D. J. Lee, Phys. Rev. B 62, 6638 (2000).
  * (12) Y. Nambu, Rev. Mod. Phys. 81, 1015 (2009).
  * (13) C. Timm, Theory of Superconductivity (Institute of theoretical Physics Dresden, 2012).
  * (14) B. V. Svistunov, E. S. Babaev, and N. Prokof’ev, Superfluid States of Matter (CRC Press, Boca Raton, 2015).
  * (15) T. Yanagisawa, Commun. Comput. Phys. 23, 459 (2017).
  * (16) Z. Y. Sun, M. M. Fogler, D. N. Basov, and Andrew J. Millis, Phys. Rev. Research 2, 023413 (2020).
  * (17) A. F. Volkov and S. M. Kogan, Zh. Eksp. Teor. Fiz 65, 2038 (1974) [Sov. Phys. JETP 38, 1018 (1974)].
  * (18) E. A. Yuzbashyan and M. Dzero, Phys. Rev. Lett 96, 230404 (2006).
  * (19) V. Gurarie, Phys. Rev. Lett. 103, 075301 (2009).
  * (20) A. Moor, P. A. Volkov, A. F. Volkov, and K. B. Efetov, Phys. Rev. B 90, 024511 (2014).
  * (21) D. Pekker and C. Varma, Annu. Rev. Condens. Matter Phys. 6, 269 (2015).
  * (22) N. Tsuji, Y. Murakami, and H. Aoki, Phys. Rev. B 94, 224519 (2016).
  * (23) F. Englert and R. Brout, Phys. Rev. Lett. 13, 321 (1964).
  * (24) P. W. Higgs, Phys. Lett. 12, 132 (1964); Phys. Rev. Lett. 13, 508 (1964).
  * (25) G. S. Guralnik, C. R. Hagen, and T. W. B. Kibble, Phys. Rev. Lett. 13, 585 (1964).
  * (26) R. Matsunaga and R. Shimano, Phys. Rev. Lett. 109, 187002 (2012).
  * (27) R. Matsunaga, Y. I. Hamada, K. Makise, Y. Uzawa, H. Terai, Z. Wang, and R. Shimano, Phys. Rev. Lett. 111, 057002 (2013).
  * (28) 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).
  * (29) R. Matsunaga, N. Tsuji, K. Makise, H. Terai, H. Aoki, and R. Shimano, Phys. Rev. B 96, 020505 (2017).
  * (30) 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).
  * (31) 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).
  * (32) K. Katsumi, Z. Z. Li, H. Raffy, Y. Gallais, R. Shimano, Phys. Rev. B 102, 054510 (2020).
  * (33) K. Isoyama, N. Yoshikawa, K. Katsumi, J. Wong, N. Shikama, Y. Sakishita, F. Nabeshima, A. Maeda, and R. Shimano, Commun. Phys. 4, 160 (2021).
  * (34) R. A. Barankov, L. S. Levitov, and B. Z. Spivak, Phys. Rev. Lett. 93, 160401 (2004).
  * (35) R. A. Barankov and L. S. Levitov, Phys. Rev. Lett. 96, 230403 (2006).
  * (36) N. Tsuji and H. Aoki, Phys. Rev. B 92, 064508 (2015).
  * (37) M. Dzero, M. Khodas, and A. Levchenko, Phys. Rev. B 91, 214505 (2015).
  * (38) M. Lu, H. W. Liu, P. Wang, and X. C. Xie, Phys. Rev. B 93, 064516 (2016).
  * (39) Y. Murotani, N. Tsuji, and H. Aoki, Phys. Rev. B 95, 104503 (2017).
  * (40) L. Schwarz and D. Manske, Phys. Rev. B 101, 184519 (2020).
  * (41) T. Papenkort, V. M. Axt, and T. Kuhn, Phys. Rev. B 76, 224522 (2007).
  * (42) T. Papenkort, T. Kuhn, and V. M. Axt, Phys. Rev. B 78, 132505 (2008).
  * (43) A. F. Kemper, M. A. Sentef, B. Moritz, J. K. Freericks, and T. P. Devereaux, Phys. Rev. B 92, 224517 (2015).
  * (44) H. Krull, N. Bittner, G. S. Uhrig, D. Manske, and A. P. Schnyder, Nat. Commun. 7, 11921 (2016).
  * (45) G. Seibold, M. Udina, C. Castellani, and L. Benfatto, Phys. Rev. B 103, 014512 (2017).
  * (46) P. W. Anderson, Phys. Rev. 112, 1900 (1958).
  * (47) S. Tsuchiya, D. Yamamoto, R. Yoshii, and M. Nitta, Phys. Rev. B 98, 094503 (2018).
  * (48) F. Yang and M. W. Wu, Phys. Rev. B 100, 104513 (2019).
  * (49) T. Cea, C. Castellani, and L. Benfatto Phys. Rev. B 93, 180507(R) (2016).
  * (50) T. Cea and L. Benfatto, Phys. Rev. B 94, 064512 (2016).
  * (51) T. Cea, P. Barone, C. Castellani, and L. Benfatto, Phys. Rev. B 97, 094516 (2018).
  * (52) A. A. Abrikosov, L. P. Gor’kov, and I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics (Prentice Hall, Englewood Cliffs, 1963).
  * (53) F. Yang and M. W. Wu, Phys. Rev. B 98, 094507 (2018).
  * (54) Y. Murotani and R. Shimano, Phys. Rev. B 99, 224510 (2019).
  * (55) M. Silaev, Phys. Rev. B 99, 224511 (2019).
  * (56) It is noted that Ref. Silaev,  also used the diagrammatic formalism in order to justify/recover the result from the Eilenberger equation. Nevertheless, a special procedure, 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 (refer to Sec. IV A in Ref. Silaev, ). This procedure is not the standard diagram formulationG1 , but follows the exact step to derive the Eilenberger equationEilen1 . Consequently, the applied diagrammatic formalism approach with this procedure is exactly same as the Eilenberger equation rather than a separate nontrivial approach.
  * (57) N. Tsuji and Y. Nomura, Phys. Rev. Research 2, 043029 (2020).
  * (58) R. Haenel, P. Froese, D. Manske, and L. Schwarz, Phys. Rev. B 104, 134504 (2021).
  * (59) G. Eilenberger, Z. Phys. 214, 195 (1968).
  * (60) T. Kita, Statistical Mechanics of Superconductivity (Springer, Berlin, 2015).
  * (61) H. G. Hugdal, J. Linder, and S. H. Jacobsen, Phys. Rev. B 95, 235403 (2017).
  * (62) F. Yang and M. W. Wu, Phys. Rev. B 102, 144508 (2020).
  * (63) F. Yang and M. W. Wu, Phys. Rev. B 102, 014511 (2020).
  * (64) A. Tagliacozzo and F. Ventriglia, Il Nuovo Cimento D 11, 141 (1989).
  * (65) It is noted that the Lagrangian in Eq. (1) is in fact a basic model (complex scalar field $\Psi$ coupled to an electromagnetic potential $A_{\mu}$) in the standard-model theory [L. Hoddeson, L. Brown, M. Riordan, and M. Dresden, The Rise of the Standard Model (Cambridge, UK, 1997)], but can be applied in superconductors due to the systematical similarityAm0 . While this Lagrangian in superconductors was revealed in Ref. Am6,  through symmetry analysis. In principle, with the well-established Ginzburg-Landau free energy, one can also use the Lorentz invariance of the local Abelian $U(1)$ model to reach the Lagrangian $L=D_{\mu}\Psi^{*}D^{\mu}\Psi-\mu|\Psi|^{2}-\lambda|\Psi|^{4}$, with $D_{\mu}=\partial_{\mu}-ieA_{\mu}$.
  * (66) P. W. Anderson, Phys. Rev. 130, 439 (1963).
  * (67) H. Haug and A. P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors (Springer, Berlin, 1996).
  * (68) M. W. Wu, J. H. Jiang, and M. Q. Weng, Phys. Rep. 493, 61 (2010).
  * (69) M. E. Peskin and D. V. Schroeder, An Introduction to Quantum Field Theory (Addison-Wesley, New York, 1995).
  * (70) I. Khavkine, H. Y. Kee, and K. Maki, Phys. Rev. B 70, 184521 (2004).
  * (71) G. Tkachov and V. I. Fal’ko, Phys. Rev. B 69, 092503 (2004).
  * (72) F. Rohlfing, G. Tkachov, F. Otto, K. Richter, D. Weiss, G. Borghs, and C. Strunk, Phys. Rev. B 80, 220507(R) (2009).
  * (73) F. Yang and M. W. Wu, Phys. Rev. B 95, 075304 (2017).
  * (74) F. Yang and M. W. Wu, J. Low Temp. Phys. 192, 241 (2018).
  * (75) The vanishing influence of the vector potential on the superconducting gap at zero temperature is expected. Specifically, with the vector potential alone at low frequency, from the Hamiltonian in Eq. (3), the vector potential drives the Doppler shift to influence the gap equation $1=g\sum\frac{f(E_{\bf k}^{-})-f(E_{\bf k}^{+})}{2E_{\bf k}}$. Whereas this influence, following the idea of the Fulde-Ferrell-Larkin-Ovchinnikov stateFF1 ; FF2 , vanishes at zero temperature for weak field, since $f(E_{\bf k}^{-})-f(E_{\bf k}^{+})\equiv 1$. Consequently, one finds the vanishing influence of the vector potential on the superconducting gap (i.e., vanishing contribution of the vector potential to the Higgs-mode generation) at zero temperature.
  * (76) J. Rammer and H. Smith, Rev. Mod. Phys. 58, 323 (1986).
  * (77) P. Fulde and R. A. Ferrell, Phys. Rev. 135, 550 (1964).
  * (78) A. I. Larkin and Y. N. Ovchinnikov, Zh. Eksp. Teor. Fiz. 47, 1136 (1964) [Sov. Phys. JETP 20, 762 (1965)].