# Collisionless dynamics of superconducting gap excited by spin-splitting
field

V. Plastovets  University of Bordeaux, LOMA UMR-CNRS 5798, F-33405 Talence
Cedex, France  Institute for Physics of Microstructures, Russian Academy of
Sciences, 603950 Nizhny Novgorod, GSP-105, Russia  A. S. Mel’nikov  Institute
for Physics of Microstructures, Russian Academy of Sciences, 603950 Nizhny
Novgorod, GSP-105, Russia  Moscow Institute of Physics and Technology, 141701
Dolgoprudny, Russia  Lobachevsky State University of Nizhny Novgorod, 603950
Nizhni Novgorod, Russia  A. I. Buzdin  University of Bordeaux, LOMA UMR-CNRS
5798, F-33405 Talence Cedex, France

###### Abstract

We study the coherent dynamic interaction of a time-dependent spin-splitting
field with the homogeneous superconducting order parameter $\Delta(t)$
mediated by spin-orbit coupling using the time-dependent Bogoliubov-de Gennes
theory. In the first part of the work we show that linear response of the
superconductor is strongly affected by the Zeeman field and spin-flip
processes, giving rise to multiple resonant frequencies of the superconducting
Higgs modes. These modes can be excited either by a quench, or by an
additional non-stationary component of the spin-splitting field, which couples
linearly to the Higgs modes. In the second part, we analyze the nonadiabatic
dynamics of quasiparticle states arising from the intersection of spectral
branches from different spin subbands, which can be provoked by a linearly
growing Zeeman field. We provide insights into the dependence of the order
parameter $\Delta(t)$ on this field and interference effects caused by
tunneling of states at the avoided crossing points. We also show that since
the nonadiabatic tunneling is related to spin-flip processes, the
quasiparticle gas experiences a dynamic magnetization that contributes to its
spin susceptibility.

##  I Introduction

Extensive studies of nonequilibrium states of superconductors [1, 2] pay
considerable attention to the so-called collisionless dynamics of a
superconducting condensate, described by the complex-valued pairing potential
$\Delta(t)$. At timescales shorter than the typical inelastic relaxation time
$t\ll\tau_{\varepsilon}$ the dynamics of Cooper pairs is in coherent regime
and is described by the Keldysh technique for Green’s functions or its
quasiclassical approximation [3, 4]. The collisionless regime manifests itself
most clearly in the existence of oscillations of the amplitude of the order
parameter $\Delta(t)=\Delta_{0}+\delta\Delta(t)$ near the equilibrium gap
value $\Delta_{0}$. This mode comes from excited interference interaction
between the wavefunctions of the quasiparticles (QP) from broken Cooper
pairs.Due to the QP dispersion the summation over all interference
contributions results in an inhomogeneous broadening of the total gap mode,
which is equivalent to a weak damping with a typical time evolution
$\delta\Delta(t)\propto\cos(2\Delta_{0}t)/\sqrt{\Delta_{0}t}$ [3]. By analogy
with electro-weak particle theory this amplitude mode is called Higgs mode [5,
6]. Since the Higgs mode is a scalar excitation, it can not be coupled to the
electromagnetic field ${\bf A}(t)$ linearly and several indirect mechanisms
have been studied, such as linear excitation by the THz radiation in the
presence of dc supercurrent [7, 8]. Also it possible to realize a nonlinear
and coherent (or incoherent [9]) Higgs mode excitation using high-intensity
THz light with frequency just above the equilibrium superconducting gap
$\Delta_{0}$, which can be detected by ultrafast pump-probe spectroscopy and
third harmonic generation measurements [10, 11, 12, 13, 14].

It is known that, in addition to electromagnetic fields, superconductors also
respond to nonstationary spin-splitting fields ${\bf h}(t)$. Typically, this
field is produced by an external magnetic field ${\bf h}=\mu_{B}{\bf H}$ or by
the exchange field of an adjacent ferromagnetic layer ${\bf h}\sim
J_{\text{ex}}{\bf m}$, which is induced by proximity to the superconductor.
Spin-split systems serve as a good platforms for spintronic applications and
extensive study of various non-equilibrium processes has been done over the
last few decades [15, 16, 17]. In particular, by inducing magnetic moment
dynamics in S/F junctions an effective spin-triplet component of the
superconducting gap is generated resulting in long-range proximity effects
[18, 19, 20, 21]. On the other hand experimental observations indicate that
the superconducting subsystem has a direct impact on the ferromagnetic
resonance in hybrid S/F structures [22, 23, 24].

In recent years there has been a growing interest in studying of the Higgs
modes in the proximized superconducting systems [25, 26] as well as the
interaction of collective modes in S/F systems. For instance, it was recently
shown that in a superconductor in the helical phase, which can be achieved in
the presence of a strong spin-orbit coupling (SOC) and an exchange field, the
Higgs mode can be linearly coupled to the electromagnetic field through the
nonzero superconducting phase gradient in the ground state [27]. Also, it was
revealed that the coupling of the Higgs mode $\delta\Delta(t)$ in a
superconductor to external light ${\bf A}(t)$ and magnetic dynamics ${\bf
m}(t)$ in the F layer allows the generation of time-dependent spin currents
[28]. These currents can themselves excite the Higgs mode in the
superconductor through the resonance of the ferromagnet due to the reciprocal
effect [28]. Another example is an interplay between the superconducting Higgs
mode and a magnon mode in the adjacent F layer in the presence of a SOC and
static proximity effect [29]. Interestingly, the Higgs mode here is coupled to
the Zeeman field ${\bf h}(t)$ linearly due to the presence of both the spin-
orbit interaction and some preferred direction given by wave-vector of the
magnetic mode.

According to the aforementioned works, the SOC is critical for interaction of
different spin subbands of the QP spectrum, which directly leads to the gap
dynamics $\Delta(t)$. Some preconditions for this can be taken from the
elementary analysis of the equilibrium state. The equilibrium superconducting
gap does not depend on the Zeeman field below the so-called paramagnetic limit
$h_{\text{cr}}=\Delta_{0}/\sqrt{2}$, so that
$\Delta(h<h_{\text{cr}},T=0)=\Delta_{0}$; and above this limit the
superconductivity is completely suppressed with
$\Delta(h>h_{\text{cr}},T=0)=0$ [30, 31]. The SOC drastically changes the
dependence $\Delta(h)$ and promotes a generation of triplet component
superconducting correlations, leading to the survival of the gap at
$h>h_{\text{cr}}$ [32]. This effect is associated with mixing of the different
spin states of QP, and the appearance of such mixing is naturally expected in
the dynamic regime.

The general description of dynamics of the superconducting condensate in the
presence of both time-dependent spin-flipping field and SOC is rather
difficult problem [28, 29]. In this paper we will focus on a fairly simple and
specific system configuration, which allows us to explicitly trace the
temporal evolution of QP states and its contribution to the order parameter
$\Delta(t)$. For the sake of simplicity we consider an uniform superconductor
at zero temperature $T=0$ and consider short timescale
$t\ll\tau_{\varepsilon}$ at which the collisionless regime holds, so one can
treat the system with the pure quantum-mechanical approach within the time-
dependent Bogoliubov-de Gennes (TDBdG) equations [33]. We confine ourselves to
addressing the homogeneous spin-splitting field with only one component ${\bf
h}(t)=h(t){\bf z}_{0}$. A simple approach based on the expansion of the QP
wave function in terms of the eigenstates of the BdG Hamiltonian
$\psi(t)=\sum_{n}C_{n}(t)\Psi_{n}$ can be developed. The behavior of the QPs
and related self-consistent gap function $\Delta(t)$ are determined by the
coefficients $C_{n}(t)$, which describe how the states with a specific spin
quantum number and momentum are refilled due to nonstationary transitions.
After introducing the TDBdG equations in Sec. II, we examine two different
regimes of coherent evolution of the order parameter.

In Sec. III we analyze linearized gap dynamics where the temporal evolution of
the Higgs modes $\delta\Delta(t)$ is traced in the presence of a stationary
spin-splitting field $h_{0}$. Since the SOC allows the transitions between the
QP states with different spins, the induced perturbation of the gap
$\delta\Delta(t)$ acquires three eigenfrequencies including the standard
$2\Delta_{0}$ and two additional frequencies $2(\Delta_{0}\pm h_{0})$. These
modes define both the free oscillation of the perturbed gap at $h_{0}^{-1}\ll
t\ll\tau_{\varepsilon}$, and resonant peaks in the case of driven
oscillations. It was also shown that in the specific configuration, the linear
coupling of the Higgs mode and the perturbation of the Zeeman field
$\delta{\bf h}(t)$ co-directional with ${\bf h}_{0}$ is possible. Note that
the frequencies shifted by the spin-splitting field have been observed in the
numerical simulation of the dynamics of the one-dimensional Fermi superfluid
exposed to the nonstationary Zeeman field and strong SOC in the Ref. [34].

In Sec. IV we consider the dynamics of the gap $\Delta(t)$ driven by linearly
growing field ${\bf h}(t)$. At some point the field becomes larger then the
equilibrium gap value $|{\bf h}(t)|>\Delta_{0}$ and thus provokes the crossing
of the branches from different spin subbands of the QP spectrum. The
appearance of non-adiabatic transitions between the states at the intersection
point is equivalent to the dynamical spin-flip process and can be described
with the Landau-Zener-Stückelberg-Majorana (LZSM) tunneling problem [35].
Corresponding refilling of the amplitudes $C_{n}(t)$ contributes to the gap
function $\Delta(t)$ and can drastically change its behavior depending on the
field growth rate. In Sec. IV.3-IV.4 we derive an analytical expression for
$\Delta(t)$ from the self-consistency equation which contains two different
terms: (i) smooth dependence $\Delta_{h}[h(t)]$ arising directly from spin-
flip tunneling and depending on the occupation probability of QP states
$|C_{n}(t)|$; (ii) oscillating part $\delta\Delta(t)$ originated from the
interference between redistributed states with terms of type
$C_{n}^{*}(t)C_{n^{\prime}}(t)$. In addition, in Sec. IV.5-IV.6 we discussed
the spin imbalance generated due to LZSM tunneling and corresponding dynamical
magnetization of the QP gas. Discussion and some experimental proposals are
presented in Sec. V.

##  II Time-dependent Bogoliubov – de Gennes equations

We consider a homogeneous s-wave superconductor in the presence of the uniform
time-dependent Zeeman field $h(t)$ and Rashba spin-orbit coupling (RSOC). The
coherent QP dynamics is governed by the TDBdG equations [33]

$\displaystyle i\frac{\partial}{\partial
t}\check{\psi}_{k}=\check{\mathcal{H}}(k,t)\check{\psi}_{k},$ (1)

where the Hamiltonian

$\displaystyle\check{\mathcal{H}}(k,t)=\begin{pmatrix}\hat{H}(k,t)&i\hat{\sigma}_{y}\Delta(t)\\\
-i\hat{\sigma}_{y}\Delta(t)&-\hat{H}^{*}(-k,t)\\\ \end{pmatrix}$ (2)

is the $4\times 4$ matrix in the Nambu$\times$Spin space with the Pauli
matrices $\hat{\sigma}_{i}$ acting on the four-component wave function
$\check{\psi}_{k}(t)$. The single particle matrix Hamiltonian in the spin
space $\hat{H}(k,t)=\xi_{k}\hat{\sigma}_{0}-h(t){\bf
z}_{0}\hat{{\bm{\sigma}}}+\alpha(\hat{{\bm{\sigma}}}\times{\bf k}){\bf z}_{0}$
depends on the modulus $k=|{\bf k}|$ and the relative phase
$\theta_{k}=\text{arg}(k_{x}+ik_{y})$ of the momentum. Here
$\xi_{k}=k^{2}/2m-E_{F}$ is a free particle spectrum measured from the Fermi
level and $\alpha$ is a strength of RSOC. Hereafter we put $\hbar=1$. For
simplicity we consider here the motion of QPs only in the $x-y$ plane
neglecting their dispersion along the ${\bf z}_{0}$ axis, so that ${\bf
k}=(k_{x},k_{y})$.

The pairing potential $\Delta(t)$ should satisfy the self-consistency
equation, which at zero temperature $T=0$ can be written as follows

$\displaystyle\Delta(t)=-\frac{\lambda}{2}\sum_{\text{i.c.}}\check{\psi}^{\dagger}_{k}(t)\check{\tau}_{\Delta}\check{\psi}_{k}(t),$
(3)

where $\lambda$ is the pairing constant,
$\check{\tau}_{\Delta}=(\hat{\tau}_{x}+i\hat{\tau}_{y})\otimes
i\hat{\sigma}_{y}/2$ and the independence of $\Delta$ on $\theta_{k}$ is taken
into account. The summation here is performed over all solutions of Eq. (1)
for different initial conditions (i.c.) at $t=0$. The information about the
dynamics as well as the distribution function of the QP excitations is
contained in the functions $\check{\psi}_{k}(t)$, which self-consistently
define the temporal evolution of the gap. In the homogeneous problem, the
initial conditions are numbered by the momentum $k$, which, in the case of a
spin-split superconductor, must be supplemented by the spin quantum number.
All possible initial configurations of the QP states are defined by an
equilibrium distribution function. The pairing potential $\Delta(t)$ can be
chosen as a real function of time, and this choice will be justified below.

Generally speaking, the concept of an energy spectrum for a dynamical system
is not clearly defined. However, in the case of adiabatic evolution one can
introduce the eikonal approximation for the QP wavefunctions
$\check{\psi}_{k}(t)\propto\check{\Psi}_{k}(t)e^{iS_{k}(t)}$, from which the
adiabatic spectrum $E_{k}(t)=-\partial_{t}S_{k}$ can be extracted. The
functions $\check{\Psi}_{k}(t)$ are the instantaneous eigenstates of the
Hamiltonian $\check{\mathcal{H}}(t)$ from Eq. (2). The resulting spectrum is

$\displaystyle E_{kn}(t)=$ (4)
$\displaystyle\pm\sqrt{E_{0}^{2}+\alpha^{2}k^{2}+h^{2}(t)\mp\text{sgn}(\sigma)2\sqrt{\xi_{k}^{2}\alpha^{2}k^{2}+h^{2}(t)E_{0}^{2}}}$

where $E_{0}=\sqrt{\xi_{k}^{2}+\Delta^{2}}$. We use the index
$n\equiv\sigma\pm=\\{\uparrow+,\downarrow+,\uparrow-,\downarrow-\\}$ which
refers to different spin subbands and positive/negative energy (these
notations will be used in the text below). There are four corresponding
instantaneous eigenstates which can be written as
$\check{\Psi}_{kn}(t)=(u_{k\uparrow n},u_{k\downarrow n},v_{k\uparrow
n},v_{k\downarrow n})^{T}$. The detailed structure of the vectors is given in
Appendix A. The functions $\check{\Psi}_{kn}(t)$ form an orthonormal basis
with he normalization condition
$\check{\Psi}^{\dagger}_{kn}\check{\Psi}_{kn^{\prime}}=\delta_{nn^{\prime}}$
and the completness relation
$\sum_{kn}\check{\Psi}_{kn}\check{\Psi}^{\dagger}_{kn}=\check{1}$. Obviously,
in the limit of the stationary Zeeman field, $\check{\Psi}_{kn}$ becomes an
exact solution of stationary problem (1).

It is important to keep in mind that in the presence of both RSOC and spin-
splitting field the equilibrium gap value depends of the values of these
fields $\Delta_{\text{eq}}=\Delta_{\text{eq}}(h,\alpha)$. In what follows, the
RSOC strength $\alpha$ will be considered as a small parameter, and the
dependence $\Delta(h,\alpha)$ will be neglected for simplicity. Thus, the
equilibrium gap value is defined as follows

$\Delta_{\text{eq}}=\Delta_{0}=2\hbar\omega_{D}e^{-\frac{1}{\lambda N(0)}},$

where $\omega_{D}$ is Debye frequency and $N(0)$ is the density of states at
Fermi energy.

##  III Linearized gap dynamics

In this section we want to address the temporal evolution of a small
fluctuation of the gap $\Delta_{0}+\delta\Delta(t)$ in the presence of the
static spin-splitting field ${\bf h}=h_{0}{\bf z}_{0}$. The gap dynamics can
be excited by some external pulse at $t=0$ or can be driven, for instance, by
time-dependent spin-splitting field $\delta{\bf h}(t)=\delta h(t){\bf z}_{0}$.
In linear order in small perturbations $\delta\Delta(t),\delta h(t)\ll
h_{0}<\Delta_{0}$, the TDBdG equations for the QP wave functions read

$\displaystyle i\frac{\partial}{\partial
t}\check{\psi}_{k}(t)=\Big{[}\check{\mathcal{H}_{0}}+\check{\mathcal{V}}(t)\Big{]}\check{\psi}_{k}(t),$
(5)

where the operators in the Nambu$\times$Spin space are

$\displaystyle\check{\mathcal{H}_{0}}=\begin{pmatrix}\hat{H}_{0}(k)&i\hat{\sigma}_{y}\Delta_{0}\\\
-i\hat{\sigma}_{y}\Delta_{0}&-\hat{H}^{*}_{0}(-k)\\\ \end{pmatrix},$ (6)
$\displaystyle\check{\mathcal{V}}(t)=\begin{pmatrix}-\delta
h(t)\hat{\sigma}_{z}&i\hat{\sigma}_{y}\delta\Delta(t)\\\
-i\hat{\sigma}_{y}\delta\Delta(t)&-\delta h(t)\hat{\sigma}_{z}\end{pmatrix},$

and single particle Hamiltonian is
$\hat{H}_{0}(k)=\xi_{k}\hat{\sigma}_{0}-h_{0}\hat{\sigma}_{z}+\alpha(k_{y}\hat{\sigma}_{x}-k_{x}\hat{\sigma}_{y})$.

Time-dependent equation (5) can be written in the adiabatic basis using
stationary eigenfunctions $\check{\Psi}_{kn}$ of the operator
$\check{\mathcal{H}}_{0}$. Additionally, the RSOC energy $\alpha
k\approx\alpha k_{F}$ is considered a perturbative parameter. By approximating
the eigenvectors up to first order in $\alpha k_{F}/\Delta_{0}$ (see Appendix
A), we can infer from equation (3) that the fluctuation in the gap will have
an order up to $\mathcal{O}(\alpha^{2}k_{F}^{2}/\Delta_{0}^{2})$. However, in
the general case, the gap $\Delta$ should not be affected by the direction of
the SOC. Therefore, the first-order change in the gap
$\delta\Delta\propto\mathcal{O}(\alpha k_{F}/\Delta_{0})$ must vanish.

Instead of the general eikonal theory, we use the perturbative approach with
the ansatz written in terms of the dynamical phase

$\displaystyle\check{\psi}_{k}(t)=\sum_{n}\check{\Psi}_{kn}C_{kn}(t)e^{-iE_{kn}t}.$
(7)

The index $n=\\{\uparrow+,\downarrow+,\uparrow-,\downarrow-\\}$ denotes the
spectral branches and all negative/positive energy terms are involved into the
dynamics of QPs. Substituting the function (7) into Eq. (5) we obtain the
equation for the dynamics of the coefficients

$\displaystyle i\frac{\partial}{\partial
t}C_{km}(t)=\sum_{n}\check{\Psi}_{m}^{\dagger}\check{\mathcal{V}}(t)\check{\Psi}_{n}e^{-i(E_{n}-E_{m})t}C_{kn}(t),$
(8)

which completely determine the temporal evolution of gap $\Delta(t)$ through
the self-consistency equation

$\displaystyle\Delta_{0}+\delta\Delta(t)=$ (9)
$\displaystyle-\frac{\lambda}{2}\sum_{\text{i.c.}}\sum_{n,n^{\prime}}C^{*}_{kn}(t)C_{kn^{\prime}}(t)e^{-i(E_{n^{\prime}}-E_{n})t}\check{\Psi}_{kn}^{\dagger}\check{\tau}_{\Delta}\check{\Psi}_{kn^{\prime}}.$

In the case of zero temperature $T=0$ there are two possible initial
configurations at $t=-\infty$: All QP states with energies below Fermi level
in the first(second) spin subband with $\sigma=\uparrow(\downarrow)$ are fully
occupied for all momenta with $\xi_{k}\in(-\omega_{D},\omega_{D})$. The
corresponding initial conditions can be written as

$\displaystyle C_{kn}(t=-\infty)=\delta_{n,l}$ (10)

with Kronecker delta $\delta_{n,n^{\prime}}$ and indices
$l=\\{\uparrow-,\downarrow-\\}$. Therefore, it is natural to linearize the
equation (8) as follows

$\displaystyle C_{kn}(t)=C_{kn}(-\infty)+\delta C_{kn}(t).$ (11)

Performing Laplace transform in the complex plane $s=i\omega+\zeta$ for the
linearized equations (8, 9, 11) (see Appendix B) we get the following
expression for the gap perturbation

$\displaystyle\delta\Delta(s)=\Big{[}\mathcal{K}_{0}(s)+\mathcal{K}_{+}(s)+\mathcal{K}_{-}(s)\Big{]}\delta\Delta(s)$
(12) $\displaystyle+\delta
h(s)\Big{[}\mathcal{F}_{+}(s)-\mathcal{F}_{-}(s)\Big{]}+\mathcal{I}(s).$

Here $\mathcal{K}_{0,\pm}(s)$ represents kernels of the self-consistency
equation; $\mathcal{F}_{\pm}(s)$ defines the dynamical structure of the
"force" term (in analogy with a mechanical oscillator) related with $\delta
h(t)$; and $\mathcal{I}(s)$ includes all terms related to perturbations at the
moment $t=0$. Due to the absence of particle-hole asymmetry, which couples the
phase and amplitude fluctuations [6], the imaginary part of $\delta\Delta(s)$
naturally vanishes and we consider only amplitude (or Higgs) modes of the
superconducting gap. Knowing the function $\mathcal{K}(s)$ one can find
eigenfrequencies and free dynamics of the system, while $\mathcal{F}_{\pm}(s)$
induces the driven dynamics. We will conduct a thorough examination of these
terms below.

###  III.1 Spin-split Higgs modes

Figure 1:  (a) Branch points
$\omega_{0}=2\Delta_{0},\omega_{\pm}=2(\Delta_{0}\pm h_{0})$ (red dots)
corresponding to the kernels $\mathcal{K}_{0}(s)$ and
$\mathcal{K}_{\pm}(s),\mathcal{F}_{\pm}(s)$ [Eq. (12)] in the complex plane
$s=\zeta+i\omega$. Red lines show the chosen branch cuts. Black crosses
correspond to the poles of the external force $\delta h(s)$. (b-d)
Illustration of physical mechanism behind the appearance of three
eigenfrequencies $\omega_{+}$ (b), $\omega_{0}$ (c), $\omega_{-}$ (d).

It is known that in the absence of a spin-splitting field and RSOC the Higgs
mode has a singular behavior in the vicinity of the eigenfrequency
$\omega=2\Delta_{0}$, which defines the free evolution of the gap perturbation
$\delta\Delta(t)\propto\cos(2\Delta_{0}t)/\sqrt{t}$ [14]. Since the energy of
the Higgs mode lies at the lower bound of the QP spectrum, the oscillatory
behavior here can be represented as a coherent decay and formation of a Cooper
pair into two QPs with opposite spins and energies $\Delta_{0}$ at $k\approx
k_{F}$. The contribution from the pairs of QPs with other momenta leads to the
inhomogeneous broadening of the mode with the corresponding damping law. The
presence of Zeeman field and RSOC makes the dynamics more complicated. To
analyze the eigenmodes of the superconductor one can set $\delta h(t)=0$ and
write the self-consistency equation as follows

$\chi_{\Delta\Delta}^{-1}(s)\delta\Delta(s)=0,$

where we define the bare pair susceptibility

$\displaystyle\chi_{\Delta\Delta}(s)=\frac{1}{1-\mathcal{K}_{0}(s)-\mathcal{K}_{+}(s)-\mathcal{K}_{-}(s)}.$
(13)

The corresponding kernels read (see Appendix B)

$\displaystyle\mathcal{K}_{0}(s)=\Big{\langle}\frac{2\xi^{2}}{E_{0}}\frac{1}{s^{2}+4E^{2}_{0}}\Big{\rangle}\propto\mathcal{O}\Big{(}\frac{\alpha^{0}k^{0}_{F}}{\Delta_{0}^{0}}\Big{)},$
(14)
$\displaystyle\mathcal{K}_{\pm}(s)=\Big{\langle}\mathcal{A}^{2}(\xi)\frac{E_{0}\pm
h_{0}}{s^{2}+4(E_{0}\pm
h_{0})^{2}}\Big{\rangle}\propto\mathcal{O}\Big{(}\frac{\alpha^{2}k^{2}_{F}}{\Delta_{0}^{2}}\Big{)},$

where the notation $\big{\langle}\dots\big{\rangle}=\lambda
N(0)\int_{-\omega_{D}}^{\omega_{D}}d\xi$ is used. The function
$\mathcal{A}(\xi)\propto\check{\Psi}^{0\dagger}_{kn}\check{\tau}_{\Delta}\check{\Psi}^{0}_{kn}\propto\alpha
k_{F}/\Delta_{0}$ is proportional to nonzero triplet component of the wave
function, therefore the kernels $\mathcal{K}_{\pm}$ are of the second order in
the RSOC parameter.

The frequencies of the eigenmodes of the superconducting condensate can be
traced out from the condition $|\chi_{\Delta\Delta}^{-1}(\omega)|=0$, which
reflects the singular points of the kernels (14). Consider these points in
more detail. Instead of straightforward integrating, we are going to implement
the analysis in the spirit of the work [3] and analytically obtain the limit
${\zeta\rightarrow 0}$. The functions
${\mathcal{K}_{0,\pm}(s\rightarrow\omega)}$ can be represented as
${\mathcal{K}(s)=\mathcal{K}^{\prime}(\omega)+i\text{sgn}(\omega\zeta)\mathcal{K}^{\prime\prime}(\omega)}$.
The real parts of the kernels

$\displaystyle\frac{\mathcal{K}_{0}^{\prime}(\omega)}{\lambda N(0)}$
$\displaystyle=\fint_{-\omega_{D}}^{\omega_{D}}\frac{2\xi^{2}}{\sqrt{\xi^{2}+\Delta_{0}^{2}}(4\xi^{2}+4\Delta_{0}^{2}-\omega^{2})}d\xi,$
(15) $\displaystyle\frac{\mathcal{K}_{\pm}^{\prime}(\omega)}{\lambda N(0)}$
$\displaystyle=\fint_{-\omega_{D}}^{\omega_{D}}\frac{\mathcal{A}^{2}(\xi)(E_{0}\pm
h_{0})}{4(E_{0}-h_{0})^{2}\pm|\omega|^{2}}d\xi,$ (16)

are regular on the imaginary axis $s=i\omega$. The imaginary parts are

$\displaystyle\frac{\mathcal{K}_{0}^{\prime\prime}(\omega)}{\lambda N(0)}$
$\displaystyle=-\frac{\pi}{2}\frac{\sqrt{\omega^{2}-\omega_{0}^{2}}}{|\omega|}\Theta[\omega^{2}-\omega_{0}^{2}],$
(17) $\displaystyle\frac{\mathcal{K}_{\pm}^{\prime\prime}(\omega)}{\lambda
N(0)}$ $\displaystyle=-\frac{\pi}{8}\frac{|\omega|\mp
2h_{0}}{\xi_{\pm}}\mathcal{A}^{2}(\xi_{\pm})\Theta[\omega^{2}-\omega_{\pm}^{2}],$
(18)

where
$\xi_{\pm}=\frac{1}{2}\sqrt{\big{(}|\omega|-\omega_{\pm}\big{)}^{2}+4\Delta_{0}(|\omega|-\omega_{\pm})}$.
The discontinuities at the real axis $\zeta$ mean the existence of the branch
points

$\displaystyle\omega_{0}=2\Delta_{0},$ (19)
$\displaystyle\omega_{+}=2(\Delta_{0}+h_{0}),$
$\displaystyle\omega_{-}=2(\Delta_{0}-h_{0}),$

and corresponding cuts in the complex plane [Fig. 1(a)].

The analysis of the general linear response of the order parameter can be
significantly simplified by expanding the susceptibility
$|\chi_{\Delta\Delta}(\omega)|$ in the powers of the small parameter $\alpha
k_{F}/\Delta_{0}$, since the kernels
$\mathcal{K}_{\pm}\propto\mathcal{O}(\alpha^{2}k^{2}_{F}/\Delta_{0}^{2})$. As
mentioned before, the maximum order we can take into account is
$|\chi_{\Delta\Delta}|\propto\mathcal{O}(\alpha^{2}k_{F}^{2}/\Delta_{0}^{2})$.
The resonance condition $|\chi_{\Delta\Delta}^{-1}(\omega)|=0$ is satisfied at
$\omega=\omega_{\pm}$ where the kernels $\mathcal{K}_{\pm}^{\prime\prime}(s)$
have a singularity (note that $\mathcal{A}$ is regular at $\xi=\xi_{\pm}$),
and at $\omega=\omega_{0}$, where the function
$\mathcal{K}_{0}^{\prime\prime}(s)$ goes to zero. Thus, the branch points (19)
define new eigenmodes of the superconductor in the presence of spin-splitting
field and weak RSOC. By taking some constant initial condition
$\mathcal{I}(s)=const$ in the RHS of Eq. (12) and using inverse Laplace
transform one can consider an impulse response of the superconductor. It can
be shown that the peculiarities in the vicinities of the eigenfrequencies lead
to three partial contribution to the long-time gap dynamics

$\displaystyle\delta\Delta(t)\propto\frac{\cos(\omega_{0}t+\pi/4)}{\sqrt{\Delta_{0}t}}$
(20) $\displaystyle+\Big{(}\frac{\alpha
k_{F}}{\Delta_{0}}\Big{)}^{2}A_{+}(h_{0})\frac{\cos(\omega_{+}t+\pi/4)}{\sqrt{\Delta_{0}t}}$
$\displaystyle+\Big{(}\frac{\alpha
k_{F}}{\Delta_{0}}\Big{)}^{2}A_{-}(h_{0})\frac{\cos(\omega_{-}t+\pi/4)}{\sqrt{\Delta_{0}t}}$

with some amplitudes $A_{\pm}(h_{0})$, which can be identified as spin-split
Higgs modes. Details of the derivation of $\delta\Delta(t)$ from Eq. (12) will
be provided in the next subsection together with Appendix C.

Appearance of the frequencies (19) and corresponding oscillations (20) in the
spin-split superconductor can be explained qualitatively. Coherent decay of
the Cooper pairs from the Fermi level can occur into two different spin
subbands of the QP spectrum. When two electrons with opposite spins from a
pair dissociate into two QP at $k\approx k_{F}$ with the energies
$\Delta_{0}\pm h_{0}$ without spin-flipping, then the total decay energy is
equal to QP threshold $\approx 2\Delta_{0}$. This process corresponds to the
mode $2\Delta_{0}$ and shown in Fig. 1(c). A decay into two QPs with the same
spins is possible in the presence of RSOC due to the effective spin-flip
scattering. The energies of such two QPs are either $\Delta_{0}+h_{0}$ or
$\Delta_{0}-h_{0}$. This process leads to the modes $2(\Delta_{0}\pm h_{0})$
correspondingly [Fig. 1(b,d)]. Note that this naive interpretation of the
complicated QP dynamics is valid for the sufficiently small RSOC $\alpha
k_{F}\ll\Delta_{0}$.

Numerically calculated susceptibility $|\chi_{\Delta\Delta}(\omega)|$ from
Eqs. (13-14) is shown in Fig. (2)(a). The observed resonances have a different
parametric order of smallness. The Higgs mode with the frequency $\omega_{0}$
which exists in the absent the RCOS becomes dominating with more pronounced
peak
$|\chi_{\Delta\Delta}(\omega\approx\omega_{0})|\propto\alpha^{0}k_{F}^{0}$,
whereas two other modes at shifted frequencies $\omega_{\pm}$ are of the order
of
$|\chi_{\Delta\Delta}(\omega\approx\omega_{\pm})|\propto\alpha^{2}k^{2}_{F}/\Delta_{0}^{2}$.
These modes merge with $\omega_{0}$ at $h_{0}\rightarrow 0$ and disappear for
$\alpha\rightarrow 0$. It is expected that the excitation of the bare response
of the superconductor can be implemented with the standard THz laser pump-
probe techniques. The electric field of the pump pulse produces a quench of
the spin-split superconductor and subsequent probe pulse detects the
multifrequency Higgs oscillations.

Note, that a similar dynamics of the order parameter has been studied in the
spin-orbit coupled Fermi gases [36, 37, 38]. In particular, the existence of
the Higgs modes modified by the Zeeman field in the presence of strong SOC
with $\alpha k_{F}\sim h(t)\sim E_{F}$ was discussed in Ref. [34]. The authors
performed a numerical simulation of the one-dimensional Fermi superfluid and
examined the excitation of the gap oscillations with few frequencies by abrupt
change of the Zeeman field. Despite the significant differences between the
models, there is a general tendency for the influence of the shift of spectral
QP branches on the behavior of the order parameter modes.

Figure 2:  (a) The susceptibility $|\chi_{\Delta\Delta}(\omega)|$ of free
dynamics of the gap perturbation $\delta\Delta(\omega)$. Features at the
frequencies $\omega_{0}=2\Delta_{0}$ and $\omega_{\pm}=2(\Delta_{0}\pm h_{0})$
correspond to Higgs modes resonances. The finite width of the resonances is
due to the small imaginary part of the frequency $\omega+i\beta$. (b)
Frequency dependence of the response functions $|\chi_{\Delta h}(\omega)|$ of
the driven gap oscillations $\delta\Delta(t)$ excited by the Zeeman field
$\delta h(t)=\delta h_{0}e^{-\beta t}\cos(\omega t)$ (see Eq. 25). Both plots
are symmetrical with respect to $\omega\rightarrow-\omega$ and have the
parameters $h_{0}=0.8\Delta_{0}$, $\alpha k_{F}=0.05\Delta_{0}$ and
$\beta=0.002\Delta_{0}$.

###  III.2 Coupling of Higgs modes and Zeeman field

We found that, in addition to an electromagnetic field, the gap dynamics in a
spin-split superconductor can be excited by a nonstationary component of
Zeeman field ${\bf h}(t)=(h_{0}+\delta h(t)){\bf z}_{0}$. In this particular
configuration the perturbation of the spin-splitting field $\delta h(t)$
appears in the self-consistency equation (12) in the first order, which is the
trace of a dot product $({\bf h}_{0}\cdot\delta{\bf h})$. Note that the field
$\delta h(s)$ is weighted by the functions

$\displaystyle\mathcal{F}_{\pm}(s)=\Big{\langle}\mathcal{A}(\xi)\mathcal{B}(\xi)\frac{(E_{0}\pm
h_{0})}{s^{2}+4(E_{0}\pm
h_{0})^{2}}\Big{\rangle}\propto\mathcal{O}\Big{(}\frac{\alpha^{2}k^{2}_{F}}{\Delta_{0}^{2}}\Big{)},$
(21)

which can be written as
$\mathcal{F}(s)=\mathcal{F}^{\prime}(\omega)+i\text{sgn}(\omega\zeta)\mathcal{F}^{\prime\prime}(\omega)$
and have the same order in $\alpha k_{F}$ and the same analytical properties
as the kernels $\mathcal{K}_{\pm}(s)$ in (16,18), because both functions
$\mathcal{A}^{2}(\xi)$ and $\mathcal{A}(\xi)\mathcal{B}(\xi)$ are regular for
$\xi\in(-\omega_{D},\omega_{D})$. The presence of the singular points in the
force term makes the analysis of Eq. (12) more sophisticated, despite the fact
that these points are shared with other kernels.

Consider the general case of forced oscillations of the order parameter driven
by some field $\delta h(t)$ which is abruptly turned on at $t=0$. Since we
want to consider dynamical effects related only to the external force, we
neglect the initial conditions $\mathcal{I}(s)$ in Eq. (12), e.g. assume the
equilibrium system with $\delta\Delta(t)=0$. It is convenient to introduce the
linear response function

$\displaystyle\chi_{\Delta
h}(s)=\frac{\mathcal{F}_{+}(s)-\mathcal{F}_{-}(s)}{1-\mathcal{K}_{0}(s)-\mathcal{K}_{+}(s)-\mathcal{K}_{-}(s)}$
(22)

and write the self-consistency equation as ${\delta\Delta(s)=\chi_{\Delta
h}(s)\delta h(s)}$. One can obtain the expression for $\delta\Delta(t)$ in the
interval $t\in[0,\infty)$ using inverse Laplace method:

$\displaystyle\delta\Delta(t)=\frac{1}{2\pi
i}\int_{-i\infty+\epsilon}^{i\infty+\epsilon}\chi_{\Delta h}(s)\delta
h(s)e^{st}ds,$ (23)

where $\epsilon$ should be larger then the real part of the poles of $\delta
h(s)$. The integral can be evaluated using closed contour shown in Fig. 1(a).
Making sure that all integrals on infinitely large and small arcs vanish and
applying residue theorem we get

$\displaystyle\delta\Delta(t)=\sum_{p}\chi_{\Delta
h}(s_{p})e^{s_{p}t}\underset{s=s_{p}}{\text{Res}}\Big{[}\delta h(s)\Big{]}$
(24)
$\displaystyle+\frac{2}{\pi}\int_{\omega_{-}}^{\infty}\text{Im}\chi_{\Delta
h}(s)\big{|}_{\zeta\rightarrow+0}\text{Im}\Big{[}e^{i\omega t}\delta
h(i\omega)\Big{]}d\omega.$

The first term represents the contribution from the poles of the external
field $\delta h(t)$, while the second term is the contribution from the
integrals along the branch cuts.

The imaginary part of the susceptibility (22) can be expanded up to the second
order in $\alpha k_{F}/\Delta_{0}$, since all the kernels
$\mathcal{F}_{\pm},\mathcal{K}_{\pm}\propto\mathcal{O}(\alpha^{2}k^{2}_{F}/\Delta_{0}^{2})$.
That allows one to distinguish different strongly dominant terms of the
function $\text{Im}\chi_{\Delta h}(\omega)$ in (24) in the vicinity of
different branch points (19) and estimate their contribution to an asymptotic
expression for $\delta\Delta(t)$. The detailed calculations are provided in
Appendix C. Here we write the result for the superconducting gap oscillations,
which at large times $h_{0}^{-1}\ll t$ reads

$\displaystyle\delta\Delta(t)\approx\sum_{p}\chi_{\Delta
h}(s_{p})e^{s_{p}t}\underset{s=s_{p}}{\text{Res}}\big{[}\delta
h(s)\big{]}+\frac{4\Delta_{0}}{\pi\sqrt{\pi}}\frac{\big{[}\mathcal{F}_{+}^{\prime}(\omega_{0})-\mathcal{F}_{-}^{\prime}(\omega_{0})\big{]}}{\lambda
N(0)}\frac{\text{Im}\Big{[}\delta
h(i\omega_{0})e^{i(\omega_{0}t+\pi/4)}\Big{]}}{\sqrt{\Delta_{0}t}}$ (25)
$\displaystyle+\frac{\sqrt{\pi}}{8\sqrt{2}}\frac{(\alpha
k_{F})^{2}\Delta_{0}}{(\Delta_{0}-h_{0})^{2}}\sum_{j=\pm}\frac{\lambda
N(0)\big{[}1-\mathcal{K}_{0}^{\prime}(\omega_{j})\big{]}}{\big{[}1-\mathcal{K}_{0}^{\prime}(\omega_{j})\big{]}^{2}+\big{[}\mathcal{K}_{0}^{\prime\prime}(\omega_{j})\big{]}^{2}}\frac{\text{Im}\Big{[}\delta
h(i\omega_{j})e^{i(\omega_{j}t+\pi/4)}\Big{]}}{\sqrt{\Delta_{0}t}}.$

The first term here is related to the forced oscillations of the gap, caused
by the Zeeman field $\delta h(t)$. The last three terms correspond to the free
oscillations triggered by $\delta h(t)$ at $t=0$ in the long time asymptote,
with three characteristic frequencies (19) and square-root damping law. The
latter can be interpreted as partial contribution from the Higgs modes in the
spin-splitting field $h_{0}$.

The eigenmodes decay at $t\rightarrow\infty$ and in the long-time asymptote
the forced oscillations prevail. Consider the steady-state behavior of
$\delta\Delta(t)$ (the first term in Eq. (25)) in the time interval restricted
by the inelastic relaxation processes where the presented description of the
coherent gap dynamics is valid. Assume the general harmonic perturbation
$\delta h(t)=\delta h_{0}e^{-\beta t}\cos(\omega t)$ at $t\geq 0$ with small
finite damping factor $\beta\rightarrow 0$.The amplitude of the driven gap
perturbation $\delta\Delta(t)$ is defined by the Zeeman field $h_{0}$ and
susceptibility $\chi_{\Delta h}(s_{0})$ taken at the pole of the force
$s_{0}=-\beta+i\omega$. The numerically integrated shape of $|\chi_{\Delta
h}(\omega)|$ is shown in Fig. (2)(a). The response of the superconductor, as
expected, has three resonance peaks at the frequencies $\omega_{0,\pm}$.
However, since the the external field $\delta{\bf h}(t)$ couples to the gap
through the RSOC, the amplitude of the susceptibility in the vicinity of the
resonances has the same order of smallness $|\chi_{\Delta
h}(\omega_{0,\pm})|\propto\mathcal{O}(\alpha^{2}k_{F}^{2}/\Delta_{0}^{2})$,
which differs from the bare response (13).

In this section, we have solely focused on the longitudinal component of the
field perturbation $\delta h(t){\bf z}_{0}$ with respect to the stationary
field $h_{0}{\bf z}_{0}$. However, it is also possible to introduce the time-
dependent transversal component $\delta{\bf h}\bot(t)$ and examine its dynamic
interaction with the superconducting system in Eq. (5). This component
generates triplet correlations, but these do not contribute to the order
parameter since only singlet pairing in (3) is considered. Consequently, in
the second-order perturbation theory with respect to $\alpha
k_{F}/\Delta_{0}$, there is no linear coupling between the field $\delta{\bf
h}\bot(t)$ and the gap $\delta\Delta(t).$ This outcome is unsurprising since
the only true scalar in this regime $(\delta{\bf h}_{\bot}\cdot{\bf h}_{0})$
is zero.

##  IV Evolution of QP states in strong Zeeman field

In this section we address the case of a linearly growing spin-splitting field
$h(t)=\gamma t$, which can exceed the equilibrium value of the superconducting
gap $\Delta_{\text{eq}}\equiv\Delta_{0}$ and thus provide the crossing of the
two QP spectral branches $E_{\uparrow+}(\xi_{k})$ and
$E_{\downarrow-}(\xi_{k})$ from different spin subbands [Fig. 3(a, c)]. In the
collisionless regime and in the absence of RSOC the intersecting spectral
branches do not interact, so that the occupation of the quasiparticle states
defined at $t=-\infty$ does not change in time. This means that the self-
consistent gap function will not change even above the paramagnetic limit
$h(t)>\Delta_{0}$ and will be defined by the initial condition
$\Delta(t)=\Delta_{0}$. It is clear from the general assumptions that the
spin-orbit coupling is capable of provoking the interplay between QP states
with different spins, and we investigate the mechanism of such an interaction
and the effect on the superconducting order parameter $\Delta(t)$. As
mentioned in the Section II, we will treat the RSOC energy as a small
parameter $\alpha k_{F}/\Delta\ll 1$. Therefore, we neglect the dependence of
the equilibrium gap $\Delta_{\text{eq}}$ on $\alpha$ and assume
$\Delta_{\text{eq}}\equiv\Delta_{0}$.

###  IV.1 Adiabatic evolution of QP states

The evolution of QP wavefunction of the TDBdG equations (1) can be regarded
with the help of general adiabatic ansatz

$\displaystyle\check{\psi}_{k}(t)=\sum_{n}C_{kn}(t)\check{\Psi}_{kn}(t),$ (26)

where $\check{\Psi}_{kn}(t)$ are the instantaneous eigenstates of the
Hamiltonian (2). Here all negative/positive energy terms with the indices
$n=\\{\uparrow+,\downarrow+,\uparrow-,\downarrow-\\}$ are taken into account.
The coefficients $C_{kn}(t)$ define the occupation of QP states and its
temporal evolution. The initial conditions for $C(t)$ should be fixed by the
equilibrium distribution at $t=0$. In the case of spin-split superconductor at
zero temperature $T=0$ there are two possible initial configurations: All QP
states with energies below Fermi level in the first(second) spin subband with
$\sigma=\uparrow(\downarrow)$ are fully occupied for all momenta with
$\xi_{k}\in(-\omega_{D},\omega_{D})$. With short notations one can write this
as $C_{kn}(t=0)=\delta_{n,l}$, where $\delta_{n,n^{\prime}}$ is Kronecker
delta and $l=\\{\uparrow-,\downarrow-\\}$. This means that for the given $l$
we have $C_{kl}(t=0)=1$, and all other $C_{k(n\neq l)}(t=0)=0$.

The vector

$\displaystyle\hat{C}_{k}(t)=(C_{k\uparrow+},C_{k\downarrow+},C_{k\uparrow-},C_{k\downarrow-})^{T}$
(27)

contains all the information about the dynamics of the QP states.
Corresponding adiabatic temporal evolution can be described with the help of
the unitary operator
$\hat{C}_{k}(t_{2})=\hat{U}_{k}(t_{2},t_{1})\hat{C}_{k}(t_{1})$ where
$\hat{U}_{k}=\text{diag}(U_{k\uparrow+},U_{k\downarrow+},U_{k\uparrow-},U_{k\downarrow-})$
and

$\displaystyle
U_{kn}(t_{2},t_{1})=\exp\Big{(}-i\int_{t_{1}}^{t_{2}}E_{kn}(t)dt\Big{)}.$ (28)

The interaction of the branches $E_{k\uparrow+}(t)$ and $E_{k\downarrow-}(t)$
leads to avoided crossing of the QP levels at fixed energy $\xi_{k}$ with the
splitting proportional to $\alpha k_{F}$. Thus the adiabatic approximation is
justified only for the levels with $E_{k}(t)\gg\alpha k_{F}$, e.g. far enough
from the crossing points. Therefore, for the Zeeman field
$h(t)\lesssim\Delta_{0}$ all nonadiabatic transitions are suppressed and the
gap function defined by the self-consistency equation (3) is equal to the
equilibrium value $\Delta(t)=\Delta_{0}$.

Figure 3:  (a,b) QP spectrum $E_{k}$ from (4) for two values of Zeeman field
$h(t)$ before (a) and after (b) avoided crossing. Colored/empty circles
correspond to filled/empty states. The parameters are chosen as follows:
$\Delta_{0}/E_{F}=0.01$, $\alpha/E_{F}=0.0025$. In (e) the schematic temporal
evolution of the filling probabilities $|C_{k}|^{2}(t)$ for two states at
fixed $\xi_{k}$ is shown. The gray lines show a tunneling process similar to
the real one in the vicinity of the avoided transition point
$t_{0}(\xi_{k})=\sqrt{\xi_{k}^{2}+\Delta_{0}^{2}}/\gamma$, while red and blue
lines refer to transition matrix approximation of the LZSM tunneling with the
probability $p_{k}$. (b, d) The QP distribution function for one spin
projection $f_{\uparrow}(E,t)$ from Eq. (43) before and after crossing of
spectral branches at $\delta_{\text{LZ}}=0.5$.

###  IV.2 Transition evolution matrix

The avoided crossing between the spectral terms at $h(t)\gtrsim\Delta_{0}$
should be described in terms of nonadiabatic dynamics. For this we consider
the branch intersection as consecutive avoided crossing of pairs of the QP
states with fixed energy $\xi_{k}$ at the time instant
$t_{0}(\xi_{k})=\sqrt{\xi_{k}^{2}+\Delta^{2}}/\gamma$ [Fig. 3]. For each
crossing at $\xi_{k}\in(-\omega_{D},\omega_{D})$ it is possible to formulate
the time-dependent Landau-Zener-Stückelberg-Majorana (LZSM) problem [35],
which describes the transitions between two QP states with different spins
during their temporal evolution. Note that resulting nonadiabatic tunneling is
equivalent to dynamical spin-flip process.

In general, the description of such a tunneling (or LZSM problem) requires
joint solution of TDBdG equation (1) and self-consistency equation (3).
However, some important results can be obtained analytically using certain
approximations:

(i) If the time variation of the gap function $\Delta(t)$ is small on the
typical tunneling time scale $\tau_{\text{LZ}}$ (see Appendix D), then the
tunneling of QP states is not affected by the dynamics of the order parameter.

(ii) On the other hand, the gap $\Delta(t)$ is defined by all states in range
$\xi_{k}\in(-\omega_{D},\omega_{D})$, and a time-dependent perturbation of the
states caused by the dynamical LZSM transition makes a small contribution to
the sum over all $\xi_{k}$. Thus, one can neglect the transient dynamics of
the coefficients $\hat{C}_{k}(t)$ in the vicinity of a transition point for
each $\xi_{k}-$th mode. This also means that one can investigate the tunneling
problem with the help of so-called transition evolution matrix [35] connecting
two adiabatic regimes before ($t<t_{0}-$) and after ($t>t_{0}+$) avoided
crossing [Fig. 3(e)]. These conditions allow one to effectively decouple the
LZSM problem from the self-consistency equation and solve them without self-
consistency.

Taking into account all these assumptions, the time evolution of the vector
$\hat{C}_{k}(t)$ from the adiabatic ansatz (26) is described as follows

$\hat{C}_{k}(t)=\left\\{\begin{aligned}
\hat{U}_{k}(t,t_{0}+)\hat{S}_{\text{LZ}}\hat{U}_{k}(t_{0}-,0)\hat{C}_{k}(0),&\quad
t>t_{0}(\xi_{k})\\\ \hat{U}_{k}(t,0)\hat{C}_{k}(0),&\quad
t<t_{0}(\xi_{k})\end{aligned}\right..$ (29)

Here the nonadiabatic transitions between QP states are included into
transition matrix $\hat{S}_{\text{LZ}}$, which acts on the state vector
$\hat{C}_{k}(t)$ at the time instant $t=t_{0}(\xi_{k})$. The matrix
$\hat{S}_{\text{LZ}}$ can be obtained by considering the interaction of two
intersecting energy branches $E_{\uparrow+}$ and $E_{\downarrow-}$ in the
TDBdG equation (1). Using so-called diabatic basis (basis of Hamiltonian (2)
in the absence of RSOC) one gets a system of dynamical equations, the
asymptotic solution of which forms a transition matrix describing the passage
through the avoided intersection point. Then we go to the original adiabatic
basis (26) and get the matrix $\hat{S}_{\text{LZ}}$. The complete derivation
of $\hat{S}_{\text{LZ}}$ is presented in Appendix D and it reads

$\displaystyle\hat{S}_{\text{LZ}}=\begin{pmatrix}\sqrt{p_{k}}&0&0&\sqrt{1-p_{k}}e^{i(\dots)}\\\
0&1&0&0\\\ 0&0&1&0\\\
-\sqrt{1-p_{k}}e^{-i(\dots)}&0&0&\sqrt{p_{k}}\end{pmatrix},$ (30)

where $(\dots)=\chi_{k}-\theta_{k}-\frac{\pi}{2}\text{sgn}(\alpha)$. The
coefficient

$p_{k}=\exp\Big{[}-\delta_{\text{LZ}}\frac{\Delta^{2}}{\xi_{k}^{2}+\Delta^{2}}\Big{]}$

is expressed through the dimensionless LZSM parameter
$\delta_{\text{LZ}}=\pi\alpha^{2}k_{F}^{2}/\gamma$ and determines the
probability of tunneling between QP states with different spins. The
transition is accompanied by the appearance of the Stokes phase $\chi_{k}$
(see Appendix D) and the phase $\theta_{k}=\arg\big{(}k_{x}+ik_{y}\big{)}$.

To avoid confusion, we use the same notations for the spectral branches (4)
before ($t<t_{0}-$) and after ($t>t_{0}+$) QP transitions, as shown in Fig. 3.
Thereby, we do not need to keep track of the indices of the eigenvectors
$\check{\Psi}_{kn}(t)$ and the evolution operators $U_{kn}(t)$ from (28). It
is sufficient that these functions take into account the permutation of the
branches of the spectrum (4), so that all QP levels change their indices after
the transition in accordance with the chosen notation.

###  IV.3 Time dependence of superconducting gap

The time-dependent order parameter subjected to the field
$h(t)\gtrsim\Delta_{0}$ depends on both the adiabatic wave function (26) and
nonadiabatic LZSM tunneling (29). The calculation of $\Delta(t)$ can be
accomplished using the self-consistent equation (3), which gets the following
form

$\displaystyle\Delta(t)=-\frac{\lambda}{2}\sum_{l}\sum_{k}\sum_{n,n^{\prime}}C_{kn}^{*}(t)C_{kn^{\prime}}(t)\check{\Psi}_{kn}^{\dagger}\check{\tau}_{\Delta}\check{\Psi}_{kn^{\prime}},$
(31)

where index $l$ means different initial configurations of the occupation of
the QP spectrum at $t=0$ (see Section IV.1). The first configuration with
$C_{kn}(t=0)=\delta_{n,\uparrow-}$ corresponds to occupation of all QP states
belonging to the spectral branch $E_{k,\uparrow-}$ for all momenta with
$\xi_{k}\in(-\omega_{D},\omega_{D})$. The evolution of the coefficients
$C_{kn}(t)$ is determined by the Eq. (29) together with Eq. (27-28). Since the
branch $E_{k,\uparrow-}$ does not cross with other branches, the coefficients
$C_{kn}(t)$ have a trivial adiabatic dynamics, which can be written as follows

$\displaystyle\hat{C}_{k}(t)=\begin{pmatrix}0\\\ 0\\\
U_{\uparrow-}\big{(}t,0\big{)}\\\ 0\end{pmatrix}.$ (32)

The second initial configuration with $C_{kn}(t=0)=\delta_{n,\downarrow-}$
leads to the intersection of the filled branch $E_{k,\downarrow-}$ and empty
branch $E_{k,\uparrow+}$. Using equation (29) we obtain a nontrivial dynamics
of the states with LZSM tunneling, which reads

$\displaystyle\hat{C}_{k}(t)=$ (33)
$\displaystyle\begin{pmatrix}\sqrt{1-p_{k}}e^{i(\dots)}U_{\uparrow+}\big{(}t,t_{0}+\big{)}U_{\downarrow-}\big{(}t_{0}-,0\big{)}\Theta\big{[}t-t_{0}\big{]}\\\
0\\\ 0\\\
\Big{(}\sqrt{p_{k}}\Theta\big{[}t-t_{0}\big{]}+\Theta\big{[}t_{0}-t\big{]}\Big{)}U_{\downarrow-}\big{(}t,0\big{)}\end{pmatrix}.$

Here $(\dots)=\chi_{k}-\theta_{k}-\frac{\pi}{2}\text{sgn}(\alpha)$ and
$\Theta(t)$ is the Heaviside function.

Substituting coefficients (32) and (33) obtained from different initial
conditions together with the QP wavefunctions $\hat{\Psi}_{kn}$ from (50) into
the self-consistency equation (31) we get

$\displaystyle\Delta(t)=\lambda\sum_{|\xi_{k}|>\sqrt{h^{2}-\Delta^{2}}}u_{0}v_{0}$
(34)
$\displaystyle+\lambda\sum_{|\xi_{k}|<\sqrt{h^{2}-\Delta^{2}}}\Bigg{[}\Big{(}\overset{=1}{|C_{k\uparrow-}|^{2}}+\overset{\propto
p_{k}}{|C_{k\downarrow-}|^{2}}-\overset{\propto
1-p_{k}}{|C_{k\uparrow+}|^{2}}\Big{)}\frac{u_{0}v_{0}}{2}$
$\displaystyle+u_{0}u_{1}ie^{-i\theta_{k}}C_{k\uparrow+}^{*}C_{k\downarrow-}+v_{0}v_{1}(-i)e^{i\theta_{k}}C_{k\downarrow-}^{*}C_{k\uparrow+}\Bigg{]}.$

The last two terms are of the order of $\mathcal{O}(\alpha k_{F}/\Delta)$, so
it is convenient to write the gap function as follows

$\displaystyle\Delta(t)=\Delta_{h}[h(t)]+\delta\Delta(t).$ (35)

We have identified two contributions that have significantly different
origins: $\Delta_{h}$ is defined by the amplitude of the LZSM tunneling and
depends on time only through the Zeeman field $h(t)$;
$\delta\Delta(t)\propto\mathcal{O}(\alpha k_{F}/\Delta)$ is defined by cross-
terms and reflects interference effects between QP wavefunctions caused by
LZSM transitions and depends on time explicitly.

If one neglects the small perturbation $\delta\Delta(t)$ in (35) then it
becomes possible to get a simplified self-consistency equation for
$\Delta_{h}[h(t)]$ from Eq. (34). In general form it reads

$\displaystyle\Delta_{h}=\Delta_{0}\exp\Bigg{(}\int_{1}^{h/\Delta_{h}}\frac{e^{-\delta_{\text{LZ}}/s^{2}}-1}{\sqrt{s^{2}-1}}ds\Bigg{)},$
(36)

where $\delta_{\text{LZ}}=\pi\alpha^{2}k_{F}^{2}/\gamma$. The numerically
integrated function $\Delta_{h}[h(t)]$ is shown in Fig. 4 and below we discuss
its behavior for different tunneling regimes.

Figure 4:  Functional dependence of the superconducting gap $\Delta_{h}$ on
the spin-splitting field from expression (36) for different values of
$\delta_{\text{LZ}}$. The dashed-dotted line separates two regions where
$\Delta_{h}\lessgtr h$. The dashed line shows the critical value of the field
$h(t)=\Delta_{0}$ and the red circle marks the point of change in the behavior
of the gap $\Delta_{h}$ in the region $\Delta_{h}>h$. After this point the
equilibrium solution $\Delta_{h}=\Delta_{0}$ should jump to one of the
solutions fixed by the parameter $\delta_{\text{LZ}}$.

(i) The value $\delta_{\text{LZ}}=0$ corresponds to zero RSOC effects
($\alpha=0$), so that the spectral branches do not change after crossing and
the trivial solution for the gap $\Delta_{h}=\Delta_{0}$ holds.

(ii) The limit of $\delta_{\text{LZ}}\ll 1$ means that
$\gamma\gtrsim\alpha^{2}k_{F}^{2}$ and the spectral branches intersect
nonadiabatically, or so rapidly that they do not feel the RSOC. The Landau-
Zener tunneling is suppressed and the gap has a weak dependence on the Zeeman
field at $h(t)>\Delta_{0}$:

$\displaystyle\Delta_{h}\approx\Delta_{0}\exp\Big{(}-\delta_{\text{LZ}}\frac{\sqrt{h^{2}(t)-\Delta_{h}^{2}}}{h(t)}\Big{)}.$

(iii) In the opposite limit of $\delta_{\text{LZ}}\gg 1$ with
$\gamma\ll\alpha^{2}k_{F}^{2}$ the QPs undergo strong spin-flip tunneling
during an almost adiabatic avoided crossing. This leads to the effective
formation of the triplet superconducting correlations (or related triplet
component of the anomalous Green function [39]) even for the small RSOC energy
$\alpha k_{F}/\Delta\ll 1$. Such dynamically generated correlations are
determined by the rate of field change $\gamma$ and their effect on the gap
can significantly exceed the static mixing of singlet-triplet pairs for
$\alpha\neq 0$ [32]. As a result, the singlet gap function (3) is suppressed
and the self-consistency equation reads

$\displaystyle\Delta_{h}\approx\begin{cases}\sqrt{\Delta_{0}(2h(t)-\Delta_{0})}\quad\text{for}\quad\Delta_{h}>h/\delta_{\text{LZ}},\\\
\Delta_{0}\exp\Big{(}-\delta_{\text{LZ}}\frac{\sqrt{h^{2}(t)-\Delta_{h}^{2}}}{h(t)}\Big{)}\frac{\exp(\sqrt{\delta_{\text{LZ}}}\sqrt{\delta_{\text{LZ}}-1})}{\sqrt{\delta_{\text{LZ}}}+\sqrt{\delta_{\text{LZ}}-1}}\\\
\quad\text{for}\quad\Delta_{h}<h/\delta_{\text{LZ}}.\end{cases}$

(iv) The critical value $\delta_{\text{LZ}}\rightarrow\infty$ corresponds to
the complete Landau-Zener spin-flip tunneling, so that there are no QPs at the
energies $E>0$. In this case we have restored the thermodynamically metastable
branch $\Delta_{h}\approx\sqrt{\Delta_{0}(2h(t)-\Delta_{0})}$ from well-known
static case [30].

The actual behavior of the gap in time must be determined by switching between
different branches of $\Delta_{h}[h]$ as the Zeeman field $h(t)$ increases.
The first solution, which is fixed by the initial condition
$\Delta_{h}(t=0)=\Delta_{0}$ holds until $h(t)=\Delta_{0}$, where $\Delta_{h}$
goes to another unique possible solution $\Delta_{h}[h]$ for a given
$\delta_{\text{LZ}}$ (see the red point and black dashed line in Fig. 4). The
question of the exact dynamics of the gap in the jump region is difficult,
because due to the rapid change in the $\Delta_{h}$, the decoupling of the
LZSM problem and self-consistency equation may not be guaranteed [Section
IV.2]. It is qualitatively expected that the jump at
$t\approx\Delta_{0}/\gamma$ should be smeared both by non-zero static
contribution of SOC to the gap (since the equilibrium gap value depends on
$\alpha$) and by the QP tunneling dynamics. At large times $h(t)\gg\Delta_{0}$
there are no transitions between the QP states ($p_{k}\rightarrow 1$), since
the splitting between the spectral branches becomes zero and therefore the gap
tends to the constant asymptotics $\Delta_{h}(\infty)$.

###  IV.4 QP interference effects

In addition to the dominating term $\Delta_{h}$, the gap equation (34) also
contains small rapidly oscillating term

$\displaystyle\delta\Delta(t)=\lambda\sum_{|\xi_{k}|<\sqrt{h^{2}-\Delta^{2}}}u_{0}u_{1}ie^{-i\theta}C_{k\uparrow+}^{*}C_{k\downarrow-}$
(37)
$\displaystyle+v_{0}v_{1}(-i)e^{i\theta}C_{k\downarrow-}^{*}C_{k\uparrow+},$

arising from the interference of the QP states which have experienced LZSM
transitions. It is obvious that in its structure this function resembles the
collective Higgs mode, which is excited in a natural way during the
redistribution of states in the QP spectrum. Let us look at it in more
details. Using the time-dependent coefficients (32-33) we obtain

$\displaystyle\frac{\delta\Delta(t)}{\lambda
N(0)}=\int_{-\sqrt{h^{2}-\Delta^{2}}}^{\sqrt{h^{2}-\Delta^{2}}}\sqrt{p_{k}}\sqrt{1-p_{k}}G(\xi,t)\cos(D_{k}(t))d\xi,$
(38)

where we introduce the dynamical phase
$D_{k}(t)=2\int_{t_{0}}^{t}(E_{0}-h(t))dt+\chi_{k}+\pi$ and the function
$G(\xi,t)=\text{sgn}(\alpha)(u_{0}u_{1}+v_{0}v_{1})$. The function $G(\xi)$ is
proportional to $\alpha k_{F}/\Delta$, which means that $\delta\Delta(t)$ is
parametrically small and can be considered against the background of the main
change in the gap $\Delta_{h}$ from the equation (36).

For the integral (38), it is easy to estimate the asymptotic behavior at large
times $\Delta_{0}/\gamma\ll t$. The dynamical phase is written as
$D_{k}(t)=-(E_{0}^{2}+\gamma^{2}t^{2})/\gamma+\chi_{k}+\pi+2E_{0}t$ for the
spin-splitting field $h(t)=\gamma t$. Here $2E_{0}t$ is a fast oscillating
term at $t\rightarrow\infty$ and one can use a stationary phase approximation
for the $\xi$-integration in Eq. (38) with the stationary phase point $\xi=0$.
Using Eq. (52) for $G(\xi=0,t)$, we find the asymptotic behavior of
$\delta\Delta(t)$:

$\displaystyle\delta\Delta(t)\approx\lambda
N(0)e^{-\delta_{\text{LZ}}/2}\sqrt{1-e^{-\delta_{\text{LZ}}}}\frac{|\alpha|k_{F}}{2(\gamma
t-\Delta_{h})}\sqrt{\frac{\pi\Delta_{h}}{t}}$ (39)
$\displaystyle\times\cos\Big{(}\frac{(\gamma
t-\Delta_{h})^{2}}{\gamma}+\frac{3\pi}{4}-\chi_{0}\Big{)},$

where $\Delta_{h}[h(t\rightarrow\infty)]$ from Eq. (36) is a constant
determined by $\delta_{\text{LZ}}$. The result obtained means that the
collective interference between the two QP states at each $\xi_{k}$ after LZSM
crossing behaves at large times as a modified Higgs mode. Due to linear
dependence $h(t)$, this mode has a modulated frequency and polynomial damping
law $\propto t^{-3/2}$ arising from the inhomogeneous broadening of the mode.
Note that for large times only the contribution from the point $\xi=0$
survives, so the amplitude of $\delta\Delta(t)$ does not depend on the number
of redistributed states in the QP spectrum.

If the linear growth of the spin-splitting field $h(t)$ stops at a certain
value $h_{f}>\Delta_{0}$ after the redistribution of some of the QP states,
then the accumulated dynamic phase $D_{k}(t)$ and the gap fluctuation will
depend only on this value $h_{f}$

$\displaystyle\delta\Delta(t)\approx\lambda
N(0)e^{-\delta_{\text{LZ}}/2}\sqrt{1-e^{-\delta_{\text{LZ}}}}\frac{|\alpha|k_{F}}{2(h_{f}-\Delta_{h})}$
(40)
$\displaystyle\times\sqrt{\frac{\pi\Delta_{h}}{t}}\cos\Big{(}2(h_{f}-\Delta_{h})t-\frac{3\pi}{4}-\frac{\Delta_{h}^{2}-h_{f}^{2}}{\gamma}+\chi_{0}\Big{)}.$

The specific spectral distortion occurring between two brancher
$E_{k\uparrow+}$ and $E_{k\downarrow-}$ during the Landau-Zener dynamics at
$h(t)<h_{f}$ acts as an initial perturbation for the gap function at
$t=h_{f}/\gamma$. The free gap dynamics at $t>h_{f}/\gamma$ resembles the
Higgs mode with $\delta\Delta(t)\propto\cos(2(h_{f}-\Delta_{h})t)/\sqrt{t}$ at
the frequency $\omega=2|\Delta_{h}-h_{f}|$ (or $\omega=\omega_{-}$ in our
previous notations) with the standard damping law. It is interesting, that the
amplitude of this mode proportional to $|\alpha|k_{F}$ instead of
$\alpha^{2}k_{F}^{2}$ as it is expected in the case of small perturbations
[Section (III)]. Such amplification is a direct consequence of the
intersection of two specific spectral branches and the subsequent non-
adiabatic dynamics. Thus, this mode turns out to be leading in comparison with
other nonadiabatic corrections arising due to the interaction of all QP
spectral branches. Note, that the method for calculating the self-consistency
equation developed in Section III can be combined with the Landau-Zener
problem (72) and all corrections can be computed within the perturbation
theory.

###  IV.5 Density of states and distribution function

Rearrangment of the spectrum as a result of the intersection of spectral
branches naturally leads to a change of the structure of the density of states
(DOS), that has become time dependent. Since the temporal evolution of the
spectrum is adiabatic except the small region where the crossing occurs one
can use the quasistatic description of the DOS. For the small RSOC the DOS for
one spin projection can be written in terms of Bogoliubov-de Gennes functions

$\displaystyle
N_{\uparrow}(E,t)\approx\sum_{k}\sum_{n=\uparrow+,\uparrow-}|u_{0}|^{2}\delta\Big{(}E-E_{kn}[h(t)]\Big{)}$
(41) $\displaystyle+|v_{0}|^{2}\delta\Big{(}E+E_{kn}[h(t)]\Big{)}.$

Here we use static QP amplitudes $u_{0}$ and $v_{0}$ (see Eq. 51) to
distinguish the particle/hole contributions and $E_{k\uparrow\pm}$ are defined
in Eq. (4). The calculation of $N_{\uparrow}$ is cumbersome, because the RSOC
shifts the spectral branches and opens a minigap $\propto\alpha k_{F}$ at
$E=0$ [Appendix E]. For the small RSOC parameter these changes are negligible
and one can use a standard expression for the DOS, which now depends on time
through the spin-splitting field

$\displaystyle\frac{N_{\uparrow}(E,t)}{N(0)}\approx\frac{|E+h(t)|}{\sqrt{(E+h(t))^{2}-\Delta_{h}[h(t)]^{2}}}.$
(42)

Here the gap function $\Delta_{h}$ is taken from (36) and two coherence peaks
are present at $E=\pm\Delta_{h}[h(t)]-h(t)$.

The amplitude of the QP wavefunction $\psi_{k}(t)$ from (26) contains the
information about filling (or occupation) of the $\xi_{k}-$th state. More
precisely the coefficients $|C_{k\uparrow\pm}(t)|^{2}$ and
$|C_{k\downarrow\pm}(t)|^{2}$ can serve as an effective distribution functions
$f_{\uparrow\downarrow}(E)$ for QPs with different spin projections. As it was
discussed in the section IV.3, the temporal evolution of these coefficients is
defined by the LZSM problem, and for spin-up states one has

$\displaystyle|C_{k\uparrow-}(t)|^{2}$ $\displaystyle=1,$
$\displaystyle|C_{k\uparrow+}(t)|^{2}$
$\displaystyle=(1-p_{k})\Theta\big{[}\sqrt{h^{2}(t)-\Delta_{h}^{2}}-\xi_{k}\big{]},$

which can be rewritten as a distribution function

$\displaystyle f_{\uparrow}(E,t)\approx\begin{cases}0,\quad E>0\\\
1-\exp\Big{[}-\frac{\delta_{\text{LZ}}\Delta_{h}^{2}[h(t)]}{(E+h(t))^{2}}\Big{]},\quad\Delta_{h}-h<E<0\\\
1,\quad E<\Delta_{h}-h\end{cases}$ (43)

The dependence $f_{\uparrow}(E)$ is shown in Fig. 4 for
$\delta_{\text{LZ}}=0.5$. The most pronounced change of the distribution
function occurs at $E\approx\Delta_{h}-h$, since for large QP energies the
LZSM tunneling is suppressed. For the opposite spin projection the DOS
$N_{\downarrow}(E)$ has the similar structure (42) with $h\rightarrow-h$,
while the corresponfing distribution function $f_{\downarrow}(E)$ is different
and is given by the Eqs. (32-33). The DOS structure and effective distribution
function enable the calculation of a system’s optical or transport response,
which can be experimentally measured.

###  IV.6 Dynamical magnetization of QP gas

Nonadiabatic LZSM tunneling of QP states causes a spin imbalance in the
spectrum, which results in the appearance of nonzero dynamical magnetization.
Using the notations from the previous section we get an expression for the
z-component of the magnetization per unit volume

$\displaystyle
m_{z}(t)=\mu_{B}\sum_{\text{i.c.}}\check{\psi}^{\dagger}_{k}(t)\check{\tau}_{m}\check{\psi}_{k}(t),$
(44)

where
$\check{\tau}_{m}=(\hat{\tau}_{0}+\hat{\tau}_{z})\otimes\hat{\sigma}_{z}/2$;
the vector $\check{\psi}_{k}(t)$ is a solution of the TDBdG problem (26) and
"i.c." here means the summation over all initial conditions (see Eq. (3)). Due
to symmetry and homogenuty of the problem for the field ${\bf h}(t)=h(t){\bf
z}_{0}$ the transversal components of the magnetization $m_{x,y}(t)$ are zero.

Taking the dynamical amplitudes $C_{n}(t)$ from (32-33) and implementing the
same procedure as for the self-consistency equation (34-35) we found that the
magnetization can be written as follows

$\displaystyle m_{z}(t)=m_{h}[h(t)]+\delta m(t).$ (45)

As in the case of the gap equation (35) we have two contributions:
$m_{h}[h(t)]$ which is a slow function of time arising from the redistribution
of the quaiparticle states, and $\delta m(t)\propto\alpha k_{F}$ which is
small oscillatory term originated from the interference of the redistributed
states. The first term can be easily calculated with the help of the
quasiparticle density

$\displaystyle n_{\sigma}(t)=\int N_{\sigma}(E,t)f_{\sigma}(E,t)dE,$ (46)

where $\sigma=\uparrow\downarrow$ and the DOS $N_{\sigma}(E)$ and distribution
function $f_{\sigma}(E)$ are defined in the previous subsection. Corresponding
spin imbalance results in the dynamical magnetization
$m_{h}[h(t)]=\mu_{B}(n_{\uparrow}-n_{\downarrow})$, which is shown in Fig.
5(a).

For $h(t)<\Delta_{0}$ there is no crossing of the QP spectral branches and
according to our model there is no tunneling between QP states, therefore the
dynamical magnetisation is zero. Once the intersection has occured at
$h(t)=\Delta_{0}$, the distribution functions $f_{\uparrow,\downarrow}(E)$
transform and nonzero spin imbalance $n_{\uparrow}-n_{\downarrow}$ is
generated. Due to the jump of $\Delta_{h}$ function at $h(t)=\Delta_{0}$ [Fig.
4] the magnetization at this point also has a sharp discontinuity. At large
times the tunneling of QP states is suppressed therefore the magnetization is
saturated to a constant value determined by the parameter
$\delta_{\text{LZ}}$. Obviously, an increase in $\delta_{\text{LZ}}$ makes the
spin-flip tunneling more efficient and thereby increases the maximum value of
$m_{h}$.

Figure 5:  (a) Dynamical magnetization $m_{h}$ per unit volume induced by the
nonadiabatic tunneling of QP states and (b) corresponding spin susceptibility
$\chi^{\text{sp}}$ versus time-dependent spin-splitting field $h(t)$ for
different values of $\delta_{\text{LZ}}$.

The second therm in (45) resembles the Higgs mode term (37) and gives
negligible contribution to $m_{z}(t)$, therefore it can be discarded.

In addition one can compute the dynamical susceptibility of the QP gas in the
Zeeman field of the general form $h(t)=\mu_{\text{B}}H(t)$. It is known that
an orbital and a spin parts of the magnetic susceptibility can be splitted in
the case of small spin-orbital effects [39]. Since we consider a homogeneous
system and neglect all orbitals effects only the spin part plays a role, which
can be written as follows

$\displaystyle\chi^{\text{sp}}[h(t)]=\mu_{B}\frac{\partial m_{h}}{\partial
h}.$ (47)

The ratio of the numerically calculated susceptibility
$\chi^{\text{sp}}[h(t)]$ and the normal susceptibility
$\chi_{N}^{\text{sp}}=2\mu_{B}^{2}N(0)$ [40] is shown if Fig. 5(b). It is seen
that spin-flip tunneling in the QP spectrum provokes a paramagnetic response
of the superconducting condensate. The function (47) should have a singularity
$\chi^{\text{sp}}\propto(h(t)-\Delta_{0})^{-1/2}$ in the vicinity of
$h(t)\approx\Delta_{0}$, which is defined by the shape of the QP spectrum at
$k\approx k_{F}$ and has the same origin as the coherence peak in the DOS
(42). However, due to the jump of the order parameter $\Delta_{h}$ at this
point we observe shifted peaks, which have to be smeared out near
$h(t)=\Delta_{0}$ if more realistic model of LZSM tunneling [Section IV.2] is
taken into account. We note again that we discuss only the dynamic
contribution to the susceptibility, which, generally speaking, have to be
added to the static one, which is not equal to zero at $T=0$ in the presence
of SOC [39, 40].

##  V Discussion and experimental perspectives

We analyzed the coherent dynamics of the superconducting condensate in the
presence of Zeeman field and SOC in collisionless regime. First, it was
established that the Higgs mode of the superconducting gap is sensitive to the
spin-splitting field $h_{0}$ and can be directly triggered by either its
harmonic perturbation $\delta h(t)$ or by an external laser pulse. Second, it
was shown, that the field $h(t)\sim t$ can provoke an avoided crossing of the
QP spectral branches and adiabatic spin-flip tunneling of the QPs between the
different branches occurs. Corresponding redistribution of the QPs in the
spectrum leads to the appearance of the dependence $\Delta[h(t)]$ and
generation of the interference effects. Emerging spin imbalance reveals itself
in the effective dynamical distribution function and in a generation of a weak
magnetization of the QP gas.

We propose superconductor-ferromagnet hybrid structures as an experimental
platform for detecting the described effects. The ferromagnetic layer can
serve as a source of both Rashba spin-orbit coupling and an exchange field.
Since it is important to remove orbital effects from the system, the most
suitable geometry for superconductor is either thin film or one-dimensional
nanowire [41].

The excitation and observation of Higgs modes in superconductors requires
frequencies of the order of $\Delta_{0}/\hbar$, which vary from the far
infrared to the terahertz range. The laser excitation of modes seems to be the
most practical and feasible, and the detection can be implemented using the
THz light source with ultrafast pump-probe spectroscopy or third harmonic
generation measurements [10, 12, 11]. Generation of a fast oscillating
component of the homogeneous Zeeman field $\delta h(t)$ inside a
superconductor is a difficult task especially for the THz range. Some
proposals can be made amid encouraging progress in the ultrafast optical
control of magnetization in various materials [42, 43, 44]. Ferromagnetic
resonance induces the time-dependent stray field which in combination with
geometric constraints may serve as a Zeeman field inside a thin
superconducting film, as it was discussed for a two-dimensional electron gas
in Ref. [45]. Upon the excitation of the Higgs modes by the Zeeman field the
THz spectroscopy measurements can be implemented again.

It is possible to make basic parameter estimates for the experimental
observation of LZSM transitions in the QP spectrum. For example, consider
$\Delta_{0}=0.1$meV (for $T_{c}\approx 1$K) and $\alpha k_{F}\sim 10^{-3}$meV
$\ll\Delta_{0}$. Then the constraint for the small tunneling rate is
$\hbar\gamma\gtrsim\alpha^{2}k_{F}^{2}$ in dimensional units, which is
equivalent to $\gamma\gtrsim 10^{-3}$meV/ns. Consider inelastic relaxation of
QP with a typical time $\tau_{\text{ph}}\sim 100$ns in the case of the
electron-phonon scattering at low temperatures [46, 47] . The collisionless
regime is maintained at $t\ll\tau_{\text{ph}}$, which corresponds to times
$t\lesssim 10$ns. Under such conditions the field $h(t)=\gamma
t\sim\Delta_{0}$ is achievable only for $\gamma\sim 10^{-2}$meV/ns. The
measurements of the various properties of the superconducting condensate above
Pauli limit at short times can be implemented with the help of ultrafast THz
techniques, such as pump-probe [12] for an optical conductivity.

###### Acknowledgements.

This work has been supported by ANR OPTOFLUXONICS, ANR SUPERFAST, the LIGHT
S$\&$T Graduate Program and the Russian Science Foundation (Grant No.
21-72-10161). A.S.M. acknowledges support from the State Contract of Ministry
of Science and Higher Education of Russian Federation No. 075-03-2022-106
(project FSMG-2023-0011) of Moscow Institute of Physics and Technology.

##  Appendix A Eigenvectors

The instantaneous eigenvectors of the Hamiltonian $\check{\mathcal{H}}(k,t)$
from Eq. (2) can be written as follows

$\displaystyle\check{\Psi}_{kn}(t)=\frac{1}{\sqrt{1+a_{1n}^{2}+a_{2n}^{2}+a_{3n}^{2}}}\begin{pmatrix}1\\\
-ia_{1n}e^{i\theta_{k}}\\\ -ia_{2n}e^{i\theta_{k}}\\\ a_{3n}\end{pmatrix},$
(48)

where we have defined the phase $\theta_{k}=\arg\big{(}k_{x}+ik_{y}\big{)}$
and real coefficients

$\displaystyle a_{1n}=\frac{(h+E_{kn})^{2}-E_{0}^{2}-\alpha^{2}k^{2}}{2\alpha
k(\xi_{k}+h)},$ (49) $\displaystyle a_{2n}=\frac{\alpha
k}{\Delta}-\frac{E_{kn}-\xi_{k}-h}{\Delta}a_{1n},$ $\displaystyle
a_{3n}=\frac{E_{kn}-\xi_{k}+h}{\Delta}-\frac{\alpha k}{\Delta}a_{1n}.$

The instantaneous eigenvalues of $\check{\mathcal{H}}(k,t)$ are

$\displaystyle E_{kn}(t)\equiv E_{k\sigma\pm}(t)=$
$\displaystyle\pm\sqrt{E_{0}^{2}+\alpha^{2}k^{2}+h^{2}(t)\mp\text{sgn}(\sigma)2\sqrt{\xi_{k}^{2}\alpha^{2}k^{2}+h^{2}(t)E_{0}^{2}}},$

where $E_{0}=\sqrt{\xi_{k}^{2}+\Delta^{2}}$; the subscript $\pm$ refers to
spectral branch above/below the Fermi level and
$\sigma=\\{\uparrow,\downarrow\\}$ denotes a spin subband. Note, that the
Hamiltonian (2) implies the symmetry relations between the energies
$E_{k\uparrow+}=-E_{k\downarrow-}$ and $E_{k\downarrow+}=-E_{k\uparrow-}$, and
between the corresponding eigenvectors
$\check{\Psi}_{k\uparrow+}=i\hat{\tau}_{y}\otimes\hat{\sigma}_{z}\check{\Psi}_{k\downarrow-}^{*}$
and
$\check{\Psi}_{k\downarrow+}=i\hat{\tau}_{y}\otimes\hat{\sigma}_{z}\check{\Psi}_{k\uparrow-}^{*}$,
where $\hat{\tau}_{i}(\hat{\sigma}_{i})$ is the Pauli matrix in the
Nambu(spin) space.

For the case of weak SOC $\alpha k_{F}\ll\\{E_{F},h(t),\Delta(t)\\}$ the
eigenvectors (48) can be expanded up to the first order in $\alpha
k_{F}/\Delta$ as follows

$\displaystyle\check{\Psi}_{k\uparrow+}\approx\begin{pmatrix}u_{0}\\\
-iu_{1}e^{i\theta_{k}}\\\ -iv_{1}e^{i\theta_{k}}\\\
v_{0}\end{pmatrix},\quad\check{\Psi}_{k\downarrow+}\approx\begin{pmatrix}iu_{1}e^{-i\theta_{k}}\\\
-u_{0}\\\ v_{0}\\\ -iv_{1}e^{-i\theta_{k}}\end{pmatrix},$ (50)
$\displaystyle\check{\Psi}_{k\uparrow-}\approx\begin{pmatrix}-v_{0}\\\
iv_{1}e^{i\theta_{k}}\\\ -iu_{1}e^{i\theta_{k}}\\\ u_{0}\\\
\end{pmatrix},\quad\check{\Psi}_{k\downarrow-}\approx\begin{pmatrix}-iv_{1}e^{-i\theta_{k}}\\\
v_{0}\\\ u_{0}\\\ -iu_{1}e^{-i\theta_{k}}\\\ \end{pmatrix}.$

Here we define equilibrium QP amplitudes

$\displaystyle
u_{0}=\frac{1}{\sqrt{2}}\sqrt{1+\frac{\xi_{k}}{\sqrt{\xi_{k}^{2}+\Delta^{2}}}},\quad
v_{0}=\frac{1}{\sqrt{2}}\sqrt{1-\frac{\xi_{k}}{\sqrt{\xi_{k}^{2}+\Delta^{2}}}},$
(51)

and $u_{1},v_{1}\propto\mathcal{O}(\alpha k_{F}/\Delta)$ correspond to the
triplet component of the QP wavefunctions.

The function $G(\xi,t)=\text{sgn}(\alpha)(u_{0}u_{1}+v_{0}v_{1})$ from Eq.
(38) can be found at the stationary phase point $\xi=0$ by expanding the
coefficients (49) in the small parameter $\alpha k_{F}/\Delta$ and putting
$\Delta\approx\Delta_{h}$ (see Eq. (36)). The result reads

$\displaystyle G(0,t)\approx\frac{|\alpha|k_{F}}{2(h(t)-\Delta_{h})}.$ (52)

##  Appendix B Derivation of linearized self-consistency equation

We start with the linearized (8, 11) dynamical equations

$\displaystyle i\frac{\partial}{\partial t}\delta
C_{km}(t)=\sum_{n}\check{\Psi}_{m}^{\dagger}\check{\mathcal{V}}(t)\check{\Psi}_{n}e^{-i(E_{n}-E_{m})t}\big{(}\delta_{n,l}+\delta
C_{kn}(t)\big{)},$ (53)

where the indices $n,m=\\{\uparrow+,\downarrow+,\uparrow-,\downarrow-\\}$ and
two possible initial configurations (10) are marked as
${l=\\{\uparrow-,\downarrow-\\}}$. The compact form of the self-consistency
equation for the gap (9) is

$\displaystyle\Delta_{\text{eq}}+\delta\Delta(t)=$ (54)
$\displaystyle-\frac{\lambda}{2}\sum_{l}\sum_{k}\sum_{n,n^{\prime}}\big{(}\delta_{n,l}+\delta
C_{kn}(t)^{*}\big{)}\big{(}\delta_{n^{\prime},l}+\delta
C_{kn^{\prime}}(t)\big{)}$ $\displaystyle\times
e^{-i(E_{n^{\prime}}-E_{n})t}\check{\Psi}_{kn}^{\dagger}\check{\tau}_{\Delta}\check{\Psi}_{kn^{\prime}}.$

As was mentioned in Section II, we neglect the effect of RSOC on the
equilibrium value of the gap, which can be taken as
$\Delta_{\text{eq}}=\Delta_{0}$. It also makes sense to omit the negligibly
small corrections from the RSOC to the energy spectrum, so one can put
$E_{n}\equiv E_{k\sigma\pm}\approx\pm E_{0}-\text{sgn}(\sigma)h_{0}$.

The equations (53, 54) can be simplified and written as follows

$\displaystyle\frac{\partial f_{1}}{\partial
t}=ie^{i(2(E_{0}-h_{0})t)}[\mathcal{A}\delta\Delta(t)-\mathcal{B}\delta
h(t)],$ (55) $\displaystyle\frac{\partial f_{2}}{\partial
t}=ie^{i(2(E_{0}+h_{0})t)}[\mathcal{A}\delta\Delta(t)+\mathcal{B}\delta
h(t)],$ $\displaystyle\frac{\partial g}{\partial
t}=i\frac{\xi}{E_{0}}e^{i(2E_{0}t)}\delta\Delta(t),$
$\displaystyle\delta\Delta(t)=\Big{\langle}\frac{\mathcal{A}}{2}\text{Re}f_{1}(t)e^{-i(2(E_{0}-h_{0})t)}\Big{\rangle}$
$\displaystyle+\Big{\langle}\frac{\mathcal{A}}{2}\text{Re}f_{2}(t)e^{-i(2(E_{0}+h_{0})t)}\Big{\rangle}+\Big{\langle}\frac{\xi}{E_{0}}\text{Re}g(t)e^{-i(2E_{0}t)}\Big{\rangle}.$

We have used the notation
$\big{\langle}\dots\big{\rangle}=\lambda\sum_{k}\approx\lambda
N(0)\int_{-\omega_{D}}^{\omega_{D}}d\xi$ and introduced new complex-valued
functions

$\displaystyle f_{1}\equiv-ie^{i\theta}\delta C_{\uparrow+},\quad
g_{1}\equiv-\delta C_{\downarrow+},$ (56) $\displaystyle f_{2}\equiv-
ie^{-i\theta}\delta C_{\downarrow+},\quad g_{2}\equiv-\delta C_{\uparrow+},$
$\displaystyle g=\frac{g_{1}+g_{2}}{2},$

where the subscript corresponds to the two possible initial conditions. The
functions

$\displaystyle\mathcal{A}(\xi)$
$\displaystyle=2(u_{0}u_{1}+v_{0}v_{1})\propto\mathcal{O}(\alpha
k_{F}/\Delta_{0}),$ (57) $\displaystyle\mathcal{B}(\xi)$
$\displaystyle=2(u_{0}v_{1}+u_{1}v_{0})\propto\mathcal{O}(\alpha
k_{F}/\Delta_{0})$

have the lowest order in $\alpha k_{F}$ parameter [Appendix A] and are even in
$\xi$. All terms odd in $\xi$ in Eq. (55) are related to the imaginary part of
$\delta\Delta(t)$ and vanish due to the approximate electron-hole symmetry of
BdG Hamiltonian (2), due to which the density of states is approximated as
$N(\xi)\approx N(0)$ in the $\langle...\rangle$-integration [48].

Applying the Laplace transform $f(s)=\int_{0}^{\infty}e^{-st}f(t)dt$ with
$s=i\omega+\zeta$ (where $\zeta\rightarrow 0$) for Eq. (55) we obtain the gap
equation in the complex plane, which is found to be

$\displaystyle\delta\Delta(s)=\delta\Delta(s)\Big{\langle}\frac{2\xi^{2}}{E_{0}}\frac{1}{s^{2}+4E^{2}_{0}}\Big{\rangle}+\delta\Delta(s)\Big{\langle}\mathcal{A}^{2}(\xi)\frac{(E_{0}+h_{0})}{s^{2}+4(E_{0}+h_{0})^{2}}\Big{\rangle}+\delta\Delta(s)\Big{\langle}\mathcal{A}^{2}(\xi)\frac{(E_{0}-h_{0})}{s^{2}+4(E_{0}-h_{0})^{2}}\Big{\rangle}$
(58) $\displaystyle+\delta
h(s)\Big{\langle}\mathcal{A}(\xi)\mathcal{B}(\xi)\frac{(E_{0}+h_{0})}{s^{2}+4(E_{0}+h_{0})^{2}}\Big{\rangle}-\delta
h(s)\Big{\langle}\mathcal{A}(\xi)\mathcal{B}(\xi)\frac{(E_{0}-h_{0})}{s^{2}+4(E_{0}-h_{0})^{2}}\Big{\rangle}$
$\displaystyle+\Big{\langle}f_{1}^{\prime}(0)\frac{\mathcal{A}(\xi)}{2}\frac{s}{s^{2}+4(E_{0}-h_{0})^{2}}\Big{\rangle}+\Big{\langle}f_{1}^{\prime\prime}(0)\frac{\mathcal{A}(\xi)(E_{0}-h_{0})}{s^{2}+4(E_{0}-h_{0})^{2}}\Big{\rangle}$
$\displaystyle+\Big{\langle}f_{2}^{\prime}(0)\frac{\mathcal{A}(\xi)}{2}\frac{s}{s^{2}+4(E_{0}+h_{0})^{2}}\Big{\rangle}+\Big{\langle}f_{2}^{\prime\prime}(0)\frac{\mathcal{A}(\xi)(E_{0}+h_{0})}{s^{2}+4(E_{0}+h_{0})^{2}}\Big{\rangle}$
$\displaystyle+\Big{\langle}g^{\prime}(0)\frac{\xi}{E_{0}}\frac{s}{s^{2}+4E^{2}_{0}}\Big{\rangle}+\Big{\langle}g^{\prime\prime}(0)\frac{2\xi
E_{0}}{s^{2}+4E^{2}_{0}}\Big{\rangle}.$

Here $f=f^{\prime}+if^{\prime\prime}$ and the initial conditions
$f_{1,2}(0)=f_{1,2}(t=0)$, $g(0)=g(t=0)$ implicitly contain the initial value
of the gap perturbation $\delta\Delta(t=0)$. Note that terms with initial
conditions will be discarded when calculating the superconductor response (see
Section III). Now we can single out functions of $s$ with different
singularities in the complex plane and denote them using short notations

$\displaystyle\mathcal{K}_{0}(s)=\Big{\langle}\frac{2\xi^{2}}{E_{0}}\frac{1}{s^{2}+4E^{2}_{0}}\Big{\rangle},$
(59)
$\displaystyle\mathcal{K}_{\pm}(s)=\Big{\langle}\mathcal{A}^{2}(\xi)\frac{(E_{0}\pm
h_{0})}{s^{2}+4(E_{0}\pm h_{0})^{2}}\Big{\rangle},$
$\displaystyle\mathcal{F}_{\pm}(s)=\Big{\langle}\mathcal{A}(\xi)\mathcal{B}(\xi)\frac{(E_{0}\pm
h_{0})}{s^{2}+4(E_{0}\pm h_{0})^{2}}\Big{\rangle}$

and $\mathcal{I}(s)$, which consists of all initial perturbations at $t=0$.
The functions $\mathcal{A}$ and $\mathcal{B}$ are of the first order in the
small parameter $\alpha k_{F}/\Delta$, therefore we have

$\displaystyle\mathcal{K}_{0}(s)\propto\mathcal{O}\Big{(}\frac{\alpha^{0}k^{0}_{F}}{\Delta_{0}^{0}}\Big{)},$
$\displaystyle\mathcal{K}_{\pm}(s),\mathcal{F}_{\pm}(s)\propto\mathcal{O}\Big{(}\frac{\alpha^{2}k^{2}_{F}}{\Delta_{0}^{2}}\Big{)}.$

It can be shown that the the difference
$[\mathcal{F}_{+}(s)-\mathcal{F}_{-}(s)]$ is proportional to $h_{0}$. This
allows one to write the terms with $\delta h(t)$ in (58) as $\delta
h(s)[\mathcal{F}_{+}(s)-\mathcal{F}_{-}(s)]$ or $({\bf h}_{0}\cdot\delta{\bf
h}(s))[\mathcal{F}_{+}(s)-\mathcal{F}_{-}(s)]/h_{0}$, where both vectors are
oriented along ${\bf z}_{0}$ axis. By rewriting the equation (58) with the new
introduced functions (59) we get the self-consistency equation (12).

##  Appendix C Long-time behavior of $\delta\Delta(t)$

The susceptibility $\text{Im}\chi_{\Delta h}(s)|_{\zeta\rightarrow
0}=\text{Im}\chi_{\Delta h}(\omega)$ in Eq. (24) has strongly dominant terms
in the vicinity of different branch points in the interval
$\omega\in[\omega_{-},\infty)$. In order to demonstrate this, the function
$\text{Im}\chi_{\Delta h}(s)$ can be expanded in a series up to the second
order in the parameter $\alpha k_{F}/\Delta$, and this expansion must be
carried out accurately near the branch points and may differ in different
regions of $\omega$. Therefore, we assume that the value of the integral is
determined by these dominant contributions of $\text{Im}\chi_{\Delta
h}(\omega)$ and can be evaluated sequentially as
$\int_{\omega_{-}}^{\infty}=\int_{\omega_{-}}^{\omega_{0}}+\int_{\omega_{0}}^{\omega_{+}}+\int_{\omega_{+}}^{\infty}.$
Let us consider the small regions $\Omega\ll\omega_{0,\pm}$ in the vicinity of
these points separately.

1: $\omega\approx\omega_{-}+\Omega$

Close to the point $\omega=\omega_{-}$ the term
$\mathcal{F}_{-}^{\prime\prime}(\omega)$ dominates:

$\displaystyle\mathcal{F}_{-}^{\prime\prime}(\Omega)\approx-\lambda
N(0)\frac{\pi\Delta_{0}\mathcal{A}(0)\mathcal{B}(0)}{4\sqrt{2\Delta_{0}\Omega}}\propto\frac{1}{\sqrt{\Omega}}.$
(60)

Despite the kernel $1-\mathcal{K}_{0}^{\prime}(\omega)$ goes to zero at
$\omega\rightarrow\omega_{0}$ there is no singularity in $\chi_{\Delta
h}(\omega)$ at this point due to the small terms of the order of $(\alpha
k_{F})^{2}$ in the denominator. Therefore, the region in the vicinity of
$\omega_{0}$ will not contrubute to the intergal. Thus, the behavior of the
first integral for $\omega\in[\omega_{-},\omega_{0})$ at large time $h_{0}t\gg
1$ can be estimated as follows

$\displaystyle\int_{\omega_{-}}^{\omega_{0}}\approx\text{Im}\Bigg{[}\frac{\delta
h(i\omega_{-})e^{i\omega_{-}t}}{\big{[}1-\mathcal{K}_{0}^{\prime}(\omega_{-})\big{]}}\int_{0}^{\omega_{0}-\omega_{-}}\mathcal{F}_{-}^{\prime\prime}(\Omega)e^{i\Omega
t}d\Omega\Bigg{]}$ (61) $\displaystyle\approx-\lambda
N(0)\frac{\pi^{3/2}\Delta_{0}\mathcal{A}(0)\mathcal{B}(0)}{4\sqrt{2\Delta_{0}t}}\frac{\text{Im}\Big{[}\delta
h(i\omega_{-})e^{i(\omega_{-}t+\pi/4)}\Big{]}}{\big{[}1-\mathcal{K}_{0}^{\prime}(\omega_{-})\big{]}}.$

2: $\omega\approx\omega_{0}+\Omega$

In the vicinity of the branch point $\omega=\omega_{0}$ the main contribution
is defined by

$\displaystyle\mathcal{K}_{0}^{\prime\prime}(\Omega)\approx-\lambda
N(0)\frac{\pi}{2\Delta_{0}}\sqrt{\Delta_{0}\Omega}\propto\sqrt{\Omega}.$ (62)

Thus at large times $h_{0}t\gg 1$ we get

$\displaystyle\int_{\omega_{0}}^{\omega_{+}}=\int_{\omega_{0}}^{\omega_{+}}\frac{\big{[}\mathcal{F}_{+}^{\prime}(\omega)-\mathcal{F}_{-}^{\prime}(\omega)\big{]}}{\mathcal{K}_{0}^{\prime\prime}(\omega)}\text{Im}\Big{[}e^{i\omega
t}\delta h(i\omega)\Big{]}d\omega$ (63)
$\displaystyle\approx-\frac{2\sqrt{\Delta_{0}}}{\sqrt{\pi
t}}\frac{\big{[}\mathcal{F}_{+}^{\prime}(\omega_{0})-\mathcal{F}_{-}^{\prime}(\omega_{0})\big{]}}{\lambda
N(0)}\text{Im}\Big{[}\delta h(i\omega_{0})e^{i(\omega_{0}t+\pi/4)}\Big{]}.$

3: $\omega\approx\omega_{+}+\Omega$

For the last branch point $\omega=\omega_{+}$ the kernel
$\mathcal{F}_{+}^{\prime\prime}(\omega)$ dominates:

$\displaystyle\mathcal{F}_{+}^{\prime\prime}(\Omega)\approx\lambda
N(0)\frac{\pi\Delta_{0}\mathcal{A}(0)\mathcal{B}(0)}{4h_{0}\sqrt{2\Delta_{0}\Omega}}\propto\frac{1}{\sqrt{\Omega}}.$
(64)

At large times $h_{0}t\gg 1$ we get

$\displaystyle\int_{\omega_{+}}^{\infty}\approx\int_{\omega_{+}}^{\infty}-\frac{\big{[}1-\mathcal{K}_{0}^{\prime}(\omega_{+})\big{]}\mathcal{F}_{+}^{\prime\prime}(\omega)}{\big{[}1-\mathcal{K}_{0}^{\prime}(\omega_{+})\big{]}^{2}+\big{[}\mathcal{K}_{0}^{\prime\prime}(\omega_{+})\big{]}^{2}}$
(65) $\displaystyle\times\text{Im}\Big{[}e^{i\omega t}\delta
h(i\omega)\big{)}\Big{]}d\omega$ $\displaystyle\approx-\lambda
N(0)\frac{\pi^{3/2}\Delta_{0}\mathcal{A}(0)\mathcal{B}(0)}{4\sqrt{2\Delta_{0}t}}\frac{\big{[}1-\mathcal{K}_{0}^{\prime}(\omega_{+})\big{]}}{\big{[}1-\mathcal{K}_{0}^{\prime}(\omega_{+})\big{]}^{2}+\big{[}\mathcal{K}_{0}^{\prime\prime}(\omega_{+})\big{]}^{2}}$
$\displaystyle\times\text{Im}\Big{[}\delta
h(i\omega_{+})\big{)}e^{i(\omega_{+}t+\pi/4)}\Big{]}.$

4: total integral

By combining all three contribution (61,63,65) we will get the equation (25)
in the main text. Note that discussed approximations work for
$0<h_{0}<\Delta_{0}$. The functions $\mathcal{A},\mathcal{B}$ from (57) at the
point $\xi=0$ can be calculated using the wavefunctions (50). By expanding the
coefficients (49) we obtain

$\displaystyle\mathcal{A}(0)\mathcal{B}(0)=\mathcal{A}^{2}(0)\approx\frac{(\alpha
k_{F})^{2}}{(\Delta_{0}-h_{0})^{2}}.$ (66)

Also, the analytical expressions for the kernel $\mathcal{K}_{0}(\omega)$ at
$\omega>0$ can be found:

$\displaystyle\frac{1-\mathcal{K}_{0}^{\prime}(\omega)}{\lambda N(0)}=$ (67)
$\displaystyle\begin{cases}\frac{\sqrt{4\Delta_{0}^{2}-\omega^{2}}}{\omega}\arctan\Big{(}\frac{\omega}{\sqrt{4\Delta_{0}^{2}-\omega^{2}}}\Big{)}\quad\text{for}\quad\omega<2\Delta_{0}\\\
-\frac{\sqrt{\omega^{2}-4\Delta_{0}^{2}}}{\omega}\frac{1}{2}\ln\Big{(}\frac{\omega-\sqrt{\omega^{2}-4\Delta_{0}^{2}}}{\omega+\sqrt{\omega^{2}-4\Delta_{0}^{2}}}\Big{)}\quad\text{for}\quad\omega>2\Delta_{0}\end{cases},$
$\displaystyle\frac{\mathcal{K}_{0}^{\prime\prime}(\omega)}{\lambda
N(0)}=-\frac{\pi}{2}\frac{\sqrt{\omega^{2}-4\Delta_{0}^{2}}}{\omega}\Theta[\omega-2\Delta_{0}].$
(68)

Finally, the expression with the kernels $\mathcal{F}_{\pm}(\omega)$ from (63)
can be calculated numerically for small $\alpha k_{F}\ll\Delta_{0}$:

$\displaystyle\frac{\big{[}\mathcal{F}_{+}^{\prime}(\omega_{0})-\mathcal{F}_{-}^{\prime}(\omega_{0})\big{]}}{\lambda
N(0)}$ (69)
$\displaystyle=h_{0}\fint_{0}^{\omega_{D}}\mathcal{A}(\xi)\mathcal{B}(\xi)\frac{h_{0}^{2}-\xi^{2}-2\Delta_{0}^{2}}{(\xi^{2}-h_{0}^{2})^{2}-4\Delta_{0}^{2}h_{0}^{2}}d\xi.$

##  Appendix D Derivation and solution of LZSM problem

The dynamics of two levels with avoided crossing can be simply described with
the help of so-called diabatic basis formed by the instantaneous eigenfunction
$\check{\Phi}^{0}_{kn}(t)$ of the time-dependent Hamiltonian (1) at
$\alpha=0$. Note here that for $\alpha=0$ the eigenstates do not depend of
$h(t)$ at all and consist only of the Bogoliubov’s amplitudes $u_{0}$ and
$v_{0}$ (one can use (50) and put $\alpha=0$ there). The complete solution of
the time-dependent hamiltonian can be written as follows

$\displaystyle\check{\Psi}_{k}(t)=\sum_{n}C^{d}_{kn}(t)\check{\Phi}^{0}_{kn}(t),$
(70)

where $n=\\{\uparrow+,\downarrow+,\uparrow-,\downarrow-\\}$. In order to avoid
confusion with adiabatic basis in (26) the superscript "d" is used to denote
the diabatic basis. The time-dependent coefficients obey the following
equation derived from (1):

$\displaystyle i\frac{\partial}{\partial
t}C^{d}_{m}=\sum_{n}C^{d}_{n}\check{\Phi}^{0\dagger}_{km}\Big{[}\check{\mathcal{H}}(t)-i\frac{\partial}{\partial
t}\Big{]}\check{\Phi}^{0}_{kn}.$ (71)

Note that here $\check{\mathcal{H}}(t)\check{\Phi}^{0}_{kn}\neq
E_{n}(t)\check{\Phi}^{0}_{kn}$. By keeping in mind that
$\hat{\Phi}^{0}_{kn}(t)$ depends on time only through $\Delta(t)$, one can
rewrite (71) as follows

$\displaystyle i\frac{\partial}{\partial t}\begin{pmatrix}C^{d}_{\uparrow+}\\\
C^{d}_{\downarrow+}\\\ C^{d}_{\uparrow-}\\\
C^{d}_{\downarrow-}\end{pmatrix}=\begin{pmatrix}{E_{0}-h(t)}&-\frac{\xi_{k}}{E_{0}}i\alpha
ke^{-i\theta_{k}}&i\frac{\xi_{k}}{2E_{0}^{2}}\frac{\partial\Delta}{\partial
t}&{\frac{\Delta}{E_{0}}i\alpha ke^{-i\theta_{k}}}\\\
\frac{\xi_{k}}{E_{0}}i\alpha
ke^{i\theta_{k}}&E_{0}+h(t)&-\frac{\Delta}{E_{0}}i\alpha
ke^{i\theta_{k}}&i\frac{\xi_{k}}{2E_{0}^{2}}\frac{\partial\Delta}{\partial
t}\\\ -i\frac{\xi_{k}}{2E_{0}^{2}}\frac{\partial\Delta}{\partial
t}&\frac{\Delta}{E_{0}}i\alpha
ke^{-i\theta_{k}}&-E_{0}-h(t)&\frac{\xi_{k}}{E_{0}}i\alpha
ke^{-i\theta_{k}}\\\ {-\frac{\Delta}{E_{0}}i\alpha
ke^{i\theta_{k}}}&-i\frac{\xi_{k}}{2E_{0}^{2}}\frac{\partial\Delta}{\partial
t}&-\frac{\xi_{k}}{E_{0}}i\alpha
ke^{i\theta_{k}}&{-E_{0}+h(t)}\end{pmatrix}\begin{pmatrix}C^{d}_{\uparrow+}\\\
C^{d}_{\downarrow+}\\\ C^{d}_{\uparrow-}\\\ C^{d}_{\downarrow-}\end{pmatrix},$
(72)

where $E_{0}=\sqrt{\xi_{k}^{2}+\Delta^{2}}$. One can remove the phase
$\theta_{k}=\arg\big{(}k_{x}+ik_{y}\big{)}$ from (72) by the unitary operator

$\displaystyle\hat{U}_{\theta}=\begin{pmatrix}e^{i(\frac{\pi}{4}-\frac{\theta_{k}}{2})\hat{\sigma}_{z}}&0\\\
0&e^{i(\frac{\pi}{4}-\frac{\theta_{k}}{2})\hat{\sigma}_{z}}\end{pmatrix},$
(73)

so that in the new basis we have

$\displaystyle i\frac{\partial}{\partial
t}\begin{pmatrix}\tilde{C}^{d}_{\uparrow+}\\\ \tilde{C}^{d}_{\downarrow+}\\\
\tilde{C}^{d}_{\uparrow-}\\\
\tilde{C}^{d}_{\downarrow-}\end{pmatrix}=\begin{pmatrix}E_{0}-h(t)&-\frac{\xi}{E_{0}}\alpha
k&i\frac{\xi}{E_{0}^{2}}\frac{\partial\Delta}{\partial
t}&\frac{\Delta}{E_{0}}\alpha k\\\ -\frac{\xi}{E_{0}}\alpha
k&E_{0}+h(t)&\frac{\Delta}{E_{0}}\alpha
k&i\frac{\xi}{E_{0}^{2}}\frac{\partial\Delta}{\partial t}\\\
-i\frac{\xi}{E_{0}^{2}}\frac{\partial\Delta}{\partial
t}&\frac{\Delta}{E_{0}}\alpha k&-E_{0}-h(t)&\frac{\xi}{E_{0}}\alpha k\\\
\frac{\Delta}{E_{0}}\alpha
k&-i\frac{\xi}{E_{0}^{2}}\frac{\partial\Delta}{\partial
t}&\frac{\xi}{E_{0}}\alpha
k&-E_{0}+h(t)\end{pmatrix}\begin{pmatrix}\tilde{C}^{d}_{\uparrow+}\\\
\tilde{C}^{d}_{\downarrow+}\\\ \tilde{C}^{d}_{\uparrow-}\\\
\tilde{C}^{d}_{\downarrow-}\end{pmatrix}.$ (74)

We assume that the time evolution of the gap function $\Delta(t)$ is adiabatic
on the timescale of the problem (74). Therefore one can assume $\Delta$ to be
constant during the transition with the typical time $\sim\tau_{\text{LZ}}$.
Since the most emphasized dynamics occurs between two crossing branches, it is
convenient to consider the interaction of only the corresponding terms
$C_{k\uparrow+}$ and $C_{k\downarrow-}$ [Fig. 3]. Hence, one can extract an
effective two-level problem for the crossing levels:

$\displaystyle i\frac{\partial}{\partial
t}\begin{pmatrix}\tilde{C}^{d}_{k\uparrow+}\\\
\tilde{C}^{d}_{k\downarrow-}\end{pmatrix}=\begin{pmatrix}E_{0}-\gamma
t&\frac{\Delta}{E_{0}}\alpha k\\\ \frac{\Delta}{E_{0}}\alpha k&-E_{0}+\gamma
t\end{pmatrix}\begin{pmatrix}\tilde{C}^{d}_{k\uparrow+}\\\
\tilde{C}^{d}_{k\downarrow-}\end{pmatrix}.$ (75)

This system can be viewed as the LZSM problem, which allows an exact solution
[35]. However, as discussed in IV.2, one can neglect the transient dynamics of
the $C_{k}(t)$ coefficients in the gap equation (3) and use the transition
matrix approach instead. Thus, we need to obtain the relation between the
long-time asymptotes of the functions
$\tilde{C}^{d}_{k\uparrow\downarrow+}(t)$ before ($t_{0}-$) and after
($t_{0}+$) transition at the point
$t_{0}(\xi_{k})=\sqrt{\xi_{k}^{2}+\Delta^{2}}/\gamma$. Here we use short
notations $(t_{0}\mp)\approx t_{0}\mp\tau_{\text{LZ}}/2$. The asymptotic
solution of the problem (75) is well-known [35] and reads

$\displaystyle\begin{pmatrix}\tilde{C}^{d}_{k\uparrow+}(t_{0}+)\\\
\tilde{C}^{d}_{k\downarrow-}(t_{0}+)\end{pmatrix}=$ (76)
$\displaystyle\begin{pmatrix}\sqrt{p_{k}}&-\text{sgn}(\alpha)\sqrt{1-p_{k}}e^{i\chi_{k}}\\\
\text{sgn}(\alpha)\sqrt{1-p_{k}}e^{-i\chi_{k}}&\sqrt{p_{k}}\end{pmatrix}$
$\displaystyle\times\begin{pmatrix}\tilde{C}^{d}_{k\uparrow+}(t_{0}-)\\\
\tilde{C}^{d}_{k\downarrow-}(t_{0}-)\end{pmatrix},$

where the coefficient

$p_{k}=\exp\Big{[}-\delta_{\text{LZ}}\frac{\Delta^{2}}{\xi_{k}^{2}+\Delta^{2}}\Big{]}$

with
$\delta_{\text{LZ}}=\pi\alpha^{2}k^{2}/\gamma\approx\pi\alpha^{2}k_{F}^{2}/\gamma$
defines the probability of tunneling. Here
$\chi_{k}=\pi/4+\text{arg}\Gamma(1+i\ln p_{k}/2\pi)-\ln p_{k}(\ln(-\ln
p_{k}/2\pi)-1)/2\pi$ is the Stokes phase with the Gamma function $\Gamma$.

For small energies $\xi_{k}\lesssim\Delta$ two different tunneling regimes are
possible:

$\displaystyle(\text{weak})\quad\gamma\gtrsim\alpha^{2}k_{F}^{2}\quad\rightarrow\quad\delta_{\text{LZ}}\approx
0\quad\rightarrow\quad p_{k}\approx 1,$
$\displaystyle(\text{strong})\quad\gamma\ll\alpha^{2}k_{F}^{2}\quad\rightarrow\quad\delta_{\text{LZ}}\gg
1\quad\rightarrow\quad p_{k}\approx 0.$

When $\xi_{k}\gg\Delta$, tunneling is suppressed ($p_{k}\rightarrow 1$) as the
quasiparticle spectrum resembles that of a normal metal with no splitting
between crossing spectral branches.

The typical transient time $\tau_{\text{LZ}}$ for the LZSM tunneling can be
estimated as follows [35]

$\tau_{\text{LZ}}\sim\sqrt{\frac{\hbar}{\gamma}}\text{max}\Bigg{\\{}1,\frac{\alpha
k_{F}}{\sqrt{2\gamma}}\frac{\Delta}{\sqrt{\xi_{k}^{2}+\Delta^{2}}}\Bigg{\\}}.$

If the intersection of the branches of the QP spectrum occurs at some
$\xi_{k}$, then it is possible to determine the interval $\Delta\xi_{k}$ in
which all QP states experience transient dynamics. The size of $\Delta\xi_{k}$
depends on transient time, however, it can be shown, that the upper limit for
this interval is $\Delta\xi_{k}\sim\alpha k_{F}\ll\Delta$. The smalness of
$\Delta\xi_{k}$ and the fact that the gap function $\Delta(t)$ is determined
by all QP states in $(-\omega_{D},\omega_{D})$ confirm the validity of the
approximations made in Section IV.2.

Combining all the results we write the asymptotic transition matrix
$\hat{S}^{d}_{\text{LZ}}$ in diabatic basis as

$\displaystyle\begin{pmatrix}C^{d}_{k\uparrow+}(t_{0}+)\\\
C^{d}_{k\downarrow+}(t_{0}+)\\\ C^{d}_{k\uparrow-}(t_{0}+)\\\
C^{d}_{k\downarrow-}(t_{0}+)\end{pmatrix}=\begin{pmatrix}\sqrt{p_{k}}&0&0&\sqrt{1-p_{k}}e^{i\chi_{k}-i\theta_{k}-i\frac{\pi}{2}\text{sgn}(\alpha)}\\\
0&1&0&0\\\ 0&0&1&0\\\
-\sqrt{1-p_{k}}e^{-i\chi_{k}+i\theta_{k}+i\frac{\pi}{2}\text{sgn}(\alpha)}&0&0&\sqrt{p_{k}}\end{pmatrix}\begin{pmatrix}C^{d}_{k\uparrow+}(t_{0}-)\\\
C^{d}_{k\downarrow+}(t_{0}-)\\\ C^{d}_{k\uparrow-}(t_{0}-)\\\
C^{d}_{k\downarrow-}(t_{0}-)\end{pmatrix}.$ (77)

The LZSM transition matrix in the adiabatic basis (26) has the form
$\hat{S}_{\text{LZ}}=\hat{R}^{-1}(t_{0}+)\hat{S}^{d}_{\text{LZ}}\hat{R}(t_{0}-)$,
where we use the relationship between the two basis (70) and (26) written in
general form as a time-dependent matrix $\hat{R}(t)$. Using the perturbation
theory with respect to the small parameter $\alpha k_{F}/\Delta$ and
considering points $t_{0}\pm$ far from the nonadiabatic region, one can show
that the matrix $\hat{R}(t_{0}\pm)$ can be approximated with an identity
matrix. The corrections proportional to $\alpha k_{F}/\Delta$ in all elements
of the matrix $\hat{R}(t_{0}\pm)$ as well as $\hat{S}_{\text{LZ}}$ can be
neglected, since in all equations of Section IV we consider the minimum
possible order of the perturbation theory with respect to the parameter
$\alpha k_{F}/\Delta$. With these approximations the matrices
$\hat{S}_{\text{LZ}}$ and $\hat{S}^{d}_{\text{LZ}}$ actually coincide and the
LZSM transition matrix in the adiabatic basis can be taken taken from (77).
Thus, we get the Eq (30).

##  Appendix E Calculation of spin-split DOS

The DOS for one spin projection can be written as follows

$\displaystyle
N_{\uparrow}(E,t)\approx\sum_{k}\sum_{n=\uparrow+,\uparrow-}|u_{0}|^{2}\delta\Big{(}E-E_{kn}[h(t)]\Big{)}+|v_{0}|^{2}\delta\Big{(}E+E_{kn}[h(t)]\Big{)}.$
(78)

Here we use static QP amplitudes $u_{0}$ and $v_{0}$ to distinguish the
particle/hole contributions and $E_{k\uparrow\pm}$ are defined in Eq. (4).
Note that the function $N_{\uparrow}(E,t)$ depends on time only through the
Zeeman field $h(t)$. The straightforward calculations for $h(t)<\Delta_{0}$
yield

$\displaystyle\frac{N_{\uparrow}(E,t)}{N(0)}\approx\frac{|E|}{\xi_{0}}\Bigg{|}1-\text{sgn}(E)\frac{\alpha^{2}k_{F}^{2}+h^{2}(t)}{\sqrt{\xi_{0}^{2}\alpha^{2}k^{2}_{F}+h^{2}(t)(\xi_{0}^{2}+\Delta_{h}^{2})}}\Bigg{|}^{-1},$
(79)

where

$\displaystyle\xi_{0}(E,t)\approx\sqrt{E^{2}+h^{2}(t)-\Delta_{h}^{2}+\alpha^{2}k_{F}^{2}+\text{sgn}(E)2\sqrt{E^{2}(h^{2}(t)+\alpha^{2}k_{F}^{2})-\Delta_{h}^{2}\alpha^{2}k_{F}^{2}}}$
(80)

and we have assumed $\alpha k\approx\alpha k_{F}$ due to the vicinity to the
Fermi energy. The time-dependent gap function $\Delta_{h}[h(t)]$ is defined in
(36). Two standard coherence peaks at the energies
$E=-\sqrt{(\Delta_{h}+h(t))^{2}+\alpha^{2}k_{F}^{2}}$ and
$E=\sqrt{(\Delta_{h}-h(t))^{2}+\alpha^{2}k_{F}^{2}}$ appear [Fig. 6(a)].

For the case of large Zeeman fields $h(t)>\Delta_{0}$ one obtains

$\displaystyle\frac{N_{\uparrow}(E,t)}{N(0)}\approx\begin{cases}\frac{|E|}{\xi_{0}}\Bigg{|}1-\text{sgn}(E)\frac{\alpha^{2}k_{F}^{2}+h^{2}(t)}{\sqrt{\xi_{1}^{2}\alpha^{2}k^{2}_{F}+h^{2}(t)(\xi_{1}^{2}+\Delta_{h}^{2})}}\Bigg{|}^{-1},\quad
E>\Delta_{m}\quad\text{and}\quad
E<-\sqrt{(\Delta_{h}+h(t))^{2}+\alpha^{2}k_{F}^{2}}\\\
\frac{|E|}{\xi_{0}}\Bigg{|}1-\frac{\alpha^{2}k_{F}^{2}+h^{2}(t)}{\sqrt{\xi_{2}^{2}\alpha^{2}k^{2}_{F}+h^{2}(t)(\xi_{2}^{2}+\Delta_{h}^{2})}}\Bigg{|}^{-1},\quad-\sqrt{(\Delta_{h}-h(t))^{2}+\alpha^{2}k_{F}^{2}}<E<-\Delta_{m}\end{cases}.$
(81)

where

$\Delta_{m}\approx\frac{\Delta_{h}\alpha
k_{F}}{\sqrt{h^{2}(t)+\alpha^{2}k_{F}^{2}}}.$

The splitting of the energy spectrum in the vicinity of $E=0$ leads to the
appearance of the two additional coherence peaks and corresponging minigap at
the energies $E=\pm\Delta_{m}$, which are shown in Fig. 6(b).

Figure 6:  Spectrum (4) and density of states (79\- 81) for the QPs with
$\uparrow$ spin for two different values of $h(t)$. The value $\Delta_{m}$
represents a minigap. Colored areas in DOS indicate the filling of the states
in the corresponding energy intervals according to Eq. (43). The parameters
are chosen as in Fig. (3): $\Delta_{0}/E_{F}=0.01$, $\alpha/E_{F}=0.0025$.

## References

  * Langenberg  _et al._ [1986] D. Langenberg, A. Larkin, and A. Larkin,  _Nonequilibrium Superconductivity_, Modern problems in condensed matter sciences (North-Holland, 1986).
  * Kopnin [2001] N. Kopnin,  _Theory of Nonequilibrium Superconductivity_, International Series of Monographs on Physics (Clarendon Press, 2001).
  * Volkov and Kogan [1974] A. F. Volkov and S. M. Kogan, Collisionless relaxation of the energy gap in superconductors, Soviet Journal of Experimental and Theoretical Physics 38, 1018 (1974).
  * Kulik  _et al._ [1981] I. O. Kulik, O. Entin-Wohlman, and R. Orbach, Pair susceptibility and mode propagation in superconductors: A microscopic approach, Journal of Low Temperature Physics 43, 591 (1981).
  * Higgs [1964] P. W. Higgs, Broken symmetries and the masses of gauge bosons, Phys. Rev. Lett. 13, 508 (1964).
  * Pekker and Varma [2015] D. Pekker and C. Varma, Amplitude/Higgs modes in condensed matter physics, Annual Review of Condensed Matter Physics 6, 269 (2015).
  * Moor  _et al._ [2017] A. Moor, A. F. Volkov, and K. B. Efetov, Amplitude Higgs mode and admittance in superconductors with a moving condensate, Phys. Rev. Lett. 118, 047001 (2017).
  * Nakamura  _et al._ [2019] S. Nakamura, Y. Iida, Y. Murotani, R. Matsunaga, H. Terai, and R. Shimano, Infrared activation of the Higgs mode by supercurrent injection in superconducting NbN, Phys. Rev. Lett. 122, 257001 (2019).
  * Bellitti  _et al._ [2022] M. Bellitti, C. R. Laumann, and B. Z. Spivak, Incoherent excitation of coherent Higgs oscillations in superconductors, Phys. Rev. B 105, 104513 (2022).
  * Papenkort  _et al._ [2007] T. Papenkort, V. M. Axt, and T. Kuhn, Coherent dynamics and pump-probe spectra of BCS superconductors, Phys. Rev. B 76, 224522 (2007).
  * Matsunaga  _et al._ [2014] 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 345, 1145 (2014).
  * Matsunaga  _et al._ [2013] R. Matsunaga, Y. I. Hamada, K. Makise, Y. Uzawa, H. Terai, Z. Wang, and R. Shimano, Higgs amplitude mode in the bcs superconductors $\text{Nb}_{1-x}\text{Ti}_{x}\text{N}$ induced by terahertz pulse excitation, Phys. Rev. Lett. 111, 057002 (2013).
  * Kemper  _et al._ [2015] 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 92, 224517 (2015).
  * Shimano and Tsuji [2020] R. Shimano and N. Tsuji, Higgs mode in superconductors, Annual Review of Condensed Matter Physics 11, 103 (2020).
  * Buzdin [2005] A. I. Buzdin, Proximity effects in superconductor-ferromagnet heterostructures, Reviews of Modern Physics 77, 935 (2005).
  * Eschrig [2015] M. Eschrig, Spin-polarized supercurrents for spintronics: a review of current progress, Reports on Progress in Physics 78, 104501 (2015).
  * Heikkilä  _et al._ [2019] T. T. Heikkilä, M. Silaev, P. Virtanen, and F. S. Bergeret, Thermal, electric and spin transport in superconductor/ferromagnetic-insulator structures, Progress in Surface Science 94, 100540 (2019).
  * Houzet [2008] M. Houzet, Ferromagnetic josephson junction with precessing magnetization, Physical Review Letters 101, 057009 (2008).
  * Barnes  _et al._ [2011] S. E. Barnes, M. Aprili, I. Petković, and S. Maekawa, Ferromagnetic resonance with a magnetic josephson junction, Superconductor Science and Technology 24, 024020 (2011).
  * Petković  _et al._ [2009] I. Petković, M. Aprili, S. E. Barnes, F. Beuneu, and S. Maekawa, Direct dynamical coupling of spin modes and singlet josephson supercurrent in ferromagnetic josephson junctions, Phys. Rev. B 80, 220502 (2009).
  * Takahashi  _et al._ [2007] S. Takahashi, S. Hikino, M. Mori, J. Martinek, and S. Maekawa, Supercurrent pumping in josephson junctions with a half-metallic ferromagnet, Phys. Rev. Lett. 99, 057003 (2007).
  * Li  _et al._ [2018] L.-L. Li, Y.-L. Zhao, X.-X. Zhang, and Y. Sun, Possible evidence for spin-transfer torque induced by spin-triplet supercurrents*, Chinese Physics Letters 35, 077401 (2018).
  * Golovchanskiy  _et al._ [2023] I. Golovchanskiy, N. Abramov, O. Emelyanova, I. Shchetinin, V. Ryazanov, A. Golubov, and V. Stolyarov, Magnetization dynamics in proximity-coupled superconductor-ferromagnet-superconductor multilayers. ii. thickness dependence of the superconducting torque, Phys. Rev. Appl. 19, 034025 (2023).
  * Silaev [2022] M. Silaev, Anderson-Higgs mass of magnons in superconductor-ferromagnet-superconductor systems, Phys. Rev. Appl. 18, L061004 (2022).
  * Vadimov  _et al._ [2019] V. L. Vadimov, I. M. Khaymovich, and A. S. Mel’nikov, Higgs modes in proximized superconducting systems, Physical Review B 100, 104515 (2019).
  * Tang  _et al._ [2020] G. Tang, W. Belzig, U. Zülicke, and C. Bruder, Signatures of the Higgs mode in transport through a normal-metal–superconductor junction, Phys. Rev. Res. 2, 022068 (2020).
  * Lu  _et al._ [2022a] Y. Lu, S. Ilić, R. Ojajärvi, T. T. Heikkilä, and F. S. Bergeret, Reducing the frequency of the Higgs mode in a helical superconductor coupled to an LC-circuit, arXiv:2212.11615 (2022a).
  * Silaev  _et al._ [2020] M. A. Silaev, R. Ojajärvi, and T. T. Heikkilä, Spin and charge currents driven by the higgs mode in high-field superconductors, Physical Review Research 2, 033416 (2020).
  * Lu  _et al._ [2022b] Y. Lu, R. Ojajärvi, P. Virtanen, M. A. Silaev, and T. T. Heikkilä, Coupling the Higgs mode and ferromagnetic resonance in spin-split superconductors with Rashba spin-orbit coupling, Physical Review B 106, 024514 (2022b).
  * Sarma [1963] G. Sarma, On the influence of a uniform exchange field acting on the spins of the conduction electrons in a superconductor, Journal of Physics and Chemistry of Solids 24, 1029 (1963).
  * Abrikosov [2017] A. Abrikosov,  _Fundamentals of the Theory of Metals_ (Dover Publications, 2017).
  * Tewari  _et al._ [2011] S. Tewari, T. D. Stanescu, J. D. Sau, and S. D. Sarma, Topologically non-trivial superconductivity in spin–orbit-coupled systems: bulk phases and quantum phase transitions, New Journal of Physics 13, 065004 (2011).
  * Ketterson  _et al._ [1999] J. Ketterson, J. Ketterson, S. Song, and K. B,  _Superconductivity_ (Cambridge University Press, 1999).
  * Genwang Fan [2022] P. Z. Genwang Fan, Xiao-Long Chen, Probing two higgs oscillations in a one-dimensional fermi superfluid with raman-type spin-orbit coupling, Frontiers of Physics 17, 52502 (2022).
  * Ivakhnenko  _et al._ [2023] O. V. Ivakhnenko, S. N. Shevchenko, and F. Nori, Nonadiabatic Landau-Zener-Stückelberg-Majorana transitions, dynamics, and interference, Physics Reports 995, 1 (2023).
  * Behrle  _et al._ [2018] A. Behrle, T. Harrison, J. Kombe, K. Gao, M. Link, J.-S. Bernier, C. Kollath, and M. Köhl, Higgs mode in a strongly interacting fermionic superfluid, Nature Physics 14, 781 (2018).
  * Wang  _et al._ [2015] P. Wang, W. Yi, and G. Xianlong, Topological phase transition in the quench dynamics of a one-dimensional fermi gas with spin-orbit coupling, New Journal of Physics 17, 013029 (2015).
  * Dong  _et al._ [2015] Y. Dong, L. Dong, M. Gong, and H. Pu, Dynamical phases in quenched spin–orbit-coupled degenerate fermi gas, Nature communications 6, 6103 (2015).
  * Gor’kov and Rashba [2001] L. P. Gor’kov and E. I. Rashba, Superconducting 2d system with lifted spin degeneracy: Mixed singlet-triplet state, Physical Review Letters 87, 037004 (2001).
  * Frigeri  _et al._ [2004] P. A. Frigeri, D. F. Agterberg, and M. Sigrist, Spin susceptibility in superconductors without inversion symmetry, New Journal of Physics 6, 115 (2004).
  * Szombati  _et al._ [2016] D. B. Szombati, S. Nadj-Perge, D. Car, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Josephson $\phi$-junction in nanowire quantum dots, Nature Physics 12, 568 (2016).
  * Kirilyuk  _et al._ [2010] A. Kirilyuk, A. V. Kimel, and T. Rasing, Ultrafast optical manipulation of magnetic order, Reviews of Modern Physics 82, 2731 (2010).
  * Kirilyuk  _et al._ [2013] A. Kirilyuk, A. V. Kimel, and T. Rasing, Laser-induced magnetization dynamics and reversal in ferrimagnetic alloys, Reports on Progress in Physics 76, 026501 (2013).
  * El-Ghazaly  _et al._ [2020] A. El-Ghazaly, J. Gorchon, R. B. Wilson, A. Pattabi, and J. Bokor, Progress towards ultrafast spintronics applications, Journal of Magnetism and Magnetic Materials 502, 166478 (2020).
  * Plekhanov  _et al._ [2019] K. Plekhanov, M. Thakurathi, D. Loss, and J. Klinovaja, Floquet second-order topological superconductor driven via ferromagnetic resonance, Physical Review Research 1, 032013 (2019).
  * Gershenzon  _et al._ [1990] E. Gershenzon, M. Gershenzon, G. Gol’tsman, A. Lyul’kin, A. Semenov, and A. Sergeev, Electron-phonon interaction in ultrathin Nb films, Sov. Phys. JETP 70, 505 (1990).
  * Kardakova  _et al._ [2013] A. Kardakova, M. Finkel, D. Morozov, V. Kovalyuk, P. An, C. Dunscombe, M. Tarkhov, P. Mauskopf, T. M. Klapwijk, and G. Goltsman, The electron-phonon relaxation time in thin superconducting titanium nitride films, Applied Physics Letters 103, 252602 (2013).
  * Tsuchiya  _et al._ [2018] S. Tsuchiya, D. Yamamoto, R. Yoshii, and M. Nitta, Hidden charge-conjugation, parity, and time-reversal symmetries and massive Goldstone (Higgs) modes in superconductors, Physical Review B 98, 094503 (2018).