# Inverse Faraday effect in superconductors with potential impurities

Maxim Dzero  Department of Physics, Kent State University, Kent, OH 44242, USA

###### Abstract

I consider a nonlinear response of conventional superconductors contaminated
with nonmagnetic impurities to a circular polarized light. I focus on dc
contributions to the induced current density which describe the emergence of
the static magnetization in a superconductor. This effect is known as inverse
Faraday effect. By employing quasiclassical theory of superconductivity I
derive an expression for the induced static magnetization as a function of
frequency of external ac field and disorder scattering rate. The scattering of
electrons off potential impurities is taken into account within the framework
of the self-consistent Born approximation. It is found that the magnitude of
the inverse Faraday effect decreases with an increase in disorder scattering
rate. I have also discovered that the value of the induced magnetization has a
characteristic minimum at a frequency which approximately equals twice the
value of the pairing gap in a clean superconductor and shifts to higher values
with an increase in disorder scattering rate.

##  I Introduction

In the last decade there have been significant advances in state-of-the-art
optical measurements and techniques that have spurred a lot of interest in the
physical phenomena related to nonlinear responses in conventional and
unconventional superconductors THz1 ; THz2 ; Shimano2012 ; Shimano2013 ;
Shimano2014 ; THz3 ; Millis2020 ; Armi2014 ; Armitage2023 . These experimental
and theoretical studies usually uncover and address a number of interesting
and important problems. Some of these problems are concerned with the
fundamental aspects of superconductivity such as various proposals to probe
the physics associated with the excitation the amplitude Higgs mode Varma2014
; Cea2016 ; Shimano2020 ; Haenel2021-Disorder ; Sherman2015-Disorder ; Spivak-
Higgs  and Eliashberg effect Eliashberg1970 ; Eremin2023  to name a few.
Another set of problems focuses on development of potential applications of
various nonlinear effects in superconductor-based devices for it offers one a
possibility to manipulate the physical properties of superconductors such as
magnetic and transport response functions on extremely short time scales
Chen2017 ; Croitoru2022 ; SpivakGiant2020 .

Inverse Faraday effect (IFE) in dispersive medium was discovered by L. P.
Pitaevskii Pit1961  who demonstrated that static magnetization can be produced
by subjecting a medium to high-intensity circular polarized light. In
subsequent years this effect has been and continues to be actively studied in
various electronic systems Pershan1966 ; Battiato2014 ; yang2022inverse ;
yang2022inverse ; mou2023reversed ; mou2023chiral ; IFE-Mott2022 ;
gao2020topological ; parchenko2023plasmonenhanced ; Han_2023 . For example,
Majedi Majedi2021 , using phenomenological set of arguments, has recently
proposed that the IFE can be induced in superconductors by subjecting them to
an external microwave radiation, Fig. 1. Almost immediately after the
publication of that work, several groups have further developed this idea by
studying the IFE in superfluid condensates Mironov2021-IFESC  and in various
superconductor-based devices Parafilo2022Fl ; Putilov2023-IFESC ; Buzdin2023 ;
Croitoru2023 . Most recently, Sharma and Balatsky Balatsky2023  has utilized
the quasiclassical approach to formulate a microscopic theory of the IFE in
superconductors.

Figure 1: Schematic representation for the realization of inverse Faraday
effect in superconductors. When a conventional superconductor is subjected to
a circular polarized light with frequency $\omega$, it induces a static
component in the current density and, as a consequence, a static magnetization
${\mathbf{M}}_{\textrm{ind}}$. The effect is of the second order in the powers
of the amplitude $E_{0}$ of an external electric field.

In the context of the present work it must be mentioned that in the earlier
works Majedi2021 ; Balatsky2023  scattering due to potential disorder, which
is almost inevitably present in superconducting samples, has been taken into
account on purely phenomenological level by formally assigning finite
imaginary part to the single particle energies,
$\varepsilon\to\varepsilon+i\Gamma$. It is well known that weak potential
disorder does not affect the physical properties of conventional
superconductors in equilibrium (Anderson theorem) AndersonTheorem ; AG1961 ;
Balatsky-RMP  and for that reason the approximations taken in the above
mentioned works may seem well justified. However, as it has been established
independently by several groups, potential impurities do affect the
superconducting response functions provided that a superconductor has been
driven out-of-equilibrium Silaev2019-Disorder ; Yang2022-Disorder ;
Yang2022-Disorder2 ; Dzero2023-Disorder .

Within the theoretical framework of the quasiclassical theory of
superconductivity LO  this can be understood as follows. The self-energy due
to scattering on potential impurities $\check{\Sigma}$ enters in the
quasiclassical equations in a commutator with the single-particle matrix
correlation function $\check{g}$. Since in Nambu representation the self-
energy part due to disorder scattering is proportional to $\check{g}$ the
corresponding commutator vanishes identically. Qualitatively, this means that
Cooper pairs scatter coherently off the weak potential impurities and pair
breaking processes do not occur. When a superconductor is driven out of
equilibrium, the commutator will include the derivatives of both
$\check{\Sigma}$ and $\check{g}$ with respect to time and, as a result, no
such cancellation occurs. In other words, the monochromatic radiation breaks
time-reversal symmetry and, as a result, scattering on potential impurities
induces pair-breaking. Therefore, in the context of the IFE in
superconductors, it is important to elucidate the role of nonmagnetic disorder
and how it affects the ac-field induced dc-magnetization.

In this paper I formulate a microscopic theory of inverse Faraday effect in
disordered superconductors in which the effects of potential disorder are
treated self-consistently. In particular I demonstrate that under the action
of external ac electric field the self-energy correction due to disorder
scattering leads to the renormalization of the quasiparticle energies and
pairing amplitude similar to the renormalization effects which appear due to
the presence of paramagnetic impurities. I analyze the IFE by calculating the
static contribution to the charge current density to the second order in
external electric field. I have found that the magnitude of the induced dc
magnetization decreases with the increase in disorder scattering rate.
Intriguingly, I also find that the magnetization as a function of disorder
scattering rate is suppressed most strongly when the frequency of electric
field coincides with the resonant frequency of the amplitude Higgs mode
$\omega_{\textrm{res}}$. Furthermore, magnetization as a function of frequency
for a fixed disorder strength exhibits a minimum at
$\omega=\omega_{\textrm{res}}$. Therefore, based on these observations I
suggest that induced dc magnetization via IFE in conventional superconductors
contains a ’fingerprint’ of the amplitude Higgs mode.

##  II Qualitative discussion

When a conventional superconductor is subjected to a monochromatic external
electric field

${\mathbf{E}}_{\textrm{ext}}(\mathbf{r},t)={\mathbf{E}}_{0}e^{i\mathbf{k}\mathbf{r}-i\omega
t}+{\mathbf{E}}_{0}^{*}e^{-i\mathbf{k}\mathbf{r}+i\omega t}$ (1)

with ${\mathbf{E}}_{0}=E_{0}(\hat{\mathbf{x}}+i\hat{\mathbf{y}})$ (here
$\hat{\mathbf{x}}$ and $\hat{\mathbf{y}}$ are the unit vectors), Fig. 1, the
superconducting order parameter, defined by a scalar quantity, must couple to
the electric field nonlinearly. As a result of this nonlinear coupling two
collective modes will be excited: Clarkson-Goldman (phase) mode and amplitude
Higgs mode. In the context of neutral superfluids, the former is referred to
as Anderson-Bogoliubov mode Kamenev2011 . Usually, these modes are analyzed
separately because in the Nambu representation they are described by the
correlation functions with different matrix structure (see discussion in the
Section III below). In addition, the phase mode is gapless while the amplitude
mode is gapped at twice the value of the pairing amplitude in equilibrium,
$\omega_{\textrm{res}}=2\Delta_{0}$.

As it turns out there are several physical effects, such as inverse Faraday
and ponderomotive force effects, for which an interplay between phase and
amplitude modes becomes crucial and so it could be important for deeper
understanding of these effects to find out how this interplay can be
identified in the physical observables. In order to succeed in achieving this
goal, one needs to treat these modes on equal footing. For example, in the IFE
the phase fluctuations produce the correction to the particle density
distribution which, due to the fact that phase mode is gapless, must be
linearly proportional to the gradient of an external electric field, $\delta
n\sim({\mathbf{v}}_{F}\cdot{\mbox{\boldmath$\nabla$}})(\mathbf{n}\cdot{\mathbf{E}}_{\textrm{ext}}(\mathbf{r},t))$
(here ${\mathbf{v}}_{F}=v_{F}\mathbf{n}$ is the Fermi velocity). On the other
hand, the amplitude mode is gapped and as a result fluctuation associated with
it produce the correction to particle velocity
$\delta{\mathbf{v}}\sim{\mathbf{v}}_{F}(\mathbf{n}\cdot{\mathbf{E}}_{\textrm{ext}}(\mathbf{r},t))$
and gradient corrections can only appear in the second order, i.e. the
amplitude mode is diffusive, $\omega_{\textrm{amp}}(q)\approx
2\Delta_{0}+Dq^{2}$ (here $D=v_{F}^{2}\tau/2$ is the diffusion coefficient for
a two-dimensional superconductor). As a result, the second order correction to
the current in powers of electric field can be computed using the hydrodynamic
expression ${\mathbf{j}}^{(2)}=\langle\langle\delta
n\delta{\mathbf{v}}\rangle_{\mathbf{n}}\rangle_{2\pi/\omega}$, where the
averaging is performed over the directions of the Fermi velocity and period of
oscillating electric field. Among several contributions, an expression for
${\mathbf{j}}^{(2)}$ will contain the terms describing the gyration of
electric field which gives rise to the _static_ magnetization
${\mathbf{M}}_{\textrm{ind}}(\omega)=if_{\omega}{\mathbf{E}}_{\textrm{ext}}\times{\mathbf{E}}_{\textrm{ext}}^{*}$
(see Fig. 1). From quite general considerations one may immediately discover
that ${\mathbf{M}}_{\textrm{ind}}(\omega)\propto e^{2}v_{F}^{4}/\omega^{3}$
and we can represent the resulting expression for the induced magnetization as
${\mathbf{M}}_{\textrm{ind}}(\omega)=ie^{2}v_{F}^{4}\gamma_{\omega}({\mathbf{E}}_{\textrm{ext}}\times{\mathbf{E}}_{\textrm{ext}}^{*})$,
where function $\gamma_{\omega}\propto 1/\omega^{3}$ is determined by the
microscopic details such as disorder, interactions between the constituent
quasiparticles etc. Balatsky2023  By virtue of the logic above, one then
concludes that the induced magnetization must also exhibit a non-monotonic
frequency dependence - a minimum - around $\omega\approx 2\Delta_{0}$ which
corresponds to the resonant frequency for the excitation of the amplitude
Higgs mode (see above). However, in order to check how prominent this feature
will be it requires a microscopic calculation. In passing we note that since
the amplitude fluctuations can only produce second order gradient corrections,
these fluctuations alone cannot contribute to the IFE and so at this point it
is not immediately clear how to capture this effect if one chooses to work in
the dirty limit $\tau\Delta_{0}\ll 1$ by employing the Usadel equation
Usadel1970 . In the following Section I provide the microscopic derivation of
the expression for the function $\gamma_{\omega}$ using the quasiclassical
theory of superconductivity LO ; Usadel1970 ; Kita2001 .

##  III Main equations

In the problems which are concerned with response of superconductors to the
application of external electromagnetic radiation, it proves important to
write down the equations for the gauge invariant single-particle propagators
Kopnin1994 ; Kita2001

$\check{g}(1,2)=\left(\begin{matrix}\hat{g}^{R}(1,2)&\hat{g}^{K}(1,2)\\\
0&\hat{g}^{A}(1,2)\end{matrix}\right)$ (2)

within the framework of the quantum field theory for non-equilibrium systems
Kamenev2009 ; Kamenev2011 . Each component of the matrix function
$\check{g}(1,2)$ is defined in Keldysh and Nambu spaces and the argument of
the propagator should be understood as
$(1,2)=(\mathbf{r}_{1},t_{1};\mathbf{r}_{2},t_{2})$. One can follow the avenue
of Refs. Kita2001 ; Kita-Review  to derive the quasiclassical equation for the
Wigner transformed single-particle propagator $\check{g}$:

$\begin{split}&[\varepsilon\check{\tau}_{3},\check{g}]+\frac{e}{c}{\mathbf{v}}_{F}{\mathbf{A}}(\mathbf{r},t)[\check{\tau}_{3},\check{g}]+i{\mathbf{v}}_{F}{\partial_{\mathbf{r}}}\check{g}-[\check{\Sigma}\check{\tau}_{3}\stackrel{{\scriptstyle\circ}}{{,}}\check{g}]+\frac{ie}{2}\\\
&\times\left[\partial_{\varepsilon}\left(\check{g}\check{\tau}_{3}\right){\cal
E}\left(i\stackrel{{\scriptstyle\leftrightarrow}}{{\partial}}_{\varepsilon
t}\right)+\partial_{\varepsilon}\left(\check{\tau}_{3}\check{g}\right){\cal
E}\left(-i\stackrel{{\scriptstyle\leftrightarrow}}{{\partial}}_{\varepsilon
t}\right)\right]{\mathbf{v}}_{F}{\mathbf{E}}(\mathbf{r},t)=0.\end{split}$ (3)

In this equation
$\stackrel{{\scriptstyle\leftrightarrow}}{{\partial}}_{\varepsilon
t}=\frac{1}{2}\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{\varepsilon}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{t}$,
$\check{g}=\check{g}(\mathbf{n}\varepsilon;\mathbf{r}t)$,
$\mathbf{r}=(\mathbf{r}_{1}+\mathbf{r}_{2})/2$, $t=(t_{1}+t_{2})/2$,
${\mathbf{v}}_{F}=v_{F}\mathbf{n}$ is the Fermi velocity,
${\mathbf{A}}(\mathbf{r},t)$ is a vector potential,
${\mathbf{E}}(\mathbf{r},t)=-(1/c)\partial_{t}{\mathbf{A}}(\mathbf{r},t)$ is
an external electric field, $\check{\tau}_{3}$ is the third Pauli matrix which
is diagonal in Keldysh subspace, function ${\cal E}(u)=(e^{u}-1)/u$,
$\check{\Sigma}$ is the self-energy part which encodes the effects of
superconducting pairing and disorder:

$\check{\Sigma}\check{\tau}_{3}=-\frac{i}{2\tau}\langle\check{g}\rangle_{\mathbf{n}}+i\check{\tau}_{2}\Delta^{\prime}+i\check{\tau}_{1}\Delta^{\prime\prime},$
(4)

$\check{\tau}_{1(2)}=\hat{\mathbbm{1}}_{2\times 2}\otimes\hat{\tau}_{1(2)}$
are first and second Pauli matrices, which act in Nambu and Keldysh spaces and
are diagonal in Keldysh space. Lastly, the commutator is defined according to

$[A\stackrel{{\scriptstyle\circ}}{{,}}B]=Ae^{\frac{i}{2}\left(\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{\varepsilon}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{t}-\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{t}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{\varepsilon}\right)}B-Be^{\frac{i}{2}\left(\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{\varepsilon}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{t}-\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{t}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{\varepsilon}\right)}A.$
(5)

Given the expression for the electric field (1), I will look for the solution
of (3) in the form

$\begin{split}\check{g}(\mathbf{n}\varepsilon;\mathbf{r}t)=\check{g}(\mathbf{n}\varepsilon;\mathbf{k}\omega)e^{i(\mathbf{k}\mathbf{r}-\omega
t)}+\check{g}(\mathbf{n}\varepsilon;-\mathbf{k},-\omega)e^{-i(\mathbf{k}\mathbf{r}-\omega
t)}.\end{split}$ (6)

It is important to emphasize here that equation (3) does not preserve the norm
of $\check{g}$, i.e. $\check{g}^{2}=\hat{\mathbbm{1}}$ will be satisfied only
when either ${\mathbf{E}}(\mathbf{r},t)=0$ or in the limit
$|\varepsilon|\to\infty$ Kita2001 . The similar situation happens in a problem
of non-equilibrium dynamics of the a superconductor with paramagnetic
impurities Dzero2023-Disorder  which has one common feature with our problem
at hand: the lack of time-reversal symmetry.

In the ground state, the Keldysh propagator is just a parametrization

$\hat{g}_{\varepsilon}^{K}=(\hat{g}_{\varepsilon}^{R}-\hat{g}_{\varepsilon}^{A})\tanh\left(\frac{\varepsilon}{2T}\right),$
(7)

where $T$ is temperature. In the ground state we can set
$\Delta^{\prime}=\Delta_{0}$, $\Delta^{\prime\prime}=0$ and choose
$\hat{g}_{\varepsilon}^{R(A)}=g_{\varepsilon}^{R(A)}\hat{\tau}_{3}-if_{\varepsilon}^{R(A)}\hat{\tau}_{2}$.
By employing the normalization condition
$\hat{g}_{\varepsilon}^{R(A)}\hat{g}_{\varepsilon}^{R(A)}=\hat{\mathbbm{1}}$ I
find $g_{\varepsilon}^{R(A)}=(\varepsilon\pm i0)/\zeta_{\varepsilon}^{R(A)}$,
$f_{\varepsilon}^{R(A)}=\Delta_{0}/\zeta_{\varepsilon}^{R(A)}$ and

$\zeta_{\varepsilon}^{R(A)}=\left\\{\begin{matrix}\pm\textrm{sign}(\varepsilon)\sqrt{(\varepsilon\pm
i0)^{2}-\Delta_{0}^{2}},\quad|\varepsilon|>\Delta_{0},\\\
i\sqrt{\Delta_{0}^{2}-\varepsilon^{2}},\quad|\varepsilon|<\Delta_{0}.\end{matrix}\right.$
(8)

The value of the superconducting order parameter in equilibrium is then
determined self-consistently from

$\Delta_{0}=\frac{\lambda}{2}\int\limits_{-\omega_{D}}^{\omega_{D}}f_{\varepsilon}^{K}d\varepsilon,$
(9)

where $\lambda$ is the dimensionless interaction strength and $\omega_{D}$ is
the Debye frequency. Naturally, the disorder scattering rate does not enter
into this expressions since the commutator in (3) with the first term in (4)
vanishes in equilibrium.

The expressions for the particle density and particle velocity for a given
${\mathbf{v}}_{F}$ are expressed in terms of the Keldysh propagator as
follows:

$\begin{split}n(\mathbf{n};\mathbf{r}t)&=\frac{\pi\nu_{F}}{2}\int\textrm{Tr}\left\\{\hat{g}^{K}(\mathbf{n}\varepsilon;\mathbf{r}t)\right\\}d\varepsilon,\\\
{\mathbf{v}}(\mathbf{n};\mathbf{r}t)&=\frac{\pi\nu_{F}}{2}\int{\mathbf{v}}_{F}\textrm{Tr}\left\\{\hat{\tau}_{3}\hat{g}^{K}(\mathbf{n}\varepsilon;\mathbf{r}t)\right\\}d\varepsilon.\end{split}$
(10)

Here $\nu_{F}$ is the single particle density of states at the Fermi level. In
terms of these quantities, the macroscopic current density is defined as

${\mathbf{j}}(\mathbf{r},t)=\langle
n(\mathbf{n};\mathbf{r},t){\mathbf{v}}(\mathbf{n};\mathbf{r},t)\rangle_{\mathbf{n}}.$
(11)

Although the expressions above are valid in any spacial dimenions, for the
calculation of the averages over the Fermi surface, it will be assumed without
loss of generality that we have a two-dimensional superconductor, Fig. 1. As I
have discussed above, our goal is to compute linear corrections for the
Keldysh Green’s function $\hat{g}^{K}$ in electric field, which will allow us
to evaluate nonlinear in electric field correction to the current density
(11).

###  III.1 Linear analysis: retarded and advanced functions

We start by computing the linear corrections to $\hat{g}^{R(A)}$. I will look
for the solution of equation (3) with

$\hat{g}^{R(A)}=\hat{g}_{\varepsilon}^{R(A)}+\hat{g}_{1}^{R(A)},$ (12)

where $\hat{g}_{1}^{R(A)}\propto\mathbf{n}{\mathbf{E}}_{0}$. Inserting this
expression into (3) and keeping only terms to the linear order in
${\mathbf{E}}_{0}$ I obtain the following equation for the matrix function
$\hat{g}_{1}^{R(A)}(\mathbf{n}\varepsilon;\mathbf{k}\omega)$:

$\begin{split}&[\varepsilon\hat{\tau}_{3},\hat{g}_{1}]+\omega\\{\hat{\tau}_{3},\hat{g}_{1}\\}-v_{F}(\mathbf{n}\mathbf{k})\hat{g}_{1}-i\Delta_{0}[\hat{\tau}_{2},\hat{g}_{1}]\\\
&-\frac{i}{2\tau}\left(\hat{g}_{1}-\langle\hat{g}_{1}\rangle_{\mathbf{n}}\right)\hat{g}_{\varepsilon-\frac{\omega}{2}}+\frac{i}{2\tau}\hat{g}_{\varepsilon+\frac{\omega}{2}}\left(\hat{g}_{1}-\langle\hat{g}_{1}\rangle_{\mathbf{n}}\right)\\\
&=\frac{\mathbf{n}{\mathbf{Q}}_{0}}{\omega}\\{\hat{\tau}_{3},\hat{g}_{\varepsilon-\omega/2}-\hat{g}_{\varepsilon+\omega/2}\\}+i\delta\overline{\Delta}_{\mathbf{k}\omega}\hat{\tau}_{1}\hat{g}_{\varepsilon-\omega/2}\\\
&-i\delta\overline{\Delta}_{\mathbf{k}\omega}\hat{g}_{\varepsilon+\omega/2}\hat{\tau}_{1}+\frac{\mathbf{n}{\mathbf{Q}}_{0}}{\omega}[\hat{\tau}_{3},\hat{g}_{\varepsilon-\omega/2}+\hat{g}_{\varepsilon+\omega/2}],\end{split}$
(13)

where ${\mathbf{Q}}_{0}=iev_{F}{\mathbf{E}}_{0}$, function
$\delta\overline{\Delta}_{\mathbf{k}\omega}$ is the Fourier component of the
linear correction to the superconducting order parameter
$\delta\Delta(\mathbf{r},t)=\delta\overline{\Delta}_{\mathbf{k}\omega}e^{i(\mathbf{k}\mathbf{r}-\omega
t)}$ and we suppressed the superscripts $R(A)$ for brevity. In deriving this
equation we have used the auxiliary expressions listed in Appendix A.

Let me first average both sides of this equation over the Fermi surface.
Denoting
$\langle\hat{g}_{1}(\mathbf{n}\varepsilon;\mathbf{k}\omega)\rangle_{\mathbf{n}}=\hat{\cal
J}_{\mathbf{k}\omega}(\varepsilon)$ I have

$\begin{split}&[\hat{\varepsilon},\hat{\cal
J}_{\mathbf{k}\omega}]+\omega\\{\hat{\tau}_{3},\hat{\cal
J}_{\mathbf{k}\omega}\\}=i\delta\overline{\Delta}_{\mathbf{k}\omega}[\hat{\tau}_{1},\hat{g}_{\varepsilon}]+v_{F}\langle\mathbf{n}\mathbf{k}\hat{g}_{1}\rangle_{\mathbf{n}},\end{split}$
(14)

where a new notation
$\hat{\varepsilon}=\varepsilon\hat{\tau}_{3}-i\Delta_{0}\hat{\tau}_{2}$ has
been introduced for brevity. Two comments are in order: (1) note that the
terms which contain disorder scattering rate do not appear in this equation;
(2) here I have neglected the higher order derivatives of the bare propagator
by replacing the commutator
$[\delta\hat{\Delta}(\mathbf{r},t)\stackrel{{\scriptstyle\circ}}{{,}}\hat{g}]$
with its equilibrium expression $[\delta\hat{\Delta}(\mathbf{r},t),\hat{g}]$
since the terms which contain $\delta\Delta_{\mathbf{k}\omega}$ will be
irrelevant for the IFE. Because $\delta\overline{\Delta}_{\mathbf{k}\omega}$
must be proportional to the electric field, it will be convenient to use
function $\delta\Delta_{\omega}$ instead, which is defined by
$\delta\overline{\Delta}_{\mathbf{k}\omega}=(\mathbf{k}{\mathbf{Q}}_{0}/2)\delta\Delta_{\omega}$.
By the same token I will re-define $\hat{\cal
J}_{\mathbf{k}\omega}(\varepsilon)=(\mathbf{k}{\mathbf{Q}}_{0}/2)\delta\hat{\cal
G}_{\varepsilon}(\omega)$. In what follows I can also neglect the anomalous
part of the third term on the right hand side of equation (14) since
$v_{F}k=(v_{F}/c)\omega\ll\omega$. The same approximation has already been
used in deriving equation (3) by ignoring the contributions from the magnetic
field Balatsky2023 ; Kita2001 .

Next, I multiply both sides of equation (13) by $\mathbf{n}\mathbf{k}$ and
average the resulting expression over the Fermi surface. This yields:

$\begin{split}&\left(\tilde{\varepsilon}_{\omega}^{R(A)}\hat{\tau}_{3}-i\tilde{\Delta}_{\omega}^{R(A)}\hat{\tau}_{2}\right)\hat{\cal
G}_{\varepsilon}^{R(A)}-\hat{\cal
G}_{\varepsilon}^{R(A)}\left(\tilde{\varepsilon}_{-\omega}^{R(A)}\hat{\tau}_{3}-i\tilde{\Delta}_{-\omega}^{R(A)}\hat{\tau}_{2}\right)\\\
&=a_{\varepsilon}^{R(A)}(\omega)\hat{\tau}_{0}-b_{\varepsilon}^{R(A)}(\omega)\hat{\tau}_{1},\end{split}$
(15)

where
$\langle(\mathbf{n}\mathbf{k})\hat{g}_{1}(\mathbf{n}\varepsilon;\mathbf{k}\omega)\rangle_{\mathbf{n}}=({\mathbf{k}{\mathbf{Q}}_{0}}/2)\hat{\cal
G}_{\varepsilon}^{R(A)}(\omega)$ and I have introduced functions

$\begin{split}\tilde{\varepsilon}_{\omega}^{R(A)}&=\omega+\varepsilon+\frac{i}{2\tau}g_{\varepsilon+\frac{\omega}{2}}^{R(A)},\\\
\tilde{\Delta}_{\omega}^{R(A)}&=\Delta_{0}+\frac{i}{2\tau}f_{\varepsilon+\frac{\omega}{2}}^{R(A)}\end{split}$
(16)

and
$a_{\varepsilon}^{R(A)}(\omega)=(g_{\varepsilon-\omega/2}^{R(A)}-g_{\varepsilon+\omega/2}^{R(A)})/\omega$,
$b_{\varepsilon}^{R(A)}(\omega)=(f_{\varepsilon-\omega/2}^{R(A)}+f_{\varepsilon+\omega/2}^{R(A)})/\omega$.
For $\omega=0$ formulas (16) coincide with the corresponding expressions in
AG1961 . By analyzing the matrix combinations which enter into both equations
(14,15), it is straightforward to realize that both of these equations can be
solved using the following ansatz

$\begin{split}&\hat{\cal G}_{\varepsilon}^{R(A)}(\omega)={\cal
G}_{\varepsilon}^{R(A)}(\omega)\hat{\tau}_{3}-i{\cal
F}_{\varepsilon}^{R(A)}(\omega)\hat{\tau}_{2},\\\ &\delta\hat{\cal
G}_{\varepsilon}^{R(A)}(\omega)=\delta{\cal
G}_{\varepsilon}^{R(A)}(\omega)\hat{\tau}_{0}+\delta{\cal
F}_{\varepsilon}^{R(A)}(\omega)\hat{\tau}_{1}.\end{split}$ (17)

To summarize, the expression for the linear correction to the retarded and
advanced parts of $\check{g}$ is formally given by

$\hat{g}_{1}^{R(A)}(\mathbf{n}\varepsilon;\mathbf{k}\omega)=\hat{\cal
G}_{\varepsilon}^{R(A)}(\omega)\mathbf{n}{\mathbf{Q}}_{0}+\delta\hat{\cal
G}_{\varepsilon}^{R(A)}(\omega)(\mathbf{n}\mathbf{k})(\mathbf{n}{\mathbf{Q}}_{0}).$
(18)

Expressions for the components of the functions $\hat{\cal
G}_{\varepsilon}^{R(A)}(\omega)$ are listed in Appendix B. Components of
$\delta\hat{\cal G}_{\varepsilon}^{R(A)}(\omega)$ can be easily found:

$\delta{\cal G}_{\varepsilon}^{R(A)}(\omega)=\frac{v_{F}}{2\omega}{\cal
G}_{\varepsilon}^{R(A)}(\omega),\quad\delta{\cal
F}_{\varepsilon}^{R(A)}(\omega)=-\frac{\delta\Delta_{\omega}}{\zeta_{\varepsilon}^{R(A)}}.$
(19)

Note that there are contributions to $\delta{\cal
G}_{\varepsilon}^{R(A)}(\omega)$ which cancel out and this is the reason why
we have had to retain the normal component of
$\langle\mathbf{n}\mathbf{k}\hat{g}_{1}(\mathbf{n}\varepsilon;\mathbf{k}\omega)\rangle_{\mathbf{n}}$
in equation (14). Furthermore, as I will show below this correction plays a
crucial role in determining the magnitude of the IFE.

###  III.2 Linear analysis: Keldysh function

The Keldysh components of the quasiclassical function $\check{g}_{1}$ can be
computed in full analogy with the calculation of $\hat{g}_{1}^{R(A)}$. In
particular, since the equation for
$\langle\mathbf{n}\mathbf{k}\hat{g}_{1}^{K}\rangle_{\mathbf{n}}$ is identical
to (14) for the normal component of $\delta\hat{\cal
G}_{\varepsilon}^{K}(\omega)$ I have

$\delta{\cal G}^{K}(\varepsilon,\omega)=\frac{v_{F}}{2\omega}{\cal
G}^{K}(\varepsilon,\omega).$ (20)

We proceed with the equation for the Keldysh component of $\hat{\cal
G}_{\varepsilon}^{K}$ which reads:

$\begin{split}&\left(\tilde{\varepsilon}_{\omega}^{R}\hat{\tau}_{3}-i\tilde{\Delta}_{\omega}^{R}\hat{\tau}_{2}\right)\hat{\cal
G}^{K}-\hat{\cal
G}^{K}\left(\tilde{\varepsilon}_{-\omega}^{A}\hat{\tau}_{3}-i\tilde{\Delta}_{-\omega}^{A}\hat{\tau}_{2}\right)\\\
&=\tilde{a}_{\varepsilon}^{K}(\omega)\hat{\tau}_{0}-\tilde{b}_{\varepsilon}^{K}(\omega)\hat{\tau}_{1}.\end{split}$
(21)

This equation is, of course, quite similar to the equation (15) for the
retarded and advanced functions above. Functions appearing in the right hand
side of this equation are defined as:

$\begin{split}\tilde{a}_{\varepsilon}^{K}(\omega)=a_{\varepsilon}^{K}(\omega)+\frac{i}{2\tau}&\left({\cal
G}_{\varepsilon}^{R}g_{\varepsilon-\frac{\omega}{2}}^{K}-g_{\varepsilon+\frac{\omega}{2}}^{K}{\cal
G}_{\varepsilon}^{A}\right.\\\
&\left.+f_{\varepsilon+\frac{\omega}{2}}^{K}{\cal F}_{\varepsilon}^{A}-{\cal
F}_{\varepsilon}^{R}f_{\varepsilon-\frac{\omega}{2}}^{K}\right),\\\
\tilde{b}_{\varepsilon}^{K}(\omega)=b_{\varepsilon}^{K}(\omega)+\frac{i}{2\tau}&\left(f_{\varepsilon+\frac{\omega}{2}}^{K}{\cal
G}_{\varepsilon}^{A}-g_{\varepsilon+\frac{\omega}{2}}^{K}{\cal
F}_{\varepsilon}^{A}\right.\\\ &\left.+{\cal
G}_{\varepsilon}^{R}f_{\varepsilon-\frac{\omega}{2}}^{K}-{\cal
F}_{\varepsilon}^{R}g_{\varepsilon-\frac{\omega}{2}}^{K}\right),\end{split}$
(22)

where
$a_{\varepsilon}^{K}(\omega)=(g_{\varepsilon-\omega/2}^{K}-g_{\varepsilon+\omega/2}^{K})/\omega$
and
$b_{\varepsilon}^{K}(\omega)=(f_{\varepsilon-\omega/2}^{K}+f_{\varepsilon+\omega/2}^{K})/\omega$.
In the case of strong disorder when $\tau\Delta_{0}\ll 1$, we see that both
functions (22) are linearly dependent on $\tau^{-1}$. We should emphasize here
that no such terms $\propto\tau^{-1}$ will appear in the expression for the
Keldysh function if one chooses to work in the dirty limit from the outset by
using the Usadel equation Moore2017 ; Eremin2023 ; Yantao2023 .

The solution of this matrix equation can be found in the form $\hat{\cal
G}_{\varepsilon}^{K}(\omega)={\cal
G}_{\varepsilon}^{K}(\omega)\hat{\tau}_{3}-i{\cal
F}_{\varepsilon}^{K}(\omega)\hat{\tau}_{2}$ with

$\begin{split}{\cal
G}^{K}(\varepsilon,\omega)&=\frac{\tilde{a}_{\varepsilon}^{K}(\omega)(\tilde{\varepsilon}_{\omega}^{R}+\tilde{\varepsilon}_{-\omega}^{A})+\tilde{b}_{\varepsilon}^{K}(\omega)(\tilde{\Delta}_{\omega}^{R}-\tilde{\Delta}_{-\omega}^{A})}{[\tilde{\varepsilon}_{\omega}^{R}]^{2}-[\tilde{\varepsilon}_{-\omega}^{A}]^{2}-[\tilde{\Delta}_{\omega}^{R}]^{2}+[\tilde{\Delta}_{-\omega}^{A}]^{2}},\\\
{\cal
F}^{K}(\varepsilon,\omega)&=\frac{\tilde{b}_{\varepsilon}^{K}(\omega)(\tilde{\varepsilon}_{\omega}^{R}-\tilde{\varepsilon}_{-\omega}^{A})+\tilde{a}_{\varepsilon}^{K}(\omega)(\tilde{\Delta}_{\omega}^{R}+\tilde{\Delta}_{-\omega}^{A})}{[\tilde{\varepsilon}_{\omega}^{R}]^{2}-[\tilde{\varepsilon}_{-\omega}^{A}]^{2}-[\tilde{\Delta}_{\omega}^{R}]^{2}+[\tilde{\Delta}_{-\omega}^{A}]^{2}}.\end{split}$
(23)

We are providing the expression for ${\cal F}_{\varepsilon}^{K}(\omega)$ for
completeness since this function does not contribute to the IFE. We note that
temperature enters into these expressions only through the relation (7) and
the temperature dependent value of the pairing gap $\Delta_{0}(T)$.

###  III.3 Induced static magnetization

Having computed the linear-in-electric-field to the Keldysh Green’s function,
we go back to the hydrodynamic expression for the current density (11). Since
linear correction to both particle density and velocity will be determined by
the normal components of $\check{g}_{1}^{K}$, it will be convenient to
represent the normal component of $\hat{g}_{1}^{K}$ as
${g}_{1}^{K}(\mathbf{n}\varepsilon;\mathbf{r}t)={\cal
G}_{1}^{K}(\mathbf{n}\varepsilon;\mathbf{r}t)+\delta{\cal
G}_{1}^{K}(\mathbf{n}\varepsilon;\mathbf{r}t)$ with

$\begin{split}{\cal
G}_{1}^{K}(\mathbf{n}\varepsilon;\mathbf{r}t)&=iev_{F}{\cal
G}^{K}(\varepsilon,\omega)\mathbf{n}{\mathbf{E}}(\mathbf{r},t)\\\
&+iev_{F}{\cal
G}^{K}(\varepsilon,-\omega)\mathbf{n}{\mathbf{E}}^{*}(\mathbf{r},t),\\\
\delta{\cal
G}_{1}^{K}(\mathbf{n}\varepsilon;\mathbf{r}t)&=\frac{ev_{F}^{2}}{2\omega}{\cal
G}^{K}(\varepsilon,\omega)(\mathbf{n}{\mbox{\boldmath$\nabla$}})\mathbf{n}{\mathbf{E}}(\mathbf{r},t)\\\
&-\frac{ev_{F}^{2}}{2\omega}{\cal
G}^{K}(\varepsilon,-\omega)(\mathbf{n}{\mbox{\boldmath$\nabla$}})\mathbf{n}{\mathbf{E}}^{*}(\mathbf{r},t).\end{split}$
(24)

Here the electric field
${\mathbf{E}}(\mathbf{r},t)={\mathbf{E}}_{0}e^{i(\mathbf{k}\mathbf{r}-\omega
t)}$ and I used equation (20). After I insert these expressions into equation
for the currnt, there will be several contributions, but we single out two
contributions which contain the following combination of ${\mathbf{E}}$ and
${\mathbf{E}}^{*}$:

$\begin{split}{\mathbf{j}}_{\textrm{dc}}(\omega)&=ie^{2}v_{F}^{4}\gamma_{\omega}\langle\mathbf{n}(\mathbf{n}{\mathbf{E}})(\mathbf{n}{\mbox{\boldmath$\nabla$}})(\mathbf{n}{\mathbf{E}}^{*})\rangle_{\mathbf{n}}\\\
&-ie^{2}v_{F}^{4}\gamma_{\omega}\langle\mathbf{n}(\mathbf{n}{\mathbf{E}}^{*})(\mathbf{n}{\mbox{\boldmath$\nabla$}})(\mathbf{n}{\mathbf{E}}\rangle_{\mathbf{n}})\rangle_{\mathbf{n}}.\end{split}$
(25)

Alternatively, one can derive (25) by averaging the full expression for
current over a period of oscillations $2\pi/\omega$. Clearly, upon averaging
over the Fermi surface there is a contribution to the current
$\propto{\mbox{\boldmath$\nabla$}}\times({\mathbf{E}}\times{\mathbf{E}}^{*})$
which can be represented as
${\mbox{\boldmath$\nabla$}}\times{\mathbf{M}}_{\textrm{ind}}$. Function
$\gamma_{\omega}$ in expression (25) is formally given by the following
expression

$\begin{split}\gamma_{\omega}&=\frac{\nu_{F}^{2}}{8\omega}\int\limits_{-\infty}^{\infty}{\cal
G}^{K}(\varepsilon,\omega)d\varepsilon\int\limits_{-\infty}^{\infty}{\cal
G}^{K}(\varepsilon^{\prime},-\omega)d\varepsilon^{\prime}\equiv\frac{I(\omega,\tau)}{\omega^{3}}.\end{split}$
(26)

Here on the last step we took into account that coefficients
$a_{\varepsilon}(\omega)$ and $b_{\varepsilon}(\omega)$ are proportional to
$1/\omega$ and, as a consequence, ${\cal G}^{K}(\varepsilon,\omega)\propto
1/\omega$, so that function $\gamma_{\omega}\propto 1/\omega^{3}$, while
function $I(\omega,\tau)$ encodes the effects of disorder on the induced
magnetization.

Figure 2: Dependence of the function $I(\omega,\tau)$, Eq. (26), on frequency
of external electric field for various values of the disorder scattering rate
$\tau^{-1}$. Since function $I(\omega,\tau)$ determines the magnitude of the
effect, from these results we conclude that disorder suppresses the value of
the induced static magnetization. It is worthwhile to note that the
suppression is strongest at frequencies $\omega_{\textrm{min}}\approx
2\Delta_{0}$, which coincides with the resonant frequency of the amplitude
Higgs mode. The value of the temperature has been set to $T=0.1\Delta_{0}$.
Figure 3: Dependence of the function $I(\omega,\tau)$, Eq. (26), on disorder
scattering rate for various values of frequency of an external electric field
$\omega$. Note that suppression is much stronger at the resonant frequency of
the amplitude Higgs mode $\omega=2\Delta_{0}$. The value of the temperature
has been set to $T=0.1\Delta_{0}$.

Function $I(\omega,\tau)$ appearing in (26) can be evaluated numerically for
different values of the disorder scattering rate and temperature. In the dirty
limit $\tau\Delta_{0}\ll 1$ this function will become weakly dependent on the
disorder scattering rate. In order to see this, we first note that in this
limit both numerator and denominator in the expression for ${\cal
G}_{\varepsilon}^{K}$ will be varying linearly with $\tau^{-1}$. This is
because the terms which are proportional to $\tau^{-2}$ in denominator cancel
out due to normalization condition (see Appendix B). In passing I note that
this fact - the dependence of the Keldysh Green function on disorder
scattering rate - is yet another manifestation of the observation made earlier
that Anderson theorem does not apply to disordered nonequilibrium
superconductors, i.e. and scattering on potential impurities does affect
various out-of-equilibrium response functions, which conceptually is in
agreement with the earlier results Silaev2019-Disorder ; Yang2022-Disorder ;
Dzero2023-Disorder .

In Fig. 2 I show the result of the calculation of the frequency dependence of
function $I(\omega,\tau)$ for different values of $1/\tau\Delta_{0}$. As we
can see this dependence is nonmonotonic. However, this function clearly
acquires smaller values with an increase in $\tau^{-1}$. Note that
$I(\omega,\tau)$ has a minimum at frequencies $\simeq 2\Delta_{0}$ which
matches the value of the resonant frequency for the excitation of the
amplitude Higgs mode. Notably, the minimum moves to higher frequencies with
increase in the values of $\tau^{-1}$, which seems to be in agreement with the
results of Ref. Dzero2023-Disorder . Since such a shift is not expected to
appear in the strongly disordered superconductors Moore2017 ; Yantao2023 , we
attribute this effect to the perturbative nature of the self-consistent Born
approximation.

In Fig. 3 the dependence of $I(\omega,\tau)$ on scattering rate $\tau^{-1}$ is
shown. This dependence turns out to be perfectly monotonic and reflects the
suppression of ${\mathbf{M}}_{\textrm{ind}}$ for four different values of
$\omega$. Notably, the suppressions is by far the strongest when
$\omega=2\Delta_{0}$, which is consistent with the results presented in Fig.
2. I am inclined to interpret this suppression as being caused by the resonant
excitation of the amplitude Higgs mode.

##  IV Conclusions

In this paper using the quasiclassical theory of superconductivity I have
investigated how potential impurities may affect the emergence of the inverse
Faraday effect in conventional superconductors. Scattering on potential
impurities have been taken into account within the self-consistent Born
approximation. The normal contribution to the current must still be determined
by the diffusion coefficient $D\propto\tau$, while the superfluid component of
the current is suppressed by the pair breaking processes due to non-
equilibrium Dzero2023-Disorder . In agreement these qualitative observations,
I have found that in the fairly broad range of frequencies potential
impurities suppress the magnitude of the IFE. In addition, it has been found
that the static magnetization ${\mathbf{M}}_{\textrm{ind}}(\omega)$ has a
minimum for $\omega\approx 2\Delta_{0}$. This frequency also corresponds to
the minimal excitation frequency for the amplitude Higgs mode. This
observation provides another example of manifestation of the amplitude Higgs
mode in response functions of a superconductor driven out-of-equilibrium
Manske2014 ; Aoki2015 ; Kemper2015 ; Cea2016 ; Moore2017 ; Eremin2023 ;
Yantao2023 . We certainly hope that our theoretical prediction of the minimum
in the dc magnetization as a function of frequency can be verified
experimentally. Performing such a measurement could provide yet another
important insight into the physics of the collective excitations and their
impact on the IFE in conventional superconductors.

##  V Acknowledgments

I would like to thank A. V. Balatsky for bringing this problem to my attention
and useful discussion. I also would like to thank Yantao Li for bringing
several important references related to the main subject of this work to my
attention. This work was financially supported by the National Science
Foundation grant NSF-DMR-2002795.

##  Appendix A Auxiliary expressions

In this Appendix I provide some auxiliary expressions which are used to derive
the quasiclassical equation (15) in the main text.

#### Time derivatives.

First simplification arises due to our choice for the time dependence of the
electric field:

$\begin{split}&\left[\partial_{\varepsilon}\left(\check{g}_{\varepsilon}\check{\tau}_{3}\right){\cal
E}\left(\frac{i}{2}\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{\varepsilon}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{t}\right)+\partial_{\varepsilon}\left(\check{\tau}_{3}\check{g}_{\varepsilon}\right){\cal
E}\left(-\frac{i}{2}\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{\varepsilon}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{t}\right)\right]e^{-i\omega
t}\\\ &=\frac{e^{-i\omega t}}{2}\left[{\cal
E}\left(-\frac{\omega}{2}\partial_{\varepsilon}\right)+{\cal
E}\left(\frac{\omega}{2}\partial_{\varepsilon}\right)\right]\left\\{\check{\tau}_{3},\partial_{\varepsilon}\check{g}_{\varepsilon}\right\\}\\\
&+\frac{e^{-i\omega t}}{2}\left[{\cal
E}\left(-\frac{\omega}{2}\partial_{\varepsilon}\right)-{\cal
E}\left(\frac{\omega}{2}\partial_{\varepsilon}\right)\right]\left[\check{\tau}_{3},\partial_{\varepsilon}\check{g}_{\varepsilon}\right].\end{split}$
(27)

Here I use the bare propagators since this terms always enter in the
combination with external electric field, which is assumed to be small. This
expression can be further simplified if we take into account the definition of
function ${\cal E}$. For example, for the combination in the second line in
(27) it obtains:

$\begin{split}&\left[{\cal
E}\left(-\frac{\omega}{2}\partial_{\varepsilon}\right)+{\cal
E}\left(\frac{\omega}{2}\partial_{\varepsilon}\right)\right]\partial_{\varepsilon}\check{g}_{\varepsilon}=\int\limits_{-1}^{1}e^{\eta(\omega/2)\partial_{\varepsilon}}\partial_{\varepsilon}\check{g}_{\varepsilon}d\eta\\\
&=\partial_{\varepsilon}\int\limits_{-1}^{1}\check{g}_{\varepsilon+\eta\frac{\omega}{2}}d\eta=\frac{2}{\omega}\left(\check{g}_{\varepsilon+\omega/2}-\check{g}_{\varepsilon-\omega/2}\right),\\\
\end{split}$ (28)

Similarly, for the expression in the third line in (27) I find

$\begin{split}&\left[{\cal
E}\left(\frac{\omega}{2}\partial_{\varepsilon}\right)-{\cal
E}\left(-\frac{\omega}{2}\partial_{\varepsilon}\right)\right]\partial_{\varepsilon}\check{g}_{\varepsilon}=\partial_{\varepsilon}\left(\int\limits_{0}^{1}-\int\limits_{-1}^{0}\right)\check{g}_{\varepsilon+\eta\omega/2}d\eta\\\
&=\frac{2}{\omega}\left(\check{g}_{\varepsilon+\omega/2}+\check{g}_{\varepsilon-\omega/2}-2\check{g}_{\varepsilon}\right).\end{split}$
(29)

Lastly, we note that the last term in the right hand side of this expression
is canceled out by the term which contains the vector potential in the main
equation, Eq. (3). Thus, collecting all these contributions together and
taking into account that the first order correction to the superconducting
order parameter vanishes after the averaging over the Fermi surface yields
equation (13) in the main text.

#### Self energy part.

For the commutator which involves the self energy part
$\check{\sigma}=-\frac{i}{2\tau}\langle\check{g}\rangle_{\mathbf{n}}$ after
retaining the linear order terms in powers of electric field I have

$\begin{split}&[\check{\sigma}_{\textrm{u}}\stackrel{{\scriptstyle\circ}}{{,}}\check{g}_{\varepsilon}+\check{g}_{1}]\approx-\frac{i}{2\tau}[\check{g}_{\varepsilon}\stackrel{{\scriptstyle\circ}}{{,}}\check{g}_{1}]-\frac{i}{2\tau}[\langle\check{g}_{1}\rangle_{\mathbf{n}}\stackrel{{\scriptstyle\circ}}{{,}}\check{g}_{\varepsilon}]\\\
&=\frac{i}{2\tau}\check{g}_{1}e^{-\frac{i}{2}\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{t}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{\varepsilon}}\check{g}_{\varepsilon}-\frac{i}{2\tau}\check{g}_{\varepsilon}e^{\frac{i}{2}\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{\varepsilon}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{t}}\check{g}_{1}+\frac{i}{2\tau}\check{g}_{\varepsilon}e^{\frac{i}{2}\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{\varepsilon}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{t}}\langle\check{g}_{1}\rangle_{\mathbf{n}}\\\
&-\frac{i}{2\tau}\langle\check{g}_{1}\rangle_{\mathbf{n}}e^{-\frac{i}{2}\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{t}\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{\varepsilon}}\check{g}_{\varepsilon}=\frac{i}{2\tau}\left(\check{g}_{1}-\langle\check{g}_{1}\rangle_{\mathbf{n}}\right)\check{g}_{\varepsilon-\frac{\omega}{2}}\\\
&-\frac{i}{2\tau}\check{g}_{\varepsilon+\frac{\omega}{2}}\left(\check{g}_{1}-\langle\check{g}_{1}\rangle_{\mathbf{n}}\right).\end{split}$
(30)

Note that the averages over the Fermi surface are nonzero only for the
contributions to $\check{g}_{1}$ which are directly proportional to the dot
product $\mathbf{n}\mathbf{k}$ (see e.g. (18) in the main text), which means
that averaging both sides of this equation over the Fermi surface will give
zero in this case. This fact is used in the derivation of the equation (14) in
the main text.

##  Appendix B Expressions for the linear corrections to the retarded and
advanced propagators

Components of the matrix function $\hat{\cal G}_{\varepsilon}(\omega)$, Eq.
(17), must satisfy the following system of linear equations

$\begin{split}&(\tilde{\varepsilon}_{\omega}-\tilde{\varepsilon}_{-\omega}){\cal
G}_{\varepsilon}(\omega)-(\tilde{\Delta}_{\omega}-\tilde{\Delta}_{-\omega}){\cal
F}_{\varepsilon}(\omega)=a_{\varepsilon}(\omega),\\\
&(\tilde{\Delta}_{\omega}+\tilde{\Delta}_{-\omega}){\cal
G}_{\varepsilon}(\omega)-(\tilde{\varepsilon}_{\omega}+\tilde{\varepsilon}_{-\omega}){\cal
F}_{\varepsilon}(\omega)=-b_{\varepsilon}(\omega).\end{split}$ (31)

Here I have omitted the $R(A)$ superscripts for brevity. The solution of these
equations is

$\begin{split}{\cal
G}_{\varepsilon}(\omega)&=\frac{a_{\varepsilon}(\omega)(\tilde{\varepsilon}_{\omega}+\tilde{\varepsilon}_{-\omega})+b_{\varepsilon}(\omega)(\tilde{\Delta}_{\omega}-\tilde{\Delta}_{-\omega})}{\tilde{\varepsilon}_{\omega}^{2}-\tilde{\varepsilon}_{-\omega}^{2}-\tilde{\Delta}_{\omega}^{2}+\tilde{\Delta}_{-\omega}^{2}},\\\
{\cal
F}_{\varepsilon}(\omega)&=\frac{b_{\varepsilon}(\omega)(\tilde{\varepsilon}_{\omega}-\tilde{\varepsilon}_{-\omega})+a_{\varepsilon}(\omega)(\tilde{\Delta}_{\omega}+\tilde{\Delta}_{-\omega})}{\tilde{\varepsilon}_{\omega}^{2}-\tilde{\varepsilon}_{-\omega}^{2}-\tilde{\Delta}_{\omega}^{2}+\tilde{\Delta}_{-\omega}^{2}}.\end{split}$
(32)

Note that in the dirty limit $\tau\Delta_{0}\ll 1$, the disorder scattering
rate drops out from these expressions. Specifically, in the limit when
$\tau^{-1}\gg\textrm{max}\\{\omega,\varepsilon\\}$ I have

$\begin{split}&\tilde{\varepsilon}_{\omega}^{2}-\tilde{\varepsilon}_{-\omega}^{2}-\tilde{\Delta}_{\omega}^{2}+\tilde{\Delta}_{-\omega}^{2}\approx\frac{i\omega}{\tau}(g_{\varepsilon+\frac{\omega}{2}}+g_{\varepsilon-\frac{\omega}{2}})\\\
&+\frac{i\varepsilon}{\tau}(g_{\varepsilon+\frac{\omega}{2}}-g_{\varepsilon-\frac{\omega}{2}})-\frac{i\Delta_{0}}{\tau}(f_{\varepsilon+\frac{\omega}{2}}-f_{\varepsilon-\frac{\omega}{2}})\end{split}$
(33)

and the terms $O(1/\tau^{2})$ cancel out due to normalization condition
$g_{\varepsilon}^{2}-f_{\varepsilon}^{2}=1$. Since the numerator in Eq. (32)
is also linear in $1/\tau$ the final result becomes independent of
$\tau^{-1}$. This result is actually consistent both with the Anderson theorem
and with the ones from non-linear $\sigma$ model where the exact averaging
over disorder is performed Kamenev2009 ; Kamenev2011  and then nonlinear
response is studied using the Usadel equation Eremin2023 ; Yantao2023 .
However, as it is shown in the main text, the Keldysh components retain their
dependence on $\tau^{-1}$ for these functions contain all the information
above non-equilibrium state of a superconductor.

## References

  * (1) J. A. Fülöp, L. Pálfalvi, G. Almási, and J. Hebling, “Design of high-energy terahertz sources based on optical rectification,” Opt. Express, vol. 18, pp. 12311–12327, Jun 2010.
  * (2) D. N. Basov, R. D. Averitt, D. van der Marel, M. Dressel, and K. Haule, “Electrodynamics of correlated electron materials,” Rev. Mod. Phys., vol. 83, pp. 471–541, Jun 2011.
  * (3) R. Matsunaga and R. Shimano, “Nonequilibrium BCS state dynamics induced by intense terahertz pulses in a superconducting NbN film,” Phys. Rev. Lett., vol. 109, p. 187002, 2012.
  * (4) R. Matsunaga, Y. I. Hamada, K. Makise, Y. Uzawa, H. Terai, Z. Wang, and R. Shimano, “Higgs amplitude mode in the BCS superconductors ${\mathrm{Nb}}_{1\mathrm{\text{-}}x}{\mathrm{Ti}}_{x}\mathrm{N}$ induced by terahertz pulse excitation,” Phys. Rev. Lett., vol. 111, p. 057002, Jul 2013.
  * (5) R. Matsunaga, N. Tsuji, H. Fujita, A. Sugioka, K. Makise, Y. Uzawa, H. Terai, Z. Wang, H. Aoki, and R. Shimano, “Light-induced collective pseudospin precession resonating with Higgs mode in a superconductor,” Science, vol. 345, no. 6201, pp. 1145–1149, 2014.
  * (6) M. Beck, I. Rousseau, M. Klammer, P. Leiderer, M. Mittendorff, S. Winnerl, M. Helm, G. N. Gol’tsman, and J. Demsar, “Transient increase of the energy gap of superconducting NbN thin films excited by resonant narrow-band terahertz pulses,” Phys. Rev. Lett., vol. 110, p. 267003, Jun 2013.
  * (7) Z. Sun, M. M. Fogler, D. N. Basov, and A. J. Millis, “Collective modes and terahertz near-field response of superconductors,” Phys. Rev. Res., vol. 2, p. 023413, Jun 2020.
  * (8) Y. Lubashevsky, L. Pan, T. Kirzhner, G. Koren, and N. P. Armitage, “Optical birefringence and dichroism of cuprate superconductors in the THz regime,” Phys. Rev. Lett., vol. 112, p. 147001, Apr 2014.
  * (9) K. Katsumi, J. Fiore, M. Udina, R. R. I. au2, D. Barbalas, J. Jesudasan, P. Raychaudhuri, G. Seibold, L. Benfatto, and N. P. Armitage, “Revealing novel aspects of light-matter coupling in terahertz two-dimensional coherent spectroscopy: the case of the amplitude mode in superconductors,” 2023.
  * (10) D. Pekker and C. Varma, “Amplitude/Higgs modes in condensed matter physics,” Annual Review of Condensed Matter Physics, vol. 6, no. 1, pp. 269–297, 2015\.
  * (11) T. Cea, C. Castellani, and L. Benfatto, “Nonlinear optical effects and third-harmonic generation in superconductors: Cooper pairs versus Higgs mode contribution,” Phys. Rev. B, vol. 93, p. 180507, May 2016.
  * (12) R. Shimano and N. Tsuji, “Higgs mode in superconductors,” Annual Review of Condensed Matter Physics, vol. 11, no. 1, pp. 103–124, 2020.
  * (13) R. Haenel, P. Froese, D. Manske, and L. Schwarz, “Time-resolved optical conductivity and Higgs oscillations in two-band dirty superconductors,” Phys. Rev. B, vol. 104, p. 134504, Oct 2021.
  * (14) D. Sherman, U. S. Pracht, B. Gorshunov, S. Poran, J. Jesudasan, M. Chand, P. Raychaudhuri, M. Swanson, N. Trivedi, A. Auerbach, M. Scheffler, A. Frydman, and M. Dressel, “The Higgs mode in disordered superconductors close to a quantum phase transition,” Nature Physics, vol. 11, no. 2, pp. 188–192, 2015.
  * (15) M. Bellitti, C. R. Laumann, and B. Z. Spivak, “Incoherent excitation of coherent Higgs oscillations in superconductors,” Phys. Rev. B, vol. 105, p. 104513, Mar 2022.
  * (16) G. Eliashberg, “Film superconductivity stimulated by a high-frequency field,” Sov. Phys. - JETP Lett., vol. 11, p. 114, 1970.
  * (17) P. Derendorf, A. F. Volkov, and I. Eremin, “Nonlinear response of diffusive superconductors to ac-electromagnetic fields,” arXiv:2308.00838, 2023.
  * (18) X.-J. Chen, “Fundamental mechanism for all-optical helicity-dependent switching of magnetization,” Scientific Reports, vol. 7, no. 1, p. 41294, 2017.
  * (19) M. D. Croitoru, S. V. Mironov, B. Lounis, and A. I. Buzdin, “Toward the light-operated superconducting devices: Circularly polarized radiation manipulates the current-carrying states in superconducting rings,” Advanced Quantum Technologies, vol. 5, no. 10, p. 2200054, 2022.
  * (20) M. Smith, A. V. Andreev, and B. Z. Spivak, “Debye mechanism of giant microwave absorption in superconductors,” Phys. Rev. B, vol. 101, p. 134508, Apr 2020\.
  * (21) L. P. Pitaevskii, “Electric forces in a transparent dispersive medium,” Sov. Phys. - JETP, vol. 12, p. 1008, 1961.
  * (22) P. S. Pershan, J. P. van der Ziel, and L. D. Malmstrom, “Theoretical discussion of the inverse Faraday effect, Raman scattering, and related phenomena,” Phys. Rev., vol. 143, pp. 574–583, Mar 1966.
  * (23) M. Battiato, G. Barbalinardo, and P. M. Oppeneer, “Quantum theory of the inverse Faraday effect,” Phys. Rev. B, vol. 89, p. 014413, Jan 2014.
  * (24) X. Yang, Y. Mou, H. Zapata, B. Reynier, B. Gallas, and M. Mivelle, “An inverse Faraday effect through linear polarized light,” 2022.
  * (25) Y. Mou, X. Yang, B. Gallas, and M. Mivelle, “A reversed inverse Faraday effect,” 2023.
  * (26) Y. Mou, X. Yang, B. Gallas, and M. Mivelle, “A chiral inverse Faraday effect mediated by an inversely designed plasmonic antenna,” 2023.
  * (27) S. Banerjee, U. Kumar, and S.-Z. Lin, “Inverse Faraday effect in Mott insulators,” Phys. Rev. B, vol. 105, p. L180414, May 2022.
  * (28) Y. Gao, C. Wang, and D. Xiao, “Topological inverse Faraday effect in Weyl semimetals,” 2020.
  * (29) S. Parchenko, K. Hofhuis, A. Ciuciulkaite, V. Kapaklis, V. Scagnoli, L. Heyderman, and A. Kleibert, “Plasmon-enhanced optical control of magnetism at the nanoscale via the inverse Faraday effect,” 2023.
  * (30) J. W. Han, P. Sai, D. B. But, E. Uykur, S. Winnerl, G. Kumar, M. L. Chin, R. L. Myers-Ward, M. T. Dejarld, K. M. Daniels, T. E. Murphy, W. Knap, and M. Mittendorff, “Strong transient magnetic fields induced by THz-driven plasmons in graphene disks,” Nature Communications, vol. 14, Nov. 2023\.
  * (31) A. H. Majedi, “Microwave-induced inverse Faraday effect in superconductors,” Phys. Rev. Lett., vol. 127, p. 187001, 2021.
  * (32) S. V. Mironov, A. S. Mel’nikov, I. D. Tokman, V. Vadimov, B. Lounis, and A. I. Buzdin, “Inverse Faraday effect for superconducting condensates,” Phys. Rev. Lett., vol. 126, p. 137002, Apr 2021.
  * (33) A. V. Parafilo, M. V. Boev, V. M. Kovalev, and I. G. Savenko, “Photogalvanic transport in fluctuating Ising superconductors,” Phys. Rev. B, vol. 106, p. 144502, Oct 2022.
  * (34) A. V. Putilov, S. V. Mironov, A. S. Mel’nikov, and A. A. Bespalov, “Inverse Faraday effect in superconductors with a finite gap in the excitation spectrum,” JETP Letters, vol. 117, no. 11, pp. 827–833, 2023.
  * (35) V. Plastovets and A. Buzdin, “Fluctuation-mediated inverse Faraday effect in superconducting rings,” Physics Letters A, vol. 481, p. 129001, 2023.
  * (36) M. D. Croitoru, B. Lounis, and A. I. Buzdin, “Helicity-controlled switching of superconducting states by radiation pulse,” Applied Physics Letters, vol. 123, p. 122601, 09 2023.
  * (37) P. Sharma and A. V. Balatsky, “Light induced magnetism in metals via inverse Faraday effect,” 2023.
  * (38) P. W. Anderson, “Knight shift in superconductors,” Phys. Rev. Lett., vol. 3, pp. 325–326, Oct 1959.
  * (39) A. A. Abrikosov and L. P. Gor’kov, “Contribution to the theory of superconducting alloys with paramagnetic impurities,” Sov. Phys. - JETP, vol. 12, p. 1243, 1961.
  * (40) A. V. Balatsky, I. Vekhter, and J.-X. Zhu, “Impurity-induced states in conventional and unconventional superconductors,” Rev. Mod. Phys., vol. 78, pp. 373–433, May 2006.
  * (41) M. Silaev, “Nonlinear electromagnetic response and Higgs-mode excitation in BCS superconductors with impurities,” Phys. Rev. B, vol. 99, p. 224511, Jun 2019.
  * (42) F. Yang and M. W. Wu, “Impurity scattering in superconductors revisited: Diagrammatic formulation of the supercurrent-supercurrent correlation and Higgs-mode damping,” Phys. Rev. B, vol. 106, p. 144509, Oct 2022.
  * (43) F. Yang and M. W. Wu, “Influence of scattering on the optical response of superconductors,” Phys. Rev. B, vol. 102, p. 144508, Oct 2020.
  * (44) M. Dzero, “Collisionless dynamics in disordered superconductors,” arXiv:2303.06750, 2023.
  * (45) A. I. Larkin and Y. Ovchinnikov, “Quasiclassical method in the theory of superconductivity,” Sov. Phys. - JETP, vol. 28, p. 1200, 1969.
  * (46) K. D. Usadel, “Generalized diffusion equation for superconducting alloys,” Phys. Rev. Lett., vol. 25, pp. 507–509, Aug 1970.
  * (47) T. Kita, “Gauge invariance and Hall terms in the quasiclassical equations of superconductivity,” Phys. Rev. B, vol. 64, p. 054503, Jun 2001.
  * (48) N. B. Kopnin, “Kinetic equations for clean superconductors: Application to the flux flow Hall effect,” Journal of Low Temperature Physics, vol. 97, no. 1, pp. 157–179, 1994.
  * (49) A. Kamenev and A. Levchenko, “Keldysh technique and non-linear $\sigma$-model: basic principles and applications,” Advances in Physics, vol. 58, no. 3, pp. 197–319, 2009.
  * (50) A. Kamenev, Field Theory of Non-Equilibrium Systems.  Cambridge University Press, 2011.
  * (51) T. Kita, “Introduction to Nonequilibrium Statistical Mechanics with Quantum Field Theory,” Progress of Theoretical Physics, vol. 123, pp. 581–658, 04 2010.
  * (52) Y. Li and M. Dzero, “Amplitude Higgs mode in superconductors with magnetic impurities,” arXiv:2311.09310, 2023.
  * (53) A. Moor, A. F. Volkov, and K. B. Efetov, “Amplitude Higgs mode and admittance in superconductors with a moving condensate,” Phys. Rev. Lett., vol. 118, p. 047001, Jan 2017.
  * (54) H. Krull, D. Manske, G. S. Uhrig, and A. P. Schnyder, “Signatures of nonadiabatic BCS state dynamics in pump-probe conductivity,” Phys. Rev. B, vol. 90, p. 014515, Jul 2014.
  * (55) N. Tsuji and H. Aoki, “Theory of Anderson pseudospin resonance with Higgs mode in superconductors,” Phys. Rev. B, vol. 92, p. 064508, Aug 2015.
  * (56) A. F. Kemper, M. A. Sentef, B. Moritz, J. K. Freericks, and T. P. Devereaux, “Direct observation of Higgs mode oscillations in the pump-probe photoemission spectra of electron-phonon mediated superconductors,” Phys. Rev. B, vol. 92, p. 224517, Dec 2015.