# Nonreciprocal optical response in parity-breaking superconductors

Hikaru Watanabe  RIKEN Center for Emergent Matter Science (CEMS), Wako
351-0198, Japan  Department of Physics, Graduate School of Science, Kyoto
University, Kyoto 606-8502, Japan  Akito Daido  Department of Physics,
Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan  Youichi
Yanase  Department of Physics, Graduate School of Science, Kyoto University,
Kyoto 606-8502, Japan  Institute for Molecular Science, Okazaki 444-8585,
Japan

###### Abstract

Superconductivity, an emergence of a macroscopic quantum state, gives rise to
unique electromagnetic responses leading to perfect shielding of magnetic
field and zero electrical resistance. In this paper, we propose that
superconductors with the space-inversion symmetry breaking host giant
nonreciprocal optical phenomena, such as photocurrent generation and second
harmonic generation. The nonreciprocal optical responses show divergent
behaviors in the low-frequency regime and originate from two-fold indicators
unique to parity-breaking superconductors, namely, the nonreciprocal
superfluid density and the Berry curvature derivative. Furthermore, the
mechanism and frequency dependence are closely tied to the preserved temporal
symmetry in the superconductor. The relation is useful for probing the space-
time structure of the superconducting symmetry. The indicators characterize
the low-frequency property of nonreciprocal optical responses and hence
quantify the performance of superconductors in nonreciprocal optics.

##  I Introduction

Nonlinear optics is one of the central fields in physics and has been offering
extensive interests from fundamental science to engineering. When the system
has no space-inversion ($\mathcal{P}$) symmetry, the nonlinear optical
response acquires the nonreciprocal property, that is, asymmetry in the
response to mutually antiparallel external fields. A prototypical example of
the nonreciprocal response is the diode effect where the electric conductivity
is allowed or almost forbidden depending on the direction of voltage.

The nonreciprocal optical (NRO) phenomena have been explored in various
systems Fiebig2005 ; Orenstein2021 . The leading NRO response is represented
by the second-order nonlinear response. The formula is given by the current
response induced by the double photo-electric field

$\displaystyle\Braket{\mathcal{J}^{\alpha}_{\text{NRO}}(\omega)}=\int\frac{d\Omega}{2\pi}\sigma^{\alpha;\beta\lambda}\,(\omega;\Omega,\omega-\Omega)E^{\beta}(\Omega)E^{\lambda}(\omega-\Omega).$
(1)

The NRO response covers two important phenomena, namely, second harmonic and
photocurrent generations which are denoted by $\sigma(2\omega;\omega,\omega)$
and $\sigma(0;\omega,-\omega)$, respectively.

The second harmonic generation, where irradiating light with frequency
$\omega$ induces the light oscillating with the doubled frequency $2\omega$,
is a useful tool for probing the microscopic $\mathcal{P}$ breaking in
materials. In the field of multiferroic material science, the second harmonic
generation led to the successful optical imaging of magnetic domains
Fiebig2005 . Furthermore, recent second harmonic generation experiments
implied exotic parity-violating ordering in quantum materials such as a family
of high-temperature copper-based superconductors Zhao2017 ; De_la_Torre2020 .

On the one hand, the photo-induced direct current, namely, the photocurrent
generation (photogalvanic effect), is also of interest. The conventional
mechanism for the photocurrent response is attributed to spatially-
inhomogeneous and asymmetric structures found in photo-diode devices, whereas
recent studies shed light on the bulk photocurrent response originating from
microscopic parity breaking Sturman1992Book . Since the bulk photocurrent
performance is influenced by the electronic property of materials, the
realization of a novel type of photo-electric converter is one of the strong
motivations for materials science. Potential candidate systems are topological
materials Liu2020SemimetalReview ; Orenstein2021  and parity-violating magnets
Ogawa2016 ; Burger2020 . In the former systems, the photocurrent response is
enhanced by the nontrivial quantum geometry in the electronic band structure,
while the latter systems enable a tunable photocurrent response due to the
good controllability of magnetic order.

Figure 1:  Nonreciprocal optics in superconductors. (a) Sketch of
nonreciprocal optical responses in a parity-breaking superconductor. The
Cooper pairs (paired electrons) are colored in green. The injected light
generates light with doubled frequency (second harmonic generation) and direct
current (photocurrent generation). (b) $s+p$-wave superconductivity is an
example of a time-reversal symmetric parity-breaking superconductor. Its
characteristic NRO response is the diverging imaginary component of the NRO
conductivity ($\Im{[\,\sigma_{\text{NRO}}]}$), which is determined by the
Berry curvature derivative $\hat{B_{d}}$. (c) $s+ip$-wave superconductivity is
a parity-time-reversal symmetric parity-breaking superconductor. The
nonreciprocal superfluid density $\hat{f}$ contributes to the diverging NRO
conductivity in the real part ($\Re{[\,\sigma_{\text{NRO}}]}$).

Stimulated by the progress mentioned above, we investigate the NRO responses
in another class of quantum materials, superconductors. Superconductors host
striking electromagnetic properties stemming from the Cooper pairs’ quantum
condensation. Well-known examples are the perfect shielding of magnetic flux
(Meissner effect) and zero resistivity phenomenon. The optical conductivity
also shows a unique spectrum in the low-frequency region and has been
exploited to evaluate the superfluidity of Cooper pairs
Tinkham2004introduction . The nontrivial quantum nature of superconductors is
also attracting from the viewpoint of nonreciprocal optics. The dissipation-
free electric conductivity of superconductors is favorable for suppressing
undesirable energy loss in the photo-electric conversion process. It is
therefore expected that superconductivity paves a new way for optoelectronics.
Supporting this prospect, enhanced nonreciprocal direct conductivity has been
reported in superconductors Tokura2018 ; Ideue2020  and the “perfect” diode
effect has recently been observed in a parity-breaking superconductor Ando2020
. Despite the potential applications to these transport responses, only a few
works were conducted for exploring NRO responses of superconductors Zhao2017 ;
Xu2019 ; Yang2019 ; Nakamura2020 ; Vaswani2020 ; Lim2020 . In particular,
unique features of the NRO response in superconductors remain to be elusive.
To motivate further interest in superconductor-based optoelectronics, in this
paper, we predict a giant NRO response of parity-violating superconductors.

The key features of the NRO response are two-fold. One is the absence of the
definite $\mathcal{P}$ parity in the superconducting state. The other is the
temporal symmetry, such as the time-reversal ($\mathcal{T}$) symmetry and the
symmetry of the combined $\mathcal{P}$ and $\mathcal{T}$ operations, called
the parity-time-reversal ($\mathcal{PT}$) symmetry. Based on the symmetry
analysis and microscopic calculations, in this paper, we demonstrate that the
parity-breaking superconductor shows the diverging NRO response in the low-
frequency regime. This newly discovered response is useful for sensitive
photo-detection and energy harvesting. We also show that the divergent
behavior shows the characteristic power law associated with two indicators,
namely, the Berry curvature derivative $\hat{B_{d}}$ and the nonreciprocal
superfluid density $\hat{f}$, which are closely tied to the temporal symmetry.
Figure 1 overviews this work with a typical example of parity-breaking
superconductors comprising the even-parity $s$-wave and odd-parity $p$-wave
pairings.

The outline of the paper is given below. In Sec. II, we present the symmetry
analysis to identify superconductors offering the NRO responses and classify
potential candidates. Section III briefly explains the formulation of the NRO
response in superconductors. Starting from the general expression obtained
with the Bogoloubov-de Gennes formalism, we propose the superconductivity-
induced NRO response, which we call the anomalous NRO response. The analytical
formulas are in agreement with those obtained by two-fold methods, namely,
density matrix and Green function methods (the derivations are detailed in
Appendices D and E). The results are summarized in Table 1. Furthermore,
making use of the adiabatic time-evolution process, we obtain formulas for the
anomalous nonreciprocal responses which can apply to many-body Hamiltonians
(Sec. IV). The analytical results are corroborated by the numerical
calculations of the NRO responses such as photocurrent and second harmonic
generations in Sec. V. The contents are concluded in Sec. VI.

Throughout this paper, we present formulas with the unit $\hbar=1$ (Dirac
constant), $q=1$ (electron charge).

##  II Symmetry and candidate superconductors

First, we present the symmetry analysis of candidate superconductors for the
nonreciprocal optics. Broken space-inversion symmetry is an essential
ingredient for the NRO response in the superconducting state as well. On the
other hand, the emergence of superconductivity gives the stiffness of the
U(1)-gauge transformation. Therefore, it is relevant to take into account the
combined operations of the usual symmetry operations and U(1)-gauge
transformation. In the case of the odd-parity superconductors, the combined
inversion operation $\mathcal{P}$$\times$U(1) is still preserved although the
usual $\mathcal{P}$ symmetry is broken. The surviving “inversion” symmetry
also prohibits the NRO conductivity. Thus, we cannot expect any NRO responses
in purely odd-parity superconductors, e.g., centrosymmetric ferromagnetic
superconductors such as UGe2, URhGe, and UCoGe Aoki_review2019 . Multiple
superconducting transitions in UPt3 are also expected to be caused by purely
odd-parity pairings Sauls_review1994 ; Joynt_review2002 . This symmetry
requirement is different from that of the diagonal long-range order invoking
parity violation such as structural transitions and magnetic orderings
Watanabe2021 .

According to the above-mentioned symmetry considerations, superconductors
hosting the NRO response have no definite parity under the $\mathcal{P}$
operation. To our best knowledge, the parity violation realizes in the
following setups; (I) superconductivity occurs in the system where
$\mathcal{P}$ symmetry has already been broken in the normal state, (II) an
external field induces parity breaking unique to superconductors, and (III)
multiple superconducting transitions are caused by coexisting even-parity and
odd-parity Cooper pairings.

Case (I) may exist in a broad range of materials since there are a series of
noncentrosymmetric bulk superconductors such as CePt3Si and Li2(Pd,Pt)3B
NCSC_book . In addition, parity-breaking effects have recently been
demonstrated Saito2016 ; Lu2015 ; STO_parity  in artificial two-dimensional
superconductors such as ultrathin MoS2 Ye2012  and SrTiO3/LaAlO3 interfaces
Reyren2007 . Other situations with broken $\mathcal{P}$ symmetry can be found
when superconductivity coexists with the structural transition Hiroi2018  and
magnetic ordering Zhao2017 ; Sumita2017 .

Corresponding to Case (II), an electric current breaks the $\mathcal{P}$ and
$\mathcal{T}$ symmetry. Different from the metal state, the superconducting
state allows the electric current to flow in thermal equilibrium with the aid
of the superfluidity. Under the electric current, Cooper pairs acquire finite
momentum and the Fulde-Ferrell superconductivity is realized FF1964 , which is
classified as the $\mathcal{PT}$ symmetric parity-breaking superconducting
state. Thus, it is expected that superconductors exhibit NRO responses under
the supercurrent flow. The experimental observations have been reported in
Refs. Yang2019 ; Nakamura2020 ; Vaswani2020 .

Finally, we consider Case (III). This class of the parity-breaking
superconductors is intrinsic, nontrivial, and exotic because the
superconductivity itself breaks both of the $\mathcal{P}$ symmetry and its
combination with the U(1) gauge transformation. The resulting parity violation
is solely induced by superconducting instability, and thus it is in sharp
contrast to the known parity-breaking superconductors classified into Cases
(I) and (II). The spontaneously parity-mixed superconductivity may be
stabilized near the topological phase transition between $s$-wave and $p$-wave
superconducting states Wang2017 . It was shown that the relative phase of the
even and odd-parity pair potentials is $\pm\pi/2$ in a coexisting phase space.
The resulting parity-mixed superconductivity preserves $\mathcal{PT}$
symmetry. Although this scenario requires a tunable parameter for the
topological phase transition, the recent discovery of the heavy fermion
superconductor UTe2 Ran2019 , which shows multiple superconducting transitions
ran2019extreme ; Braithwaite2019 ; Aoki2020 ; Knebel2020 ; Lin2020 ;
aoki2021fieldinduced , led to a theoretical prediction of the spontaneously
parity-mixed superconductivity Ishizuka2021 .

Other superconductivity-driven symmetry breaking was also studied in terms of
chiral superconductivity Kallin_2016 ; Brydon2019  and nematic
superconductivity Yonezawa2019 . The time-reversal symmetry breaking has been
discussed based on the $\mu$SR and Kerr rotation measurements Kallin_2016 ,
while the rotation symmetry breaking has been studied by various bulk
measurements Yonezawa2019 . However, experimental probes for detecting
intrinsic space-inversion symmetry breaking are awaited to verify a new phase
of matter. The NRO response may be useful for detecting parity violation
induced by superconductivity and for evaluating its parity-mixing effect. In
the numerical calculations, we demonstrate the NRO responses arising from the
parity violation of Case (III) as an example, while our formulation applies to
all the cases (I)-(III).

##  III Formulation of NRO responses in superconductors

###  III.1 General formula for NRO responses

Our formulation is based on the Bogoliubov-de Gennes (BdG) Hamiltonian

$\mathcal{H}_{\text{BdG}}=\frac{1}{2}\bm{\Psi}^{\dagger}H_{\text{BdG}}\bm{\Psi}+\text{const.},$
(2)

where Nambu spinor $\bm{\Psi}^{\dagger}=(\bm{c}^{\dagger},\bm{c}^{T})$ with
the creation ($\bm{c}^{\dagger}$) and annihilation ($\bm{c}$) operators of
electrons. The BdG Hamiltonian $H_{\text{BdG}}$ consists of the normal-state
Hamiltonian $H_{\text{N}}$ and the pair potential $\Delta$

$H_{\text{BdG}}(\bm{A})=\begin{pmatrix}H_{\text{N}}(\bm{A})&\Delta\\\
\Delta^{\dagger}&-[H_{\text{N}}(\bm{A})]^{T}\end{pmatrix}.$ (3)

The vector potential ${\bm{A}}$ dependence is taken into account for the
electromagnetic field perturbation, while we choose ${\bm{A}}=\bm{0}$ in the
equilibrium. Following the standard perturbative treatments (Appendix A), we
evaluate the expectation value of the electric current density as

$\Braket{\mathcal{J}^{\alpha}(\omega)}=\sum_{n=1}\Braket{\mathcal{J}^{\alpha}(\omega)}_{(n)},$
(4)

where $\Braket{\mathcal{J}^{\alpha}(\omega)}_{(n)}$ is the electric current of
the $n$-th order in $\bm{A}$. With the velocity gauge
$\bm{E}=-\partial_{t}\bm{A}(t)$, we obtain the NRO response function in Eq.
(1) corresponding to the second-order component
$\Braket{\mathcal{J}^{\alpha}(\omega)}_{(n=2)}$. The formula is given by

$\displaystyle\sigma^{\alpha;\beta\gamma}(\omega;\omega_{1},\omega_{2})$
$\displaystyle=\frac{1}{2\left(i\omega_{1}-\eta\right)\left(i\omega_{2}-\eta\right)}\left[\sum_{a}\frac{1}{2}{J}^{\alpha\beta\gamma}_{aa}f_{a}\right.$
(5)
$\displaystyle\left.+\sum_{a,b}\frac{1}{2}\left(\frac{{J}^{\alpha\beta}_{ab}{J}^{\gamma}_{ba}f_{ab}}{\omega_{2}+i\eta-
E_{ba}}+\frac{{J}^{\alpha\gamma}_{ab}{J}^{\beta}_{ba}f_{ab}}{\omega_{1}+i\eta-
E_{ba}}\right)\right.$ (6)
$\displaystyle\left.+\sum_{a,b}\frac{1}{2}\frac{{J}^{\alpha}_{ab}{J}^{\beta\gamma}_{ba}f_{ab}}{\omega+2i\eta-
E_{ba}}\right.$ (7)
$\displaystyle\left.+\sum_{a,b,c}\frac{1}{2}\frac{{J}^{\alpha}_{ab}}{\omega+2i\eta-
E_{ba}}\left(\frac{{J}^{\beta}_{bc}{J}^{\gamma}_{ca}f_{ac}}{\omega_{2}+i\eta-
E_{ca}}-\frac{{J}^{\beta}_{ca}{J}^{\gamma}_{bc}f_{cb}}{\omega_{2}+i\eta-
E_{bc}}\right)\right.$ (8)
$\displaystyle\left.+\sum_{a,b,c}\frac{1}{2}\frac{{J}^{\alpha}_{ab}}{\omega+2i\eta-
E_{ba}}\left(\frac{{J}^{\gamma}_{bc}{J}^{\beta}_{ca}f_{ac}}{\omega_{1}+i\eta-
E_{ca}}-\frac{{J}^{\gamma}_{ca}{J}^{\beta}_{bc}f_{cb}}{\omega_{1}+i\eta-
E_{bc}}\right)\right],$ (9)

where indices $a,b,c$ are spanned by the energy eigenstates of the unperturbed
BdG Hamiltonian $H_{\text{BdG}}(\bm{0})$ and the eigenenergy difference is
$E_{ab}=E_{a}-E_{b}$. We introduced the Fermi-Dirac distribution function
$f_{a}=(e^{\beta E_{a}}+1)^{-1}$ and accordingly defined $f_{ab}=f_{a}-f_{b}$.
The infinitesimal positive parameter $\eta$ appears due to the adiabatic
application of the external field. The electric current operators are defined
by

${J}^{\alpha_{1}\cdots\alpha_{n}}=(-1)^{n}\frac{\partial^{n}H_{\text{BdG}}({\bm{A}})}{\partial
A^{\alpha_{1}}\cdots\partial A^{\alpha_{n}}}\Bigr{|}_{{\bm{A}}=0}.$ (10)

The first- and second-order ones ($n=1,2$) are called paramagnetic and
diamagnetic current density operators, respectively. The detail of the
derivation is shown in Appendix A.

When the pair potential vanishes, the formula for the NRO conductivity agrees
with the prior results for the normal state Passos2018 ; Michishita2020 . We
can check the consistency by replacing ${\bm{A}}$ with ${\bm{k}}$ in the
derivation. The replacement is justified by the minimal coupling
$\bm{p}\rightarrow\bm{p}-q{\bm{A}}$ in the normal state, whereas it fails in
the superconducting state because particles with opposite charges, that is
electron and hole, are treated on equal footing in the BdG formalism.

The failure of the replacement also implies that some useful relations do not
hold for the BdG Hamiltonian. For instance, according to the Hellmann-Feynman
theorem ensured in the normal state, we obtain the relation between the
paramagnetic current operator $J^{\alpha}$ and Berry connection embedded in
the band structure as

$J^{\alpha}_{ab}({\bm{k}})=\frac{\partial}{\partial
k_{\alpha}}E_{a}\delta_{ab}+iE_{ab}\,\xi^{\alpha}_{ab}({\bm{k}}),$ (11)

in the Bloch representation Vanderbilt2018berry . The Berry connection is
defined as

$\xi^{\alpha}_{ab}({\bm{k}})=i\Braket{u_{a}({\bm{k}})}{\frac{\partial
u_{b}({\bm{k}})}{\partial k_{\alpha}}},$ (12)

with the periodic part of the Bloch states $\\{\ket{u_{a}({\bm{k}})}\\}$. This
relation is essential for a remarkable simplification of formulas for
transport and optical phenomena in the normal state Von_Baltz1981 ; Sipe2000 ;
De_Juan2020 ; Watanabe2021 . On the other hand, it is unclear whether the
nonlinear conductivity of superconductors can be computed similarly to the
optoelectronic phenomena in the normal state.

The failure of the Hellmann-Feynman relation in Eq. (11) may make it hard to
formulate the optical conductivity in a physically transparent way. The
difficulty, however, is eliminated by vector potential parametrization. Taking
into account the minimal coupling, we introduce the variational parameter
$\bm{\lambda}$ by

$\mathcal{H}[\hat{\bm{p}}]\rightarrow\mathcal{H}_{\bm{\lambda}}=\mathcal{H}[\hat{\bm{p}}-\bm{\lambda}].$
(13)

Here and hereafter in this section, we drop “BdG” of
$\mathcal{H}_{\mathrm{BdG}}$ and $H_{\mathrm{BdG}}$ for simplicity. Since the
parameter $\bm{\lambda}$ plays the same role as the (spatially-uniform and
time-independent) vector potential, we obtain the relation

${J}^{\alpha}_{ab}=-\lim_{\bm{\lambda}\rightarrow\bm{0}}\Braket{a_{\bm{\lambda}}}{\frac{\partial
H_{\bm{\lambda}}}{\partial\lambda_{\alpha}}}{b_{\bm{\lambda}}},$ (14)

for the paramagnetic current operator. Note that this relation holds even in
the BdG formulation. The Hellmann-Feynman relation for the vector potential is
obtained as

$\Braket{a_{\bm{\lambda}}}{\frac{\partial
H_{\bm{\lambda}}}{\partial\lambda_{\alpha}}}{b_{\lambda}}\left[=-\left({J}_{\bm{\lambda}}^{\alpha}\right)_{ab}\right]=\frac{\partial
E_{a}({\bm{\lambda}})}{\partial\lambda_{\alpha}}\delta_{ab}-i[\xi^{\lambda_{\alpha}},H_{\bm{\lambda}}]_{ab},$
(15)

where we define the connection $\xi^{\lambda_{\alpha}}$ with the replacement
$\partial_{k_{\alpha}}\to\partial_{\lambda_{\alpha}}$ in Eq. (12):

$\xi^{\lambda_{\alpha}}_{ab}=i\Braket{a_{\bm{\lambda}}}{\frac{\partial
b_{\bm{\lambda}}}{\partial\lambda_{\alpha}}}.$ (16)

The eigenstates $\ket{a_{\bm{\lambda}}},\ket{b_{\bm{\lambda}}}$, eigenenergy
$E_{a}({\bm{\lambda}})$, and modified velocity operator
${J}_{\bm{\lambda}}^{\alpha}$ are defined on the basis of $H_{\bm{\lambda}}$.
The quantities defined with the parameter ${\bm{\lambda}}$ become those with
${\bm{\lambda}}=\bm{0}$ in the limit ${\bm{\lambda}}\rightarrow\bm{0}$. We can
apply the Hellmann-Feynman relation to the simplification of the optical
conductivity formulas. For instance, the linear optical conductivity and
photocurrent response are decomposed into the normal and anomalous
contributions with the aid of the ${\bm{\lambda}}$-parametrization (Appendices
C and D).

Although the introduction of the variational parameter reminds us of the
Kohn’s classical work Resta2018 ; Kohn1964 , we note that there is a slight
difference. In Ref. Kohn1964 , the variational parameter is continuous and can
be taken smaller than $2\pi/L$ ($L$ is the linear dimension of the system),
whereas we adopt the $2\pi/L$-discretized values for ${\bm{\lambda}}$ and
finally take the thermodynamic limit $L\to\infty$. In other words, in a finite
system, the variational parameter ${\bm{\lambda}}$ is defined for the
$2\pi/L$-discretized values and the ${\bm{\lambda}}$ derivative is instead
evaluated by the neighboring values for ${\bm{\lambda}}$. The finite size
effect on Eq. (15) vanishes after the thermodynamic limit. Note that all
physical observables have $2\pi/L$-periodic ${\bm{\lambda}}$ dependence in the
normal state, while generally they do not in the superconducting state based
on the BdG formalism with a fixed pair potential. The reason is the following:
In the normal state, the introduction of ${\bm{\lambda}}$ just shifts the wave
number, and ${\bm{\lambda}}$ dependence of physical quantities vanishes after
the $\bm{k}$ summation. By contrast, in the superconducting state,
${\bm{\lambda}}$ corresponds to the center-of-mass momentum of the Cooper
pairs and, thus, the states with different ${\bm{\lambda}}$ should be
distinguished from each other. While the value of ${\bm{\lambda}}$ minimizing
the free energy (denoted by ${\bm{\lambda}}=\bm{0}$ in this paper) is realized
in equilibrium, the derivatives of ${\bm{\lambda}}$ have physical implications
such as the superfluid weight as we see below. Thus, the vector potential
parametrization offers us a good means to extract physical properties
intrinsic to superconductors.

###  III.2 Anomalous photocurrent responses

Emergence of the superconductivity leads to unique optical phenomena which we
call anomalous optical responses. For example, we investigate the photocurrent
response given by the NRO conductivity
$\sigma^{\alpha;\beta\gamma}(0;\Omega,-\Omega)$, and frequency dependence will
be suppressed in this section.

The real and imaginary parts denote the photocurrent conductivity induced by
linearly-polarized and circularly-polarized-lights, respectively
Sturman1992Book . According to the prior studies investigating the normal
state, the mechanism for the photocurrent generation and its dependence on the
polarization state of light are closely related to the temporal symmetry such
as $\mathcal{T}$ and $\mathcal{PT}$ symmetries Von_Baltz1981 ; Sipe2000 ;
Zhang2019 ; De_Juan2020 ; Holder2020 ; Watanabe2021 . The same classification
applies to the superconducting state as revealed in our derivation.

Leaving the derivation to Appendix D, we show the formulas for the NRO
response used in the following calculation. In the gapful superconductors at
low temperatures, the total photocurrent conductivity is given by the two
components;

$\sigma=\sigma_{\text{n}}+\sigma_{\text{a}}.$ (17)

The former ($\sigma_{\text{n}}$) is a normal photocurrent which is a
counterpart of the known photocurrent in the normal state, while the latter
new contribution ($\sigma_{\text{a}}$) represents an anomalous photocurrent
unique to the superconducting state.

The normal part consists of the four contributions

$\sigma_{\text{n}}=\sigma_{\text{Einj}}+\sigma_{\text{Minj}}+\sigma_{\text{shift}}+\sigma_{\text{gyro}},$
(18)

which are termed electric injection current, magnetic injection current, shift
current, and gyration current, respectively. The formulas read

$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{Einj}}=-\frac{i\pi}{8\eta}\sum_{a\neq
b}\left({J}^{\alpha}_{aa}-{J}^{\alpha}_{bb}\right)\Omega_{ba}^{\lambda_{\beta}\lambda_{\gamma}}F_{ab},$
(19)
$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{Minj}}=\frac{\pi}{4\eta}\sum_{a\neq
b}\left({J}^{\alpha}_{aa}-{J}^{\alpha}_{bb}\right)g_{ba}^{\lambda_{\beta}\lambda_{\gamma}}F_{ab},$
(20)
$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{shift}}=-\frac{\pi}{4}\sum_{a\neq
b}\Im{\Biggl{[}[D_{\lambda_{\alpha}}\xi^{\lambda_{\beta}}]_{ab}\xi^{\lambda_{\gamma}}_{ba}+[D_{\lambda_{\alpha}}\xi^{\lambda_{\gamma}}]_{ab}\xi^{\lambda_{\beta}}_{ba}\Biggr{]}}F_{ab},$
(21)
$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{gyro}}=-\frac{i\pi}{4}\sum_{a\neq
b}\Re{\Biggl{[}[D_{\lambda_{\alpha}}\xi^{\lambda_{\beta}}]_{ab}\xi^{\lambda_{\gamma}}_{ba}-[D_{\lambda_{\alpha}}\xi^{\lambda_{\gamma}}]_{ab}\xi^{\lambda_{\beta}}_{ba}\Biggr{]}}F_{ab}.$
(22)

$F_{ab}=f_{ab}\,\delta(\Omega-E_{ba})$ means the Pauli exclusion principle at
the optical transition. We also defined the geometric quantities such as Berry
curvature ($\Omega_{ab}^{\lambda_{\alpha}\lambda_{\beta}}$), quantum metric
($g_{ab}^{\lambda_{\alpha}\lambda_{\beta}}$), and the covariant derivative
($D_{\lambda_{\alpha}}$). The covariant derivative is associated with the
equi-energy Berry connection where bra- and ket-states satisfy $E_{a}=E_{b}$.
Note that we finally take the limit ${\bm{\lambda}}\rightarrow\bm{0}$. The
shift current formula is in agreement with that derived in Ref. Xu2019 . We
also note that the known formulas for the normal photocurrent response
Watanabe2021  are reproduced by the replacement of ${\bm{\lambda}}$ with the
crystal momentum ${\bm{k}}$. The replacement is justified only in the normal
state where the pair potential $\Delta$ is zero.

Because of the inequivalence of ${\bm{\lambda}}$ and ${\bm{k}}$, the anomalous
NRO conductivity appears in the superconducting state. The anomalous part
consists of two contributions

$\sigma_{\text{a}}=\sigma_{\text{NRSF}}+\sigma_{\text{CD}}.$ (23)

The expressions are given by

$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{NRSF}}=\lim_{\bm{\lambda}\rightarrow\bm{0}}-\frac{1}{2\Omega^{2}}\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}\partial_{\lambda_{\gamma}}F_{\bm{\lambda}},$
(24)
$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{CD}}=\lim_{\bm{\lambda}\rightarrow\bm{0}}\frac{1}{4\Omega^{2}}\partial_{\lambda_{\alpha}}\left[\sum_{a\neq
b}{J}^{\beta}_{ab}{J}^{\gamma}_{ba}f_{ab}\left(\frac{1}{\Omega-
E_{ab}}+\frac{1}{E_{ab}}\right)\right].$ (25)

In Eq. (24), $F_{\bm{\lambda}}$ is the free energy of the BdG Hamiltonian. The
second-order derivative of the free energy is the superfluid density
$\rho_{\text{s}}^{\alpha\beta}=\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}F_{\bm{\lambda}}$,
which supports the macroscopic quantum coherence of the superconducting state
Tinkham2004introduction . Thus, $\sigma^{\alpha;\beta\gamma}_{\text{NRSF}}$ is
determined by the nonreciprocal correction to the superfluid density
$\rho_{\text{s}}^{\alpha\beta}$, which we call nonreciprocal superfluid
density
$f^{\alpha\beta\gamma}=\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}\partial_{\lambda_{\gamma}}F_{\bm{\lambda}}$.
The anomalous term in Eq. (24) is rewritten by

$\sigma^{\alpha;\beta\gamma}_{\text{NRSF}}=-\frac{1}{2\Omega^{2}}f^{\alpha\beta\gamma}.$
(26)

The nonreciprocal superfluid density $\hat{f}$ is a totally-symmetric rank-3
tensor. It shows the odd-parity under the $\mathcal{P}$ and $\mathcal{T}$
symmetry operations, whereas the even-parity under the $\mathcal{PT}$
operation.

Another contribution in Eq. (25) is named conductivity derivative effect,
because the formula roughly corresponds to the ${\bm{\lambda}}$-derivative of
the linear optical conductivity. Supporting this argument, the formula in the
low-frequency limit ($\Omega\rightarrow 0$) leads to the
${\bm{\lambda}}$-derivative of the total Berry curvature
$\sum_{a}\epsilon_{\beta\gamma\delta}\Omega^{\lambda_{\delta}}_{a}f_{a}$,
which gives the anomalous Hall conductivity Niu1984 . Thus, we can recast Eq.
(25) as

$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{CD}}\rightarrow\sigma^{\alpha;\beta\gamma}_{\text{sCD}}$
$\displaystyle=\lim_{\bm{\lambda}\rightarrow\bm{0}}\frac{i}{4\Omega}\epsilon_{\beta\gamma\delta}\partial_{\lambda_{\alpha}}\left(\sum_{a}\Omega^{\lambda_{\delta}}_{a}f_{a}\right),$
(27)
$\displaystyle\equiv\frac{i}{4\Omega}\epsilon_{\beta\gamma\delta}B_{d}^{\,\alpha\delta},$
(28)

in the low-frequency regime. Here, we represented the conductivity derivative
effect in the static limit ($\sigma^{\alpha;\beta\gamma}_{\text{sCD}}$) by
Berry curvature derivative $\hat{B}_{d}$ and suppressed $O(\Omega^{0})$ terms.
The temporal symmetry of $\hat{B}_{d}$ is contrasting with the nonreciprocal
superfluid density $\hat{f}$; $\hat{B}_{d}$ is $\mathcal{T}$-even and
$\mathcal{PT}$-odd.

The anomalous photocurrent responses are particularly of interest due to the
giant photo-electric conversion for the low-frequency light. The anomalous
photocurrent response shows high performance for the light with frequency up
to the energy scale of the superconductivity $(\Omega\lesssim\text{0.1\,-\,1
THz}$). The $\mathcal{PT}$ symmetric parity-breaking superconductors show the
divergent response as large as $O(\Omega^{-2})$ due to the nonreciprocal
superfluid density, whereas the $\mathcal{T}$ symmetric superconductors host
the $O(\Omega^{-1})$-diverging response arising from the Berry curvature
derivative. Furthermore, the photocurrent converted from the low-frequency
light may be carried by Cooper pairs and hence flow without undesirable Joule
heating. The energy-saving property originates from the quantum mechanical
nature of the superconductors and cannot be achieved by existing photo-
electric converters such as the p-n junction of semiconductors.

###  III.3 Anomalous contributions to general NRO responses

In Sec. III.2, the anomalous and normal photocurrent responses are detailed on
the basis of the gapped parity-breaking superconducting state. The low-
frequency divergent behavior is similarly observed in other NRO responses such
as the second harmonic generation. In this subsection, the anomalous
contribution to the NRO responses is obtained by the Green function method.
Although the derivation is left to Appendix E, the low-frequency anomalous
contribution reads

$\displaystyle\sigma^{\alpha;\beta\gamma}\left(\omega_{1}+\omega_{2};\omega_{1},\omega_{2}\right)$
$\displaystyle=\sigma^{\alpha;\beta\gamma}_{\text{NRSF}}\left(\omega_{1}+\omega_{2};\omega_{1},\omega_{2}\right)+\sigma^{\alpha;\beta\gamma}_{\text{sCD}}\left(\omega_{1}+\omega_{2};\omega_{1},\omega_{2}\right),$
(29)

where we suppressed $O(\omega_{1}^{a}\omega_{2}^{b})$ terms ($a+b\geq 0$). The
formula is given by the nonreciprocal superfluid density and the static
conductivity derivative effects

$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{NRSF}}=\frac{1}{2\omega_{1}\omega_{2}}f^{\alpha\beta\gamma},$
(30)
$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{sCD}}=-\frac{i}{4}\lim_{\bm{\lambda}\rightarrow\bm{0}}\left(\frac{1}{\omega_{2}}\partial_{\lambda_{\beta}}\sigma_{\alpha\gamma}^{({\bm{\lambda}})}+\frac{1}{\omega_{1}}\partial_{\lambda_{\gamma}}\sigma_{\alpha\beta}^{({\bm{\lambda}})}\right).$
(31)

The tensor $\hat{\sigma}^{({\bm{\lambda}})}$ represents the regular linear
conductivity. Since we formulated the NRO conductivity by the Green function
method with neglecting vertex corrections, the formulas also hold in the
presence of a self-energy correction such as that originates from disorder
scattering. The regular linear conductivity consists of the Drude and Berry
curvature contributions. Thus, when the conductivity tensor
$\sigma_{\alpha\beta}^{({\bm{\lambda}})}$ is decomposed into the symmetric and
antisymmetric components under the permutation $\alpha\leftrightarrow\beta$,
the static conductivity derivative is further decomposed into the Drude
derivative $\hat{D}_{d}$ and Berry curvature derivative $\hat{B}_{d}$

$\displaystyle\partial_{\lambda_{\gamma}}\sigma_{\alpha\beta}^{({\bm{\lambda}})}$
$\displaystyle=\partial_{\lambda_{\gamma}}\left[\frac{1}{2}\left(\sigma_{\alpha\beta}^{({\bm{\lambda}})}+\sigma_{\beta\alpha}^{({\bm{\lambda}})}\right)+\frac{1}{2}\left(\sigma_{\alpha\beta}^{({\bm{\lambda}})}-\sigma_{\beta\alpha}^{({\bm{\lambda}})}\right)\right],$
(32)
$\displaystyle=D_{d}^{\gamma;\alpha\beta}+\epsilon_{\alpha\beta\delta}B_{d}^{\gamma\delta}.$
(33)

The Drude derivative defined by Eq. (237) vanishes without quasiparticle
excitation. In contrast to the Berry curvature derivative ($\mathcal{T}$-even
and $\mathcal{PT}$-odd), the Drude derivative is forbidden by the
$\mathcal{T}$ symmetry but allowed in the $\mathcal{PT}$ symmetric systems. As
a result, the static conductivity derivative effect is rewritten by

$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{sCD}}$
$\displaystyle=-\frac{i}{4}\left[\frac{1}{\omega_{2}}\left(D_{d}^{\beta;\alpha\gamma}+\epsilon_{\alpha\gamma\delta}B_{d}^{\beta\delta}\right)+\frac{1}{\omega_{1}}\left(D_{d}^{\gamma;\alpha\beta}+\epsilon_{\alpha\beta\delta}B_{d}^{\gamma\delta}\right)\right].$
(34)

The obtained expressions reveal that the nonreciprocal superfluid density and
Berry curvature derivative contribute to the general NRO conductivity as well
as the photocurrent generation. Furthermore, the equation (29) includes the
Fermi surface effect such as the Drude derivative contribution, while the
gapped superconductor is assumed in Sec. III.2. If the Fermi-surface
contribution is neglected, we obtain the relation of the Berry curvature
derivative

$\epsilon_{\beta\gamma\delta}B_{d}^{\alpha\delta}=\epsilon_{\alpha\gamma\delta}B_{d}^{\beta\delta}-\epsilon_{\alpha\beta\delta}B_{d}^{\gamma\delta}.$
(35)

Taking $\omega_{1}=-\omega_{2}=\Omega$ in Eq. (29) and assuming the gapped
superconducting state, we reproduce the Berry curvature derivative effect in
Eq. (28).

We comment on the tensor symmetry of the above-mentioned quantities induced by
the parity-breaking superconductivity. First, the nonreciprocal superfluid
density $\hat{f}$ is a totally-symmetric rank-3 polar tensor and thus allowed
in the noncentrosymmetric crystal point groups except for the several high-
symmetric chiral groups where $C_{n}\,(n\geq 4)$ rotation symmetry holds
[$432\,(O)$, $422$ ($D_{4}$), $622$ $(D_{6})$]. Next, we consider the static
conductivity derivative. The Berry curvature derivative $\hat{B}_{d}$ is the
axial rank-2 tensor and hence allowed in the gyrotropic groups
Halasyamani1998noncentrosymmetric . Here, we denote the gyrotropic groups as
the noncentrosymmetric crystal point groups other than the three crystal point
groups $\bar{6}\,(C_{3h}),~{}\bar{6}m2\,(D_{3h}),~{}\bar{4}3m\,(T_{d})$. The
Drude derivative $\hat{D}_{d}$ is the direct product of the rank-1 and rank-2
polar tensors. Thus, it has the same symmetry as that of the piezoelectric
response coefficient and allowed in all the noncentrosymmetric crystal point
groups except for the cubic and chiral group $432\,(O)$.

Finally, we provide a classification table of photocurrent generation (Table
1) to summarize the results obtained by the density matrix (Sec. III.2) and
Green function (Sec. III.3) methods. Mechanism of normal and anomalous
photocurrent is classified based on the temporal symmetry of the system and
the polarization state of light.

Table 1:  Classification of normal and anomalous photocurrent responses. Parities under the $\mathcal{T}$ and $\mathcal{PT}$ symmetry operations are denoted by $\pm$. In the column ‘$\bm{e}$’ we specify the linearly-polarized-light-induced ($\updownarrow$) or circularly-polarized-light-induced ($\circlearrowleft$) photocurrent. The check mark ($\checkmark$) in the column ‘FS (Fermi Surface effect)’ indicates the photocurrent which vanishes in the absence of the Fermi surface. As for the normal photocurrent, we reproduce the classification performed in Ref. Watanabe2021  where the normal state is assumed. On the other hand, the anomalous photocurrent is unique to superconductors. Mechanism | $\mathcal{T}$ | $\mathcal{PT}$ | $\bm{e}$ | FS | Ref.
---|---|---|---|---|---
(Normal photocurrent) |  |  |  |  |
Drude | - | + | $\updownarrow$ | $\checkmark$ | Holder2020
Berry curvature dipole | + | - | $\circlearrowleft$ | $\checkmark$ | Moore2010
Intrinsic Fermi surface (electric) | + | - | $\circlearrowleft$ | $\checkmark$ | De_Juan2020
Intrinsic Fermi surface (magnetic) | - | + | $\updownarrow$ | $\checkmark$ | Watanabe2021
Injection current (electric) | + | - | $\circlearrowleft$ |  | Sipe2000
Injection current (magnetic) | - | + | $\updownarrow$ |  | Zhang2019
Shift current | + | - | $\updownarrow$ |  | Von_Baltz1981 ; Sipe2000
Gyration current | - | + | $\circlearrowleft$ |  | Watanabe2021 ; Ahn2020
(Anomalous photocurrent) |  |  |  |  |
Nonreciprocal superfluid density | - | + | $\updownarrow$ |  | This work
Berry curvature derivative | + | - | $\circlearrowleft$ |  | This work
Drude derivative | - | + | $\circlearrowleft$ | $\checkmark$ | This work

##  IV Generalization to many-body systems

In Sec. III, we derived the anomalous NRO responses from the second-order
perturbation theory. In this section, we derive a more general formula for
anomalous responses unique to superconductors, which applies to interacting
systems.

We consider the many-body Hamiltonian comprising parameters $\bm{\lambda}$
which adiabatically evolve in time. We assume that the ground state is
separated from the excited states by a finite energy gap. We also assume the
zero temperature. The secular equation of the Hamitonian
$H=H({\bm{\lambda}}(t))$ reads

$H({\bm{\lambda}}(t))\ket{a_{\bm{\lambda}}}=\varepsilon_{a}({\bm{\lambda}})\ket{a_{\bm{\lambda}}},$
(36)

where the instantaneous eigenstates and eigenvalues are introduced. If the
many-body wavefunction is initially in the ground state $\ket{0}$ and
adiabatically evolves in time, the wavefunction is obtained as Thouless1983

$\displaystyle\ket{\psi(t)}$ $\displaystyle=\ket{0_{\bm{\lambda}}}+\sum_{a\neq
0}\frac{-i}{\varepsilon_{0}({\bm{\lambda}})-\varepsilon_{a}({\bm{\lambda}})}\ket{a_{\bm{\lambda}}}\Braket{a_{\bm{\lambda}}}{\partial_{t}0_{\bm{\lambda}}},$
(37) $\displaystyle=\ket{0_{\bm{\lambda}}}-\sum_{a\neq
0}\frac{1}{\varepsilon_{0}({\bm{\lambda}})-\varepsilon_{a}({\bm{\lambda}})}\ket{a_{\bm{\lambda}}}\Xi_{a0}^{\lambda_{\alpha}}\dot{\lambda_{\alpha}}.$
(38)

with the many-body Berry connection $\Xi_{a0}^{\lambda_{\alpha}}\equiv
i\braket{a_{\bm{\lambda}}}{\partial_{\lambda_{\alpha}}0_{\bm{\lambda}}}$. The
operator $\partial_{\lambda_{\alpha}}H({\bm{\lambda}}(t))$ is evaluated by the
perturbed wavefunction $\ket{\psi(t)}$ as

$\displaystyle\Braket{\frac{\partial
H({\bm{\lambda}}(t))}{\partial\lambda_{\alpha}}}$
$\displaystyle=\Braket{\psi(t)}{\frac{\partial
H({\bm{\lambda}}(t))}{\partial\lambda_{\alpha}}}{\psi(t)},$ (39)
$\displaystyle=\frac{\partial\varepsilon_{0}({\bm{\lambda}})}{\partial\lambda_{\alpha}}+\Omega^{\lambda_{\alpha}\lambda_{\beta}}\dot{\lambda_{\beta}}.$
(40)

The many-body Berry curvature is

$\displaystyle\Omega^{\lambda_{\alpha}\lambda_{\beta}}$
$\displaystyle=2\Im{\,[\Braket{\partial_{\lambda_{\alpha}}0_{\bm{\lambda}}}{\partial_{\lambda_{\beta}}0_{\bm{\lambda}}}]},$
(41) $\displaystyle=\sum_{a\neq
0}i\left(\Xi^{\lambda_{\alpha}}_{0a}\Xi^{\lambda_{\beta}}_{a0}-\Xi^{\lambda_{\alpha}}_{a0}\Xi^{\lambda_{\beta}}_{0a}\right).$
(42)

When we take the vector potential ${\bm{A}}$ as the adiabatic parameter, the
expectation value of the electric current is obtained as

$\displaystyle\Braket{\mathcal{J}^{\alpha}(t)}$
$\displaystyle=-\Braket{\frac{\partial H({\bm{A}}(t))}{\partial A_{\alpha}}},$
(43) $\displaystyle=-\frac{\partial\varepsilon_{0}({\bm{A}})}{\partial
A_{\alpha}}-\Omega^{A_{\alpha}A_{\beta}}\dot{A_{\beta}}.$ (44)

The first term is related to the generalized Drude weight evaluated at the
zero temperature Drudeweight . When the ground state energy
$\varepsilon_{0}({\bm{A}})$ is expanded by the vector potential, we obtain

$\displaystyle\varepsilon_{0}({\bm{A}})$
$\displaystyle=\varepsilon_{0}(\bm{0})+\frac{1}{2}\left.\frac{\partial^{2}\varepsilon_{0}({\bm{A}})}{\partial
A_{\alpha}\partial A_{\beta}}\right|_{{\bm{A}}=\bm{0}}A_{\alpha}A_{\beta}$
$\displaystyle+\frac{1}{6}\left.\frac{\partial^{3}\varepsilon_{0}({\bm{A}})}{\partial
A_{\alpha}\partial A_{\beta}\partial
A_{\gamma}}\right|_{{\bm{A}}=\bm{0}}A_{\alpha}A_{\beta}A_{\gamma}+O({\bm{A}}^{4}).$
(45)

The first-order term is dropped because the electric current is absent in the
unperturbed ground state. $\partial_{A}^{2}\varepsilon_{0}({\bm{A}})$ and
$\partial_{A}^{3}\varepsilon_{0}({\bm{A}})$ are derivatives of the free energy
in the thermodynamic limit which means linear and nonreciprocal superfluidity.
The agreement between the Drude weight and superfluid density (Meissner
weight) is a property of the gapful superconductor at the zero temperature
Scalapino1993 . Converting the vector potential into the electric field
$\bm{E}=-\dot{{\bm{A}}}$, we obtain the anomalous linear and nonreciprocal
optical responses in the many-body representation. For instance, the response
formula for the nonreciprocal part is written by

$\displaystyle\sigma_{\text{NRSF}}^{\alpha;\beta\gamma}(\omega;\omega_{1},\omega_{2})$
$\displaystyle=\frac{1}{\omega_{1}\omega_{2}}\left.\frac{\partial^{3}\varepsilon_{0}({\bm{A}})}{\partial
A_{\alpha}\partial A_{\beta}\partial A_{\gamma}}\right|_{{\bm{A}}=\bm{0}},$
(46)

where $\omega=\omega_{1}+\omega_{2}$ and the adiabaticity parameter $\eta$ is
suppressed. When we apply the formula to one-body Hamiltonian, Eqs. (30) and
(104) are reproduced.

Similarly, we formulate the Berry curvature derivative effect on the basis of
the many-body Hamilonian. The many-body Berry curvature
$\Omega^{A_{\alpha}A_{\beta}}$ is parametrized by the vector potential, and
the Berry curvature derivative appears in the first-order expansion
coefficient.

$\Omega^{A_{\alpha}A_{\beta}}=\left.\Omega^{A_{\alpha}A_{\beta}}\right|_{{\bm{A}}=\bm{0}}+\left.\frac{\partial\Omega^{A_{\alpha}A_{\beta}}}{\partial
A_{\gamma}}\right|_{{\bm{A}}=\bm{0}}A_{\gamma}+O({\bm{A}}^{2}).$ (47)

When a one-body Hamiltonian is assumed, the many-body Berry curvature leads to
the one-body expression. Accordingly, the first-order expansion coefficient
$\partial_{A_{\gamma}}\Omega^{A_{\alpha}A_{\beta}}$ reproduces the Berry
curvature derivative. The Berry curvature derivative contributes to the
nonreciprocal conductivitity as noteresta .

$\displaystyle\sigma_{\text{sCD}}^{\alpha;\beta\gamma}(\omega;\omega_{1},\omega_{2})$
$\displaystyle=\frac{-i}{2\omega_{2}}\left.\frac{\partial\Omega^{A_{\alpha}A_{\beta}}}{\partial
A_{\gamma}}\right|_{{\bm{A}}=\bm{0}}+\frac{-i}{2\omega_{1}}\left.\frac{\partial\Omega^{A_{\alpha}A_{\gamma}}}{\partial
A_{\beta}}\right|_{{\bm{A}}=\bm{0}}.$ (48)

Thus, the obtained formulas (46) and (48) are consistent with the low-
frequency limit of the anomalous NRO response formulas shown in the previous
section.

Although we here consider the vector potential effect on the electric current
up to the second-order, higher-order anomalous contributions are similarly
obtained as $O(|{\bm{A}}|^{n})$ or $O(|{\bm{A}}|^{n-1}|\dot{{\bm{A}}}|)$ terms
($n\geq 3$), which are higher-order ${\bm{A}}$-derivative of the superfluid
density and Berry curvature, respectively. These contributions may dominate
the nonlinear optical responses such as the high-order harmonic generation in
the low-frequency regime, while the resonant contributions have been
intensively studied Shimano2020 .

Finally, we give a few comments. We generalized the formulas of the anomalous
NRO responses to those for the many-body Hamiltonian [Eqs. (46), (48)]. The
nonreciprocal contributions in Eq. (44) can cover a broad range of
electromagnetic responses, while only the photo-electric field is treated in
Sec. III. For instance, with a spatially-nonuniform vector potential, the
anomalous nonreciprocal contribution leads to the nonreciprocal Meissner
response HDY_meissner .

##  V Model study of anomalous NRO responses

Here, we demonstrate the anomalous NRO responses with numerical calculations.
In this section, we consider the parity-breaking superconductor classified
into Case (III) (see Sec. II), where superconductivity itself breaks the
$\mathcal{P}$ symmetry, by referring to the recent proposal of the spontaneous
parity-mixed superconductivity in a heavy fermion superconductor UTe2
Ishizuka2021 .

With the model Hamiltonian introduced in Sec. V.1, we investigate the normal
and anomalous contributions to the photocurrent response and corroborate our
classification based on the preserved temporal symmetry (Sec. V.2). The
dominant role of the anomalous contribution appears in the numerical
calculations of the photocurrent and second harmonic generations (Sec. V.3).
Finally, we investigate scattering rate dependence of the anomalous NRO
response (Sec. V.4).

###  V.1 Model Hamiltonian

The model Hamiltonian consists of spinful fermions on the two-dimensional
square lattice. In the tight-binding approximation, the normal part
$H_{\text{N}}$ of $H_{\text{BdG}}$ reads

$H^{{\bm{k}}}_{\text{N}}=\left(\varepsilon_{0}-\mu\right)+V\rho_{x}+\bm{g}\cdot\bm{\sigma}\rho_{z},$
(49)

at each crystal momentum ${\bm{k}}$. Pauli matrices $\bm{\sigma}$ and
$\bm{\rho}$ respectively denote the spin and crystal sublattice degrees of
freedom. The model is constructed on the square network where two crystal
sublattices are placed in the checkerboard pattern. The Hamiltonian comprises
the chemical potential $\mu$ and the intra- and intersublattice hoppings
defined as

$\displaystyle\varepsilon_{0}=-4t_{1}\cos{k_{x}}\cos{k_{y}},$ (50)
$\displaystyle V=-2t_{2}\left(\cos{k_{x}}+\cos{k_{y}}\right).$ (51)

Assuming the locally noncentrosymmetric structure, the staggered spin-orbit
coupling is given by the local Rashba type KaneMele2005 ; Yanase2014 ;
Zelezny2014

$\bm{g}=4\alpha\left(\cos{k_{x}}\sin{k_{y}},-\sin{k_{x}}\cos{k_{y}},0\right).$
(52)

The parity-mixed pair potential consists of the $s$-wave component
$\psi_{\bm{k}}=\Delta_{\text{e}}$ and $p$-wave component
$\bm{d}_{\bm{k}}=4\Delta_{\text{o}}\sin{k_{x}}\cos{k_{y}}\hat{z}$. For a
$\mathcal{T}$ symmetric parity-breaking superconductor, the pair potential is
taken in the $s+p$-wave form

$\Delta_{\bm{k}}i\sigma_{y}=\psi_{\bm{k}}+\bm{d}_{\bm{k}}\cdot\bm{\sigma},$
(53)

at each crystal momentum ${\bm{k}}$, while the $s+ip$-wave type pair potential
is given by

$\Delta_{\bm{k}}i\sigma_{y}=\psi_{\bm{k}}+i\bm{d}_{\bm{k}}\cdot\bm{\sigma},$
(54)

for a $\mathcal{PT}$ symmetric superconductor. The superconducting state is
the admixture of the $A_{g}$\- and $E_{u}$-type pairings in terms of the
irreducible representation of the point group $D_{4h}$ ($4/mmm$). The
preserved unitary symmetry operations in the parity-breaking superconducting
state are characterized by the $C_{2v}$ ($m2m$) symmetry, whose two-fold
rotation axis is the $y$-axis. This symmetry allows the NRO conductivity,
$\sigma^{y;yy}$, $\sigma^{y;xx}$, $\sigma^{x;xy}$, $\sigma^{x;yx}$,
$\sigma^{z;yz}$, $\sigma^{z;zy}$, and $\sigma^{y;zz}$. The last three
components having the index $z$ are not considered because the two-dimensional
model is adopted.

The parameters are $\mu=-4$, $\alpha=0.3$, $t_{1}=0.6$, $t_{2}=1.0$,
$\Delta_{\text{e}}+4\Delta_{\text{o}}=0.1$, and
${\Delta_{\text{e}}}/{\Delta_{\text{o}}}=4$. For a quantitative estimate, we
set $t_{2}=1\,{\mathrm{\,eV\,}}$ and calculate the response coefficients in
the SI unit. Thus, the numerical results of the NRO conductivity are given in
the unit A$\cdot$V-2. Although the model Hamiltonian is two-dimensional, we
convert the obtained NRO conductivity into the three-dimensional value by
using the thickness of the two-dimensional net $l=1\,{\mathrm{\,nm\,}}$. The
numerical calculations are performed on the $N^{2}$-discretized Brillouin
zone. For numerical convergence, we introduce a phenomenological scattering
rate $\gamma$ (Appendix B) and a finite temperature $T$ for the Fermi-Dirac
distribution function.

###  V.2 Normal and anomalous photocurrent responses

Figure 2:  Frequency dependence of the photocurrent conductivity
$\sigma^{x;xy}$ in (a,b) a $\mathcal{T}$ symmetric parity-breaking
superconductor and (c, d) a $\mathcal{PT}$ symmetric parity-breaking
superconductor. (a) and (c) plot the real part while (b) and (d) are the
imaginary part. The total photocurrent conductivity ($\sigma$) is shown in
red, while the normal and anomalous photocurrent conductivities are shown in
blue and green, respectively. The calculation is performed with $N=4000$ and
$T=10^{-4}$. We also adopted $\gamma=10^{-4}$ for the $\mathcal{T}$ symmetric
case and $\gamma=3\times 10^{-3}$ for the $\mathcal{PT}$ symmetric case.

Here, numerical analysis of the photocurrent conductivity
$\sigma^{x;yx}(0;\Omega,-\Omega)$ is presented. The frequency dependence of
the conductivity is implicit unless otherwise mentioned. Figure 2 plots the
frequency dependence of the photocurrent conductivity $\sigma^{x;yx}$ under
the linearly-polarized and circularly-polarized lights. To distinguish the
normal and anomalous photocurrent responses, we calculate each contribution as
well as the total photocurrent conductivity.

In the case of the $\mathcal{T}$ symmetric system, the real part
$\Re{\,[\sigma^{x;yx}]}$ for the linearly-polarized light in Fig. 2(a) stems
from the normal photocurrent mechanism, that is, the shift current
$(\sigma_{\text{shift}})$ Von_Baltz1981 ; Sipe2000 . Since the normal
photocurrent requires quasiparticle excitations, the contribution vanishes in
the low-frequency regime $(\Omega\lesssim 0.01)$ due to the superconducting
gap. On the other hand, the imaginary part $\Im{\,[\sigma^{x;yx}]}$ for the
circularly-polarized light is solely determined by the anomalous contribution,
namely the conductivity derivative $(\sigma_{\text{CD}})$ [Fig. 2(b)]. In
agreement with Eq. (28), we clearly observe the low-frequency divergence
proportional to $\Omega^{-1}$ in the circular photocurrent conductivity
$\Im{\,[\sigma^{x;yx}]}$. Although the normal photocurrent called the electric
injection current $(\sigma_{\text{Einj}})$ also contributes to the circular
photocurrent conductivity, it gives a negligible contribution compared to the
anomalous contribution. This result implies that the anomalous mechanism
provides distinctive optoelectronic properties of parity-breaking
superconductors.

Similarly, it is shown that the anomalous photocurrent gives rise to the
dominant photocurrent creation in the $\mathcal{PT}$ symmetric parity-breaking
superconductors [Fig. 2(c, d)]. Owing to the $\mathcal{PT}$ symmetry, the
situation is in contrast to the $\mathcal{T}$ symmetric superconductors. The
circular photocurrent conductivity $\Im{\,[\sigma^{x;yx}]}$ is determined by
the normal photocurrent called gyration current $(\sigma_{\text{gyro}})$
Ahn2020 ; Watanabe2021  [Fig. 2(d)], while the linear photocurrent
$\Re{\,[\sigma^{x;yx}]}$ is given by the anomalous contribution, that is due
to the nonreciprocal superfluid density effect $(\sigma_{\text{NRSF}})$, with
a negligible normal photocurrent originating from the magnetic injection
current $(\sigma_{\text{Minj}})$ Zhang2019  [Fig. 2(c)]. The low-frequency
divergence of $\Re{\,[\sigma^{x;yx}]}$ is inversely proportional to
$\Omega^{2}$ in agreement with the formula (26).

Here, we compare the photocurrent conductivity and the numerically calculated
nonreciprocal superfluid density to make the analytical result in Eq. (26)
more compelling. Figure 3 shows the comparison between
$\Re{[\,\Omega^{2}\sigma^{x;yx}]}$ and the nonreciprocal superfluid density
$f^{xyx}$. In agreement with the analytical result, the NRO conductivity
multiplied by the square of frequency $\Omega^{2}$ asymptotically approaches
to the nonreciprocal superfluid density in the low-frequency limit. Since
$\Re{[\,\sigma^{x;yx}]}$ includes a non-divergent part of the conductivity
derivative effect in Eq. (25), the dispersion in
$\Re{[\,\Omega^{2}\sigma^{x;yx}]}\propto\Omega^{2}$ is observed. The agreement
between the analytical and numerical results indicates that the nonreciprocal
superfluid density is an indicator of the low-frequency NRO response.
Conversely, the NRO measurements such as the second harmonic generation and
photocurrent response may provide a quantitative estimate of the nonreciprocal
superfluid density.

Figure 3:  Comparison between the photocurrent conductivity and numerically
calculated nonreciprocal superfluid density. The photocurrent conductivity
$\Re{[\,\sigma^{x;yx}]}$ (blue solid line) and the nonreciprocal superfluid
density $f^{xyx}$ (red dashed line) are calculated by Eqs. (17) and (173),
respectively. We take $T=10^{-4}$ and $N=3000$. The calculation of
$\Re{[\,\sigma^{x;yx}]}$ is performed with $\gamma=10^{-3}$.

###  V.3 Anomalous contributions to NRO responses

The divergent NRO responses are similarly observed in the other photocurrent
conductivity components and in other NRO responses such as the second harmonic
generation. In this section, we show the frequency dependence of all the
allowed components of the photocurrent and second harmonic generation
coefficients. We distinguish the two response functions by the second harmonic
generation coefficient
$\sigma^{\alpha;\beta\gamma}_{\text{SHG}}=\sigma^{\alpha;\beta\gamma}(2\Omega;\Omega,\Omega)$
and the photocurrent conductivity
$\sigma^{\alpha;\beta\gamma}_{\text{PC}}=\sigma^{\alpha;\beta\gamma}(0;\Omega,-\Omega)$.
The numerical calculations are performed with the total NRO conductivity
formula [Eqs. (5) to (9)] where frequencies are taken by
$\omega_{1}=\pm\omega_{2}=\Omega$.

First, we consider the $\mathcal{T}$ symmetric parity-mixed superconductor.
The low-frequency divergence manifests in the imaginary part of the NRO
responses [Figs. 4(b) and 4(d)]. The divergence originates from the Berry
curvature derivative $\hat{B}_{d}$. In agreement with the analytical
expression in Eq. (34),
$\Im{[\,\sigma_{\text{SHG}}^{x;xy}]}=\Im{[\,\sigma_{\text{SHG}}^{x;yx}]}$ and
$\Im{[\,\sigma_{\text{PC}}^{x;xy}]}=-\Im{[\,\sigma_{\text{PC}}^{x;yx}]}$
converge to the same value in the low-frequency limit [see green-colored plots
in Figs. 4(b) and 4(d)]. $\Im{[\,\sigma_{\text{SHG}}^{y;yy}]}$ does not show
divergence in the low-frequency regime, since the vanishing totally-
antisymmetric tensor $\epsilon_{yy\delta}=0$ forbids the Berry curvature
derivative effect. Note that the divergence does not appear in the real part
of the NRO conductivity except for an artificial contribution arising from the
phenomenological scattering rate as discussed in Sec. V.4.

Figure 4:  Frequency dependence of the second harmonic and photocurrent
generation in a $\mathcal{T}$ symmetric superconducting state. (a)
$\Re{[\,\sigma_{\text{SHG}}^{\alpha;\beta\gamma}]}$ and (b)
$\Im{[\,\sigma_{\text{SHG}}^{\alpha;\beta\gamma}]}$ for the second harmonic
generation. (c) $\Re{[\,\sigma_{\text{PC}}^{\alpha;\beta\gamma}]}$ and (d)
$\Im{[\,\sigma_{\text{PC}}^{\alpha;\beta\gamma}]}$ for the photocurrent
generation. We take $\gamma=10^{-4}$, $T=10^{-4}$, and $N=4000$.

Next, we look into the NRO response functions for the $\mathcal{PT}$ symmetric
parity-mixed superconductor. The low-frequency divergence is observed in the
real part of the NRO responses [Figs. 5(a) and 5(c)] due to the nonreciprocal
superfluid density $\hat{f}$. Since $\hat{f}$ is the totally-symmetric tensor,
$f^{xxy}=f^{yxx}$ leads to the same low-frequency divergence in
$\Re{[\,\sigma^{x;xy}]}$ and $\Re{[\,\sigma^{y;xx}]}$, while
$\Re{[\,\sigma_{\text{SHG}}^{x;xy}]}=\Re{[\,\sigma_{\text{SHG}}^{x;yx}]}$ and
$\Re{[\,\sigma_{\text{PC}}^{x;xy}]}=\Re{[\,\sigma_{\text{PC}}^{x;yx}]}$ are
satisfied by definition.

Figure 5:  Frequency dependence of the second harmonic and photocurrent
generation in a $\mathcal{PT}$ symmetric superconducting state. (a)
$\Re{[\,\sigma_{\text{SHG}}^{\alpha;\beta\gamma}]}$ and (b)
$\Im{[\,\sigma_{\text{SHG}}^{\alpha;\beta\gamma}]}$ for the second harmonic
generation. (c) $\Re{[\,\sigma_{\text{PC}}^{\alpha;\beta\gamma}]}$ and (d)
$\Im{[\,\sigma_{\text{PC}}^{\alpha;\beta\gamma}]}$ for the photocurrent
generation. We take $\gamma=3\times 10^{-3}$, $T=10^{-4}$, and $N=4000$.

###  V.4 Scattering rate dependence of anomalous NRO responses

We address the scattering rate dependence of the anomalous NRO responses by
taking the photocurrent conductivity as an example.

For the anomalous NRO responses in the low-frequency regime, numerical results
show the diverging behavior. $\Re{[\sigma^{x;yx}]}$ ($\Im{[\sigma^{x;yx}]}$)
in the low-frequency regime is, however, influenced by the scattering rate for
the $\mathcal{T}$ symmetric ($\mathcal{PT}$ symmetric) parity-breaking
superconductor. In contrast, $\Im{[\sigma^{x;yx}]}$ ($\Re{[\sigma^{x;yx}]}$)
for the $\mathcal{T}$ symmetric ($\mathcal{PT}$ symmetric) parity-breaking
superconductor is intrinsic and hence robust to scattering-rate effects. These
contrasting properties are revealed by calculating the photocurrent
conductivity with scattering rates ranging over several orders of magnitude
(Fig. 6). For instance, in Fig. 6(a) for the $\mathcal{T}$ symmetric case,
$\Re{[\sigma^{x;yx}]}$ disappears with decreasing the scattering rate, whereas
$\Im{[\sigma^{x;yx}]}$ shows the remarkable tolerance to the scattering
effects. According to Table 1, the scattering-rate-sensitive part may be
artificial or extrinsic. Thus, we did not show such part in the previous
subsections.

In the $\mathcal{T}$ symmetric system, we do not have any diverging linearly-
polarized-light-induced photocurrent response when we take into account
scattering effects by implementing the Green function method [Eq. (29)]. Thus,
the low-frequency divergence of $\Re{[\sigma^{x;yx}]}$ is an artificial effect
arising from the phenomenology of scattering effects. On the other hand, a
tiny circularly-polarized-light-induced photocurrent response in the
$\mathcal{PT}$ symmetric system [Fig. 6(d)] may be an extrinsic effect due to
the Drude derivative effect. We can identify a genuine extrinsic contribution
by treating the scattering effect in a rigorous way beyond the
phenomenological scattering rate approximation. Such rigorous calculation is
left for future study. The main finding of this numerical work is the
intrinsic low-frequency divergence of the imaginary part of the NRO
conductivity in the $\mathcal{T}$ symmetric superconductor and the real part
in the $\mathcal{PT}$ symmetric superconductor, which are robust to scattering
effects. The scattering-tolerant property ensures that the Berry curvature
derivative $\hat{B}_{d}$ and the nonreciprocal superfluid density $\hat{f}$
are promising indicators of the NRO phenomena in superconductors.

Figure 6:  Scattering rate dependence of the total photocurrent conductivity
$\sigma^{x;yx}$. Frequency dependence of the photocurrent conductivity is
plotted for (a, b) the case of $\mathcal{T}$ symmetric superconductors and for
(c, d) the case of $\mathcal{PT}$ symmetric superconductors. In each panel, a
series of scattering rates are distinguished by color. In (b) and (c), all the
plots are overlapped, indicating the tolerance to scattering effects. We take
the temperature $T=10^{-4}$ and $N=1500$.

##  VI Discussion and Summary

This paper elaborated the nonreciprocal optical responses characteristic to
parity-breaking superconductors. Based on the analytical and numerical
calculations, we showed that the anomalous NRO responses play essential roles
in superconductors. The key ingredients, namely, nonreciprocal superfluid
density and Berry curvature derivative, manifest in the superconducting state,
although they have no counterpart in the normal state. Since the anomalous NRO
response is intrinsically divergent under the low-frequency light such as sub-
terahertz light, it surpasses the NRO responses of other origins clarified in
prior studies Sturman1992Book . Our microscopic calculations confirmed that
the anomalous NRO responses dominate the photocurrent and second harmonic
generation responses. For instance, for the linearly-polarized-light-induced
photocurrent in the $\mathcal{PT}$ symmetric parity-breaking superconductor,
the normal photocurrent called injection current is negligible compared to the
nonreciprocal superfluid density effect [Fig. 2(c)].

A promising application of our work is a probe for superconducting symmetry.
The parity-mixing effect of the superconducting pairings has been drawing much
attention in superconducting science from the discovery of noncentrosymmetric
superconductors NCSC_book . The NRO responses can appear in the parity-mixed
superconductors and thus provide an estimate of the parity-mixing effect. The
estimate may be demonstrated by incorporating chemical or physical operation,
which tunes the parity-mixing effect Harada2012 ; Kitagawa2020 . Furthermore,
the anomalous NRO response is in close relation with the temporal symmetry
preserved in superconductors. Indeed, the photocurrent response induced by the
circularly-polarized (linearly-polarized) light shows the characteristic
divergence for the low-frequency light when the parity-breaking superconductor
preserves $\mathcal{T}$ ($\mathcal{PT}$) symmetry. Recent studies implied the
parity violation in exotic quantum materials such as cuprates Zhao2017 ;
De_la_Torre2020 ; Lim2020  and heavy-fermion systems Braithwaite2019 ;
Ishizuka2021 . Thus, our classification of the NRO responses will be helpful
for identifying the spatial and temporal symmetry of these materials.

It is also essential to explore novel functional devices utilizing the
remarkable properties of the anomalous NRO response in superconductors. Owing
to the enhanced NRO response for the low-frequency light, the parity-breaking
superconductor may realize a highly efficient optical apparatus for future
terahertz light technology. The anomalous NRO response arises from the
condensed Cooper pairs, and hence the photo-electric conversion is expected to
be carried out without undesirable Joule heating. The photo-diode without
energy loss possibly originates from the quantum nature of superconductivity.
Interestingly, the energy-saving NRO performance of superconductors can be
turned on by injecting the supercurrent Nakamura2020  [Case (II) in Sec. II].
The tunability is a consequence of the superfluidity and cannot be found in
systems without macroscopic quantum condensate. The supercurrent-induced
parity violation can be induced in various superconductors such as
prototypical superconducting alloys and high-temperature copper-based
superconductors. As a result, superconductivity-based optoelectronics are
available in a broad range of materials.

Our systematic study of the nonreciprocal optical response is expected to
stimulate further investigations of superconductors in the scientific and
engineering fields and make essential building blocks for superconductor-based
optoelectronics.

## Acknowledgements

We thank Kyosuke Adachi, Shota Kanasugi, Hiroto Tanaka, and Yoshihiro
Michishita for fruitful discussions. H.W. is grateful to Naoto Nagaosa for
enlightening comments. This work was supported by JSPS KAKENHI (Grants No.
JP18H05227, No. JP18H01178, and No. 20H05159), SPIRITS 2020 of Kyoto
University. H.W. is a JSPS research fellow and supported by JSPS KAKENHI
(Grant No. 18J23115 and No. 21J00453). A.D. is supported by JSPS KAKENHI
(Grant No. 21K13880).

##  Appendix A Density matrix theory based on Bogoliubov-de Gennes Hamiltonian

This section is devoted to the density matrix formulation on the basis of the
BdG Hamiltonian. We note that some expressions are overlapped with those in
the main text.

The Hamiltonian is given by the BdG form,

$\mathcal{H}_{\text{BdG}}=\frac{1}{2}\bm{\Psi}^{\dagger}H_{\text{BdG}}\bm{\Psi}+\text{const.},$
(55)

where Nambu spinor $\bm{\Psi}^{\dagger}=(\bm{c}^{\dagger},\bm{c}^{T})$ with
the creation ($\bm{c}^{\dagger}$) and annihilation ($\bm{c}$) operators of
electrons. Here, the components of $\bm{c}$ are the normal-state degrees of
freedom, namely, spatial position, sublattice, orbital, spin, and so on. The
BdG Hamiltonian $H_{\text{BdG}}$ includes the Hamiltonian in the normal state
$H_{\text{N}}$ and the pair potential of superconductivity $\Delta$

$H_{\text{BdG}}(\bm{A})=\begin{pmatrix}H_{\text{N}}(\bm{A})&\Delta\\\
\Delta^{\dagger}&-H_{\text{N}}(\bm{A})^{T}\end{pmatrix}.$ (56)

The electromagnetic perturbation is taken into account by the vector-potential
$\bm{A}$. In equilibrium, we take $\bm{A}=0$.

Owing to the particle-hole symmetry of $H_{\text{BdG}}$, the spinor satisfies
the relation

$[\mathcal{H}_{\text{BdG}},\Psi_{i}]=-\left(H_{\text{BdG}}\right)_{ij}\Psi_{j}.$
(57)

Here, the index $i,j$ run over both the normal and the Nambu degrees of the
freedom.

It is convenient to retake the basis of the Nambu spinor $\Psi_{i}$ to the
energy eigenstates of $H_{\text{BdG}}$,
$[H_{\text{BdG}},\Phi_{a}]=-E_{a}\Phi_{a}$. This commutation relation
indicates that the density-matrix formalism straightforwardly applies to the
BdG Hamiltonian as follows. We introduce the density matrix $\hat{\rho}(t)$ as
in the normal state. Without perturbation, it is written as

$\hat{\rho}^{(0)}=e^{-\beta\mathcal{H}_{\text{BdG}}}/\operatorname{Tr}[e^{-\beta\mathcal{H}_{\text{BdG}}}],$
(58)

and the time evolution of $\hat{\rho}(t)$ is described by the von-Neumann
equation

$i\partial_{t}\hat{\rho}(t)=[\mathcal{H}_{\text{BdG}}+\Delta\mathcal{H}(t),\hat{\rho}(t)].$
(59)

Here, $\Delta\mathcal{H}(t)=\bm{\Psi}^{\dagger}\Delta H(t)\bm{\Psi}/2$ is the
perturbative Hamiltonian. It is convenient to introduce the reduced density
matrix described by the Bogoliubov quasiparticles

$\rho_{ab}(t)=\operatorname{Tr}[\Phi_{b}^{\dagger}\Phi_{a}\hat{\rho}(t)],$
(60)

and accordingly obtain the von-Neumann equation

$\displaystyle i\partial_{t}\rho_{ab}=E_{ab}\rho_{ab}(t)+\left[\Delta
H(t),{\rho}(t)\right]_{ab},$ (61)

where $E_{ab}=E_{a}-E_{b}$ and $[A,B]_{ab}=A_{ac}B_{cb}-A_{cb}B_{ac}$. In
equilibrium, the density matrix $\rho_{ab}^{(0)}(t)=f_{a}\delta_{ab}$ is
reduced to the Fermi-Dirac distribution function, $f_{a}\equiv(e^{\beta
E_{a}}+1)^{-1}$. Throughout this work, we consider the photoelectric field
expressed in the velocity gauge. The perturbative Hamiltonian is given by
$\Delta H(t)=H_{\mathrm{BdG}}(\bm{A}(t))-H_{\mathrm{BdG}}(\bm{0})$, that is,

$\Delta H(t)=\sum_{n=1}\frac{1}{n!}\left(-1\right)^{n}A^{\alpha_{1}}(t)\cdots
A^{\alpha_{n}}(t){J}^{\alpha_{1}\cdots\alpha_{n}},$ (62)

where we introduce the generalized velocity operator

${J}^{\alpha_{1}\cdots\alpha_{n}}=(-1)^{n}\frac{\partial^{n}H_{\text{BdG}}({\bm{A}})}{\partial
A^{\alpha_{1}}\cdots\partial A^{\alpha_{n}}}\Bigr{|}_{{\bm{A}}=0}.$ (63)

For the BdG Hamiltonian $H_{\text{BdG}}({\bm{A}})$, the vector potential is
coupled to the diagonal part in the form of the minimal coupling [Eq. (56)].
Note that we do not take into account field-induced corrections to the pair
potential for simplicity. Integrating Eq. (61) iteratively, we obtain the
perturbed density matrix $\hat{\rho}(t)$. Then, we immediately obtain the
expectation value of the electric current as

$\Braket{\mathcal{J}^{\alpha}}(t)=\operatorname{Tr}[\hat{\rho}(t)\mathcal{J}^{\alpha}(t)],$
(64)

with the velocity operator

$\mathcal{J}^{\alpha}(t)=\frac{1}{2}\bm{\Psi}^{\dagger}{J}^{\alpha}_{\bm{A}}(t)\bm{\Psi}.$
(65)

The matrix denote the modified velocity operator

${J}_{\bm{A}}^{\alpha}(t)=\sum_{m=0}\frac{1}{m!}\left(-1\right)^{m}A^{\beta_{1}}(t)\cdots
A^{\beta_{m}}(t){J}^{\alpha\,\beta_{1}\cdots\beta_{m}}.$ (66)

Here, we give the perturbed density matrix by expanding the density matrix in
Eq. (61) as

$\rho_{ab}(t)=\sum_{n=0}\rho_{ab}^{(n)}(t).$ (67)

$\rho_{ab}^{(n)}$ denotes the $O(|\bm{A}|^{n})$ correction to the density
matrix. Accordingly, the perturbed density matrix up to the second-order are
obtained as

$\displaystyle\left(i\partial_{t}-E_{ab}\right)\rho_{ab}^{(0)}(t)=0,$ (68)
$\displaystyle\left(i\partial_{t}-E_{ab}\right)\rho_{ab}^{(1)}(t)=-A^{\alpha}(t)[{J}^{\alpha},\rho^{(0)}(t)]_{ab},$
(69)
$\displaystyle\left(i\partial_{t}-E_{ab}\right)\rho_{ab}^{(2)}(t)=-A^{\alpha}(t)[{J}^{\alpha},\rho^{(1)}(t)]_{ab}+\frac{1}{2}A^{\alpha}(t)A^{\beta}(t)[{J}^{\alpha\beta},\rho^{(0)}(t)]_{ab}.$
(70)

Taking the convention for Fourier transformation as

$X(t)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}X(\omega)e^{-i\omega t-\eta
t},$ (71)

we obtain the perturbed density matrix in the frequency domain. The
adiabaticity parameter $\eta$ is positive infinitesimal. The density matrix
are $\rho_{ab}^{(0)}(\omega)=f_{a}\delta_{ab}~{}2\pi\delta(\omega)$ and

$\displaystyle\rho_{ab}^{(1)}(\omega)$ $\displaystyle=-\frac{1}{\omega+i\eta-
E_{ab}}A^{\alpha}(\omega){J}^{\alpha}_{ab}f_{ba},$ (72)
$\displaystyle\rho_{ab}^{(2)}(\omega)$ $\displaystyle=\frac{1}{\omega+2i\eta-
E_{ab}}\int\frac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}(2\pi)\delta(\omega-\omega_{1}-\omega_{2})$
$\displaystyle\times\left(-A^{\alpha}(\omega_{1})[{J}^{\alpha},\rho^{(1)}(\omega_{2})]_{ab}+\frac{1}{2}A^{\alpha}(\omega_{1})A^{\beta}(\omega_{2}){J}^{\alpha\beta}_{ab}f_{ba}\right).$
(73)

The expectation value of the electric current is calculated by the current
operator [Eq. (63)] and the perturbed density matrix,

$2\Braket{\mathcal{J}^{\alpha}(t)}=2\operatorname{Tr}{[\mathcal{J}^{\alpha}(t)\hat{\rho}(t)]}=\sum_{a,b}{J}_{ab}^{\alpha}(t)\rho_{ba}(t).$
(74)

The perturbed electric current operator reads

$\displaystyle\mathcal{J}^{\alpha}(t)$
$\displaystyle=\mathcal{J}_{(0)}^{\alpha}(t)+\mathcal{J}_{(1)}^{\alpha}(t)+\mathcal{J}_{(2)}^{\alpha}(t)+\cdots,$
(75) $\displaystyle\mathcal{J}_{(0)}^{\alpha}(t)$
$\displaystyle=\int\frac{d\omega}{2\pi}e^{-i\omega
t}\,2\pi\delta(\omega)\mathcal{J}^{\alpha},$ (76)
$\displaystyle\mathcal{J}_{(1)}^{\alpha}(t)$
$\displaystyle=-\int\frac{d\omega}{2\pi}e^{-i\omega t+\eta
t}\mathcal{J}^{\alpha\beta}A^{\beta}(\omega),$ (77)
$\displaystyle\mathcal{J}_{(2)}^{\alpha}(t)$
$\displaystyle=\int\frac{d\omega}{2\pi}e^{-i\omega t+2\eta
t}\int\frac{d\omega_{1}d\omega_{2}}{2\pi}\delta(\omega_{1}+\omega_{2}-\omega)\frac{1}{2}\mathcal{J}^{\alpha\beta\gamma}A^{\beta}(\omega_{1})A^{\gamma}(\omega_{2}).$
(78)

Before proceeding to the second-order nonlinear optical conductivity, let us
consider the linear optical conductivity. The linear response is given by

$\displaystyle\Braket{2\mathcal{J}^{\alpha}(\omega)}_{(1)}$
$\displaystyle=\int dt\ e^{-i\omega t-\eta
t}\operatorname{Tr}{[2\mathcal{J}^{\alpha}_{(1)}(t)\hat{\rho}^{(0)}(t)+2\mathcal{J}^{\alpha}_{(0)}(t)\hat{\rho}^{(1)}(t)]},$
(79)
$\displaystyle=\sum_{a,b}\left(\frac{-{J}^{\alpha}_{ab}{J}^{\beta}_{ba}f_{ab}}{\omega+i\eta-
E_{ba}}A^{\beta}(\omega)-A^{\beta}(\omega){J}^{\alpha\beta}_{ab}f_{a}\delta_{ab}\right).$
(80)

In the last line, the former and the latter originate from the correction to
the density matrix and current operator, which are so-called paramagnetic and
diamagnetic contributions, respectively. With the velocity gauge
$\bm{E}=-\partial_{t}\bm{A}(t)$, the linear optical conductivity is obtained
as

$2\sigma^{\alpha\beta}(\omega)=\frac{1}{i\omega-\eta}\sum_{a,b}\left(\frac{-{J}^{\alpha}_{ab}{J}^{\beta}_{ba}f_{ab}}{\omega+i\eta-
E_{ba}}-{J}^{\alpha\beta}_{ab}f_{a}\delta_{ab}\right).$ (81)

Similarly, the second-order electic current is divided into three components,

$\displaystyle\Braket{\mathcal{J}^{\alpha}(\omega)}_{(2)}$ $\displaystyle=\int
dt\ e^{-i\omega t-\eta
t}\operatorname{Tr}{[\mathcal{J}^{\alpha}_{(2)}(t)\hat{\rho}^{(0)}(t)+\mathcal{J}^{\alpha}_{(1)}(t)\hat{\rho}^{(1)}(t)+\mathcal{J}^{\alpha}_{(0)}(t)\hat{\rho}^{(2)}(t)]}.$
(82)

The second-order optical conductivity is explicitly given by

$\displaystyle 2\sigma^{\alpha;\beta\gamma}(\omega;\omega_{1},\omega_{2})$
$\displaystyle=\frac{1}{i\omega_{1}-\eta}\frac{1}{i\omega_{2}-\eta}\left[\sum_{a}\frac{1}{2}{J}^{\alpha\beta\gamma}_{aa}f_{a}\right.$
(83)
$\displaystyle\left.+\sum_{a,b}\frac{1}{2}\left(\frac{{J}^{\alpha\beta}_{ab}{J}^{\gamma}_{ba}f_{ab}}{\omega_{2}+i\eta-
E_{ba}}+\frac{{J}^{\alpha\gamma}_{ab}{J}^{\beta}_{ba}f_{ab}}{\omega_{1}+i\eta-
E_{ba}}\right)\right.$ (84)
$\displaystyle\left.+\sum_{a,b}\frac{1}{2}\frac{{J}^{\alpha}_{ab}{J}^{\beta\gamma}_{ba}f_{ab}}{\omega+2i\eta-
E_{ba}}\right.$ (85)
$\displaystyle\left.+\sum_{a,b,c}\frac{1}{2}\frac{{J}^{\alpha}_{ab}}{\omega+2i\eta-
E_{ba}}\left(\frac{{J}^{\beta}_{bc}{J}^{\gamma}_{ca}f_{ac}}{\omega_{2}+i\eta-
E_{ca}}-\frac{{J}^{\beta}_{ca}{J}^{\gamma}_{bc}f_{cb}}{\omega_{2}+i\eta-
E_{bc}}\right)\right.$ (86)
$\displaystyle\left.+\sum_{a,b,c}\frac{1}{2}\frac{{J}^{\alpha}_{ab}}{\omega+2i\eta-
E_{ba}}\left(\frac{{J}^{\gamma}_{bc}{J}^{\beta}_{ca}f_{ac}}{\omega_{1}+i\eta-
E_{ca}}-\frac{{J}^{\gamma}_{ca}{J}^{\beta}_{bc}f_{cb}}{\omega_{1}+i\eta-
E_{bc}}\right)\right],$ (87)

in which the first line, second line, and third to fifth lines correspond to
the first, second, and third components in Eq. (82), respectively. Note that
only Eqs. (86) and (87) appear in electron systems in the continuum space
Von_Baltz1981 . We confirmed that the obtained formulas are consistent with
those in the normal state Passos2018 ; Michishita2020 ; Gao2020Photo  by
replacing the vector potential with the crystal momentum ${\bm{k}}$ and by
ignoring the pair potential.

##  Appendix B Phenomenological scattering rate

Here, we introduce the phenomenological treatment of the scattering rate. The
scattering rate $\gamma$ is introduced by modifying the von-Neumann equation
in Eq. (59) to

$i\partial_{t}\hat{\rho}^{(n)}=\sum_{m=0}^{n}\left[\mathcal{H}^{(n-m)}(t),\hat{\rho}^{(m)}(t)\right]-in\gamma\hat{\rho}^{(n)},$
(88)

where $\hat{\rho}^{(n)}$ and $\mathcal{H}^{(n)}$ are $O(F^{n})$ where $F$
represents the strength of the external field. The scattering rate is
therefore multiplied by the perturbation order $n$ of the density matrix. This
description is reasonable to obtain physically meaningful results.
Implementing the phenomenology in Eq. (88), Ref. Passos2018  obtained the
(third-order) nonlinear optical conductivity spectrum consistent with the
gauge symmetry.

Following Eq. (88), we obtain the formula for the nonreciprocal optical
conductivity by replacing the adiabaticity parameter $\eta$ with $\gamma$ in
Eqs. (5)-(9). Note that the prefactor in Eq. (5),
$(i\omega_{1}-\eta)^{-1}(i\omega_{2}-\eta)^{-1}$, does not include the
scattering rate $\gamma$ since it appears by converting the vector potential
into the electric field.

##  Appendix C Linear optical conductivity

We demonstrate that $\bm{\lambda}$-parametrization introduced in Sec. III.1
leads to a physically-transparent expression by investigating the linear
conductivity in Eq. (81). Similarly to Eq. (15), the diamagnetic current
operator obeys the relation

${J}^{\alpha\beta}_{aa}=\braket{a}{\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}H}{a}=-\partial_{\lambda_{\beta}}{J}^{\alpha}_{aa}+i\left[\xi^{\lambda_{\beta}},{J}^{\alpha}\right]_{aa},$
(89)

where $\bm{\lambda}$ dependence is suppressed. The linear conductivity in Eq.
(81) is therefore simplified as

$\displaystyle 2\sigma^{\alpha\beta}(\omega)$
$\displaystyle=\frac{1}{i\omega-\eta}\left[\sum_{a}-\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}E_{a}f_{a}-\sum_{a\neq
b}\left(\frac{-{J}^{\alpha}_{ab}{J}^{\beta}_{ba}f_{ab}}{\omega+i\eta-
E_{ba}}-i\xi^{\lambda_{\beta}}_{ab}{J}^{\alpha}_{ba}f_{ab}\right)\right],$
(90)
$\displaystyle=\frac{1}{i\omega-\eta}\left[\sum_{a}-\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}E_{a}f_{a}+\sum_{a\neq
b}{J}^{\alpha}_{ab}{J}^{\beta}_{ba}f_{ab}\left(\frac{1}{\omega+i\eta-
E_{ba}}-\frac{1}{E_{ab}}\right)\right],$ (91)
$\displaystyle=\frac{1}{i\omega-\eta}\left[\sum_{a}-\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}E_{a}f_{a}-(\omega+i\eta)\sum_{a\neq
b}\frac{{J}^{\alpha}_{ab}{J}^{\beta}_{ba}f_{ab}}{(\omega+i\eta-
E_{ba})E_{ab}}\right].$ (92)

From the obtained expression, we see the clear physical meaning of each
contribution to the linear optical conductivity. For instance, the second term
in the right hand side represents an interlevel transition process and is
rewritten by

$\displaystyle i\sum_{a\neq
b}\frac{{J}^{\alpha}_{ab}{J}^{\beta}_{ba}f_{ab}}{(\omega+i\eta-
E_{ba})E_{ab}}=-i\sum_{a\neq b}\frac{E_{ab}}{\omega+i\eta-
E_{ba}}\left(g^{\lambda_{\alpha}\lambda_{\beta}}_{ab}-\frac{i}{2}\Omega^{\lambda_{\alpha}\lambda_{\beta}}_{ab}\right)f_{ab},$
(93)

where we introduced the band-resolved quantum metric and Berry curvature
Gao2020 ; Watanabe2021

$g^{\lambda_{\alpha}\lambda_{\beta}}_{ab}=\Re{[\xi^{\lambda_{\alpha}}_{ab}\xi^{\lambda_{\beta}}_{ba}]},~{}\Omega^{\lambda_{\alpha}\lambda_{\beta}}_{ab}=-2\Im{[\xi^{\lambda_{\alpha}}_{ab}\xi^{\lambda_{\beta}}_{ba}]}.$
(94)

In this way, this term has a geometric meaning. In particular, taking the
static limit $\omega\rightarrow 0$, the interlevel transition process is given
by the Berry curvature

$\sum_{b(\neq
a)}\Omega^{\lambda_{\alpha}\lambda_{\beta}}_{ab}f_{a}=\epsilon_{\alpha\beta\gamma}\Omega^{\lambda_{\gamma}}_{a}f_{a}.$
(95)

We introduce the Berry curvature parametrized by the variational parameter
$\bm{\lambda}$

$\Omega^{\lambda_{\gamma}}_{a}=\frac{1}{2}\sum_{b(\neq
a)}\epsilon_{\alpha\beta\gamma}\Omega^{\lambda_{\alpha}\lambda_{\beta}}_{ab}.$
(96)

Replacing ${\bm{\lambda}}$ with ${\bm{k}}$, we can see that the obtained
expression reproduces the formula for the anomalous Hall conductivity in the
normal state. Note that the Berry curvature effect does not lead to the
quantized Hall conductivity in the superconducting state, since the
variational parameter has no periodicity in terms of $2\pi/L$ ($L$ is the
system size). In the normal state, the recovered periodicity of
${\bm{k}}({\bm{\lambda}})$ leads to the quantized Hall conductivity in the
gapped phase at the zero temperature Niu1984 . Note also that the Berry
curvature defined in Eq. (96) is different from the Berry curvature
contributing to the thermal Hall effect in superconductors Nomura2012 ;
Sumiyoshi2013 . The latter one (“thermal” Berry curvature) is defined by the
connection Eq. (12), and treats electrons and holes on equal footing, giving
rise to the thermal transport. On the other hand, Eq. (96) (“charge” Berry
curvature) is defined by Eq. (16), and includes electrons and holes with
different signs, contributing to the charge transport as is derived here.

The remaining term in the linear conductivity [Eq. (92)] is the intralevel
transition process denoted by

$\frac{1}{i\omega-\eta}\sum_{a}-\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}E_{a}f_{a}.$
(97)

Recalling the Drude conductivity in the normal state, we can deduce the Drude
part in Eq. (97). Conducting the partial integration, we obtain

$\frac{1}{i\omega-\eta}\sum_{a}\left[-\partial_{\lambda_{\beta}}\left(\partial_{\lambda_{\alpha}}E_{a}f_{a}\right)+\partial_{\lambda_{\alpha}}E_{a}\partial_{\lambda_{\beta}}f_{a}\right].$
(98)

The second term gives the Drude conductivity, whereas the first term should
vanish in the normal state due to the periodicity of the Brillouin zone. On
the other hand, the disappearance of the first term does not hold for the
superconducting state. With the free energy calculated with the Hamiltonian
$H_{\bm{\lambda}}$,

$\displaystyle\sum_{a}\partial_{\lambda_{\alpha}}E_{a}f_{a}$
$\displaystyle=\operatorname{Tr}[(\partial_{\lambda_{\alpha}}H)f(H_{\bm{\lambda}})]$
(99)
$\displaystyle=-\frac{1}{\beta}\partial_{\lambda_{\alpha}}\operatorname{Tr}\ln(1+e^{-\beta
H_{\bm{\lambda}}})$ (100)
$\displaystyle=2\partial_{\lambda_{\alpha}}F_{\bm{\lambda}}.$ (101)

The prefactor 2 comes from the doubling of particles and holes. Accordingly,
we obtain

$\displaystyle\frac{-1}{i\omega-\eta}\sum_{a}\partial_{\lambda_{\beta}}\left(\partial_{\lambda_{\alpha}}E_{a}f_{a}\right)$
$\displaystyle=\frac{-2}{i\omega-\eta}\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}F_{\bm{\lambda}},$
(102) $\displaystyle=\frac{-2}{i\omega-\eta}\rho^{\alpha\beta}_{\text{s}},$
(103)
$\displaystyle=-2\rho^{\alpha\beta}_{\text{s}}\left(\text{P}\frac{1}{i\omega}-\pi\delta(\omega)\right),$
(104)

where the superfluid density is defined by
$\rho^{\alpha\beta}_{\text{s}}=\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}F_{\bm{\lambda}}$.
In the final line, P denotes the principal integral for the frequency
$\Omega$. Finally, taking the limit $\bm{\lambda}\rightarrow\bm{0}$, we obtain
the real part,

$\pi\rho^{\alpha\beta}_{\text{s}}\delta(\omega),$ (105)

which ensures the zero electrical resistivity. As a result, we reproduced the
linear optical conductivity in the superconducting state Tinkham1959SCskin ;
Tinkham2004introduction  based on the formulation using the variational
parameter $\bm{\lambda}$.

##  Appendix D Derivation of formulas for photocurrent response

The ${\bm{\lambda}}$-parametrization works in the formulation of NRO
conductivity in superconductors as in the case of the linear optical
conductivity (Appendix C). For example, we study the photocurrent response
with the condition $\omega_{1}=-\omega_{2}=\Omega$ in Eqs. (5) to (9) and
suppress the frequency dependence in the following.

The entire expression is given by

$\sigma^{\alpha;\beta\gamma}=\sigma^{\alpha;\beta\gamma}_{a}+\sigma^{\alpha;\beta\gamma}_{b}+\sigma^{\alpha;\beta\gamma}_{c}+\sigma^{\alpha;\beta\gamma}_{d},$
(106)

where each contribution is given by

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{a}$
$\displaystyle=\frac{1}{(i\Omega+i\delta-\eta)(-i\Omega+i\delta-\eta)}\sum_{a}\frac{1}{2}{J}^{\alpha\beta\gamma}_{aa}f_{a},$
(107) $\displaystyle 2\sigma^{\alpha;\beta\gamma}_{b}$
$\displaystyle=\frac{1}{(i\Omega+i\delta-\eta)(-i\Omega+i\delta-\eta)}\sum_{a,b}\left(\frac{1}{2}\frac{{J}^{\alpha\beta}_{ab}{J}^{\gamma}_{ba}f_{ab}}{-\Omega+\delta+i\eta-
E_{ba}}+\frac{1}{2}\frac{{J}^{\alpha\gamma}_{ab}{J}^{\beta}_{ba}f_{ab}}{\Omega+\delta+i\eta-
E_{ba}}\right),$ (108) $\displaystyle 2\sigma^{\alpha;\beta\gamma}_{c}$
$\displaystyle=\frac{1}{(i\Omega+i\delta-\eta)(-i\Omega+i\delta-\eta)}\sum_{a,b}\frac{1}{2}\frac{{J}^{\alpha}_{ab}{J}^{\beta\gamma}_{ba}f_{ab}}{2\delta+2i\eta-
E_{ba}},$ (109) $\displaystyle 2\sigma^{\alpha;\beta\gamma}_{d}$
$\displaystyle=\frac{1}{(i\Omega+i\delta-\eta)(-i\Omega+i\delta-\eta)}$
$\displaystyle\times\sum_{a,b,c}\Biggl{[}\frac{1}{2}\frac{{J}^{\alpha}_{ab}}{2\delta+2i\eta-
E_{ba}}\left(\frac{{J}^{\beta}_{bc}{J}^{\gamma}_{ca}f_{ac}}{-\Omega+\delta+i\eta-
E_{ca}}-\frac{{J}^{\beta}_{ca}{J}^{\gamma}_{bc}f_{cb}}{-\Omega+\delta+i\eta-
E_{bc}}\right)$
$\displaystyle~{}~{}~{}~{}+\sum_{a,b,c}\frac{1}{2}\frac{{J}^{\alpha}_{ab}}{2\delta+2i\eta-
E_{ba}}\left(\frac{{J}^{\gamma}_{bc}{J}^{\beta}_{ca}f_{ac}}{\Omega+\delta+i\eta-
E_{ca}}-\frac{{J}^{\gamma}_{ca}{J}^{\beta}_{bc}f_{cb}}{\Omega+\delta+i\eta-
E_{bc}}\right)\Biggr{]}.$ (110)

For a technical reason, we introduce a finite sum-frequency
$2\delta=\omega_{1}+\omega_{2}$ Sipe2000 ; De_Juan2020  and will take the
limit $\delta\rightarrow 0$.

First, we take the component $\sigma_{d}$. When $E_{a}=E_{b}$, the expression
diverges in the limit $\delta\rightarrow 0$ due to $E_{ab}=0$. Thus, we
perform the expansion of what is enclosed in the parenthesis for $\delta$. For
instance, taking the first term in Eq. (110),

$\displaystyle\frac{1}{2}\frac{{J}^{\alpha}_{aa}{J}^{\beta}_{ac}{J}^{\gamma}_{ca}f_{ac}}{2\delta+2i\eta}\left(\frac{1}{-\Omega+\delta+i\eta-
E_{ca}}\right)$
$\displaystyle=\frac{1}{2}\frac{{J}^{\alpha}_{aa}{J}^{\beta}_{ac}{J}^{\gamma}_{ca}f_{ac}}{2\delta+2i\eta}\left(\frac{1}{-\Omega+i\eta-
E_{ca}}-\frac{\delta}{(-\Omega+i\eta-E_{ca})^{2}}+O(\delta^{2})\right),$ (111)
$\displaystyle\rightarrow\frac{1}{2}{J}^{\alpha}_{aa}{J}^{\beta}_{ac}{J}^{\gamma}_{ca}f_{ac}\left(\frac{1}{2i\eta(-\Omega+i\eta-
E_{ca})}-\frac{1/2}{(-\Omega+i\eta-E_{ca})^{2}}\right).$ (112)

In the last line, we take the limit $\delta\rightarrow 0$. Thus,
$O(\delta^{-1})$ and $O(\delta^{0})$ terms are non-zero in Eq. (110). The
careful treatment can be found in the prior studies where the perturbative
calculations are performed on the basis of the Bloch states Sipe2000 ;
De_Juan2020 ; Watanabe2021 . $O(\delta^{0})$ terms are given by

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{\text{intI}}$
$\displaystyle=\frac{1}{\Omega^{2}+\eta^{2}}\sum_{a,c}-\frac{1}{4}\Delta^{\alpha}_{ac}{J}^{\beta}_{ac}{J}^{\gamma}_{ca}f_{ac}\frac{1}{(\Omega-i\eta-
E_{ac})^{2}}+\left[\left(\beta,\gamma,\Omega\right)\leftrightarrow\left(\gamma,\beta,-\Omega\right)\right],$
(113)
$\displaystyle=\frac{1}{\Omega^{2}+\eta^{2}}\sum_{a,c}\frac{1}{4}{J}^{\beta}_{ac}{J}^{\gamma}_{ca}f_{ac}\partial_{\alpha}\frac{1}{\Omega-i\eta-
E_{ac}}+\left[\left(\beta,\gamma,\Omega\right)\leftrightarrow\left(\gamma,\beta,-\Omega\right)\right],$
(114)
$\displaystyle=\frac{1}{\Omega^{2}+\eta^{2}}\sum_{a,c}\frac{1}{2}{J}^{\beta}_{ac}{J}^{\gamma}_{ca}f_{ac}\partial_{\alpha}\text{P}\frac{1}{\Omega-
E_{ac}}.$ (115)

This contribution will be discussed later.

Taking the $O(\delta^{-1})$ term, we obtain

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{\text{inj}}$
$\displaystyle=\frac{1/(4i\eta)}{\Omega^{2}+\eta^{2}}\sum_{a\neq
c}-2i\pi\Delta^{\alpha}_{ac}{J}_{ca}^{\beta}{J}^{\gamma}_{ac}f_{ac}\delta(\Omega-
E_{ac}),$ (116)
$\displaystyle=\frac{1/(4i\eta)}{\Omega^{2}+\eta^{2}}\sum_{a\neq
c}-2i\pi\Delta^{\alpha}_{ac}E_{ca}^{2}\xi_{ca}^{\beta}\xi^{\gamma}_{ac}f_{ac}\delta(\Omega-
E_{ac}),$ (117)
$\displaystyle=\frac{\Omega^{2}/(4i\eta)}{\Omega^{2}+\eta^{2}}\sum_{a\neq
c}-2i\pi\Delta^{\alpha}_{ac}\xi_{ca}^{\beta}\xi^{\gamma}_{ac}f_{ac}\delta(\Omega-
E_{ac}),$ (118) $\displaystyle=\frac{-\pi}{2\eta}\sum_{a\neq
c}\Delta^{\alpha}_{ac}\left(g_{ca}^{\beta\gamma}-\frac{i}{2}\Omega_{ca}^{\beta\gamma}\right)f_{ac}\delta(\Omega-
E_{ac}).$ (119)

We defined the velocity difference matrix
$\Delta^{\alpha}_{ac}={J}^{\alpha}_{aa}-{J}^{\alpha}_{cc}=-\partial_{\alpha}E_{ac}$
and used Eq. (15) and Eq. (94) for the band-resolved quantum metric and Berry
curvature. For a simplified notation, we suppress $\lambda$ for the
connection, derivative, and geometric quantities.

The obtained component $\sigma_{\text{inj}}$ is the so-called injection
current. We can decompose the injection current into two terms; one is
determined by the quantum metric and symmetric under the interchange
$\beta\leftrightarrow\gamma$ (photocurrent induced by linearly-polarized
light, LP-photocurrent), while the other determined by the Berry curvature is
anti-symmetric (photocurrent induced by circularly-polarized light, CP-
photocurrent). Labeling these photocurrents by “Minj” and ”Einj”, we rewrite
the injection current term by

$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{inj}}=\sigma^{\alpha;\beta\gamma}_{\text{Minj}}+\sigma^{\alpha;\beta\gamma}_{\text{Einj}},$
(120) $\displaystyle
2\sigma^{\alpha;\beta\gamma}_{\text{Einj}}=\frac{i\pi}{4\eta}\sum_{a\neq
b}\Delta^{\alpha}_{ab}\Omega_{ba}^{\beta\gamma}f_{ab}\delta(\Omega-E_{ab}),$
(121) $\displaystyle
2\sigma^{\alpha;\beta\gamma}_{\text{Minj}}=\frac{-\pi}{2\eta}\sum_{a\neq
b}\Delta^{\alpha}_{ab}g_{ba}^{\beta\gamma}f_{ab}\delta(\Omega-E_{ab}).$ (122)

Because the relevant geometric quantities have the contrasting symmetry
property, the LP-(CP-)photocurrent is prohibited by the $\mathcal{T}$
($\mathcal{PT}$) symmetry Zhang2019 . Conversely, the LP-(CP-)photocurrent due
to the injection current term is allowed in systems preserving the
$\mathcal{PT}$ ($\mathcal{T}$) symmetry, and it is called magnetic (electric)
injection current. Although the injection current formula shows the diverging
behavior due to $\eta=+0$, the scattering effect suppresses this divergence
DeJuan2016 . The regularized expression is given by

$2\sigma^{\alpha;\beta\gamma}_{\text{inj}}=\frac{-\pi}{2i\gamma}\sum_{a\neq
b}\Delta^{\alpha}_{ab}\left(g_{ba}^{\beta\gamma}-\frac{i}{2}\Omega_{ba}^{\beta\gamma}\right)f_{ab}\mathcal{L}\left(\Omega-
E_{ab}\right).$ (123)

We defined the Lorentzian function parametrized by the scattering rate
$\gamma$,

$\mathcal{L}(x)=\frac{\gamma}{\pi}\frac{1}{\gamma^{2}+x^{2}}.$ (124)

Finally, we obtain the two components from Eq. (110) with the intralevel
transition condition $E_{a}=E_{b}$ as

$\sigma^{\alpha;\beta\gamma}_{d,{\rm
intra}}=\sigma^{\alpha;\beta\gamma}_{\text{inj}}+\sigma^{\alpha;\beta\gamma}_{\text{intI}}.$
(125)

Next, we consider the interlevel transition process ($E_{ab}\neq 0$) for
$\sigma_{d}$. In this case, we safely take the limit $\delta\rightarrow 0$.
The expression is given by

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{d,{\rm inter}}$
$\displaystyle=\frac{1}{\Omega^{2}+\eta^{2}}\Biggl{[}\sum_{a,b}^{\prime}\sum_{c}\frac{1}{2}\frac{{J}^{\alpha}_{ab}}{2i\eta-
E_{ba}}\left(\frac{{J}^{\gamma}_{bc}{J}^{\beta}_{ca}f_{ac}}{\Omega+i\eta-
E_{ca}}-\frac{{J}^{\gamma}_{ca}{J}^{\beta}_{bc}f_{cb}}{\Omega+i\eta-
E_{bc}}\right)\Biggr{]}$
$\displaystyle+\left[\left(\beta,\gamma,\Omega\right)\leftrightarrow\left(\gamma,\beta,-\Omega\right)\right].$
(126)

Here, $\Sigma^{\prime}$ indicates the summation of $a,b$ with the condition
$E_{ab}\neq 0$. Making use of Eq. (15), we rewrite the current operator by the
Berry connection,

$\frac{{J}^{\alpha}_{ab}}{E_{ab}}{J}^{\gamma}_{bc}=-i\xi^{\alpha}_{ab}{J}^{\gamma}_{bc}=\Braket{a}{\partial_{\alpha}b}\Braket{b}{{J}^{\gamma}}{c}=-\Braket{\partial_{\alpha}a}{b}\Braket{b}{{J}^{\gamma}}{c}=-\Braket{D_{\alpha}a}{b}\Braket{b}{{J}^{\gamma}}{c},$
(127)

where we use the covariant derivative

$\ket{D_{\alpha}a}=\ket{\partial_{\alpha}a}+i\sum_{d}^{E_{ad}=0}\xi^{\alpha}_{da}\ket{d}.$
(128)

The gauge degree of freedom associated with the covariant derivative is
defined by the eigenstates having the same energy. The degeneracy does not
exist in some cases, and hence the covariant derivative is characterized by
the $U(1)$ gauge Von_Baltz1981 . On the other hand, for the $\mathcal{PT}$
symmetric and spinful systems, the $\mathcal{PT}$-ensured Kramers degeneracy
leads to the $U(2)$ covariant derivative Watanabe2021  even in the presence of
the variational parameter $\bm{\lambda}$. With the covariant derivative, Eq.
(126) is transformed as

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{d,{\rm inter}}$
$\displaystyle=-\frac{1}{2\left(\Omega^{2}+\eta^{2}\right)}\Biggl{[}\sum_{a,c}\frac{\Braket{D_{\alpha}a}{{J}^{\gamma}}{c}{J}^{\beta}_{ca}f_{ac}}{\Omega+i\eta-
E_{ca}}+\sum_{b,c}\frac{\Braket{c}{{J}^{\gamma}}{D_{\beta}b}{J}^{\beta}_{bc}f_{cb}}{\Omega+i\eta-
E_{bc}}\Biggr{]}$
$\displaystyle+\left[\left(\beta,\gamma,\Omega\right)\leftrightarrow\left(\gamma,\beta,-\Omega\right)\right],$
(129)
$\displaystyle=-\frac{1}{2\left(\Omega^{2}+\eta^{2}\right)}\sum_{a,b}\Biggl{[}\Braket{D_{\alpha}a}{{J}^{\gamma}}{b}{J}^{\beta}_{ba}f_{ab}+\Braket{a}{{J}^{\gamma}}{D_{\alpha}b}{J}^{\beta}_{ba}f_{ab}\Biggr{]}$
$\displaystyle~{}~{}~{}~{}~{}~{}\times\left(\text{P}\frac{1}{\Omega-
E_{ba}}-i\pi\delta(\Omega-
E_{ba})\right)+\left[\left(\beta,\gamma,\Omega\right)\leftrightarrow\left(\gamma,\beta,-\Omega\right)\right].$
(130)

Introducing the Berry connection $\mathcal{A}^{\alpha}$ associated with the
covariant derivative $D_{\alpha}$, we symbolically define the covariant
derivative of operators by

$[D_{\alpha}O]_{ab}=\partial_{\alpha}O_{ab}-i\left[\mathcal{A}^{\alpha},O\right]_{ab}.$
(131)

Accordingly, we obtain

$\Braket{D_{\alpha}a}{{J}^{\gamma}}{b}+\Braket{a}{{J}^{\gamma}}{D_{\alpha}b}=[D_{\alpha}{J}^{\gamma}]_{ab}-(\partial_{\alpha}{J}^{\gamma})_{ab}=[D_{\alpha}{J}^{\gamma}]_{ab}+{J}^{\alpha\gamma}_{ab}.$
(132)

Meanwhile, $\sigma_{b}$ in Eq. (108) reads with the limit $\delta\rightarrow
0$,

$2\sigma_{b}=\frac{1}{2\left(\Omega^{2}+\eta^{2}\right)}\sum_{a,b}{J}^{\alpha\gamma}_{ab}{J}^{\beta}_{ba}f_{ab}\left[\text{P}\frac{1}{\Omega-
E_{ba}}-i\pi\delta(\Omega-
E_{ba})\right]+\left[\left(\beta,\gamma,\Omega\right)\leftrightarrow\left(\gamma,\beta,-\Omega\right)\right],$
(133)

Thus, by summing this term and Eq. (130), the component including
${J}^{\alpha\gamma}$ in Eq. (130) is canceled. The remaining term is given by

$\sigma^{\alpha;\beta\gamma}_{d,inter}+\sigma^{\alpha;\beta\gamma}_{b}=\sigma^{\alpha;\beta\gamma}_{\text{abs}}+\sigma^{\alpha;\beta\gamma}_{\text{rea}},$
(134)

where

$2\sigma^{\alpha;\beta\gamma}_{\text{abs}}=-\frac{1}{2\left(\Omega^{2}+\eta^{2}\right)}\sum_{a,b}-i\pi\Biggl{[}[D_{\alpha}{J}^{\gamma}]_{ab}{J}^{\beta}_{ba}-[D_{\alpha}{J}^{\beta}]_{ba}{J}^{\gamma}_{ab}\Biggr{]}f_{ab}\delta(\Omega-
E_{ba}),$ (135)

and

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{\text{rea}}$
$\displaystyle=-\frac{1}{2\left(\Omega^{2}+\eta^{2}\right)}\sum_{a,b}\Biggl{[}[D_{\alpha}{J}^{\gamma}]_{ab}{J}^{\beta}_{ba}+[D_{\alpha}{J}^{\beta}]_{ba}{J}^{\gamma}_{ab}\Biggr{]}f_{ab}\text{P}\frac{1}{\Omega-
E_{ba}}.$ (136)

$\sigma^{\alpha;\beta\gamma}_{\text{abs}}$ characterized by the resonant
contribution is distinguished from the injection current term [Eq. (123)]
because it is insensitive to the scattering rate. This term is divided into
the LP and CP-photocurrents as

$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{abs}}=\sigma^{\alpha;\beta\gamma}_{\text{shift}}+\sigma^{\alpha;\beta\gamma}_{\text{gyro}},$
(137) $\displaystyle
2\sigma^{\alpha;\beta\gamma}_{\text{shift}}=-\frac{\pi}{2\Omega^{2}}\sum_{a\neq
b}\Im{\Biggl{[}[D_{\alpha}{J}^{\beta}]_{ab}{J}^{\gamma}_{ba}+[D_{\alpha}{J}^{\gamma}]_{ab}{J}^{\beta}_{ba}\Biggr{]}}f_{ab}\delta(\Omega-
E_{ba}),$ (138) $\displaystyle
2\sigma^{\alpha;\beta\gamma}_{\text{gyro}}=-\frac{i\pi}{2\Omega^{2}}\sum_{a\neq
b}\Re{\Biggl{[}[D_{\alpha}{J}^{\beta}]_{ab}{J}^{\gamma}_{ba}-[D_{\alpha}{J}^{\gamma}]_{ab}{J}^{\beta}_{ba}\Biggr{]}}f_{ab}\delta(\Omega-
E_{ba}),$ (139)

where we use the identity

$[D_{\alpha}{J}^{\gamma}]_{ab}{J}^{\beta}_{ba}-[D_{\alpha}{J}^{\beta}]_{ba}{J}^{\gamma}_{ab}=i\Im\,\left([D_{\alpha}{J}^{\beta}]_{ab}{J}^{\gamma}_{ba}+[D_{\alpha}{J}^{\gamma}]_{ab}{J}^{\beta}_{ba}\right)-\Re\,\left([D_{\alpha}{J}^{\beta}]_{ab}{J}^{\gamma}_{ba}-[D_{\alpha}{J}^{\gamma}]_{ab}{J}^{\beta}_{ba}\right).$
(140)

The resulting terms $\sigma_{\text{shift}}$ and $\sigma_{\text{gyro}}$ are
called shift current Von_Baltz1981 ; Sipe2000 ; Sturman1992Book  and gyration
current Watanabe2021 ; Ahn2020 . Owing to the Hellmann-Feynman relation, the
velocity operators are rewritten by the Berry connections as

$\displaystyle
2\sigma^{\alpha;\beta\gamma}_{\text{shift}}=-\frac{\pi}{2}\sum_{a\neq
b}\Im{\Biggl{[}[D_{\alpha}\xi^{\beta}]_{ab}\xi^{\gamma}_{ba}+[D_{\alpha}\xi^{\gamma}]_{ab}\xi^{\beta}_{ba}\Biggr{]}}f_{ab}\delta(\Omega-
E_{ba}),$ (141) $\displaystyle
2\sigma^{\alpha;\beta\gamma}_{\text{gyro}}=-\frac{i\pi}{2}\sum_{a\neq
b}\Re{\Biggl{[}[D_{\alpha}\xi^{\beta}]_{ab}\xi^{\gamma}_{ba}-[D_{\alpha}\xi^{\gamma}]_{ab}\xi^{\beta}_{ba}\Biggr{]}}f_{ab}\delta(\Omega-
E_{ba}).$ (142)

Following the parallel discussion in Refs. Von_Baltz1981 ; Watanabe2021 , we
can show that these photocurrent responses are written by the shift vector and
the chiral shift vector, respectively. For the numerical calculations, we
approximate the delta function into the Lorentzian function as in Eq. (123).

The off-resonant term $\sigma_{\text{rea}}$ is added to Eq. (115) and gives

$\displaystyle
2\sigma^{\alpha;\beta\gamma}_{\text{intI+rea}}=2\sigma^{\alpha;\beta\gamma}_{\text{intI}}+2\sigma^{\alpha;\beta\gamma}_{\text{rea}},$
(143)
$\displaystyle=\frac{1}{2\left(\Omega^{2}+\eta^{2}\right)}\left[\sum_{a\neq
b}\partial_{\alpha}\left({J}^{\beta}_{ab}{J}^{\gamma}_{ba}\text{P}\frac{1}{\Omega-
E_{ab}}\right)f_{ab}+f_{ab}\text{P}\frac{1}{\Omega-
E_{ba}}\left(i[\mathcal{A}^{\alpha},{J}^{\beta}]_{ba}{J}^{\gamma}_{ab}+i{J}^{\beta}_{ba}[\mathcal{A}^{\alpha},{J}^{\gamma}]_{ab}\right)\right].$
(144)

Owing to the definition of the Berry connection $\mathcal{A}^{\alpha}$, the
second term vanishes in the last line. Performing the partial integration of
the first term, we obtain

$2\sigma^{\alpha;\beta\gamma}_{\text{intI+rea}}=\frac{1}{2\left(\Omega^{2}+\eta^{2}\right)}\left[-\sum_{a\neq
b}{J}^{\beta}_{ab}{J}^{\gamma}_{ba}\text{P}\frac{1}{\Omega-
E_{ab}}\partial_{\alpha}f_{ab}+\partial_{\alpha}\left(\sum_{a\neq
b}{J}^{\beta}_{ab}{J}^{\gamma}_{ba}\text{P}\frac{1}{\Omega-
E_{ab}}f_{ab}\right)\right].$ (145)

We should note that the second term enclosed by the derivative of
$\lambda_{\alpha}$ does not vanish in the BdG Hamiltonian while it is dropped
due to the periodicitity in the crystal momentum space in the normal state. We
label the first and second terms in Eq. (145) by $\sigma_{\text{FS1}}$ and
$\sigma_{\text{SC1}}$,

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{\text{FS1}}$
$\displaystyle=-\frac{1}{2\left(\Omega^{2}+\eta^{2}\right)}\sum_{a\neq
b}{J}^{\beta}_{ab}{J}^{\gamma}_{ba}\text{P}\frac{1}{\Omega-
E_{ab}}\partial_{\alpha}f_{ab},$ (146) $\displaystyle
2\sigma^{\alpha;\beta\gamma}_{\text{SC1}}$
$\displaystyle=\frac{1}{2\left(\Omega^{2}+\eta^{2}\right)}\partial_{\alpha}\left(\sum_{a\neq
b}{J}^{\beta}_{ab}{J}^{\gamma}_{ba}\text{P}\frac{1}{\Omega-
E_{ab}}f_{ab}\right).$ (147)

Next, we consider $\sigma_{a}$ [Eq. (107)]. Taking $\delta\rightarrow 0$ and
using the relation

${J}^{\alpha\beta\gamma}_{aa}=-\partial_{\alpha}{J}^{\beta\gamma}_{aa}+i[\xi^{\alpha},{J}^{\beta\gamma}]_{aa},$
(148)

we obtain

$2\sigma^{\alpha;\beta\gamma}_{a}=\frac{1}{2(i\Omega-\eta)(-i\Omega-\eta)}\left[\sum_{a}-\partial_{\alpha}{J}^{\beta\gamma}_{aa}f_{a}+\sum_{a,b}\frac{{J}^{\alpha}_{ab}{J}^{\beta\gamma}_{ba}f_{ba}}{E_{ab}}\right].$
(149)

The second component in the right hand side is canceled out by $\sigma_{c}$
[Eq. (109)]. The first component is further transformed by Eq. (89) into

$\displaystyle 2\partial_{\alpha}{J}^{\beta\gamma}_{aa}$
$\displaystyle=-\partial_{\alpha}\partial_{\beta}{J}^{\gamma}_{aa}+i\partial_{\alpha}\left([\xi^{\beta},{J}^{\gamma}]_{aa}\right)-\partial_{\alpha}\partial_{\gamma}{J}^{\beta}_{aa}+i\partial_{\alpha}\left([\xi^{\gamma},{J}^{\beta}]_{aa}\right),$
(150)
$\displaystyle=2\partial_{\alpha}\partial_{\beta}\partial_{\gamma}\epsilon_{a}+i\partial_{\alpha}\left([\xi^{\beta},{J}^{\gamma}]_{aa}\right)+i\partial_{\alpha}\left([\xi^{\gamma},{J}^{\beta}]_{aa}\right),$
(151)

and we accordingly obtain

$2\sigma^{\alpha;\beta\gamma}_{a}+2\sigma^{\alpha;\beta\gamma}_{c}=-\frac{1}{2\Omega^{2}}\sum_{a}f_{a}\partial_{\alpha}\partial_{\beta}\partial_{\gamma}\epsilon_{a}+\frac{1}{\Omega^{2}}\left[-\sum_{a\neq
b}\frac{1}{2}\frac{{J}^{\beta}_{ab}{J}^{\gamma}_{ba}}{E_{ab}}\partial_{\alpha}f_{ab}+\partial_{\alpha}\left(\sum_{a\neq
b}\frac{1}{2}\frac{{J}^{\beta}_{ab}{J}^{\gamma}_{ba}}{E_{ab}}f_{ab}\right)\right].$
(152)

To see the correspondence between the photocurrent responses in the normal
state and those in superconducting state, we perform the partial integration
for the first term,

$\frac{1}{2\Omega^{2}}\left[\partial_{\beta}\partial_{\gamma}\left(\sum_{a}f_{a}\partial_{\alpha}\epsilon_{a}\right)-\partial_{\gamma}\left(\sum_{a}\partial_{\beta}f_{a}\partial_{\alpha}\epsilon_{a}\right)-\partial_{\beta}\left(\sum_{a}\partial_{\gamma}f_{a}\partial_{\alpha}\epsilon_{a}\right)+\sum_{a}\partial_{\beta}\partial_{\gamma}f_{a}\partial_{\alpha}\epsilon_{a}\right].$
(153)

What is not differentiated by $\lambda$ is the (nonlinear) Drude term
Holder2020  denoted by

$2\sigma^{\alpha;\beta\gamma}_{\text{D}}=-\frac{1}{2\Omega^{2}}\sum_{a}\partial_{\alpha}\epsilon_{a}\partial_{\beta}\partial_{\gamma}f_{a}.$
(154)

This term can be finite in the normal state. The remaining term is finite only
in the superconducting state and labeled by

$2\sigma^{\alpha;\beta\gamma}_{\text{SC2}}=-\frac{1}{2\Omega^{2}}\left[\partial_{\beta}\partial_{\gamma}\left(\sum_{a}f_{a}\partial_{\alpha}\epsilon_{a}\right)-\partial_{\gamma}\left(\sum_{a}\partial_{\beta}f_{a}\partial_{\alpha}\epsilon_{a}\right)-\partial_{\beta}\left(\sum_{a}\partial_{\gamma}f_{a}\partial_{\alpha}\epsilon_{a}\right)\right].$
(155)

We also denote the $\lambda$-derivative term in Eq. (152) as

$2\sigma^{\alpha;\beta\gamma}_{\text{SC3}}=\frac{1}{2\Omega^{2}}\partial_{\alpha}\left(\sum_{a\neq
b}\frac{{J}^{\beta}_{ab}{J}^{\gamma}_{ba}}{E_{ab}}f_{ab}\right),$ (156)

which is also finite only in the superconducting state. Accordingly, we
decompose Eq. (152) as

$\sigma^{\alpha;\beta\gamma}_{a}+\sigma^{\alpha;\beta\gamma}_{c}=\sigma^{\alpha;\beta\gamma}_{\text{D}}+\sigma^{\alpha;\beta\gamma}_{\text{SC2}}+\sigma^{\alpha;\beta\gamma}_{\text{SC3}}+\sigma^{\alpha;\beta\gamma}_{\text{FS2}},$
(157)

where we defined

$2\sigma^{\alpha;\beta\gamma}_{\text{FS2}}=-\frac{1}{2\Omega^{2}}\sum_{a\neq
b}\frac{{J}^{\beta}_{ab}{J}^{\gamma}_{ba}}{E_{ab}}\partial_{\alpha}f_{ab}.$
(158)

Combining $\sigma_{\text{FS1}}$ and $\sigma_{\text{FS2}}$, we obtain a Fermi
surface contribution

$\displaystyle 2\sigma_{\text{FS1}}+2\sigma_{\text{FS2}}=$
$\displaystyle-\frac{1}{\Omega}\sum_{a\neq
b}\frac{1}{2}\frac{{J}^{\beta}_{ab}{J}^{\gamma}_{ba}}{E_{ab}\left(\Omega-
E_{ab}\right)}\partial_{\alpha}f_{ab},$ (159) $\displaystyle=-\sum_{a\neq
b}\frac{1}{2}\frac{{J}^{\beta}_{ab}{J}^{\gamma}_{ba}}{E_{ab}^{2}}\left(\frac{1}{\Omega-
E_{ab}}-\frac{1}{\Omega}\right)\partial_{\alpha}f_{ab}.$ (160)

In the final line, the first and second components are called intrinsic Fermi
surface term De_Juan2020 ; Watanabe2021  and Berry curvature dipole effect
Moore2010 ; Sodemann2015 . We denote them as

$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{IFS}}=-\frac{1}{2}\sum_{a\neq
b}\xi^{\beta}_{ab}\xi^{\gamma}_{ba}\frac{1}{\Omega-
E_{ab}}\partial_{\alpha}f_{ab},$ (161)
$\displaystyle\sigma^{\alpha;\beta\gamma}_{\text{BCD}}=\frac{1}{2\Omega}\sum_{a\neq
b}\xi^{\beta}_{ab}\xi^{\gamma}_{ba}\partial_{\alpha}f_{ab}=-\frac{i}{2\Omega}\sum_{a}\Omega^{\beta\gamma}_{a}\partial_{\alpha}f_{a}.$
(162)

The Berry curvature dipole effect is classified as the CP-photocurrent since
it is anti-symmetric under $\beta\leftrightarrow\gamma$. On the other hand,
the intrinsic Fermi surface effect can be decomposed into the LP and CP-
photocurrents as in the injection current [Eq. (120)]. The LP and CP parts are
called electric and magnetic intrinsic Fermi surface effects, respectively
De_Juan2020 ; Watanabe2021 .

Finally, we obtain the normal photocurrent conductivity

$\sigma_{\text{n}}^{\alpha;\beta\gamma}=\sigma_{\text{Einj}}^{\alpha;\beta\gamma}+\sigma_{\text{Minj}}^{\alpha;\beta\gamma}+\sigma_{\text{shift}}^{\alpha;\beta\gamma}+\sigma_{\text{gyro}}^{\alpha;\beta\gamma}+\sigma_{\text{D}}^{\alpha;\beta\gamma}+\sigma_{\text{IFS}}^{\alpha;\beta\gamma}+\sigma_{\text{BCD}}^{\alpha;\beta\gamma}.$
(163)

Note that we can easily reproduce the photocurrent response in the normal
state by replacing $\bm{\lambda}$ with ${\bm{k}}$ in the normal photocurrent.
The nonlinear Drude term, intrinsic Fermi surface term, and Berry curvature
dipole effect include the derivative of the Fermi distribution function and
thus vanish in a gapped system at low temperature.

Summing the contributions enclosed by the ${\bm{\lambda}}$-differentiation, we
obtain the anomalous photocurrent arising from the superconductivity

$\displaystyle 2\sigma_{\text{a}}^{\alpha;\beta\gamma}$
$\displaystyle=2\sigma_{\text{SC1}}^{\alpha;\beta\gamma}+2\sigma_{\text{SC2}}^{\alpha;\beta\gamma}+2\sigma_{\text{SC3}}^{\alpha;\beta\gamma},$
(164)
$\displaystyle=\frac{1}{2\Omega^{2}}\partial_{\alpha}\left(\sum_{a,b}{J}^{\beta}_{ab}{J}^{\gamma}_{ba}\text{P}\frac{1}{\Omega-
E_{ab}}f_{ab}\right)+\frac{1}{\Omega^{2}}\partial_{\alpha}\left(\sum_{a\neq
b}\frac{1}{2}\frac{{J}^{\beta}_{ab}{J}^{\gamma}_{ba}}{E_{ab}}f_{ab}\right)$
$\displaystyle-\frac{1}{2\Omega^{2}}\left[\partial_{\beta}\partial_{\gamma}\left(\sum_{a}f(\epsilon_{a})\partial_{\alpha}\epsilon_{a}\right)-\partial_{\beta}\left(\sum_{a}\partial_{\gamma}f(\epsilon_{a})\partial_{\alpha}\epsilon_{a}\right)-\partial_{\gamma}\left(\sum_{a}\partial_{\beta}f(\epsilon_{a})\partial_{\alpha}\epsilon_{a}\right)\right],$
(165) $\displaystyle=\frac{1}{2\Omega^{2}}\partial_{\alpha}\left[\sum_{a\neq
b}{J}^{\beta}_{ab}{J}^{\gamma}_{ba}\left(\text{P}\frac{1}{\Omega-
E_{ab}}+\frac{1}{E_{ab}}\right)f_{ab}\right]-\frac{1}{\Omega^{2}}\partial_{\alpha}\partial_{\beta}\partial_{\gamma}F_{\bm{\lambda}}$
$\displaystyle+\frac{1}{2\Omega^{2}}\left[\partial_{\beta}\left(\sum_{a}\partial_{\gamma}f(\epsilon_{a})\partial_{\alpha}\epsilon_{a}\right)+\partial_{\gamma}\left(\sum_{a}\partial_{\beta}f(\epsilon_{a})\partial_{\alpha}\epsilon_{a}\right)\right].$
(166)

This is the new contribution to the photocurrent clarified by this work.
Assuming the gapped system at low temperature, we here ignore the Fermi
surface contribution in the second line of Eq. (166). We denote the remaining
new contributions as

$\displaystyle
2\sigma_{\text{NRSF}}^{\alpha;\beta\gamma}=-\frac{1}{\Omega^{2}}\partial_{\alpha}\partial_{\beta}\partial_{\gamma}F_{\bm{\lambda}},$
(167) $\displaystyle
2\sigma_{\text{CD}}^{\alpha;\beta\gamma}=\frac{1}{2\Omega^{2}}\partial_{\alpha}\left[\sum_{a\neq
b}{J}^{\beta}_{ab}{J}^{\gamma}_{ba}\left(\text{P}\frac{1}{\Omega-
E_{ab}}+\frac{1}{E_{ab}}\right)f_{ab}\right],$ (168)

where “NRSF“ denotes nonreciprocal superfluid density while “CD” means
conductivity derivative. The nonreciprocal superfluid density is defined as
$f^{\alpha\beta\gamma}=\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\alpha}\partial_{\beta}\partial_{\gamma}F_{\bm{\lambda}}$.

It is convenient for numerical computations to rewrite the formulas by
generalized velocity operators. A straightforward example is the injection
current contributions in Eq. (119). Although Eq. (119) implicitly includes the
infinitesimal variational parameter $\bm{\lambda}$, we can obtain the formula
without $\bm{\lambda}$ by taking the limit $\bm{\lambda}\rightarrow\bm{0}$.
The resulting formula is written by the generalized velocity operators and
energy spectrum of the BdG Hamiltonian. The quantum metric and Berry curvature
are given by

$\displaystyle
g^{\alpha\beta}_{ab}=\frac{1}{E_{ab}^{2}}\Re{[{J}^{\alpha}_{ab}{J}^{\beta}_{ba}]},$
(169)
$\displaystyle\Omega^{\alpha\beta}_{ab}=-\frac{2}{E_{ab}^{2}}\Im{[{J}^{\alpha}_{ab}{J}^{\beta}_{ba}]},$
(170)

where we use Eq. (11) with $a\neq b$. Since the BdG Hamiltonian can be labeled
by the crystal momentum ${\bm{k}}$, the injection current formula is obtained
as

$2\sigma^{\alpha;\beta\gamma}_{\text{inj}}=\frac{-\pi}{2\eta}\sum_{\bm{k}}\sum_{a,b}\left[{J}^{\alpha}_{aa}({\bm{k}})-{J}^{\alpha}_{bb}({\bm{k}})\right]\frac{{J}^{\beta}_{ba}({\bm{k}}){J}^{\gamma}_{ab}({\bm{k}})}{E_{ba}^{2}}f_{ab}\delta(\Omega-
E_{ab}),$ (171)

where $a,b$ denote the band indices of Bogoliubov quasiparticles. Following
the parallel discussion, we can derive the numerically convenient formulas for
other contributions by transforming the ${\bm{\lambda}}$-parametrized formulas
into what is given by the generalized velocity operators and energy spectrum.
In the case of the shift current and gyration current contributions [Eqs.
(141) and (142)], the covariant derivative of the Berry connection is
rewritten by

$\left[D_{\alpha}\xi^{\beta}\right]_{ab}=i{J}^{\beta}_{ab}\frac{\Delta_{ab}^{\alpha}}{E_{ab}^{2}}+\frac{i}{E_{ab}}\left(-{J}_{ab}^{\alpha\beta}-\sum_{a\neq
c}\frac{{J}^{\alpha}_{ac}{J}^{\beta}_{cb}}{E_{ac}}+\sum_{b\neq
c}\frac{{J}^{\alpha}_{cb}{J}^{\beta}_{ac}}{E_{cb}}\right).$ (172)

After taking the limit ${\bm{\lambda}}\rightarrow\bm{0}$, the right hand side
becomes a numerically convenient form.

Next, we derive the anomalous photocurrent formulas in the current operator
representation. For instance, we take the nonreciprocal superfluid density
effect [Eq. (167)]. When we assume the gapped and low-temperature state for
simplicity, we can neglect the term including the derivative of the Fermi
distribution function. The obtained formula for the nonreciprocal superfluid
density is

$\displaystyle-2\partial_{\alpha}\partial_{\beta}\partial_{\gamma}F_{\bm{\lambda}}$
$\displaystyle=\sum_{a}{J}^{\alpha\beta\gamma}_{aa}f_{a}+\sum_{a\neq
b}\frac{1}{E_{ab}}\left({J}^{\alpha}_{ab}{J}^{\beta\gamma}_{ba}+{J}^{\beta}_{ab}{J}^{\gamma\alpha}_{ba}+{J}^{\gamma}_{ab}{J}^{\alpha\beta}_{ba}+c.c.\right)f_{a}$
$\displaystyle-\sum_{a\neq
b}\frac{\Delta^{\gamma}_{ab}}{E_{ab}^{2}}\left({J}^{\beta}_{ab}{J}^{\alpha}_{ba}+{J}^{\beta}_{ba}{J}^{\alpha}_{ab}\right)f_{a}$
$\displaystyle+\sum_{a\neq
b,c}\frac{1}{E_{ab}E_{ac}}f_{a}\left[{J}^{\alpha}_{ab}{J}^{\beta}_{bc}{J}^{\gamma}_{ca}+{J}^{\beta}_{ab}{J}^{\alpha}_{bc}{J}^{\gamma}_{ca}+c.c\right]$
$\displaystyle+\sum_{b\neq
a,c}\frac{1}{E_{ab}E_{bc}}f_{a}\left[{J}^{\beta}_{ab}{J}^{\gamma}_{bc}{J}^{\alpha}_{ca}+{J}^{\alpha}_{ab}{J}^{\gamma}_{bc}{J}^{\beta}_{ca}+c.c\right].$
(173)

Similarly, we can rewrite the formula for the conductivity derivative effect
[Eq. (168)] as

$\displaystyle\partial_{\alpha}\left[{J}^{\beta}_{ab}{J}^{\gamma}_{ba}\left(\text{P}\frac{1}{\Omega-
E_{ab}}+\frac{1}{E_{ab}}\right)f_{ab}\right]$ (174)
$\displaystyle=-{J}^{\beta}_{ab}{J}^{\gamma}_{ba}\Delta^{\alpha}_{ab}\left(\text{P}\frac{1}{\left(\Omega-
E_{ab}\right)^{2}}-\frac{1}{E_{ab}^{2}}\right)f_{ab}$ $\displaystyle F$
$\displaystyle+\left(-{J}^{\alpha\beta}_{ab}{J}^{\gamma}_{ba}-{J}^{\beta}_{ab}{J}^{\alpha\gamma}_{ba}\right)\left(\text{P}\frac{1}{\Omega-
E_{ab}}+\frac{1}{E_{ab}}\right)f_{ab}$ $\displaystyle+\left[\sum_{c\neq
a}\left(-\frac{{J}^{\alpha}_{ac}{J}^{\beta}_{cb}{J}^{\gamma}_{ba}}{E_{ac}}+\frac{{J}^{\beta}_{ab}{J}^{\alpha}_{ca}{J}^{\gamma}_{bc}}{E_{ca}}\right)+\sum_{c\neq
b}\left(\frac{{J}^{\alpha}_{cb}{J}^{\beta}_{ac}{J}^{\gamma}_{ba}}{E_{cb}}-\frac{{J}^{\beta}_{ab}{J}^{\alpha}_{bc}{J}^{\gamma}_{ca}}{E_{bc}}\right)\right]\left(\text{P}\frac{1}{\Omega-
E_{ab}}+\frac{1}{E_{ab}}\right)f_{ab},$ (175)

where we omit the prefactor $1/(2\Omega^{2})$ and neglect the derivative of
the Fermi distribution function.

Note that there is an arbitrariness in decomposing the total photocurrent into
the normal and anomalous contributions within the independent particle
approximation. For instance, although the Berry curvature dipole effect in Eq.
(162) is seemingly different from what has been derived in prior works
Moore2010 ; Deyo2009 ; Sodemann2015 , we can reproduce the known expression.
Since the Berry curvature is given by the derivative of the intraband Berry
connection
($\Omega_{a}^{\alpha\beta}=\partial_{\alpha}\xi^{\beta}_{aa}-\partial_{\beta}\xi^{\alpha}_{aa}$),
Eq. (162) is rewritten by

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{\text{BCD}}$
$\displaystyle=-i\frac{1}{2\Omega}\left(\sum_{a}\partial_{\alpha}\left(\Omega^{\beta\gamma}_{a}f_{a}\right)-\partial_{\alpha}\Omega^{\beta\gamma}_{a}f_{a}\right),$
(176)
$\displaystyle=-i\frac{1}{2\Omega}\sum_{a}\left(\partial_{\alpha}\left(\Omega^{\beta\gamma}_{a}f_{a}\right)-\left(\partial_{\beta}\Omega_{a}^{\alpha\gamma}-\partial_{\gamma}\Omega_{a}^{\alpha\beta}\right)f_{a}\right),$
(177)
$\displaystyle=-i\frac{1}{2\Omega}\sum_{a}\left[\partial_{\alpha}\left(\Omega^{\beta\gamma}_{a}f_{a}\right)-\partial_{\beta}\left(\Omega_{a}^{\alpha\gamma}f_{a}\right)+\partial_{\gamma}\left(\Omega_{a}^{\alpha\beta}f_{a}\right)+\left(\Omega_{a}^{\alpha\gamma}\partial_{\beta}f_{a}-\Omega_{a}^{\alpha\beta}\partial_{\gamma}f_{a}\right)\right],$
(178)

where we obtain the “conventional” Berry curvature dipole contribution

$2\sigma^{\alpha;\beta\gamma}_{\text{cBCD}}=i\frac{1}{2\Omega}\sum_{a}\left(\Omega_{a}^{\alpha\beta}\partial_{\gamma}f_{a}-\Omega_{a}^{\alpha\gamma}\partial_{\beta}f_{a}\right),$
(179)

The remaining contributions are enclosed by the $\lambda$-derivative,

$-i\frac{1}{2\Omega}\sum_{a}\left[\partial_{\alpha}\left(\Omega^{\beta\gamma}_{a}f_{a}\right)-\partial_{\beta}\left(\Omega_{a}^{\alpha\gamma}f_{a}\right)+\partial_{\gamma}\left(\Omega_{a}^{\alpha\beta}f_{a}\right)\right],$
(180)

which should vanish in the normal state. Accordingly, the two formulas for the
Berry curvature dipole effect [Eqs. (162) and (179)] are equivalent in the
normal state. On the other hand, when we adopt the formula in Eq. (179) as the
Berry curvature dipole effect, it is necessary to include Eq. (180) in the
anomalous contribution. Such arbitrariness in the normal and anomalous
photocurrents may occur because we formulate the response function without any
scattering effect. This difficulty is resolved by the Green function method in
Appendix E.

Some of the normal photocurrents diverge in the static limit, for instance,
the Drude photocurrent in Eq. (154). On the other hand, when the scattering
effect is properly introduced, the normal photocurrent is bounded to be finite
even in the zero-frequency limit ($\Omega\rightarrow 0$) Michishita2020 . On
the contrary, the anomalous photocurrent such as the nonreciprocal superfluid
density effect should diverge in the low-frequency limit even with scattering
effects since the divergence arises from rewriting the vector potential by the
electric field as $\bm{E}=i\Omega\bm{A}$.

##  Appendix E Derivation by Green function method

We formulate the nonreciprocal optical response in superconductors based on
the Green function method. The Green function formalism for the nonlinear
conductivity has been reported in prior works Michishita2020 ; Holder2020 ;
Parker2019 , and we apply it to the BdG Hamiltonian. Starting from the
Hamiltonian having the electron correlation, we take the molecular field
approximation for the Cooper channel of the electron-electron interaction.
Accordingly, the action is given by

$S=\int d\tau\int
d\bm{r}\,\frac{1}{2}\,\overline{\bm{\Psi}}\left(\partial_{\tau}+H_{\text{BdG}}\right)\bm{\Psi},$
(181)

where we suppressed a constant arising from the molecular field approximation
and neglected other channels. $\bm{\Psi}$ is the Nambu Grassman vector
consisting of the Grassmann numbers for electrons and holes,
$\bm{\Psi}^{T}=\left(\psi_{1},\cdots\psi_{n},\overline{\psi}_{1},\cdots,\overline{\psi}_{n}\right)$.
Since the action is quadratic in the Grassmann numbers, we can analytically
perform the Grassmann integration. The partition function is obtained as

$Z=\operatorname{Tr}{e^{-\beta\mathcal{H}}}=\int\mathcal{D}\overline{\psi}\mathcal{D}{\psi}e^{-S}=\int
d\tau\int d\bm{r}\frac{1}{2}\operatorname{Tr}{\log{\mathcal{G}^{-1}(\tau)}}.$
(182)

We introduced the Gor’kov Green function $\mathcal{G}$

$\mathcal{G}^{-1}=\begin{pmatrix}\partial_{\tau}+H_{\text{N}}&\Delta\\\
\Delta^{\dagger}&\partial_{\tau}-\left(H_{\text{N}}\right)^{T}\end{pmatrix},$
(183)

where the imaginary-time derivative is performed on the right hand side. It is
necessary for calculations of the conductivity to take into account the
electromagnetic fields. In this work, we introduce the minimal coupling
between the vector potential and the canonical momentum in the normal-state
Hamiltonian as

$H_{\text{N}}[\hat{\bm{p}}]\rightarrow
H_{\text{N}}[\hat{\bm{p}}-{\bm{A}}(\tau)]\equiv H_{\text{N}}({\bm{A}}(\tau)).$
(184)

We take the vector potential which yields the photo-electric field in the
real-time domain as $\bm{E}=-\partial_{t}\bm{A}(t)$. Correspondingly, the
Gor’kov Green function is replaced with $\mathcal{G}_{\bm{A}}$ whose the
normal-state Hamiltonian is $H_{\text{N}}({\bm{A}}(\tau))$.

The expectation value of the electric current is obtained in a straightforward
way once the partition function is given. The electric current is calculated
as

$\displaystyle\Braket{\mathcal{J}^{\alpha}(\tau)}$
$\displaystyle=-\frac{\delta}{\delta A_{\alpha}(\tau)}\int d\tau^{\prime}\int
d\bm{r}^{\prime}\frac{1}{2}\operatorname{Tr}{[\log{\mathcal{G}_{\bm{A}}^{-1}}(\tau^{\prime})]},$
(185) $\displaystyle=\frac{1}{2}\int
d\bm{r}\operatorname{Tr}{{J}^{\alpha}_{{\bm{A}}}(\tau)\mathcal{G}_{\bm{A}}(\tau)},$
(186)

where the coefficient $1/2$ is due to the particle hole doubling in the BdG
Hamiltonian. ${J}^{\alpha}_{{\bm{A}}}$ is the current operator for the BdG
Hamiltonian and defined as

${J}^{\alpha}_{\bm{A}}(\tau)=-\frac{\delta}{\delta
A_{\alpha}(\tau)}\begin{pmatrix}H_{\text{N}}({\bm{A}}(\tau))&0\\\
0&-\left[H_{\text{N}}({\bm{A}}(\tau))\right]^{T}\end{pmatrix}.$ (187)

The linear and nonlinear conductivities are systematically calculated by
expanding the current operator ${J}^{\alpha}_{{\bm{A}}}$ and Green function
$\mathcal{G}_{\bm{A}}$ with respect to the vector potential ${\bm{A}}$. The
derivation is similarly performed as in Ref. Michishita2020 . We show the
final result of the NRO conductivity

$\displaystyle
2\sigma^{\alpha;\beta\gamma}\left(\omega;\omega_{1},\omega_{2}\right)$
$\displaystyle=\frac{1}{\omega_{1}\omega_{2}}\int_{-\infty}^{\infty}\frac{d\Omega}{2\pi
i}f(\Omega)\frac{1}{2}\left\\{\frac{1}{2}\operatorname{Tr}\left[{J}^{\alpha\beta\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\right]\right.$
(188)
$\displaystyle+\operatorname{Tr}\left[{J}^{\alpha\beta}G^{R}\left(\Omega+\omega_{2}\right){J}^{\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)+{J}^{\alpha\beta}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}G^{A}\left(\Omega-\omega_{2}\right)\right]$
(189)
$\displaystyle+\frac{1}{2}\operatorname{Tr}\left[{J}^{\alpha}G^{R}\left(\Omega+\omega\right){J}^{\beta\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)+{J}^{\alpha}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\beta\gamma}G^{A}\left(\Omega-\omega\right)\right]$
(190)
$\displaystyle+\operatorname{Tr}\left[{J}^{\alpha}G^{R}\left(\Omega+\omega\right){J}^{\beta}G^{R}\left(\Omega+\omega_{2}\right){J}^{\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\right.$
(191)
$\displaystyle+{J}^{\alpha}G^{R}\left(\Omega+\omega_{1}\right){J}^{\beta}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}G^{A}\left(\Omega-\omega_{2}\right)$
(192)
$\displaystyle+\left.{J}^{\alpha}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\beta}G^{A}\left(\Omega-\omega_{1}\right){J}^{\gamma}G^{A}\left(\Omega-\omega\right)\right]$
(193)
$\displaystyle\left.+\left[\left(\beta,\omega_{1}\right)\leftrightarrow\left(\gamma,\omega_{2}\right)\right]\right\\},$

where $\operatorname{Tr}{}$ represents the trace over eigenstates of the
Bogoliubov quasiparticle. We introduced the retarded ($G^{R}$) and advanced
($G^{A}$) Green functions as

$G^{R}(\omega)=\frac{1}{\omega+i\eta-
H_{\text{BdG}}},~{}~{}~{}G^{A}(\omega)=\left(G^{R}(\omega)\right)^{\dagger}.$
(194)

Here, $\eta>0$ is an infinitesimal parameter building the causality into the
Green functions, $H_{\text{BdG}}$ is the BdG Hamiltonian having no
electromagnetic perturbation, and
$\left[\left(\beta,\omega_{1}\right)\leftrightarrow\left(\gamma,\omega_{2}\right)\right]$
means the symmetrization of the applied electric fields due to the intrinsic
permutation symmetry of external fields. The doubling of the conductivity
($2\sigma^{\alpha;\beta\gamma}$) respects the coefficient $1/2$ in Eq. (186).

The low-frequency behavior of the NRO conductivity is of current interest.
Thus, we take the $(\omega,\omega_{1},\omega_{2})$\- expansion from the static
limit with the frequency conservation condition,
$\omega=\omega_{1}+\omega_{2}$. The leading terms are given by

$2\sigma^{\alpha;\beta\gamma}\left(\omega_{1}+\omega_{2};\omega_{1},\omega_{2}\right)=\frac{1}{2\omega_{1}\omega_{2}}A_{\alpha\beta\gamma}+\frac{1}{2\omega_{2}}B_{\alpha\beta\gamma}+\frac{1}{2\omega_{1}}B^{\prime}_{\alpha\beta\gamma}+O(\omega^{0}).$
(195)

The coefficients are obtained as

$\displaystyle A_{\alpha\beta\gamma}$
$\displaystyle=\int_{-\infty}^{\infty}\frac{d\Omega}{2\pi
i}f(\Omega)\left\\{\frac{1}{2}\operatorname{Tr}\left[{J}^{\alpha\beta\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\right]\right.$
(196)
$\displaystyle+\operatorname{Tr}\left[{J}^{\alpha\beta}G^{R}\left(\Omega\right){J}^{\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)+{J}^{\alpha\beta}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}G^{A}\left(\Omega\right)\right]$
(197)
$\displaystyle+\frac{1}{2}\operatorname{Tr}\left[{J}^{\alpha}G^{R}\left(\Omega+\omega\right){J}^{\beta\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)+{J}^{\alpha}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\beta\gamma}G^{A}\left(\Omega-\omega\right)\right]$
(198)
$\displaystyle+\operatorname{Tr}\left[{J}^{\alpha}G^{R}\left(\Omega\right){J}^{\beta}G^{R}\left(\Omega\right){J}^{\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\right.$
(199)
$\displaystyle+{J}^{\alpha}G^{R}\left(\Omega\right){J}^{\beta}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}G^{A}\left(\Omega\right)$
(200)
$\displaystyle+\left.{J}^{\alpha}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\beta}G^{A}\left(\Omega\right){J}^{\gamma}G^{A}\left(\Omega\right)\right]$
(201) $\displaystyle\left.+\left(\beta\leftrightarrow\gamma\right)\right\\}.$

and

$\displaystyle B_{\alpha\beta\gamma}=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\Omega}{2\pi
i}f(\Omega)\Bigl{\\{}\operatorname{Tr}\left[{J}^{\alpha\beta}F^{R}(\Omega){J}^{\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)-{J}^{\alpha\beta}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}F^{A}\left(\Omega\right)\right]$
(202)
$\displaystyle+\operatorname{Tr}\left[{J}^{\alpha}F^{R}\left(\Omega\right){J}^{\beta\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)-{J}^{\alpha}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\beta\gamma}F^{A}\left(\Omega\right)\right]$
(203)
$\displaystyle+\operatorname{Tr}\left[{J}^{\alpha}F^{R}\left(\Omega\right){J}^{\beta}G^{R}\left(\Omega\right){J}^{\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)+{J}^{\alpha}G^{R}\left(\Omega\right){J}^{\beta}F^{R}\left(\Omega\right){J}^{\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\right.$
$\displaystyle~{}~{}~{}~{}~{}\left.+{J}^{\alpha}F^{R}\left(\Omega\right){J}^{\gamma}G^{R}\left(\Omega\right){J}^{\beta}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\right.$
$\displaystyle~{}~{}~{}~{}~{}\left.-{J}^{\alpha}G^{R}\left(\Omega\right){J}^{\beta}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}F^{A}\left(\Omega\right)+{J}^{\alpha}F^{R}\left(\Omega\right){J}^{\beta}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}G^{A}\left(\Omega\right)\right.$
$\displaystyle~{}~{}~{}~{}~{}\left.-{J}^{\alpha}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\beta}G^{A}\left(\Omega\right){J}^{\gamma}F^{A}\left(\Omega\right)-{J}^{\alpha}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}G^{A}\left(\Omega\right){J}^{\beta}F^{A}\left(\Omega\right)\right.$
$\displaystyle~{}~{}~{}~{}~{}\left.-{J}^{\alpha}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}F^{A}\left(\Omega\right){J}^{\beta}G^{A}\left(\Omega\right)\right]\Bigr{\\}}.$
(204)

$F^{R(A)}(\omega)=\partial_{\omega}G^{R(A)}$ is the frequency derivative of
the Green function. The coefficient $B^{\prime}_{\alpha\beta\gamma}$ is
obtained by the permutation $\beta\leftrightarrow\gamma$ in
$B_{\alpha\beta\gamma}$.

Simplification of these coefficients is performed by making use of the
variational parameter ${\bm{\lambda}}$ as in the density matrix formulation.
Let us consider the adiabatic switching of ${\bm{\lambda}}$ in the BdG
Hamiltonian as

$H_{\text{BdG}}=\lim_{\bm{\lambda}\rightarrow\bm{0}}H_{\text{BdG}}^{({\bm{\lambda}})}.$
(205)

The right hand side is written by

$H_{\text{BdG}}^{({\bm{\lambda}})}=\begin{pmatrix}H_{\text{N}}^{{\bm{k}}-{\bm{\lambda}}}&\hat{\Delta}_{\bm{k}}\\\
\hat{\Delta}_{\bm{k}}^{{\dagger}}&-(H_{\text{N}}^{-{\bm{k}}-{\bm{\lambda}}})^{T}\end{pmatrix},$
(206)

in the crystal momentum (${\bm{k}}$) basis. Normal-state Hamiltonian and pair
potentials are accordingly given in the ${\bm{k}}$ representation. Following
the definition of the current operators, we can rewrite the operators by

$\displaystyle
J^{\alpha}=\lim_{\bm{\lambda}\rightarrow\bm{0}}J^{\alpha}_{\bm{\lambda}}=\lim_{\bm{\lambda}\rightarrow\bm{0}}-\frac{\partial
H_{\text{BdG}}^{({\bm{\lambda}})}}{\partial\lambda_{\alpha}},$ (207)
$\displaystyle
J^{\alpha\beta}=\lim_{\bm{\lambda}\rightarrow\bm{0}}J^{\alpha\beta}_{\bm{\lambda}}=\lim_{\bm{\lambda}\rightarrow\bm{0}}\frac{\partial^{2}H_{\text{BdG}}^{({\bm{\lambda}})}}{\partial\lambda_{\alpha}\partial\lambda_{\beta}},$
(208) $\displaystyle
J^{\alpha\beta\gamma}=\lim_{\bm{\lambda}\rightarrow\bm{0}}J^{\alpha\beta\gamma}_{\bm{\lambda}}=\lim_{\bm{\lambda}\rightarrow\bm{0}}-\frac{\partial^{3}H_{\text{BdG}}^{({\bm{\lambda}})}}{\partial\lambda_{\alpha}\partial\lambda_{\beta}\partial\lambda_{\gamma}}.$
(209)

Owing to the correction to the BdG Hamiltonian, the Green functions are
modified. For instance, the retarded Green function is

$G^{R}_{\bm{\lambda}}(\omega)=\frac{1}{\omega+i\eta-
H_{\text{BdG}}^{({\bm{\lambda}})}}.$ (210)

Differentiating the Green functions by ${\bm{\lambda}}$, we obtain

$\frac{\partial
G^{R}_{\bm{\lambda}}(\omega)}{\partial\lambda_{\alpha}}=-G^{R}_{\bm{\lambda}}(\omega){J}^{\alpha}_{\bm{\lambda}}G^{R}_{\bm{\lambda}}(\omega).$
(211)

When we make use of the ${\bm{\lambda}}$-switching and the relation (211), we
can simplify the coefficients
$A_{\alpha\beta\gamma},~{}B_{\alpha\beta\gamma},$ and
$B^{\prime}_{\alpha\beta\gamma}$. First, we consider $A_{\alpha\beta\gamma}$.
Based on the adiabatic switching of ${\bm{\lambda}}$, the current operators
and Green functions are converted into those labeled by ${\bm{\lambda}}$. For
instance, we take the line (196). The third-order current operator is related
with the diamagnetic current operator as
${J}^{\alpha\beta\gamma}_{\bm{\lambda}}=-\partial_{\lambda_{\alpha}}{J}^{\beta\gamma}_{\bm{\lambda}}$.
Accordingly,

$\displaystyle\operatorname{Tr}\left[{J}^{\alpha\beta\gamma}_{{\bm{\lambda}}}\left(G^{R}_{\bm{\lambda}}(\Omega)-G^{A}_{\bm{\lambda}}(\Omega)\right)\right]$
$\displaystyle=-\partial_{\lambda_{\alpha}}\left(\operatorname{Tr}\left[{J}^{\beta\gamma}_{{\bm{\lambda}}}\left(G^{R}_{\bm{\lambda}}(\Omega)-G^{A}_{\bm{\lambda}}(\Omega)\right)\right]\right)+\operatorname{Tr}\left[{J}^{\beta\gamma}_{{\bm{\lambda}}}\partial_{\lambda_{\alpha}}\left(G^{R}_{\bm{\lambda}}(\Omega)-G^{A}_{\bm{\lambda}}(\Omega)\right)\right],$
(212)
$\displaystyle=-\partial_{\lambda_{\alpha}}\left(\operatorname{Tr}\left[{J}^{\beta\gamma}_{{\bm{\lambda}}}\left(G^{R}_{\bm{\lambda}}(\Omega)-G^{A}_{\bm{\lambda}}(\Omega)\right)\right]\right)-\operatorname{Tr}\left[{J}^{\beta\gamma}_{{\bm{\lambda}}}\left(G^{R}_{\bm{\lambda}}(\Omega){J}^{\alpha}_{{\bm{\lambda}}}G^{R}_{\bm{\lambda}}(\Omega)-G^{A}_{\bm{\lambda}}(\Omega){J}^{\alpha}_{{\bm{\lambda}}}G^{A}_{\bm{\lambda}}(\Omega)\right)\right],$
(213)

Using the cyclic property of trace, the second term in the last line is
canceled with Eq. (198). Performing a similar algebra, we finally rewrite
$A_{\alpha\beta\gamma}$ by

$\displaystyle A_{\alpha\beta\gamma}$ $\displaystyle=\int\frac{d\Omega}{2\pi
i}f(\Omega)\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\beta}}\partial_{\lambda_{\gamma}}\operatorname{Tr}{\left[{J}^{\alpha}_{{\bm{\lambda}}}\left(G^{R}_{\bm{\lambda}}(\Omega)-G^{A}_{\bm{\lambda}}(\Omega)\right)\right]},$
(214)
$\displaystyle=-\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\beta}}\partial_{\lambda_{\gamma}}\operatorname{Tr}{\left[{J}^{\alpha}_{{\bm{\lambda}}}f(H_{\text{BdG}}^{({\bm{\lambda}})})\right]},$
(215)