# Optical Pumping of Bardeen-Cooper-Schrieffer Superconductors

Vanessa Sulaiman<EMAIL_ADDRESS>Condensed Matter Theory, TU
Dortmund University, Otto-Hahn Straße 4, 44227 Dortmund, Germany  Götz S.
Uhrig<EMAIL_ADDRESS>Condensed Matter Theory, TU Dortmund
University, Otto-Hahn Straße 4, 44227 Dortmund, Germany

###### Abstract

Motivated by the generation by optical pulses of non-thermal distributions of
nuclear spins in quantum dots we investigate the effect of optical pulses
applied to Bardeen-Cooper-Schrieffer (BCS) superconductors. Using time-
dependent mean-field theory formulated with Anderson pseudospins, we study the
electronic configurations and the energy deposited in the system by optical
pulses. The pulses are included by Peierls substitution and we study short
rectangular pulses as well as idealized $\delta$ pulses. Already a few and
even a single pulse generates highly non-trivial distributions of electron
expectation values which we simulate numerically and explain analytically
based on the linearization of the equations of motion. These results suggest
so far unexplored experimental possibilities for the optical control of
superconducting states.

##  I Introduction

Since the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity was
first established [1], many different experiments with continuous irradiation
[2, 3, 4] were carried out to investigate whether an increase of the critical
current and of the energy gap can be reached in this way. In thin-film
bridges, an increase of the critical current could be reached [2, 3]. However,
this approach did not enhance the gap in a measurable way in similar
experiments with tunnel junctions [3].

Recent advances in the generation of ultrashort pulses have enabled further
studies using terahertz (THz) pump-probe setups [5, 6]. These pulses allow the
investigation of non-adiabatic excitations, for instance, for exciting Higgs
modes, because they are fast enough to shake the Cooper pair condensate and
their THz carrier frequency is in the same range as the typical gap energy
${\Delta/h}$ divided by the Planck constant $h$. This has renewed the
theoretical interest in the BCS model, in particular in the non-equilibrium
regime. Many studies investigated quenches, i.e., instantaneous changes of the
order parameter or the coupling constant. Some focus on the resulting dynamics
of the order parameter [7, 8, 9, 10, 11], others on the excitation of
collective modes such as the Higgs mode or the Anderson-Bogoliubov (phase)
mode [12, 13, 14]. In addition, many studies on the response induced by
external electromagnetic fields were carried out [15, 16, 17].

The mean-field dynamics of the BCS model share great similarities with those
of the central spin model [18]. The latter describes the electron spin
dynamics in singly charged quantum dots for which it has been established
experimentally [19, 20, 21] and theoretically [22, 23, 24, 25, 26] that
periodic pulses lead to quasi-stationary states favoring periodic revival of
the spin polarization. This enhances the coherence of the central spin
significantly which can be useful for applications such as quantum state
purification [27]. The appearance of commensurate structures in the
distribution of nuclear spins of quantum dots by optical pulses motivates us
to look for similar non-trivial effects in related systems. Here we
investigate which structures are induced by optical pulses applied to BCS
superconductors. Hence, this study is complementary to the numerous studies
investigating THz pulses.

We investigate the effect of optical pulses on the electronic configurations
in the BCS model. Two different pulse shapes are considered: a rectangular
pulse and an idealized $\delta$ pulse. The aim is to simulate and to analyze
the pulse-induced distribution of expectation values theoretically in order to
motivate experimental studies of this type and to predict the results to be
expected. This will extend the possibilities of actively controlling matter.

This article is organized as follows: In Sect. II, we introduce the model
including the description of the electromagnetic pulse. Subsequently, we
investigate rectangular pulses in Sect. III and in particular the influence of
the pulse duration. In Sect. IV, we examine an instantaneous $\delta$ pulses
and the influence of the delay time. We conclude our analysis in Sect. V.

##  II Model

The Hamiltonian for $s$-wave BCS superconductors is given by

$H_{\text{BCS}}=\sum_{k,\sigma}\varepsilon_{k}c_{k,\sigma}^{\dagger}c_{k,\sigma}-\frac{V}{N}\sum_{k,k^{\prime}}c_{k,\uparrow}^{\dagger}c_{-k,\downarrow}^{\dagger}c_{-k^{\prime},\downarrow}c_{k^{\prime},\uparrow},$
(1)

where $c_{k,\sigma}^{(\dagger)}$ are the fermionic annihilaton (creation)
operators with spin ${\sigma\in\\{\uparrow,\downarrow\\}}$ and
$\varepsilon_{k}$ the electronic dispersion. The second term describes the
attractive electron-electron interaction with $V>0$ between the electrons of
the Cooper pairs. It results from the elimination of the electron-phonon
interaction, for instance by unitary transformations, see, e.g., Ref. 28.
Since the relevant energies are much smaller than the Fermi energy
$\varepsilon_{\text{F}}$, we consider the energy difference

$\varepsilon_{k}=\frac{\hbar^{2}k^{2}}{2m}-\varepsilon_{\text{F}}.$ (2)

The order parameter

$\Delta=\frac{V}{N}\sum_{k^{\prime}}\langle
c_{-k^{\prime},\downarrow}c_{k^{\prime},\uparrow}\rangle$ (3)

is also the energy gap of the system. The binding energy of a single Cooper
pair is given by $2\Delta$ which is also the energy of a Higgs boson [29, 30].

The Anderson pseudospin operators [31, 18]

$\displaystyle\sigma_{k}^{z}$
$\displaystyle=\left(c_{k,\uparrow}^{\dagger}c_{k,\uparrow}+c_{-k,\downarrow}^{\dagger}c_{-k,\downarrow}-1\right)/2$
(4a) $\displaystyle\sigma_{k}^{+}$
$\displaystyle=c_{k,\uparrow}^{\dagger}c_{-k\downarrow}^{\dagger}$ (4b)
$\displaystyle\sigma_{k}^{-}$ $\displaystyle=c_{-k,\downarrow}c_{k,\uparrow}$
(4c)

have all properties of spin-$\frac{1}{2}$ operators. Thus, we can utilize the
same methods as for spin systems to analyze the Hamiltonian and its dynamics.
We rewrite the Hamiltonian with the pseudospin operators in the form

$H_{\text{BCS}}=\sum_{k}2\varepsilon_{k}\sigma_{k}^{z}-\frac{V}{N}L_{+}L_{-}\quad\text{with}\quad\bm{L}=\sum_{k}\bm{\sigma}_{k}$
(5)

where we used bold face symbols to denote vectors which we continue to do
below. Then, the self-consistency implies

$\Delta(t)=\Delta_{x}(t)-i\Delta_{y}(t)=\frac{V}{N}\sum_{k^{\prime}}\langle\sigma_{k}^{-}(t)\rangle.$
(6)

Now we apply mean-field BCS theory where the effective field seen by each
pseudospin $\bm{\sigma}_{k}$ is replaced by the quantum mechanical expectation
values

$\bm{b}_{k}(t)=(-2\Delta_{x}(t),-2\Delta_{y}(t),2\varepsilon_{k})^{\top}.$ (7)

The resulting Hamiltonian

$H_{\text{BCS}}=\sum_{k}\bm{b}_{k}\cdot\bm{\sigma}_{k}$ (8)

induces the equations of motion

$\dot{\bm{\sigma}}_{k}=i[H_{\text{BCS}},\bm{\sigma}_{k}]=\bm{b}_{k}\times\bm{\sigma}_{k}.$
(9)

Now, we can easily take the quantum mechanical expectation value

$\bm{s}_{k}(t)=2\langle\bm{\sigma}_{k}(t)\rangle$ (10)

because the equations of motion are linear. The resulting classical equations
read

$\dot{\bm{s}}_{k}=\bm{b}_{k}\times\bm{s}_{k}$ (11)

with

$\displaystyle\bm{b}_{k}$
$\displaystyle=(-VJ_{x},-VJ_{y},2\varepsilon_{k})^{\top}$ (12) and
$\displaystyle\bm{J}$
$\displaystyle=\frac{1}{N}\sum_{k^{\prime}}\bm{s}_{k^{\prime}}.$ (13)

As the complex phase of the order parameter is constant due to particle-hole
symmetry [32], we can choose $\Delta$ such that $\Delta$ is real, i.e.,
${\Delta_{y}=0}$, for simplicity. Then the self-consistency equation
simplifies to

$\Delta(t)=\frac{V}{N}\sum_{k^{\prime}}\frac{s_{k}^{-}(t)}{2}=\frac{1}{2}\frac{V}{N}\sum_{k^{\prime}}s_{k^{\prime}}^{x}(t).$
(14)

This also implies ${J_{y}=0}$ and the equations of motion reduce to

$\displaystyle\dot{s}_{k}^{x}$ $\displaystyle=-2\varepsilon_{k}s_{k}^{y}$
(15a) $\displaystyle\dot{s}_{k}^{y}$
$\displaystyle=2\varepsilon_{k}s_{k}^{x}+VJ_{x}s_{k}^{z}$ (15b)
$\displaystyle\dot{s}_{k}^{z}$ $\displaystyle=-VJ_{x}s_{k}^{y}.$ (15c)

The initial conditions are given by the ground state of the BCS model

$\langle
c_{k,\uparrow}^{\dagger}c_{k,\uparrow}\rangle=\frac{1}{2}-\frac{\varepsilon_{k}}{2E_{k}},\quad\langle
c_{-k,\downarrow}c_{k,\uparrow}\rangle=\frac{\Delta}{2E_{k}}$ (16)

with the quasiparticle energy
$E_{k}=\sqrt{\varepsilon_{k}^{2}+\lvert\Delta\rvert^{2}}$. The equilibrium
state is characterized by

$s_{k}^{x}=\frac{\Delta}{E_{k}},\quad s_{k}^{y}=0,\quad
s_{k}^{z}=-\frac{\varepsilon_{k}}{E_{k}}$ (17)

for the Anderson pseudospins. Note that these values represent a static
solution of the self-consistency equations (14).

We use the Peierls substitution [33] to include the electromagnetic field of
the pulse in the model. First, we consider the tight-binding model

$H_{0}=-\sum_{i,j=1}^{N}t_{ij}c_{i}^{\dagger}c_{j}-\mu\sum_{i=1}^{N}c_{i}^{\dagger}c_{i}.$
(18)

The electromagnetic field couples the hopping to $t_{ij}$. Using the
Hamiltonian gauge [34] we set ${\phi(\bm{r},t)=0}$ and the hopping term is
changed by phase resulting from the vector potential. This vector potential
for a sinusoidal field is given by

$\bm{A}(t)=\frac{cE}{\gamma}\cos(\gamma t),$ (19)

where $\gamma$ is the angular frequency and $E$ the amplitude of this
electromagnetic field while the speed of light is given by $c$. For
calculational simplicity, we consider a diagonal orientation of the vector
potential, i.e., $\bm{A}$ is oriented in the (1,1,1) spatial direction so that
all nearest-neighbor hopping elements are affected in the same way. Inserting
the vector potential into the Peierls substitution

$t_{ij}\to t_{ij}\exp\left(-\frac{ie}{\hbar
c}\int_{\bm{R}_{i}}^{\bm{R}_{j}}\bm{A}(\bm{r},t)\,\text{d}\bm{r}\right)$ (20)

and averaging over the oscillations [35] in the spirit of a Magnus expansion
[36], we obtain

$t_{ij}\to t_{ij}J_{0}\left(\frac{E_{0}}{\gamma}\right)\quad\text{with}\quad
E_{0}=\frac{eaE}{\hbar}$ (21)

with the zeroth Bessel function $J_{0}$ and the lattice constant $a$.

Averaging over the oscillations is well justified because we focus on optical
pulses. Even though they are short on an absolute scale of picoseconds, they
are still long relative to a single oscillation period and, expressed
reversely, the oscillations are fast relative to the pulse duration. In
addition, they are fast relative to the internal time scales of the
superconductors. Their energy
(${\hbar\gamma\approx$1\text{\,}\mathrm{e}\mathrm{V}$}$) is by about three
orders of magnitude higher than the value given by the superconducting gap
${\Delta\approx$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$}$. The corrections
to the averaged Hamiltonian can be estimated to be of the order
${\Delta/(\hbar\gamma)}$.

In the BCS model, the effect of the optical pulse on the energies
$\varepsilon_{k}$ consists in the reduction [33, 35]

$\varepsilon_{k}\to\varepsilon_{k}J_{0}\left(\frac{E_{0}}{\gamma}\right)$ (22)

during the pulse. Our estimate is that the argument ${E_{0}/\gamma}$ of the
Bessel function is much less than 1 so that $J_{0}$ is close to 1.

In experiments, a typical material [37] is NbN with a lattice constant of
${a_{\text{NbN}}=$4.39\text{\times}{10}^{-10}\text{\,}\mathrm{m}\text{/}$}$.
Usually, an optical pulse [38] has a photon energy of
${\hbar\gamma=$2.2\text{\,}\mathrm{e}\mathrm{V}$}$ and a length of
${t_{\text{puls}}=$100\text{\,}\mathrm{f}\mathrm{s}$}$ or longer. The
amplitude is approximately
${E\approx$9.5\text{\times}{10}^{8}\text{\,}\mathrm{V}\text{/}\mathrm{m}$}$.
With these parameters, we obtain approximately ${\gamma
t_{\text{puls}}/(2\pi)\approx 50}$ oscillations during one pulse. Thus, we
consider the approximation to use the averaged Hamiltonian as well justified
for pulses of at least ${t_{\text{puls}}=$100\text{\,}\mathrm{f}\mathrm{s}$}$.
If the pulse duration reaches $1\text{\,}\mathrm{p}\mathrm{s}$ it comprises
even 500 oscillations. In addition, we estimate the argument of the Bessel
function

$\frac{eaE}{\hbar\gamma}\approx$0.19$,$ (23)

yielding

$J_{0}\approx$0.99$.$ (24)

##  III Rectangular pulse

We consider a rectangular pulse as shown in Fig. 1. Using the Peierls
substitution, the energies $\varepsilon_{k}$ are modified during the pulse by
the constant factor $J_{0}$. As we consider a diagonal orientation of the
vector field, all components of the pseudospin are affected by this
substitution which simplifies the calculation. We do not expect that this
particular orientation has a big effect on the results.

Figure 1: Time dependence of the effective $\varepsilon_{k}$ due to the
Peierls substitution and the averaging for a rectangular pulse.

Unless otherwise stated, the initial value of the order parameter $\Delta_{0}$
is chosen to be $1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. The BCS
interaction is constant up to a cutoff that is given by the Debye frequency
[28] which we choose to be
${\varepsilon_{\text{c}}=$50\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$}$ for the
dispersion. The simulation was realized for ${N=10^{5}}$ equidistant discrete
energies between ${-\varepsilon_{\text{c}}}$ to $\varepsilon_{\text{c}}$
measured relative to the Fermi energy. The results will illustrate that it is
indeed necessary to choose such a fine discretization.

Figure 2: Dynamics of the superconducting order parameter. The pulse itself is
shown in orange in arbitrary units. The blue oscillations are too fast to be
resolved. The initial gap is
${\Delta_{0}=$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$}$, the pulse depth is
${J_{0}=$0.99$}$ and the pulse duration
${t_{\text{puls}}=$1\text{\,}\mathrm{p}\mathrm{s}$}$.

To illustrate the effect of a train of rectangular pulses, we choose the value
of the Bessel function to be ${J_{0}=0.99}$, the pulse duration
${t_{\text{puls}}=$1\text{\,}\mathrm{p}\mathrm{s}$}$, and the repetition rate
of the pulses ${T_{\text{rep}}=$100\text{\,}\mathrm{p}\mathrm{s}$}$. The
resulting dynamics of the order parameter $\Delta$ are shown in Fig. 2. It
consists of fast sinusoidal oscillations which are not resolved in Fig. 2 with
an envelope decaying like $\frac{1}{\sqrt{t}}$ similar to the dynamics after a
quench [39, 32, 40, 41, 9, 10]. The decay differs slightly from pulse to pulse
and the average order parameter $\Delta_{\text{ave}}$ around which the order
parameter oscillates decreases with each pulse. This is an expected behavior
because quenches yield the long-time limit ${\Delta_{\infty}<\Delta_{0}}$ [39,
10]. So it is natural that each pulse reduces the order parameter by a certain
amount because it perturbs the superconducting order and deposits energy in
the system, see also below.

Figure 3: Expectation values of the pseudospins at ${t=0}$ (thin black solid
lines) and at ${t=T_{\text{rep}}}$, i.e., just before the second pulse (thick
dashed colored lines).

In Fig. 3 the distribution of the expectation values of the pseudospin just
before the second pulse, i.e., at ${t=T_{\text{rep}}}$, is depicted. For
comparison, the initial distribution (17) is shown in black. After the pulse,
some pseudospin values are amplified while others are suppressed depending on
$\varepsilon_{k}$. This can be seen in more detail in Fig. 4. The distances
between the peaks seem to be constant while the height of the peaks decreases
with increasing energy $\varepsilon_{k}$. In addition, there is a finer
structure in the peaks that is shown in the inset. To examine this further, we
will focus on $s_{y}$ because the effects are easier to discern in the $y$
component than in the other components because the $y$ component is initially
zero for all $\varepsilon_{k}$.

Figure 4: Upper envelope of the $y$ component of the spin distribution after
the first pulse and a delay time ${t_{\text{del}}=T_{\text{rep}}}$. It is
fitted by a sinc function as given in Eq. (25).

###  III.1 Pseudospin distribution

To analyze the pseudospin distribution after the first pulse, we examine its
envelope first. The behavior strongly resembles a sinc function

$f(x)=a\left\lvert\frac{\sin(b(x+c))}{x+c}\right\rvert$ (25)

so that we use it as fit function. The resemblance of the envelope to a sinc
function provides evidence that the envelope is essentially given by the
Fourier transform of the rectangular pulse shape. This can indeed be
underpinned by an analytical approximation for the pseudospin distribution the
details of which are provided in App. A. We find that the $y$ component after
a pulse can be approximated by

$\begin{split}s_{k,1}^{y}(t)&=2\frac{\Delta_{0}}{E_{k,0}}(J_{0}-1)\varepsilon_{k,0}t_{\text{puls}}\\\
&\hphantom{=}\cdot\cos\left(2E_{k,0}\left(t-\frac{t_{\text{puls}}}{2}\right)\right)\mathrm{sinc}\left(E_{k,0}t_{\text{puls}}\right).\end{split}$
(26)

The fit as displayed in Fig. 4 is done for negative and positive energies
$\varepsilon_{k}$ separately. It is noticeable that the envelope decays
slightly faster than the sinc function. For low energies, the fit remains
slightly below the envelope, but for higher energies, the amplitude of the fit
is above the envelope. In contrast, the positions of the nodes are quite
accurately reproduced by the fit.

The finding that the sinc function in Eqs. (25) and (26), respectively,
results from the Fourier transform indicates that the distance
$\Delta\varepsilon$ between the nodes depends on the pulse duration
$t_{\text{puls}}$. Hence, we expect

$\Delta\varepsilon=\frac{\pi}{b}=\frac{\pi}{t_{\text{puls}}}.$ (27)

This relation is confirmed by the fit parameters for simulations with various
$J_{0}$ and pulse durations $t_{\text{puls}}$ in Fig. 5 where a very good
agreement between the fitted $b_{\text{sinc}}$ and $t_{\text{puls}}$ is
established.

Figure 5: Fit parameter $b_{\text{sinc}}$ found from the sinc fit plotted as
function of the pulse duration $t_{\text{puls}}$. The linear fit shows the
agreement between both quantities.

###  III.2 Energy deposited by a pulse

In addition to the distribution of the pseudospin values, we are interested in
the energy deposited by the pulse in the superconductor. To this end, we must
know the total energy for given pseudospin values. In mean-field theory, the
mean-field Hamiltonian in Eq. 8 reads

$H_{\text{MF}}=\sum_{k}\bm{b}_{k}\,\cdot\bm{\sigma}_{k}\quad\text{with}\quad\bm{b}_{k}=\begin{pmatrix}-2\Delta_{x}\\\
-2\Delta_{y}\\\ 2\varepsilon_{k}\end{pmatrix}.$

Taking its expectation value (10) leads to

$E=\langle H\rangle=\frac{1}{2}\sum_{k}\tilde{\bm{b}}_{k}\,\cdot\bm{s}_{k}.$
(28)

But this overcounts the interaction energy by a factor of two [42]. To obtain
the correct energy, $\bm{b}_{k}$ has to be modified according to

$\tilde{\bm{b}}_{k}=(-\Delta_{x},\,-\Delta_{y},\,2\varepsilon_{k}).$ (29)

Using the initial value ${\Delta_{y}=0}$ and

$\Delta_{x}=\frac{V}{2N}\sum_{k^{\prime}}s_{k^{\prime}}^{x}=\frac{V}{2}J_{x}$
(30)

we can simplify $E$ to

$E=\sum_{k}\left(\varepsilon_{k}s_{k}^{z}-\frac{V}{4}J_{x}s_{k}^{x}\right).$
(31)

With this relation, we evaluate the energy per pseudospin before and after the
pulse to determine the energy differences and sum over all of them.

Figure 6: Energy difference induced by one pulse as a function of the pulse
depth ${(1-J_{0})}$. Each color represents a different pulse duration
$t_{\text{puls}}$. Note the non-monotonic evolution of the parabolae upon
increasing pulse duration.

In Fig. 6, the total energy difference is plotted as a function of the pulse
depth $J_{0}$ for various pulse durations $t_{\text{puls}}$ encoded by the
given color scale. The lines are quadratic fits

$\Delta E(1-J_{0})=a(1-J_{0})^{2}$ (32)

which describe the dependence nicely. This is in accordance with Fermi’s
golden rule which states that the energy intake is proportional to the modulus
of the squared matrix element. The latter is proportional to the pulse depth
$(1-J_{0})$. Close inspection of Fig. 6 reveals that there is no monotonic
behavior of the prefactor $a$ as function of the pulse duration
$t_{\text{puls}}$. This observation requires further investigation.

Figure 7: Fit parameter $a$ of the quadratic fit (32) of Fig. 6. The non-
monotonic oscillatory behavior is clearly visible.

In Fig. 7 the parameter of the fits of Fig. 6 are plotted against the pulse
duration. There is a clear oscillatory dependence on the pulse duration
$t_{\text{puls}}$ discernible. The same oscillations can be seen in the
explicit energy difference as a function of the pulse duration
$t_{\text{puls}}$ for given pulse depth as is depicted in Fig. 8. Here, we
also find that the frequency of the oscillation does not depend significantly
on the pulse depth $J_{0}$. To quantify this behavior, we fit the data by the
function

$\Delta
E(t_{\text{puls}})=a\sin(bt_{\text{puls}}+\varphi)e^{-dt_{\text{puls}}}+c$
(33)

in order to extract the angular frequencies $b$; they are given in Tab. 1. The
left table illustrates the weak dependence on the pulse depth. The right table
shows that the oscillation frequency scales with the initial order parameter
$\Delta_{0}$. There is roughly a factor of 2 between both quantities ${\hbar
b/\Delta_{0}\approx 2}$. Our data (not shown) indicates that this factor
applies to weak pulses and grows upon increasing pulse depth.

Figure 8: Energy difference induced by one pulse as a function of the pulse
duration $t_{\text{puls}}$. The black dotted lines show the fits to the data
according to Eq. 33.

$1-J_{0}$ | $b\,[$\mathrm{m}\mathrm{e}\mathrm{V}$\mathbin{/}\hbar]$
---|---
$0.01$ |  $2.14$${}\pm{}$ | $0.02$
$0.02$ |  $2.21$${}\pm{}$ | $0.02$
$0.04$ |  $2.37$${}\pm{}$ | $0.01$
$0.06$ |  $2.52$${}\pm{}$ | $0.01$
$0.10$ |  $2.84$${}\pm{}$ | $0.01$

$\Delta_{0}\,[$\mathrm{m}\mathrm{e}\mathrm{V}$]$ | $b\,[$\mathrm{m}\mathrm{e}\mathrm{V}$\mathbin{/}\hbar]$
---|---
| ${1-J_{0}}=0.01$ | ${1-J_{0}}=0.04$
$0.25$ |  $0.5$${}\pm{}$ | $0.1$ |  $0.55$${}\pm{}$ | $0.01$
$0.5$ |  $1.10$${}\pm{}$ | $0.01$ |  $1.25$${}\pm{}$ | $0.01$
$1.0$ |  $2.14$${}\pm{}$ | $0.02$ |  $2.37$${}\pm{}$ | $0.01$
$2.0$ |  $4.17$${}\pm{}$ | $0.02$ |  $4.53$${}\pm{}$ | $0.02$
$4.0$ |  $8.20$${}\pm{}$ | $0.02$ |  $8.74$${}\pm{}$ | $0.02$

Table 1: Angular frequency $b$ of the deposited energy as a function of the
pulse depth ${1-J_{0}}$ (left table) for initial
${\Delta_{0}=$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$}$ and as a function of
the initial order parameter $\Delta_{0}$ (right table) for pulse depths
${{1-J_{0}}=0.01}$ and ${{1-J_{0}}=0.04}$.

For completeness, we provide results for the deposited energy for various
initial values of $\Delta_{0}$ in Fig. 9 vs. the pulse duration.

Figure 9: Energy difference induced by one pulse as a function of the pulse
duration $t_{\text{puls}}$ for various initial order parameters $\Delta_{0}$
at a pulse depth of ${1-J_{0}=$0.01$}$. The black dotted lines indicate the
fits according to Eq. 33.

These observations clearly reveal that the pulses induce to a precession of
the pseudospins which depends on the duration of the pulse and on the initial
order parameters. Naively, one would have assumed that longer pulses always
change the system more than shorter pulses, i.e., put more energy into them.
This is suggested by the naive application of Fermi’s golden rule neglecting
coherent effects. Yet, this turns out to be not appropriate. The dynamics we
are investigating here is coherent on the time scale of the pulses and the
considered delay. Note that we have not included a relaxation mechanism in the
model beyond dephasing. Thus, the pulse can put energy into the
superconductor, but it can also take it out again. We find that the angular
frequency of these oscillations corresponds roughly to ${b\approx
2\Delta_{0}/\hbar}$ which is the frequency of the Higgs mode. The ratio
${\hbar b/\Delta_{0}}$ is closest to 2 for weak pulses with a decreasing trend
for increasing gap values, see Tab. 1.

##  IV $\delta$ pulse

Until now, we varied the parameters of the pulse, but not the repetition time
of the pulses $T_{\text{rep}}$ between the pulses. In order to avoid an
overwhelming amount of structural details we want to simplify the pulses
before addressing the influence of $T_{\text{rep}}$. We achieve this by
eliminating the dependence on the pulse duration by using idealized $\delta$
pulses, i.e., instantaneous pulses. We approximate the rectangular pulse by a
$\delta$ pulse with the same area ${(1-J_{0})t_{\text{puls}}}$. Note that this
appears counterintuitive since an infinite depth ${1-J_{0}\to\infty}$ does not
make sense as renormalization of the dispersion by an optical pulse. But we
will show that this idealization of rectangular pulses works quite well and
helps simplifying further investigations.

At first, we rewrite the time-dependent $\varepsilon_{k}$ as a sum of a
constant and a time-dependent rectangular function

$\varepsilon_{k}(t)=\varepsilon_{k,0}+\Delta\varepsilon_{k}\cdot\text{rect}\left(\frac{t\text{
mod }T_{\text{rep}}}{t_{\text{puls}}}-\frac{1}{2}\right)$ (34)

with ${\Delta\varepsilon_{k}=-(1-J_{0})\varepsilon_{k,0}}$, where the
rectangular function is defined as

$\text{rect}(t)=\begin{cases}0,&\lvert t\rvert>\frac{1}{2}\\\ 1,&\lvert
t\rvert<\frac{1}{2}\end{cases}.$ (35)

This represents the rectangular pulse as shown above in Fig. 1.

If we take the limit ${t_{\text{puls}}\to 0}$ while keeping the area of the
pulse ${(1-J_{0})t_{\text{puls}}}$ constant, the prefactor
${\Delta\varepsilon_{k}}$ of the rectangular function diverges to $\infty$.
Then, this term dominates in the equations of motion. Neglecting all non-
diverging terms in the equations, we arrive at equations that describe a
simple rotation about the $z$ axis with angular frequency
${2\Delta\varepsilon_{k}=2(J_{0}-1)\varepsilon_{k,0}}$ having set ${\hbar=1}$.
Then, the effect of the pulse is fully captured by the rotation matrix

$\bm{s}_{k}(t_{\text{puls}})=\begin{pmatrix}\cos(\varphi_{k})&-\sin(\varphi_{k})&0\\\
\sin(\varphi_{k})&\cos(\varphi_{k})&0\\\ 0&0&1\end{pmatrix}\bm{s}_{k}(0)$
(36a) with the rotation angle $\displaystyle\varphi_{k}$
$\displaystyle=2\Delta\varepsilon_{k}t_{\text{puls}}$ (36b)
$\displaystyle=2(J_{0}-1)\varepsilon_{k,0}t_{\text{puls}}\eqqcolon-2\varphi\varepsilon_{k,0}.$
(36c)

The parameters $J_{0}$ and $t_{\text{puls}}$ are combined to the parameter
${\varphi=(1-J_{0})t_{\text{puls}}}$ which we keep constant in the
idealization by keeping the rectangular area constant. In this way, we reduce
the dynamics during the pulse to an instantaneous rotation around the $z$
axis. Solving the whole dynamics means that we alternate between the original
equations of motion (15) _between_ the pulses and the instantaneous rotation
(36a). This corresponds to the pulse shape of a $\delta$, significantly
simplifies the numerics, and reduces the number of parameters.

Figure 10: Comparison of the dynamics for different rectangular pulses of the
same area ${(1-J_{0})t_{\text{puls}}=$0.01\text{\,}\mathrm{p}\mathrm{s}$}$ and
the corresponding $\delta$ pulse. The convergence of the results for the
rectangular pulses with decreasing duration to the $\delta$ pulse is good.

First, we compare the results of rectangular pulses and the $\delta$ pulse to
investigate if the idea to approximate the rectangular pulse by an
instantaneous one makes sense. In Figs. 10 and 11, we compare the dynamics of
the order parameter for several rectangular pulses of the same area to the
dynamics of a $\delta$ pulse with the same area. The agreement between the
shortest rectangular pulse and the $\delta$ pulse is good enough to conclude
that the rectangular pulse converges reasonably well to the $\delta$ pulse in
the limit ${t_{\text{puls}}\to 0}$. This can also be confirmed by a comparison
of the distributions of the pseudospin expectation values after the first
pulse, i.e., at ${t=T_{\text{rep}}=$100\text{\,}\mathrm{p}\mathrm{s}$}$.

Figure 11: Comparison of the dynamics for trains of different rectangular
pulses of the same area
${(1-J_{0})t_{\text{puls}}=$0.01\text{\,}\mathrm{p}\mathrm{s}$}$ and the
corresponding train of $\delta$ pulses. Five pulses with a repetition time of
${T_{\text{rep}}=$100\text{\,}\mathrm{p}\mathrm{s}$}$ are shown.

Figure 12 shows the envelope of the $y$ component of the pseudospin
expectation values after a single pulse, but with a certain delay of
${t_{\text{del}}=$100\text{\,}\mathrm{p}\mathrm{s}$}$. As was to be expected
the nodes of the envelope vanish for the $\delta$ pulse because we established
an inverse relationship between the distance between the nodes and the pulse
duration. Hence the distance tends to infinity for vanishing pulse duration.
The convergence to the result induced by the $\delta$ pulse is satisfactory,
but far from perfect which has to be kept in mind for the subsequent
investigations.

Figure 12: Comparison of the envelopes of the pseudospin distribution after a
delay time ${t_{\text{del}}=$100\text{\,}\mathrm{p}\mathrm{s}$}$ after a
single pulse.

In order to study the effect of the time between the pulses quantitatively, we
first focus on a single pulse and compute the pseudospin values after some
delay time $t_{\text{del}}$ after this pulse in Fig. 13. It can be seen that
the envelope does hardly change while the oscillations in the fine structure
are strongly modified. To analyze these modifications, we calculate the
distances between neighboring zeros in the fine structure and plot them in
Fig. 14 vs. the energy $\varepsilon$. Close to ${\varepsilon=0}$, the node
distances are quite large, but converge quickly to a constant value for
increasing $\varepsilon$. We display the distances $\delta_{0}$ in the
converged regime at large $\varepsilon$ against the inverse delay time
${1/t_{\text{del}}}$ in Fig. 15. This unambiguously shows that the converged
distances are given by ${c/t_{\text{del}}}$ with the prefactor
${a=\frac{\pi}{2}\hbar}$; recall that we set $\hbar$ to one in the
calculations. Thus, the energy interval of a full oscillation is given by
${\hbar\pi/t_{\text{del}}}$.

Figure 13: Zoom of the pseudospin distribution at different delay times after
a pulse. Figure 14: Analysis of the oscillation period of the pseudospin
values as function of the delay time $t_{\text{del}}$ after a pulse. The
distances between the zeros are plotted vs. the energy $\varepsilon$. The
dashed lines indicate the converged limit for high energies $\varepsilon$.
Figure 15: The limits shown by dashed lines in Fig. 14 are compared to the
inverse delay time $t_{\text{del}}^{-1}$. The linear fit works perfectly.

In App. A, we show that the $y$ component at a time delay $t_{\text{del}}$
after a pulse is approximated by

$s_{k,1}^{y}(t_{\text{del}})=-2\frac{\Delta_{0}}{E_{k,0}}\varphi\varepsilon_{k,0}\cos\left(2E_{k,0}t_{\text{del}}\right).$
(37)

Thus, for a given delay time $t_{\text{del}}$, the pseudospin distribution is
proportional to

$\displaystyle\cos\left(2E_{k,0}t_{\text{del}}\right)$
$\displaystyle=\cos\left(2t_{\text{del}}\sqrt{\varepsilon_{k}^{2}+\Delta_{0}^{2}}\right)$
(38a)
$\displaystyle\approx\cos\left(2t_{\text{del}}\varepsilon_{k}\right)\text{ for
}\varepsilon_{k}\gg\Delta_{0},$ (38b)

i.e., to cosine. As we are investigating the pseudospin distribution vs.
$\varepsilon_{k}$, the prefactor of $\varepsilon_{k}$ in the argument of the
cosine is a time given by ${2t_{\text{del}}}$. Hence, the distance between two
zeros takes the value

$\delta_{0}=\frac{\pi}{2}t_{\text{del}}^{-1}.$ (39)

The same factor occurs in the approximation (26) for the rectangular pulse,
which again corroborates that the $\delta$ pulse induces the same fine
structure as the rectangular pulse except for the sinc envelope.

##  V Conclusion

The present study was motivated by the unconventional structures [23, 25, 26]
and the reduction of entropy [27] that can be achieved by periodic trains of
optical pulses applied to the electron spin in quantum dots with hyperfine
interaction to their baths of nuclear spins. These systems are conveniently
described by central spin models. In the Anderson pseudospin formalism the
time-dependent mean-field theory for BCS superconductors bears strong
similarity with the central spin model used for the spins in quantum dots.
Hence, a study of the effects of optical pulses and trains of them on
superconductors was called for and represents the main goal of the present
study. Indeed, remarkable structures in the distribution of modes appeared.

We numerically investigated the BCS model under the influence of optical
pulses. Their electromagnetic field couples to the superconducting electrons
through Peierls substitution. For optical pulses the effect of the fast
oscillating fields is to reduce the bare dispersion upon averaging over these
oscillations. We found non-equilibrium structures in the occupation of the
electronic modes stemming from such perturbations of the superconductors. The
structures depend on the amplitude and duration of the pulses. We found a
detailed rich distribution of the expectation values of the pseudospins of
which the envelope is governed by the Fourier transform of the pulse shape and
the fast oscillations by the delay time elapsed after the pulse. This is
confirmed by an analytical approximation based on the linearization of the
equations of motion. In addition, we found a periodic dependence of the energy
deposited in the system on the duration of the pulses. The periodicity can be
explained by the coherent precession of the pseudospins such that the overall
effect depends roughly on whether half revolutions fit into the pulse or not,
i.e., the angular frequency of the periodicity is given by $2\Delta_{0}$ where
$\Delta_{0}$ is the initial superconducting gap. The factor of 2 appears in
the limit of very weak pulses only.

The fascinating structures that our theoretical simulations unearthed call for
experimental verification. It would be highly interesting to investigate how
such structures appear in experiments, i.e., in conventional superconductors
which are pumped by ultrashort optical pulses. The results presented here
provide a guideline to such experiments by showing what is to be expected. For
realization, of course, one must be aware of the effects of heating. These do
not appear in our treatment where we average the effect of the Peierls phases
on the hopping. We stress that our setup does not require any particular value
of the optical frequency as long as ${f_{\text{opt}}t_{\text{puls}}\gg 1}$ and
${hf_{\text{opt}}/\Delta_{0}\gg 1}$ holds. The material system best suited for
realization should have a transparency window, i.e., a window of optical
frequencies, within which no or hardly any resonance occurs so that no
absorption takes place. Certainly, this represents quite a challenge to find
and to create.

On the theoretical side, the next steps should include explicit calculations
of further observables that can be directly measured in an experiment, for
instance, the optical conductivity. Another interesting extension consists in
the inclusion of dissipation because such effects are surely present in
experiments. They can be included either phenomenologically as in Ref. 11 in
theoretical treatments or by either to higher order diagrams inducing finite
imaginary parts in the self-energies or by resorting to Lindbladian dynamics
for the pseudospins. These intriguing issues are beyond the scope of the
present investigation and thus are left to future research.

###### Acknowledgements.

This research was supported by the Stiftung Mercator in project no.
KO-2021-0027. We acknowledge very helpful discussions with J. Althüser, A.
Böhmer, J. Bünemann, I. Eremin, and J. Stolze.

##  Appendix A Analytical approximation for the $y$ component of the
pseudospins

We start from the equations of motion

$\displaystyle\dot{s}_{k}^{x}$ $\displaystyle=-2\varepsilon_{k}s_{k}^{y}$
(40a) $\displaystyle\dot{s}_{k}^{y}$
$\displaystyle=2\varepsilon_{k}s_{k}^{x}+VJ_{x}s_{k}^{z}$ (40b)
$\displaystyle\dot{s}_{k}^{z}$ $\displaystyle=-VJ_{x}s_{k}^{y}$ (40c)

where $J_{x}=\frac{1}{N}\sum_{k}s_{k}^{x}$. In general, the effect of the
pulse on the dispersion can be written as

$\varepsilon_{k}(t)=\varepsilon_{k,0}+\lambda f_{k}(t)$ (41)

with a small parameter $\lambda$ and a time-dependent function $f_{k}(t)$, for
instance of rectangular shape in our study. As long as the pulse depth is
small, the dynamics are occuring as a small perturbation of the equilibrium
solution in linear response

$s_{k}^{\alpha}=s_{k,0}^{\alpha}+\lambda
s_{k,1}^{\alpha}\quad\text{with}\quad\alpha\in\\{x,y,z\\}$ (42a) and
$J_{x}=J_{x,0}+\lambda J_{x,1}.$ (42b)

Inserting this into the equations of motion leads to

$\displaystyle\begin{split}\dot{s}_{k,0}^{x}+\lambda\dot{s}_{k,1}^{x}&=-2\varepsilon_{k,0}s_{k,0}^{y}\\\
&\phantom{=}+\lambda(-2\varepsilon_{k,0}s_{k,1}^{y}-2f_{k}(t)s_{k,0}^{y})\\\
&\phantom{=}-\lambda^{2}2f_{k}(t)s_{k,1}^{y}\end{split}$ (43a)
$\displaystyle\begin{split}\dot{s}_{k,0}^{y}+\lambda\dot{s}_{k,1}^{y}&=2\varepsilon_{k,0}s_{k,0}^{x}+gJ_{x,0}s_{k,0}^{z}\\\
&\phantom{=}+\lambda(2\varepsilon_{k,0}s_{k,1}^{x}+2f_{k}(t)s_{k,0}^{x}\\\
&\phantom{=+\lambda(}+VJ_{x,0}s_{k,1}^{z}+VJ_{x,1}s_{k,0}^{z})\\\
&\phantom{=}+\lambda^{2}(2f_{k}(t)s_{k,1}^{x}+VJ_{x,1}s_{k,1}^{z})\end{split}$
(43b)
$\displaystyle\begin{split}\dot{s}_{k,0}^{z}+\lambda\dot{s}_{k,1}^{z}&=-gJ_{x,0}s_{k,0}^{y}\\\
&\phantom{=}-\lambda(VJ_{x,0}s_{k,1}^{y}+VJ_{x,1}s_{k,0}^{y})\\\
&\phantom{=}-\lambda^{2}VJ_{x,1}s_{k,1}^{y}.\end{split}$ (43c)

Next, we neglect all terms of order $\lambda^{2}$ and split the equations into
two parts, namely of zeroth and of first order. The zeroth-order equations are
given by

$\displaystyle\dot{s}_{k,0}^{x}$
$\displaystyle=-2\varepsilon_{k,0}s_{k,0}^{y}$ (44a)
$\displaystyle\dot{s}_{k,0}^{y}$
$\displaystyle=2\varepsilon_{k,0}s_{k,0}^{x}+VJ_{x,0}s_{k,0}^{z}$ (44b)
$\displaystyle\dot{s}_{k,0}^{z}$ $\displaystyle=-VJ_{x,0}s_{k,0}^{y}$ (44c)

and the first-order equations are given by

$\displaystyle\dot{s}_{k,1}^{x}$
$\displaystyle=-2\varepsilon_{k,0}s_{k,1}^{y}-2f_{k}(t)s_{k,0}^{y}$ (45a)
$\displaystyle\begin{split}\dot{s}_{k,1}^{y}&=2\varepsilon_{k,0}s_{k,1}^{x}+VJ_{x,0}s_{k,1}^{z}\\\
&\phantom{=}+2f_{k}(t)s_{k,0}^{x}+VJ_{x,1}s_{k,0}^{z}\end{split}$ (45b)
$\displaystyle\dot{s}_{k,1}^{z}$
$\displaystyle=-VJ_{x,0}s_{k,1}^{y}-VJ_{x,1}s_{k,0}^{y}.$ (45c)

The zeroth order reproduces the equilibrium solution which we insert into the
first-order equations. With the matrix

$\bm{M}\coloneqq\begin{pmatrix}0&-2\varepsilon_{k,0}&0\\\
2\varepsilon_{k,0}&0&VJ_{x,0}\\\ 0&-VJ_{x,0}&0\end{pmatrix}$ (46)

the resulting equations can be written as

$\dot{\bm{s}}_{k,0}(t)=\bm{M}\bm{s}_{k,0}(t)$ (47a) and
$\dot{\bm{s}}_{k,1}(t)=D_{1}+D_{2}+D_{\text{ext}}$ (47b) with $\displaystyle
D_{1}$ $\displaystyle=\bm{M}\bm{s}_{k,1}(t)$ (47c) $\displaystyle D_{2}$
$\displaystyle=VJ_{x,1}(t)\begin{pmatrix}0\\\ s_{k,0}^{z}(t)\\\
-s_{k,0}^{y}(t)\end{pmatrix}$ (47d) $\displaystyle D_{\text{ext}}$
$\displaystyle=2f_{k}(t)\begin{pmatrix}-s_{k,0}^{y}(t)\\\ s_{k,0}^{x}(t)\\\
0\end{pmatrix}.$ (47e)

Of course, the equilibrium solution is known

$s_{k,0}^{x}=\frac{\Delta_{0}}{E_{k,0}},\quad s_{k,0}^{y}=0,\quad
s_{k,0}^{z}=-\frac{\varepsilon_{k,0}}{E_{k,0}}$ (48)

with $E_{k,0}=\sqrt{\varepsilon_{k,0}^{2}+\lvert\Delta_{0}\rvert^{2}}$.

We want to neglect $D_{2}$ in Eq. (47b) to reach a simpler equation that can
be solved analytically. To verify that this is a valid approximation, we
compare $D_{2}$ to $D_{1}$. For this purpose, we need a measure for the size
of $D_{1}$, namely the projections of the spin configuration onto the
eigenvectors of $\bm{M}$, namely
${\lvert\lambda_{2/3}\rvert\sqrt{\lvert\innerproduct{\bm{s}_{k,1}|\bm{e}_{2}}{\bm{s}_{k,1}|\bm{e}_{2}}\rvert^{2}+\lvert\innerproduct{\bm{s}_{k,1}|\bm{e}_{3}}{\bm{s}_{k,1}|\bm{e}_{3}}\rvert^{2}}}$.
Given the equilibrium solution, the absolute value of $D_{2}$ reads
${V\left\lvert J_{x,1}\frac{\varepsilon_{k,0}}{E_{k,0}}\right\rvert}$. To
compare the two terms, we plot them against $\varepsilon_{k,0}$ for exemplary
pulse parameters and various times in Fig. 16 for the rectangular pulse and in
Fig. 17 for the $\delta$ pulse. Values for $\bm{s}_{k,1}$ are taken from the
numerical evaluation of the differential equations.

Figure 16: Measures of ${|D_{1}|}$ (blue lines) and ${|D_{2}|}$ (orange lines)
as defined in Eqs. (47c) and (47d) for a rectangular pulse with
${J_{0}=$0.99$}$, ${t_{\text{puls}}=$1\text{\,}\mathrm{p}\mathrm{s}$}$ at
various delay times $t_{\text{del}}$. Figure 17: Same as in Fig. 16 for a
$\delta$ pulse with
${\varphi=(J_{0}-1)t_{\text{puls}}=$0.01\text{\,}\mathrm{p}\mathrm{s}$}$ at
various delay times ${t_{\text{del}}}$.

For both kinds of pulses, the approximation is valid for large values of
${|\varepsilon_{k,0}|}$ because $D_{2}$ is significantly smaller than $D_{1}$.
For the $\delta$ pulse (Fig. 17) $D_{2}$ is almost everywhere smaller than
$D_{1}$ except in the immediate vicinity of ${\varepsilon_{k,0}=0}$. After the
rectangular pulse (lower row in Fig. 16), $D_{1}$ also dominates clearly over
$D_{2}$. But during the rectangular pulse (upper row in Fig. 16), $D_{2}$ is
of the same order of magnitude as $D_{1}$ up to moderate values of the
argument. We conclude that the approximation is justified except for small
values of ${|\varepsilon_{k,0}|}$. This is in line with our previous
comparison where the analytic approximations worked well for large
${|\varepsilon_{k,0}|}$, but showed deviations for small values.

Neglecting the term $D_{2}$, the differential equation can be solved
analytically by multi-dimensional variation of constants

$\displaystyle\bm{s}_{k,1}(t)$
$\displaystyle=\exp(\bm{M}t)\int_{0}^{t}\exp(-\bm{M}t^{\prime})\bm{v}(t^{\prime})\,\text{d}t^{\prime}$
(49a) with $\displaystyle\bm{v}(t)$
$\displaystyle=2f_{k}(t)\begin{pmatrix}-s_{k,0}^{y}(t)\\\ s_{k,0}^{x}(t)\\\
0\end{pmatrix}.$ (49b)

The eigenvalues and eigenvectors of ${-\bm{M}}$ are given by

$\displaystyle\lambda_{1}$ $\displaystyle=0$ $\displaystyle\bm{e}_{1}$
$\displaystyle=\frac{1}{2E_{k,0}}\begin{pmatrix}-VJ_{x,0}\\\ 0\\\
2\varepsilon_{k,0}\end{pmatrix}$ (50a) $\displaystyle\lambda_{2}$
$\displaystyle=-2iE_{k,0}$ $\displaystyle\bm{e}_{2}$
$\displaystyle=\frac{1}{2\sqrt{2}\,E_{k,0}}\begin{pmatrix}2\varepsilon_{k,0}\\\
-2iE_{k,0}\\\ VJ_{x,0}\end{pmatrix}$ (50b) $\displaystyle\lambda_{3}$
$\displaystyle=2iE_{k,0}$ $\displaystyle\bm{e}_{3}$
$\displaystyle=\frac{1}{2\sqrt{2}\,E_{k,0}}\begin{pmatrix}2\varepsilon_{k,0}\\\
2iE_{k,0}\\\ VJ_{x,0}\end{pmatrix}$ (50c) with $\displaystyle 2E_{k,0}$
$\displaystyle=\sqrt{(2\varepsilon_{k,0})^{2}+\left({V}J_{x,0}/N\right)^{2}}$
(50d) $\displaystyle=\sqrt{(2\varepsilon_{k,0})^{2}+(2\Delta_{0})^{2}}.$ (50e)

Then, the exponential function of $\bm{M}$ and its inverse are given by

$\exp(\mp\bm{M}t^{\prime})=\bm{U}\exp(\pm\bm{D}t^{\prime})\bm{U}^{-1}$ (51a)
with the unitary matrix $\displaystyle\bm{U}$
$\displaystyle=\begin{pmatrix}|&|&|\\\ \bm{e}_{1}&\bm{e}_{2}&\bm{e}_{3}\\\
|&|&|\end{pmatrix}$ (51b) $\displaystyle\exp(\pm\bm{D}t^{\prime})$
$\displaystyle=\begin{pmatrix}e^{\pm\lambda_{1}t^{\prime}}&0&0\\\
0&e^{\pm\lambda_{2}t^{\prime}}&0\\\
0&0&e^{\pm\lambda_{3}t^{\prime}}\end{pmatrix}.$ (51c)

This enables us to write the complete solution as

$\displaystyle\bm{s}_{k,1}(t)$
$\displaystyle=e^{\bm{M}t}\int_{0}^{t}e^{-\bm{M}t^{\prime}}\bm{v}(t^{\prime})\,\text{d}t^{\prime}$
(52a)
$\displaystyle=\bm{U}e^{-\bm{D}t}\int_{0}^{t}e^{\bm{D}t^{\prime}}\bm{U}^{-1}\bm{v}(t^{\prime})\,\text{d}t^{\prime}.$
(52b)

From these formulae, we deduce the $y$ component obtaining

$\displaystyle\begin{split}s_{k,1}^{y}(t)&=e_{1}^{y}\,e^{-\lambda_{1}t}\int_{0}^{t}e^{\lambda_{1}t^{\prime}}\bm{e}_{1}^{*}\cdot\bm{v}(t^{\prime})\,\text{d}t^{\prime}\\\
&\quad+e_{2}^{y}\,e^{-\lambda_{2}t}\int_{0}^{t}e^{\lambda_{2}t^{\prime}}\bm{e}_{2}^{*}\cdot\bm{v}(t^{\prime})\,\text{d}t^{\prime}\\\
&\quad+e_{3}^{y}\,e^{-\lambda_{3}t}\int_{0}^{t}e^{\lambda_{3}t^{\prime}}\bm{e}_{3}^{*}\cdot\bm{v}(t^{\prime})\,\text{d}t^{\prime}.\end{split}$
(53)

The last two terms are the complex conjugate of each other because
${\bm{e}_{2}=\bm{e}_{3}^{*}}$ and ${\lambda_{2}=\lambda_{3}^{*}}$. Thus, one
obtains

$\displaystyle s_{k,1}^{y}(t)$
$\displaystyle=2\operatorname{Re}\Big{(}-\frac{i}{2E_{k,0}}e^{2iE_{k,0}t}\int_{0}^{t}e^{-2iE_{k,0}t^{\prime}}f_{k}(t^{\prime})\cdot$
$\displaystyle\quad\left(-2\varepsilon_{k,0}s_{k,0}^{y}(t^{\prime})+2iE_{k,0}s_{k,0}^{x}(t^{\prime})\right)\,\text{d}t^{\prime}\Big{)}.$
(54)

Inserting the equilibrium solution, we arrive at

$\displaystyle s_{k,1}^{y}(t)$
$\displaystyle=2\operatorname{Re}\left(e^{2iE_{k,0}t}\int_{0}^{t}e^{-2iE_{k,0}t^{\prime}}f_{k}(t^{\prime})\frac{\Delta_{0}}{E_{k,0}}\,\text{d}t^{\prime}\right).$
(55)

For ${t>t_{\text{puls}}}$, the integral can be extended to ${\pm\infty}$
because ${f_{k}(t^{\prime})=0}$ outside of the support of the pulse so that we
obtain

$\displaystyle s_{k,1}^{y}(t)$
$\displaystyle=2\operatorname{Re}\left(e^{2iE_{k,0}t}\int_{-\infty}^{\infty}e^{-2iE_{k,0}t^{\prime}}f_{k}(t^{\prime})\frac{\Delta_{0}}{E_{k,0}}\,\text{d}t^{\prime}\right).$
(56)

Next, we shift the time variable by ${t_{\text{puls}}/2}$ because
${f_{k}(t^{\prime})}$ is symmetric around ${t^{\prime}=t_{\text{puls}}/2}$

$\displaystyle\begin{split}s_{k,1}^{y}(t)=2\frac{\Delta_{0}}{E_{k,0}}\operatorname{Re}\Bigg{(}&e^{2iE_{k,0}\left(t-{t_{\text{puls}}/2}\right)}\\\
&\int_{-\infty}^{\infty}e^{-2iE_{k,0}t^{\prime}}\tilde{f}_{k}(t^{\prime})\,\text{d}t^{\prime}\Bigg{)}.\end{split}$
(57)

The shifted function ${\tilde{f}_{k}(t^{\prime})}$ is a real and even function
so that its Fourier integral is real yielding

$\displaystyle s_{k,1}^{y}(t)$
$\displaystyle=2\frac{\Delta_{0}}{E_{k,0}}\operatorname{Re}\left(e^{2iE_{k,0}\left(t-{t_{\text{puls}}/2}\right)}\right)$
$\displaystyle\quad\cdot\int_{-\infty}^{\infty}e^{-2iE_{k,0}t^{\prime}}\tilde{f}_{k}(t^{\prime})\,\text{d}t^{\prime}$
(58a)
$\displaystyle=2\frac{\Delta_{0}}{E_{k,0}}\cos\left(2E_{k,0}\left(t-{t_{\text{puls}}/2}\right)\right)\mathcal{F}(\tilde{f}_{k}(t^{\prime})).$
(58b)

with ${\mathcal{F}(\tilde{f}_{k}(t^{\prime}))}$ being the Fourier transform of
${\tilde{f}_{k}(t^{\prime})}$. This provides us with an analytical
approximation for the $y$ component of the spin dynamics for small pulse
depths. One can insert various pulse shapes and compare the results to the
numerical solutions.

###  A.1 Rectangular pulse

Figure 18: Comparison of the analytical approximation (orange lines) and the
numerical result (blue lines) for a rectangular pulse of duration
${t_{\text{puls}}=$1\text{\,}\mathrm{p}\mathrm{s}$}$ at two different pulse
depths of ${J_{0}=$0.99$}$ (upper panel) and ${J_{0}=$0.999$}$ (lower panel).
We compare the $y$ component of the pseudospin distribution after a single
pulse and a delay time of
${t_{\text{del}}=$100\text{\,}\mathrm{p}\mathrm{s}$}$.

The shifted rectangular pulse is given by

$\tilde{f}_{k}(t^{\prime})=(J_{0}-1)\varepsilon_{k,0}\begin{cases}1&\lvert
t\rvert\leq{t_{\text{puls}}/2}\\\ 0&\lvert
t\rvert>{t_{\text{puls}}/2}\end{cases}$ (59)

and its Fourier transform reads

$\displaystyle\mathcal{F}(\tilde{f}_{k}(t^{\prime}))$
$\displaystyle=\int_{-\infty}^{\infty}e^{-2iE_{k,0}t^{\prime}}\tilde{f}_{k}(t^{\prime})\,\text{d}t^{\prime}$
(60a)
$\displaystyle=(J_{0}-1)\varepsilon_{k,0}t_{\text{puls}}\,\mathrm{sinc}\left(E_{k,0}t_{\text{puls}}\right)$
(60b)

with the sinc function $\mathrm{sinc}(x)=\sin(x)/x$. Thus, the $y$ component
of the spin dynamics according to Eq. 58b is given by

$\displaystyle s_{k,1}^{y}(t)$
$\displaystyle=2\frac{\Delta_{0}}{E_{k,0}}(J_{0}-1)\varepsilon_{k,0}t_{\text{puls}}\cos\left(2E_{k,0}\left(t-{t_{\text{puls}}/2}\right)\right)$
$\displaystyle\hphantom{=}\cdot\mathrm{sinc}\left(E_{k,0}t_{\text{puls}}\right).$
(61)

This analytical result is compared to the numerical solution in Fig. 18 for
two pulse depths. As expected, the approximation does not work for small
values of $\varepsilon_{k,0}$, but the agreement for higher values is good.
For ${J_{0}=$0.99$}$, the results are slightly shifted, but the overall shape
of approximate and numerical behavior is the same. The envelope of the
numerical results is captured very nicely by the $\mathrm{sinc}$ function. The
fast oscillations are reproduced by the cosine function very well.

###  A.2 $\delta$ pulse

The $\delta$ pulse

$\tilde{f}_{k}(t^{\prime})=-\varphi\varepsilon_{k,0}\delta(t^{\prime})$ (62)

implies the Fourier transform

$\displaystyle\mathcal{F}(\tilde{f}_{k}(t^{\prime}))$
$\displaystyle=\int_{-\infty}^{\infty}e^{-2iE_{k,0}t^{\prime}}\tilde{f}_{k}(t^{\prime})\,\text{d}t^{\prime}$
(63a) $\displaystyle=-\varphi\varepsilon_{k,0},$ (63b)

which leads to

$s_{k,1}^{y}(t)=-2\frac{\Delta_{0}}{E_{k,0}}\varphi\varepsilon_{k,0}\cos\left(2E_{k,0}t\right).$
(64)

We compare this approximation to the numerical result in Fig. 19. Here again,
the fast oscillations are captured well by the cosine function, but the
envelope deviates from the numerical result. For both pulse shapes, there
seems to be some additional envelope function missing. This might be caused by
neglecting $D_{2}$ in (47b). In addition, we recall that the employed
approximation is valid for small pulse depths which is not the appropriate
regime for the application to a $\delta$ pulse.

Figure 19: Comparison of the analytical approximation and the numerical result
for a $\delta$ pulse at ${t_{\text{del}}=$100\text{\,}\mathrm{p}\mathrm{s}$}$
for amplitude ${\varphi=$0.01\text{\,}\mathrm{p}\mathrm{s}$}$ corresponding to
a rectangular pulse of duration
${t_{\text{puls}}=$1\text{\,}\mathrm{p}\mathrm{s}$}$ and pulse depth
${1-J_{0}=$0.01$}$.

## References

  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957).
  * Dayem and Wiegand [1967] A. H. Dayem and J. J. Wiegand, Phys. Rev. 155, 419 (1967).
  * Dahlberg  _et al._ [1979] E. D. Dahlberg, R. L. Orbach, and I. Schuller, Journal of Low Temperature Physics 36, 367 (1979).
  * Parker and Williams [1972] W. H. Parker and W. D. Williams, Phys. Rev. Lett. 29, 924 (1972).
  * Beyer  _et al._ [2011] M. Beyer, D. Städter, M. Beck, H. Schäfer, V. V. Kabanov, G. Logvenov, I. Bozovic, G. Koren, and J. Demsar, Phys. Rev. B 83, 214515 (2011).
  * Beck  _et al._ [2013] M. Beck, I. Rousseau, M. Klammer, P. Leiderer, M. Mittendorff, S. Winnerl, M. Helm, G. N. Gol’tsman, and J. Demsar, Phys. Rev. Lett. 110, 267003 (2013).
  * Yuzbashyan  _et al._ [2005a] E. A. Yuzbashyan, B. L. Altshuler, V. B. Kuznetsov, and V. Z. Enolskii, Phys. Rev. B 72, 220503 (2005a).
  * Dzero  _et al._ [2009] M. Dzero, E. A. Yuzbashyan, and B. L. Altshuler, Europhysics Letters 85, 20004 (2009).
  * Kirmani and Dzero [2019] A. A. Kirmani and M. Dzero, Journal of Superconductivity and Novel Magnetism 32, 3473 (2019).
  * Cui  _et al._ [2019a] T. Cui, M. Schütt, P. P. Orth, and R. M. Fernandes, Phys. Rev. B 100, 144513 (2019a).
  * Cui  _et al._ [2019b] T. Cui, X. Yang, C. Vaswani, J. Wang, R. M. Fernandes, and P. P. Orth, Phys. Rev. B 100, 054504 (2019b).
  * Tsuji and Aoki [2015] N. Tsuji and H. Aoki, Phys. Rev. B 92, 064508 (2015).
  * Müller  _et al._ [2019] M. A. Müller, P. A. Volkov, I. Paul, and I. M. Eremin, Phys. Rev. B 100, 140501 (2019).
  * Derendorf  _et al._ [2022] P. Derendorf, M. A. Müller, and I. M. Eremin, Faraday Discussions 237, 186 (2022).
  * Gangopadhyay  _et al._ [2010] A. Gangopadhyay, M. Dzero, and V. Galitski, Phys. Rev. B 82, 024303 (2010).
  * Li and Dzero [2024] Y. Li and M. Dzero, Phys. Rev. B 109, 054520 (2024).
  * Derendorf  _et al._ [2024] P. Derendorf, A. F. Volkov, and I. M. Eremin, Phys. Rev. B 109, 024510 (2024).
  * Yuzbashyan  _et al._ [2005b] E. A. Yuzbashyan, B. L. Altshuler, V. B. Kuznetsov, and V. Z. Enolskii, Journal of Physics A: Mathematical and General 38, 7831 (2005b).
  * Greilich  _et al._ [2006] A. Greilich, D. R. Yakovlev, A. Shabaev, A. L. Efros, I. A. Yugova, R. Oulton, V. Stavarache, D. Reuter, A. Wieck, and M. Bayer, Science 313, 341 (2006).
  * Greilich  _et al._ [2007a] A. Greilich, A. Shabaev, D. R. Yakovlev, A. L. Efros, I. A. Yugova, D. Reuter, A. D. Wieck, and M. Bayer, Science 317, 1896 (2007a).
  * Greilich  _et al._ [2007b] A. Greilich, M. Wiemann, F. G. G. Hernandez, D. R. Yakovlev, I. A. Yugova, M. Bayer, A. Shabaev, A. L. Efros, D. Reuter, and A. D. Wieck, Phys. Rev. B 75, 233301 (2007b).
  * Beugeling  _et al._ [2016] W. Beugeling, G. S. Uhrig, and F. B. Anders, Phys. Rev. B 94, 245308 (2016).
  * Jäschke  _et al._ [2017] N. Jäschke, A. Fischer, E. Evers, V. V. Belykh, A. Greilich, M. Bayer, and F. B. Anders, Phys. Rev. B 96, 205419 (2017).
  * Beugeling  _et al._ [2017] W. Beugeling, G. S. Uhrig, and F. B. Anders, Phys. Rev. B 96, 115303 (2017).
  * Kleinjohann  _et al._ [2018] I. Kleinjohann, E. Evers, P. Schering, A. Greilich, G. S. Uhrig, M. Bayer, and F. B. Anders, Phys. Rev. B 98, 155318 (2018).
  * Schering  _et al._ [2020] P. Schering, P. W. Scherer, and G. S. Uhrig, Phys. Rev. B 102, 115301 (2020).
  * Uhrig [2020] G. S. Uhrig, SciPost Phys. 8, 040 (2020).
  * Krull  _et al._ [2012] H. Krull, N. A. Drescher, and G. S. Uhrig, Phys. Rev. B 86, 125113 (2012).
  * Pekker and Varma [2015] D. Pekker and C. M. Varma, Ann. Rev. Cond. Matt. Phys. 6, 269 (2015).
  * Krull  _et al._ [2016] H. Krull, N. Bittner, G. S. Uhrig, D. Manske, and A. P. Schnyder, Nat. Comm. 7, 11921 (2016).
  * Anderson [1958] P. W. Anderson, Phys. Rev. 112, 1900 (1958).
  * Barankov and Levitov [2006] R. A. Barankov and L. S. Levitov, Phys. Rev. Lett. 96, 230403 (2006).
  * Peierls [1933] R. Peierls, Zeitschrift für Physik 80, 763 (1933).
  * Jackson and Okun [2001] J. D. Jackson and L. B. Okun, Rev. Mod. Phys. 73, 663 (2001).
  * Kalthoff  _et al._ [2018] M. H. Kalthoff, G. S. Uhrig, and J. K. Freericks, Phys. Rev. B 98, 035138 (2018).
  * Blanes  _et al._ [2009] S. Blanes, F. Casas, J. A. Oteo, and J. Ros, Phys. Rep. 470, 151 (2009).
  * Matsunaga and Shimano [2012] R. Matsunaga and R. Shimano, Phys. Rev. Lett. 109, 187002 (2012).
  * Bossini  _et al._ [2014] D. Bossini, A. M. Kalashnikova, R. V. Pisarev, T. Rasing, and A. V. Kimel, Phys. Rev. B 89, 060405 (2014).
  * Yuzbashyan  _et al._ [2006] E. A. Yuzbashyan, O. Tsyplyatyev, and B. L. Altshuler, Phys. Rev. Lett. 96, 097005 (2006).
  * Müller  _et al._ [2018] M. A. Müller, P. Shen, M. Dzero, and I. Eremin, Phys. Rev. B 98, 024522 (2018).
  * Akbari  _et al._ [2013] A. Akbari, A. P. Schnyder, D. Manske, and I. Eremin, Europhysics Letters 101, 17002 (2013).
  * Gross  _et al._ [1991] E. K. U. Gross, E. Runge, and O. Heinonen,  _Many-Particle Theory_ (Adam Hilger, 1991).