conclusions here. Let $N$ be even. Our first step is to calculate the reduced
density matrix $\rho_{A}$ for the subsystem $A$ consisting of $\frac{N}{2}$
spins for the early-time many-body wave function (103),

$\rho_{A}=2^{-\frac{N}{2}}\sum_{K,K^{\prime}=0}^{\frac{N}{2}}\Gamma_{K}\Gamma^{*}_{K^{\prime}}e^{-\frac{\tau^{2}}{4N}(K-K^{\prime})^{2}}\lvert
K\rangle\langle K^{\prime}\rvert_{A},$ (106)

where

$\Gamma_{K}=\mbinom{\frac{N}{2}}{K}^{\frac{1}{2}}e^{-\frac{i\tau}{N}\left[K^{2}-K\left(\frac{3N}{2}+1\right)\right]}$
(107)

and $|K\rangle_{A}$ is the state of the subsystem $A$ with a definite number
$K$ of up spins that is symmetric with respect to arbitrary permutations of
spins, i.e., with the maximum possible total spin $\frac{N}{4}$ and a definite
value $K-\frac{N}{4}$ of its $z$-projection. Using this result, we evaluate
$\mathrm{Tr}\left[\rho_{A}^{n}\right]$ in the limit $N\to\infty$ by replacing
the sums with integrals and employing the multi-dimensional saddle point
method,

$\mathrm{Tr}\left[\rho_{A}^{n}\right]=\prod_{j=0}^{n-1}\left[1+\frac{\tau^{2}}{4}\sin^{2}\left(\frac{\pi
j}{n}\right)\right]^{-\frac{1}{2}}.$ (108)

Now consider the Rényi entanglement entropies for integer $n$ defined as

$\mathtt{S}^{\mathrm{R}}_{n}=\frac{\ln\mathrm{Tr}\left[\rho_{A}^{n}\right]}{1-n}.$
(109)

Treating $n$ as a replica index that can be analytically continued [85], we
obtain the von Neumann entanglement entropy in the $n\to 1$ limit as

$\begin{gathered}\mathtt{S}_{\mathrm{ent}}=\sqrt{1+\frac{\tau^{2}}{4}}\coth^{-1}\left[\sqrt{1+\frac{\tau^{2}}{4}}\right]+\ln\frac{\tau}{4}.\end{gathered}$
(110)

We see that $\mathtt{S}_{\mathrm{ent}}$ is intensive in the thermodynamic
limit. It monotonically grows from zero at $\tau=0$ $(t=t_{0})$ behaving as
$\mathtt{S}_{\mathrm{ent}}\approx\ln\tau$ at large $\tau$.

There is a simple picture of the early-time dynamics that explains Eq. (110)
and describes the saturation and subsequent oscillations of the entanglement
entropy for finite $N$ seen in Fig. 10. First of all it is not difficult to
see that $\mathtt{S}_{\mathrm{ent}}$ is bounded from above by
$\ln\left(\frac{N}{2}+1\right)$. To show this, observe that the Hamiltonian
with which the system evolves with respect to $\tau$ is
$\hat{\jmath}_{+}\hat{\jmath}_{-}$ up to a multiplicative constant, see Eq.
(102). The collective spin
$\hat{\bm{\jmath}}=\hat{\bm{\jmath}}_{A}+\hat{\bm{\jmath}}_{\bar{A}}$, where
$\hat{\bm{\jmath}}_{A}$ and $\hat{\bm{\jmath}}_{\bar{A}}$ are the total spins
of subsystems $A$ and $\bar{A}$. Initially all spins-$\frac{1}{2}$ are along
the $x$-axis and therefore the magnitudes of $\hat{\bm{\jmath}}$,
$\hat{\bm{\jmath}}_{A}$, and $\hat{\bm{\jmath}}_{\bar{A}}$ are
$\jmath=\frac{N}{2}$, $\jmath_{A}=\frac{N}{4}$, and
$\jmath_{\bar{A}}=\frac{N}{4}$. These magnitudes are conserved. As a result
the dynamics of subsystem $A$ is confined to a Hilbert space of dimension
$\mathrm{dim}\mathcal{H}_{A}=2\jmath_{A}+1=\frac{N}{2}+1$. It is well-known
[86] that the entanglement entropy is bounded from above by
$\ln(\mathrm{dim}\mathcal{H}_{A})$, i.e.,
$\mathtt{S}_{\mathrm{ent}}\leq\ln\left(\frac{N}{2}+1\right)$.

We see that the $\ln N$ scaling of the entanglement entropy in the BCS theory
(see also Sec. VII.2) is due to all to all interactions, which give rise to
the collective spin $\hat{\bm{\jmath}}$. Note also that the BCS order
parameter is related to the collective spin as
$\Delta=g\langle\hat{\jmath}_{-}\rangle$. The $\ln N$ behavior of
$\mathtt{S}_{\mathrm{ent}}$ is a general feature of late-time asymptotic
states of the BCS dynamics as well as the ground state and other low energy
stationary states unrelated to the nonautonomous character of the Hamiltonian
we study in this paper. We nevertheless discuss it for completeness.

A useful general result [87] for the von Neumann entanglement entropy in
collective spin models, such as the early-time Hamiltonian (101), is

$\mathtt{S}_{\mathrm{ent}}=\sqrt{1+\langle\hat{n}_{\mathrm{ex}}\rangle}\coth^{-1}\\!\left[\sqrt{1+\langle\hat{n}_{\mathrm{ex}}\rangle}\,\right]+\frac{1}{2}\ln\frac{\langle\hat{n}_{\mathrm{ex}}\rangle}{4},$
(111)

where
$\langle\hat{n}_{\mathrm{ex}}\rangle=\langle\hat{b}^{\dagger}\hat{b}\rangle$
is the number of excitations – the number of Holstein-Primakoff bosons for the
collective spin $\hat{\bm{\jmath}}$ bosonized via a Holstein-Primakoff
transformation around the direction of
$\langle\hat{\bm{\jmath}}(\tau)\rangle$. To gain further insight into various
features of $\mathtt{S}_{\mathrm{ent}}$, consider the semiclassical motion of
this bosonic mode. This motion is one-dimensional and Hamiltonian. Nearby
trajectories separate linearly in time, i.e., the growth of the momentum and
position with $\tau$ is linear. Since $\langle\hat{n}_{\mathrm{ex}}\rangle$ is
quadratic in the momentum and position, we expect
$\langle\hat{n}_{\mathrm{ex}}\rangle\sim\tau^{2}$. Indeed, comparing Eqs.
(111) and (110), we see that
$\langle\hat{n}_{\mathrm{ex}}\rangle=\frac{\tau^{2}}{4}$.

Bounded one-dimensional Hamiltonian motion is periodic, which explains the
post-saturation oscillations of $\mathtt{S}_{\mathrm{ent}}$ (see Fig. 10)
approximately with the period of the underlying classical trajectory [87]. The
time at which quantum fluctuations of the collective spin $\hat{\bm{\jmath}}$
become comparable to its magnitude, i.e.,
$\langle\hat{n}_{\mathrm{ex}}\rangle\sim\frac{N}{2}$, is the Ehrenfest time
$\tau_{E}$. Using $\langle\hat{n}_{\mathrm{ex}}\rangle=\frac{\tau^{2}}{4}$ and
$\tau=\eta\ln\frac{t}{t_{0}}$, we obtain

$\tau_{E}=c_{E}\sqrt{2N},\quad t_{E}=t_{0}e^{c_{E}\frac{\sqrt{2N}}{\eta}},$
(112)

where $c_{E}$ is a coefficient of order one. Recall that the entanglement
entropy in the thermodynamic limit (110) grows as $\ln\tau$ for large $\tau$.
This growth stops and finite size effects kick in at the Ehrenfest time scale
as at this point the number of bosonic excitations reaches its maximum
possible value. Note also that $t_{E}$ is the time when the argument of the
exponential function in Eq. (104) becomes of order one, i.e., the quantum
dephasing effect becomes appreciable. Numerically, we observe that Eq. (112)
provides a reasonable estimate of the time when the finite $N$ entanglement
entropy deviates from the $N\to\infty$ result (110), see, e.g., the $t=t_{E}$
line for $N=10$ and $c_{E}=1.12$ in Fig. 10 (dash-dotted vertical line). We
can also estimate the saturation value of the von Neumann entanglement entropy
for large $N$ as $\mathtt{S}_{\mathrm{ent}}\sim\ln\tau_{E}\sim\frac{1}{2}\ln
N$, cf. Eqs. (72) and (76).

##  XIII Conclusion

In this paper, we demonstrated that the far from equilibrium dynamics of BCS
superconductors is classical in the thermodynamic limit under certain
conditions. Specifically, we obtained exact solutions for the quantum and
classical (mean-field) dynamics of the BCS Hamiltonian with time-dependent
coupling, $g(t)=\frac{1}{\nu t}$ launched from the ground state at $t=t_{0}\to
0^{+}$. We explicitly determined exact quantum and mean-field wave functions
at long times, evaluated quantum averages of a generic local observable in
them, and proved that they coincide in the thermodynamic limit. It is clear
that this must remain true for a broad class of $g(t)$ and initial conditions.
Nevertheless, it is worthwhile to verify that the mean field is similarly
exact when $g(t)$ does not vanish at $t\to+\infty$. This can be done, for
example, by solving for the backward time evolution from $t=+\infty$ to
$t=0^{+}$ of our model using the method we developed in this paper.

On the other hand, the classical picture breaks down dramatically for global
quantities, such as the bipartite von Neumann entanglement entropy as we saw
above and the Loschmidt echo as shown in Ref. 39. Also noteworthy is the
behavior of anomalous averages – expectation values of operators that do not
conserve the total fermion number, such as the equal time anomalous Green’s
function. For these kinds of observables, the thermodynamic limit does not
commute with the adiabatic ($\nu\to 0^{+}$) and $t_{0}\to 0^{+}$ limits due to
quantum dephasing. Their quantum fluctuations are controlled by a new
parameter that can be much larger than the inverse particle number – the
parameter controlling the magnitude of equilibrium quantum fluctuations
(finite size corrections).

These results provide a deeper understanding of the reasons behind the success
of mean-field theories in general, beyond our focus on far from equilibrium
superconductivity. In situations where they are believed to be accurate (e.g.,
above a certain dimension), we expect the mean-field wave function to capture
the order parameter and other local observables, but not the entanglement and
other global properties of the true many-particle state. A fascinating project
is to investigate the interplay between quantumness and nonlocality in other
“mean-field” models, such as topological $p$-wave superconductors [7, 61],
infinite-dimensional Hubbard model [88], and magnon Bose-Einstein condensates
[89, 90]. In particular, an interesting question to ask in this connection is
how well the mean-field approximation describes non-trivial topological
properties, which are inherently non-local.

We saw that the unitary time evolution brings our system into a steady state
similar to one of the nonequilibrium phases in interaction quenched BCS
superfluids [35]. This state is a gapless superconductor with vanishing
superfluid density, order parameter, and pair-breaking excitation gap. As a
consequence, its superfluid properties can only be revealed in energy resolved
observables, such as the spectral supercurrent density. The steady state is
nonthermal but is described by a nontrivial emergent generalized Gibbs
ensemble, where the emergent integrals of motion are the single-particle level
occupancy operators. We only touched on the prospects of realizing our model
and testing our predictions in experiment in Introduction and Sec. III, so
this remains an important topic for future research.

We determined the dynamics of the von Neumann entanglement entropy
$\mathtt{S}_{\mathrm{ent}}$ starting from the unentangled BCS (mean-field)
ground state. Interestingly, the growth of the entanglement is entirely due to
the interaction part of the quantum BCS Hamiltonian. The entanglement entropy
is finite in the thermodynamic limit and grows monotonically as a function of
time. Note that within the mean-field treatment $\mathtt{S}_{\mathrm{ent}}=0$
at all times. In a finite system, $\mathtt{S}_{\mathrm{ent}}$ saturates at a
value that scales as $\ln V$ with the system volume $V$.

This paper paves the way to a comprehensive theory of integrability of
nonautonomous quantum Hamiltonians. We demonstrated that the off-shell Bethe
Ansatz [53, 52, 55, 56], which provides an integral representation of
solutions of the nonstationary Schrödinger equation, is a key ingredient of
this theory. Another key ingredient is the systematic method we developed here
to obtain explicit physical results, such as the ones listed above in this
section, from the integral representation. Our method should work equally well
for other models that go through the off-shell Bethe Ansatz, e.g., the problem
of molecular production in an atomic Fermi gas swept through a Feshbach
resonance and Demkov-Osherov, bow-tie, and generalized bow-tie multi-level
Landau-Zener models [52]. It is interesting to apply our approach to these
models and even more important to generalize and broaden the scope of the off-
shell Bethe Ansatz to include time-dependent models unrelated to the Gaudin
algebra.

###### Acknowledgements.

We are grateful to V. Gurarie for helpful discussions. A.W. and J.H.P. are
partially supported by NSF CAREER Grant No. DMR1941569 and the Alfred P. Sloan
Foundation through a Sloan Research Fellowship. We acknowledge the Beowulf
cluster at the Department of Physics and Astronomy of Rutgers University, and
the Rusty cluster from Flatiron Institute Research Computing, Simons
Foundation that were used in obtaining the numerical results.

##  Appendix A Solution of the classical equations of motion

As discussed in the main text, there are two ways to make Anderson pseudospins
classical. One option is to replace spin-$\frac{1}{2}$ with spin-$s$ and take
the limit $s\to\infty$. The other option is to perform mean-field decoupling,
$\langle\hat{A}_{1}\hat{A}_{2}\rangle\to\langle\hat{A}_{1}\rangle\langle\hat{A}_{2}\rangle$,
in the Heisenberg equations of motion for Anderson pseudospins. Both
approaches result in the same equations of motion (18). Here we solve these
classical equations of motion exactly and determine classical spins
$\bm{S}_{j}$ explicitly at long times. In the main text, we compare this
answer with the thermodynamic limit of the exact quantum (spin-$\frac{1}{2}$)
solution. We note that our result in this appendix also solves the generalized
SU(2) Knizhnik-Zamolodchikov equations in the limit of classical spins and the
boundary field $B\equiv\nu t\to+\infty$ [52, 55, 56]. Because we expect
classical (mean-field) dynamics to be able to reproduce quantum dynamics only
for $N\to\infty$, we primarily focus on this limit, even though our method is
suitable for finite $N$ as well.

Equation (33) is an exact integral representation for the solution of the
nonstationary Schrödinger equation for spins of arbitrary magnitudes. Let
magnitudes of all spins be $s$. Our aim is to take the classical limit
$s\to\infty$ and $\hbar\to 0$ of the solution keeping $\hbar s=S$ fixed. The
first step is to restore $\hbar$ in the Schrödinger equation,

$i\hbar\frac{d\Psi}{dt}=\Biggl{[}\sum_{j=1}^{N}2\varepsilon_{j}(\hbar\hat{s}_{j}^{z})-\frac{1}{\nu
t}\sum_{j,k=1}^{N}(\hbar\hat{s}_{j}^{+})(\hbar\hat{s}_{k}^{-})\Biggr{]}\Psi.$
(113)

Note that $\hbar$ is dimensionless in our units. Canceling one factor of
$\hbar$ on both sides and comparing to Eq. (13) with $g(t)=\frac{1}{\nu t}$,
we see that the restoration of $\hbar$ amounts to the replacement
$\nu\to\frac{\nu}{\hbar}$. Making this replacement in Eq. (33), we further
observe that the Yang-Yang action becomes,

${\cal S}(\bm{\lambda},\bm{\varepsilon},t)=2\nu
t\sum_{\alpha}\lambda_{\alpha}+2S\sum_{j}\sum_{\alpha}\ln(\varepsilon_{j}-\lambda_{\alpha})-\hbar\sum_{\alpha}\sum_{\beta\neq\alpha}\ln(\lambda_{\beta}-\lambda_{\alpha}),$
(114)

where

$S=\hbar s.$ (115)

Expressions for $\Psi(t)$, $\Xi(\bm{\lambda},\varepsilon)$, and
$\hat{L}^{+}(\lambda)$ stay the same.

Similar to the spin-$\frac{1}{2}$ case in Sec. V, the integral in Eq. (33)
localizes to the vicinities of the stationary points. The stationary point
equations $\frac{\partial\mathcal{S}}{\partial{\lambda_{\alpha}}}=0$ now are

$\nu
t+S\sum_{j}\frac{1}{\lambda_{\alpha}-\varepsilon_{j}}=\hbar\sum_{\beta\neq\alpha}\frac{1}{\lambda_{\alpha}-\lambda_{\beta}},\quad\alpha=1,\dots,N_{\uparrow}.$
(116)

The difference with Eq. (36) is in the number $N_{\uparrow}$ of rapidities
$\lambda_{\alpha}$. In Eq. (36), $N_{\uparrow}$ is the number of up spins
which is smaller than the number $N$ of $\varepsilon_{j}$. In Eq. (116),
$N_{\uparrow}$ is the amount by which the $z$-component $\jmath_{z}$ of the
total spin is raised, $\jmath_{z}=N_{\uparrow}-Ns$. Since the magnitude $s$ of
spins diverges, $N_{\uparrow}$ also diverges and there are many more
$\lambda_{\alpha}$ than $\varepsilon_{j}$.

As before, each $\lambda_{\alpha}$ must tend to one of $\varepsilon_{j}$ as
$t\to+\infty$. Suppose $n_{j}$ of $\lambda_{\alpha}$ tend to
$\varepsilon_{j}$. We denote the elements of this $j^{\mathrm{th}}$ degenerate
subset of $\lambda_{\alpha}$ as $\lambda_{k}^{j}$, where $k=1,\dots,n_{j}$.
Note that $\sum_{j=1}^{N}n_{j}=N_{\uparrow}$. Let

$\lambda_{k}^{j}=\varepsilon_{j}+\chi_{k}^{j},\quad\chi_{k}^{j}=\frac{\hbar
z_{k}^{j}}{2\nu t}.$ (117)

We will see shortly that the new variables $z_{k}^{j}$ are independent of $t$
to the leading order in $t^{-1}$. All terms in Eq. (116) that contain
$\lambda_{\alpha}-\varepsilon_{j}$ with $\lambda_{\alpha}$ not in the
degenerate subset corresponding to $\varepsilon_{j}$ or terms that contain
$\lambda_{\alpha}-\lambda_{\beta}$ with $\lambda_{\alpha}$ and
$\lambda_{\beta}$ not in the same degenerate subset are negligible. Therefore,
to the leading order in $t^{-1}$, Eq. (116) splits into $N$ decoupled sets of
equations for $z_{k}^{j}$,

$1+\frac{2s}{z_{k}^{j}}=\sum_{k^{\prime}\neq
k}\frac{2}{z_{k}^{j}-z_{k^{\prime}}^{j}},\quad k=1,\dots,n_{j}.$ (118)

This is a set of $n_{j}$ equations for $n_{j}$ variables for each
$j=1,\dots,N$. It is clear from Eq. (118) that $z_{k}^{j}$ are
$t$-independent. In fact, $z_{k}^{j}$ are the roots of the associated Laguerre
polynomial $L^{-2s-1}_{n_{j}}(z)$ [70]. We will not need this property below
as we will not evaluate integrals over $\chi_{k}^{j}$ by the steepest descent
method.

Let us see how the vector $\Xi(\bm{\lambda},\bm{\varepsilon})$ given by Eq.
(34) behaves near stationary points. Observe that
$\hat{L}^{+}(\lambda_{\alpha})\to\frac{\hat{s}_{j}^{+}}{\chi_{k}^{r}}$ when
$\lambda_{\alpha}\to\varepsilon_{r}$. Therefore,

$\Xi(\bm{\lambda},\bm{\varepsilon})\to\sum_{\\{n_{j}\\}}\left[\frac{N_{\uparrow}!}{\prod_{j}n_{j}!}\prod_{j}\frac{(\hat{s}_{j}^{+})^{n_{j}}}{\prod_{k}\chi_{k}^{j}}\right]|0\rangle.$
(119)

The summation here is over all sets $\\{n_{j}\\}$ of $N$ positive integers
such that

$\sum_{j=1}^{N}n_{j}=N_{\uparrow}.$ (120)

This corresponds to summing over all stationary points. The combinatorial
factor in Eq. (119) is the number of ways to choose $N$ groups of variables
out of $N_{\uparrow}$ variables with $n_{j}$ elements in the $j^{\mathrm{th}}$
group. The action of powers of spin-$s$ raising operators $s_{r}^{+}$ on the
state $|0\rangle$ where all spins point in the negative $z$-direction is

$\prod_{r=1}^{N}(\hat{s}_{r}^{+})^{n_{r}}|0\rangle=|n_{1}\dots
n_{N}\rangle\prod_{r=1}^{N}\sqrt{\frac{(2s)!n_{r}!}{(2s-n_{r})!}},$ (121)

where $|n_{1}\dots n_{N}\rangle$ is a normalized eigenstate of all
$\hat{s}_{r}^{z}$ with eigenvalues $m_{r}=n_{r}-s$. Eq. (119) becomes

$\Xi(\bm{\lambda},\bm{\varepsilon})\to\sum_{\\{n_{j}\\}}\prod_{j=1}^{N}\sqrt{\frac{(2s)!}{(2s-n_{j})!n_{j}!}}\,\,\frac{|n_{1}\dots
n_{N}\rangle}{\prod_{k}\chi_{k}^{j}}.$ (122)

We manipulate the Yang-Yang action (114) similarly to how we manipulated the
stationary point equations (116) neglecting $\chi_{k}^{l}$ in
$\varepsilon_{j}-\lambda_{\alpha}$ when $\lambda_{\alpha}$ does not belong to
the $j^{\mathrm{th}}$ degenerate subset and in
$\lambda_{\alpha}-\lambda_{\beta}$ when $\lambda_{\alpha}$ and
$\lambda_{\beta}$ are not in the same degenerate subset. We find

${\cal S}(\bm{\chi},\bm{\varepsilon},t)=2\nu
t\sum_{j}n_{j}\varepsilon_{j}+2S\sum_{j}\sum_{k\neq
j}n_{k}\ln(\varepsilon_{j}-\varepsilon_{k})-\hbar\sum_{j}\sum_{k\neq
j}n_{j}n_{k}\ln(\varepsilon_{j}-\varepsilon_{k})+\sum_{j}\mathcal{S}^{j}_{\mathrm{deg}},$
(123)

where

$\mathcal{S}^{j}_{\mathrm{deg}}=\hbar\sum_{k}z_{k}^{j}+2S\sum_{k}\ln(-\chi_{k}^{j})-\hbar\sum_{k}\sum_{k^{\prime}\neq
k}\ln(\chi_{k}^{j}-\chi_{k^{\prime}}^{j}).$ (124)

Recall that $\nu$ is of order $N$ and therefore the first three terms in
${\cal S}(\bm{\chi},\bm{\varepsilon},t)$ are of order $N^{2}$ because each
summation over the single-particle level index gives a factor of $N$. On the
other hand, sums over the degenerate subspace in
$\mathcal{S}^{j}_{\mathrm{deg}}$ are of order $s$ and consequently
$\sum_{j}\mathcal{S}^{j}_{\mathrm{deg}}$ is of order $Ns$. Take the limits
$N\to\infty$ and $s\to\infty$ so that $\frac{s}{N}\to 0$. The last term in Eq.
(123) is then negligible when $N\to\infty$. As for the case of
spin-$\frac{1}{2}$ in Sec. V, we need to choose the branch of the logarithm so
that $\ln(-1)=-i\pi$. It is not difficult to show that to order $N^{2}$ the
imaginary part of the third term on the right hand side of Eq. (123)
contributes only to the overall normalization constant. Together these
observations allow us to rewrite Eq. (123) as

$\hbar{\cal S}(\bm{\chi},\bm{\varepsilon},t)=2\nu
t\sum_{j}S_{j}^{z}\varepsilon_{j}-2\sum_{k>j}S_{j}^{z}S_{k}^{z}\ln|\varepsilon_{j}-\varepsilon_{k}|-2i\pi
S\sum_{j}jS_{j}^{z},$ (125)

where $S_{j}^{z}=\hbar m_{j}$ are the $z$-components of the classical spins in
the limit $\hbar\to 0$.

We also need to express the square root of the binomial coefficient in Eq.
(122) in terms of $S_{j}^{z}=\hbar m_{j}=\hbar n_{j}-S$. Using Stirling’s
approximation, we obtain

$\sqrt{\frac{(2s)!}{(2s-n_{j})!n_{j}!}}=\exp\left[-\frac{(S+S_{j}^{z})\ln(S+S_{j}^{z})+(S-S_{j}^{z})\ln(S-S_{j}^{z})}{2\hbar}\right].$
(126)

Let us summarize what we did so far. We evaluated the ingredients in the
general expression (33) for $\Psi(t)$ in the vicinity of the stationary points
of the Yang-Yang action $\mathcal{S}(\bm{\lambda},\bm{\varepsilon},t)$ since
at large $t$ the integral localizes to these vicinities. In Eq. (126), we took
advantage of the fact that in the classical limit $\hbar\to 0$ the magnitude
of the quantum spin $\hat{\bm{s}}^{z}_{j}$ diverges, $s\to\infty$, while the
magnitude $S=\hbar s$ and components, e.g., $S_{j}^{z}=\hbar m_{j}$, of the
classical spin $\bm{S}_{j}$ remain finite. Substituting Eqs. (122) and (125)
into Eq. (33) and using Eq. (126), we find

$\Psi(t)=\sum_{\\{m_{j}\\}}e^{-\frac{iEt}{\hbar}}e^{\frac{i\Omega}{\hbar\nu}}e^{-\frac{A}{2\hbar}}|m_{1}\dots
m_{N}\rangle,$ (127)

where the sum is over all $z$-projections $m_{j}$ such that
$\sum_{j}m_{j}=\jmath_{z}$ and we integrated over all $\chi_{k}^{j}$ to derive
Eq. (127) . Since we neglected $\mathcal{S}^{j}_{\mathrm{deg}}$, these
integrals are $\oint\frac{d\chi_{k}^{j}}{\chi_{k}^{j}}=2\pi i$ and contribute
only to the overall multiplicative constant. The contours of integration over
$\chi_{k}^{j}$ are guaranteed to enclose the origin as we start the time
evolution in the instantaneous ground state at $t=0^{+}$, see the discussion
below Eq. (36). Quantities $A$, $E$, and $\Omega$ are

$\displaystyle E$ $\displaystyle=\sum_{j}2\varepsilon_{j}S_{j}^{z},$ (128a)
$\displaystyle\Omega$
$\displaystyle=2\sum_{k>j}S_{j}^{z}S_{k}^{z}\ln|\varepsilon_{j}-\varepsilon_{k}|,$
(128b) $\displaystyle A$
$\displaystyle=(S+S_{j}^{z})\ln(S+S_{j}^{z})+(S-S_{j}^{z})\ln(S-S_{j}^{z})+\frac{4\pi
jSS_{j}^{z}}{\nu}.$ (128c)

In the limit $\hbar\to 0$ only terms that minimize $A$ survive in Eq. (127).
Minimizing Eq. (128c) with respect to $S_{j}^{z}$ at fixed
$\sum_{j}S_{j}^{z}$, we obtain

$S_{j}^{z}=-S\tanh(2Sa_{j}),\quad a_{j}\equiv\frac{\pi(j-\mu)}{\nu},$ (129)

where $\mu$ is the Lagrange multiplier (chemical potential) corresponding to
$\sum_{j}S_{j}^{z}$ and we used
$\tanh^{-1}z=\frac{1}{2}\ln\left(\frac{1+z}{1-z}\right)$.

The state of the system (127) is now a sum over only $m_{j}$ such that $\hbar
m_{j}\to S_{j}^{z}$ with $S_{j}^{z}$ given by Eq. (129). Note that this is
more than one value of $m_{j}$, because, for example, $m_{j}$ and $m_{j}+1$
correspond to the same $S_{j}^{z}$ in the limit $\hbar\to 0$. This allows us
to also evaluate the $x$ and $y$ components of classical spins by taking the
expectation value $\langle\hat{s}_{j}^{-}\rangle$ of the quantum spin in this
state. We have

$S_{j}^{x}-iS_{j}^{y}\equiv S_{j}^{-}=|S_{j}^{-}|e^{-\frac{i\Delta
Et}{\hbar}}e^{\frac{i\Delta\Omega}{\hbar\nu}},$ (130)

where $\Delta E=E(m_{j}+1)-E(m_{j})$ and
$\Delta\Omega=\Omega(m_{j}+1)-\Omega(m_{j})$ are the amounts by which $E$ and
$\Omega$ change when we increase $m_{j}$ by 1 or, equivalently, increase
$S_{j}^{z}$ by $\hbar$. We have

$S_{j}^{-}=\frac{Se^{-2i\varepsilon_{j}t+i\varphi_{j}}}{\cosh(2S\zeta_{j})},\quad\varphi_{j}=\frac{2\eta}{N}\sum_{k\neq
j}S_{k}^{z}\ln|\varepsilon_{k}-\varepsilon_{j}|,\quad
S_{j}^{z}=-S\tanh(2S\zeta_{j}),\quad\zeta_{j}\equiv\frac{\pi\eta(j-\mu)}{N}.$
(131)

Here we traded $\nu$ for $\frac{N}{\eta}$. This is the exact solutions of the
Hamilton’s equations of motion (18) for the classical time-dependent BCS
Hamiltonian (15) for $t\to+\infty$ and $N\to\infty$ (these two limits commute
as we discuss in Sec. IX for the quantum and at the end of Sec. XI for the
classical dynamics). As noted below Eq. (21), to compare to mean-field
dynamics starting from the BCS product state, we need to set the spin length
$S=\frac{1}{2}$ in Eq. (131).

##  Appendix B Rényi entanglement entropy

One of the key differences between quantum and classical systems is the
presence of entanglement in the quantum case, i.e., of statistical
correlations between subsystems that prevent the system from being in a state
that is a product of the states of the individual subsystems. The degree of
quantum entanglement between two subsystems $(A,\bar{A})$ is quantified by the
von Neumann entanglement entropy defined as
$\mathtt{S}_{\mathrm{ent}}=\mathrm{Tr}\rho_{A}\ln\rho_{A}$ where
$\rho_{A}=\mathrm{Tr}_{\bar{A}}\rho$ is the reduced density matrix of $A$ and
$\rho$ is the density matrix of the combined system. Rather than computing
$\mathtt{S}_{\mathrm{ent}}$ directly, it is more convenient to work with the
Rényi entanglement entropy $\mathtt{S}^{\mathrm{R}}_{n}$ defined in Eq. (109).
Its usefulness lies in the index $n$ which can be treated as a replica index
to evaluate the von Neumann entropy in the limit $n\to 1$. The Rényi entropy
also encodes other entanglement measures such as the Hartley entropy for $n=0$
and the purity for $n=2$. In what follows, we derive an expression for the
Rényi entanglement entropy for arbitrary $n$ (and the von Neumann entanglement
entropy by analytic continuation) for a bipartition of the early-time wave
function (103) for the dynamics with the quantum BCS Hamiltonian (11) [or
equivalently with Hamiltonian (13) for $s=\frac{1}{2}$] with
$g(t)=\frac{\eta}{Nt}$ starting from the infinite coupling BCS ground state
(91) at $t=t_{0}\to 0^{+}$.

###  B.1 Reduced density matrix

In order to determine the Rényi entropy $\mathtt{S}^{\mathrm{R}}_{n}$, we
require the reduced density matrix $\rho_{A}$ for a bipartition of the system.
As in Sec. VII.2, we choose $A$ to be the set of spins $\hat{\bm{s}}_{j}$ that
correspond the lower half of the single-particle energy spectrum
$\varepsilon_{j}$. We need to evaluate $\rho_{A}$ for the state (103). To do
so, we first rewrite Eq. (103) in a more convenient form by substituting Eq.
(74) into it and rearranging the coefficients by first splitting the sum into
two terms ($N_{\uparrow}\leq\frac{N}{2}$ and $N_{\uparrow}>\frac{N}{2}$) and
then recombining it,

$\Phi(t)=2^{-\frac{N}{2}}\sum_{N_{\uparrow}^{A}=0}^{\frac{N}{2}}\binom{\frac{N}{2}}{N_{\uparrow}^{A}}^{\frac{1}{2}}\left|N_{\uparrow}^{A}\right\rangle_{A}\otimes\sum_{N_{\uparrow}^{\bar{A}}=0}^{\frac{N}{2}}e^{-iE_{0}(N_{\uparrow})\tau}\binom{\frac{N}{2}}{N_{\uparrow}^{\bar{A}}}^{\frac{1}{2}}\left|N_{\uparrow}^{\bar{A}}\right\rangle_{\bar{A}},$
(132)

where $N_{\uparrow}=N_{\uparrow}^{A}+N_{\uparrow}^{\bar{A}}$ and
$\tau=\eta\ln\frac{t}{t_{0}}$. At $\tau=0$ the system is in the product state
(91) [all spins point along $x$] and we see that the two sums decouple in this
case as they should. The density matrix of the system is then

$\rho(t)=2^{-N}\sum_{K,K^{\prime}=0}^{\frac{N}{2}}\binom{\frac{N}{2}}{K}^{\frac{1}{2}}\binom{\frac{N}{2}}{K^{\prime}}^{\frac{1}{2}}\lvert
K\rangle\langle
K^{\prime}\rvert_{A}\otimes\sum_{Q,Q^{\prime}=0}^{\frac{N}{2}}e^{-i\left[E_{0}(K+Q)-E_{0}(K^{\prime}+Q^{\prime})\right]\tau}\binom{\frac{N}{2}}{Q}^{\frac{1}{2}}\binom{\frac{N}{2}}{Q^{\prime}}^{\frac{1}{2}}\lvert
Q\rangle\langle Q^{\prime}\rvert_{\bar{A}},$ (133)

where we renamed the summation indices $N_{\uparrow}^{A}\to K$ and
$N_{\uparrow}^{\bar{A}}\to Q$ for simplicity. It remains to trace over
$\bar{A}$ to find the reduced density matrix of $A$. Performing the trace
yields a factor of $\delta_{Q,Q^{\prime}}$ which consumes the sum over
$Q^{\prime}$. In the large $N$ limit, we replace the remaining sum over $Q$
with an integral over $z=\frac{2Q}{N}$ and evaluate it using the saddle point
method to arrive at

$\rho_{A}(t)=2^{-\frac{N}{2}}\sum_{K,K^{\prime}=0}^{\frac{N}{2}}\sqrt{\binom{\frac{N}{2}}{K}}e^{-\frac{i\tau}{N}\left[K^{2}-K(\frac{3N}{2}+1)\right]}\sqrt{\binom{\frac{N}{2}}{K^{\prime}}}e^{\frac{i\tau}{N}\left[K^{\prime
2}-K^{\prime}(\frac{3N}{2}+1)\right]}e^{-\frac{\tau^{2}}{4N}(K-K^{\prime})^{2}}\lvert
K\rangle\langle K^{\prime}\rvert_{A}.$ (134)

###  B.2 Rényi entropy

Having determined $\rho_{A}(t)$, we now turn to calculating the Rényi
entanglement entropy defined by Eq. (109). When taking the trace of powers of
$\rho_{A}$ the phases cancel so that

$\mathrm{Tr}\left[\rho_{A}^{n}\right]=2^{-n\frac{N}{2}}\sum_{x_{1},\ldots,x_{n}}^{\frac{N}{2}}\binom{\frac{N}{2}}{x_{1}}\ldots\binom{\frac{N}{2}}{x_{n}}\exp\left[-\frac{\tau^{2}}{4N}\left[(x_{1}-x_{2})^{2}+\ldots+(x_{n}-x_{1})^{2}\right]\right].$
(135)

As before, in the large $N$ limit we replace the sums with integrals and
evaluate them using the multi-dimensional saddle point method. The stationary
point occurs at $x_{1}=x_{2}=\ldots=x_{n}=\frac{1}{2}$ and the matrix elements
of the Hessian matrix are

$\mathrm{Hess}_{ij}=-\left(4+\frac{\tau^{2}}{2}\right)\delta_{i,j}+\frac{\tau^{2}}{4}\delta_{i+1,j}+\frac{\tau^{2}}{4}\delta_{i-1,j}$
(136)

with eigenvalues given by [91]

$\kappa_{j}=-\left(4+\frac{\tau^{2}}{2}\right)+\frac{\tau^{2}}{4}\omega^{j}+\frac{\tau^{2}}{4}\omega^{(n-1)j}$
(137)

where $\omega=e^{\frac{2\pi i}{n}}$ is a primitive $n^{\text{th}}$ root of
unity. Substituting $\prod_{j}\kappa_{j}$ for the determinant into the saddle
point formula and simplifying we have

$\mathrm{Tr}\left[\rho_{A}^{n}\right]=\prod_{j=0}^{n-1}\left[1+\frac{\tau^{2}}{16}\left(2-e^{i\frac{2\pi
j}{n}}-e^{-i\frac{2\pi j}{n}}\right)\right]^{-\frac{1}{2}}$ (138)

from which the Rényi entanglement entropy follows

$\mathtt{S}^{\mathrm{R}}_{n}(\tau)=\frac{-1}{2(1-n)}\sum_{j=0}^{n-1}\ln\left[1+\frac{\tau^{2}}{16}\left(2-e^{i\frac{2\pi
j}{n}}-e^{-i\frac{2\pi j}{n}}\right)\right].$ (139)

To evaluate the von Neumann entanglement entropy we analytically continue to
$n=1$ by writing the sum as a contour integral [85]

$\mathtt{S}^{\mathrm{R}}_{n}(\tau)=\frac{-1}{2(1-n)}\sum_{j=0}^{n-1}\oint\frac{\mathrm{d}u}{2\pi
i}\frac{\ln\left[1+\frac{\tau^{2}}{16}\left(2-u-u^{-1}\right)\right]}{u-e^{i2\pi
j/n}}=\frac{-1}{2(1-n)}\oint\frac{\mathrm{d}u}{2\pi
i}\left(\frac{nu^{n-1}}{u^{n}-1}\right)\ln\left[1+\frac{\tau^{2}}{16}\left(2-u-u^{-1}\right)\right].$
(140)

In the limit $n\to 1$ [85]

$\mathtt{S}^{\mathrm{R}}_{n\to
1}(\tau)\to-\frac{1}{2}\oint\frac{\mathrm{d}u}{2\pi
i}\left(\frac{n}{n-1}\frac{1}{1-u}+\frac{\ln
u}{(u-1)^{2}}+\mathcal{O}(n-1)\right)\ln\left[1+\frac{\tau^{2}}{16}\left(2-u-u^{-1}\right)\right],$
(141)

which becomes

$\mathtt{S}^{\mathrm{R}}_{n\to
1}(\tau)=\mathtt{S}_{\mathrm{ent}}(\tau)=\int_{1}^{\infty}\mathrm{d}\lambda\frac{\ln\left[1+\frac{\tau^{2}}{16}\left(2+\lambda+\lambda^{-1}\right)\right]}{(\lambda+1)^{2}}=\sqrt{1+\frac{\tau^{2}}{4}}\coth^{-1}\sqrt{1+\frac{\tau^{2}}{4}}+\ln\frac{\tau}{4}.$
(142)

An immediate observation is that $\mathtt{S}_{\mathrm{ent}}(\tau)$ is finite
in the thermodynamic limit. In the main text we provide a thorough qualitative
explanation of Eq. (142) and use it together with other considerations to
understand the behavior of $\mathtt{S}_{\mathrm{ent}}(\tau)$ for finite $N$ as
well, namely, its saturation at $\mathtt{S}_{\mathrm{ent}}\sim\frac{1}{2}\ln
N$ at $\tau\sim\sqrt{N}$ and subsequent oscillations as a function of $\tau$.

## References

  * [1] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory of Superconductivity, Phys. Rev. 108, 1175 (1957).
  * [2] P. W. Anderson, Random-Phase Approximation in the Theory of Superconductivity, Phys. Rev. 112, 1900 (1958).
  * [3] E. A. Yuzbashyan, B. L. Altshuler, V. B. Kuznetsov, and V. Z. Enolskii, Nonequilibrium Cooper pairing in the nonadiabatic regime, Phys. Rev. B 72, 220503(R) (2005).
  * [4] P. Coleman, Introduction to Many-Body Physics (Cambridge University Press, 1st ed., 2016).
  * [5] The mean-field ground state and excited states of the condensate correspond to equilibria of the classical spins. Removing or adding a fermion corresponds to removing a classical spin.
  * [6] A. L. DiRienzo and R. A. Young, A coupled angular momentum model for the Josephson junction, Am. J. Phys. 51, 587 (1983).
  * [7] N. Read and D. Green, Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum Hall effect, Phys. Rev. B 61, 10267 (2000).
  * [8] R. W. Richardson, Pairing in the limit of a large number of particles, J. Math. Phys. 18, 1802 (1977).
  * [9] J. Roman, G. Sierra, and J. Dukelsky, Large N limit of the exactly solvable BCS model: Analytics versus numerics, Nucl. Phys. B 634, 483 (2002).
  * [10] E. A. Yuzbashyan, A. A. Baytin, and B. L. Altshuler, Finite-size corrections for the pairing Hamiltonian, Phys. Rev. B 71, 094505 (2005).
  * [11] R. A. Barankov, L. S. Levitov, and B. Z. Spivak, Collective Rabi Oscillations and Solitons in a Time-Dependent BCS Pairing Problem, Phys. Rev. Lett. 93, 160401 (2004).
  * [12] A. V. Andreev, V. Gurarie, and L. Radzihovsky, Nonequilibrium Dynamics and Thermodynamics of a Degenerate Fermi Gas Across a Feshbach Resonance, Phys. Rev. Lett. , 130402 (2004).
  * [13] M. H. Szymanska, B. D. Simons, and K. Burnett, Dynamics of the BCS-BEC Crossover in a Degenerate Fermi Gas, Phys. Rev. Lett. 94, 170402 (2005).
  * [14] E. A. Yuzbashyan, B. L. Altshuler, V. B. Kuznetsov, and V. Z. Enolskii, Solution for the dynamics of the BCS and central spin problems, J. Phys. A 38, 7831 (2005).
  * [15] E. A. Yuzbashyan, O. Tsyplyatyev, and B. L. Altshuler, Relaxation and Persistent Oscillations of the Order Parameter in Fermionic Condensates, Phys. Rev. Lett. 96, 097005 (2006).
  * [16] R. A. Barankov and L. S. Levitov, Synchronization in the BCS Pairing Dynamics as a Critical Phenomenon, Phys. Rev. Lett. 96, 230403 (2006).
  * [17] E. A. Yuzbashyan and M. Dzero, Dynamical Vanishing of the Order Parameter in a Fermionic Condensate, Phys. Rev. Lett. 96, 230404 (2006).
  * [18] M. Dzero, E. A. Yuzbashyan, B. L. Altshuler, and P. Coleman, Spectroscopic Signatures of Nonequilibrium Pairing in Atomic Fermi Gases, Phys. Rev. Lett. 99, 160402 (2007).
  * [19] T. Papenkort, V. M. Axt, and T. Kuhn, Coherent dynamics and pump-probe spectra of BCS superconductors, Phys. Rev. B 76 (2007).
  * [20] A. Tomadin, M. Polini, M. P. Tosi, and R. Fazio, Nonequilibrium pairing instability in ultracold Fermi gases with population imbalance, Phys. Rev. A 77, 033605 (2008).
  * [21] A. Nahum and E. Bettelheim, Dissipationless BCS dynamics with large branch imbalance, Phys. Rev. B 78, 184510 (2008).
  * [22] V. Gurarie, Nonequilibrium Dynamics of Weakly and Strongly Paired Superconductors, Phys. Rev. Lett. 103, 075301 (2009).
  * [23] A. Faribault, P. Calabrese, and J.-S. Caux, Bethe ansatz approach to quench dynamics in the Richardson model, J. Math. Phys. 50, 095212 (2009).
  * [24] A. Faribault, P. Calabrese, and J.-S. Caux, Quantum quenches from integrability: the fermionic pairing model, J. Stat. Mech. P03018 (2009).
  * [25] R. Matsunaga and R. Shimano, Nonequilibrium BCS State Dynamics Induced by Intense Terahertz Pulses in a Superconducting NbN Film, Phys. Rev. Lett. 109, 187002 (2012).
  * [26] C. Sträter, O. Tsyplyatyev, A. Faribault, Nonequilibrum dynamics in the strongly excited inhomogeneous Dicke model, Phys. Rev. B 86, 195101 (2012).
  * [27] M. Beck, I. Rousseau, M. Klammer, P. Leiderer, M. Mittendorff, S. Winnerl, M. Helm, G. N. Gol’tsman, and J. Demsar, Transient Increase of the Energy Gap of Superconducting NbN Thin Films Excited by Resonant Narrow-Band Terahertz Pulses, Phys. Rev. Lett. 110, 267003 (2013).
  * [28] R. Matsunaga, Y. I. Hamada, K. Makise, Y. Uzawa, H. Terai, Z. Wang, and R. Shimano, Higgs Amplitude Mode in the BCS Superconductors Nb1-xTixN Induced by Terahertz Pulse Excitation, Phys. Rev. Lett. 111, 057002 (2013).
  * [29] M. S. Foster, M. Dzero, V. Gurarie, and E. A. Yuzbashyan, Quantum quench in a $p+ip$ superfluid: Winding numbers and topological states far from equilibrium, Phys. Rev. B 88, 104511 (2013).
  * [30] 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).
  * [31] H. Krull, D. Manske, G. S. Uhrig, and A. P. Schnyder, Signatures of nonadiabatic BCS state dynamics in pump-probe conductivity, Phys. Rev. B 90, 014515 (2014).
  * [32] M. S. Foster, V. Gurarie, M. Dzero, and E. A. Yuzbashyan, Quench-Induced Floquet Topological $p$-Wave Superfluids, Phys. Rev. Lett. 113, 076403 (2014).
  * [33] N. Tsuji and H. Aoki, Theory of Anderson pseudospin resonance with Higgs mode in superconductors, Phys. Rev. B 92, 064508 (2015).
  * [34] Y.-Z. Chou, Y. Liao, and M. S. Foster, Twisting Anderson pseudospins with light: Quench dynamics in terahertz-pumped BCS superconductors, Phys. Rev. B 95, 104507 (2017).
  * [35] E. A. Yuzbashyan, M. Dzero, V. Gurarie, and M. S. Foster, Quantum quench phase diagrams of an s-wave BCS-BEC condensate, Phys. Rev. A 91, 033628 (2015).
  * [36] F. Li, V. Y. Chernyak, and N. A. Sinitsyn, Quantum annealing and thermalization: insights from integrability, Phys. Rev. Lett. 121, 190601 (2018).
  * [37] M. Dzero, E. A. Yuzbashyan and B. L. Altshuler, Cooper pair turbulence in atomic Fermi gases, EPL 85, 20004 (2009).
  * [38] More precisely, $n$ must be such that $\frac{n}{N_{\mathrm{f}}}\to 0$ in the thermodynamic limit, where $N_{\mathrm{f}}$ is the total number of fermions.
  * [39] S. Gaur, V. Gurarie, E. A. Yuzbashyan, Singularities in the Loschmidt echo of quenched topological superconductors, arXiv:2207.08131.
  * [40] S. Smale, P. He, B. A. Olsen, K. G. Jackson, H. Sharum, S. Trotzky, J. Marino, A. M. Rey, and J. H. Thywissen, Observation of a transition between dynamical phases in a quantum degenerate Fermi gas, Science Advances 5, 10.1126/sciadv.aax1568 (2019).
  * [41] J. A. Muniz, D. Barberena, R. J. Lewis-Swan, D. J. Young, J. R. K. Cline, A. M. Rey, and J. K. Thompson, Exploring dynamical phase transitions with cold atoms in an optical cavity, Nature 580, 602 (2020).
  * [42] R. J. Lewis-Swan, D. Barberena, J. R. K. Cline, D. J. Young, J. K. Thompson, and A. M. Rey, Cavity-QED Quantum Simulator of Dynamical Phases of a Bardeen-Cooper-Schrieffer Superconductor, Phys. Rev. Lett. 126, 173601 (2021).
  * [43] A. Shankar, E. A. Yuzbashyan, V. Gurarie, P. Zoller, J. J. Bollinger, and A. M. Rey, Simulating dynamical phases of chiral $p+ip$ superconductors with a trapped ion magnet, arXiv:2204.05671.
  * [44] D. Rainer, J. A. Sauls, and D. Waxman, Current carried by bound states of a superconducting vortex, Phys. Rev. B 54, 10094 (1996).
  * [45] S.-K. Yip, Energy-resolved supercurrent between two superconductors, Phys. Rev. B 58, 5803 (1998).
  * [46] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a Completely Integrable Many-Body Quantum System: An Ab Initio Study of the Dynamics of the Highly Excited States of 1D Lattice Hard-Core Bosons Phys. Rev. Lett. 98, 050405 (2007).
  * [47] M. Rigol, A. Muramatsu, and M. Olshanii, Hard-core bosons on optical superlattices: Dynamics and relaxation in the superfluid and insulating regimes, Phys. Rev. A 74, 053616 (2006).
  * [48] P. Mehta and N. Andrei, Nonequilibrium Transport in Quantum Impurity Models: The Bethe Ansatz for Open Systems, Phys. Rev. Lett. 96, 216802 (2006).
  * [49] E. Bettelheim, Towards a Non-equilibrium Bethe ansatz for the Kondo Model, J. Phys. A 48, 165003 (2015).
  * [50] J.-S. Caux and F. H. L. Essler, Time Evolution of Local Observables After Quenching to an Integrable Model, Phys. Rev. Lett. 110, 257203 (2013).
  * [51] B. Wouters, J. De Nardis, M. Brockmann, D. Fioretto, M. Rigol, and J.-S. Caux, Quenching the Anisotropic Heisenberg Chain: Exact Solution and Generalized Gibbs Ensemble Predictions, Phys. Rev. Lett. 113, 117202 (2014).
  * [52] E. A. Yuzbashyan, Integrable time-dependent Hamiltonians, solvable Landau-Zener models and Gaudin magnets, Ann. Phys. 392, 323 (2018).
  * [53] H. M. Babujian, Off-shell Bethe ansatz equations and N-point correlators in the SU(2) WZNW theory, J. Phys. A 26, 6981 (1993).
  * [54] V. G. Knizhnik and A. B. Zamolodchikov, Current algebra and Wess-Zumino model in two dimensions, Nucl. Phys. B 247, 83 (1984).
  * [55] H. M. Babujian and A. V. Kitaev, Generalized Knizhnik-Zamolodchikov equations and isomonodromy quantization of the equations integrable via the Inverse Scattering Transform: Maxwell-Bloch system with pumping, J. Math. Phys. 39, 2499 (1998).
  * [56] D. Fioretto, J.-S. Caux, and V. Gritsev, Exact out-of-equilibrium central spin dynamics from integrability, New J. Phys. 16, 043024 (2014).
  * [57] A. Patra and E. A. Yuzbashyan, Quantum integrability in the multistate Landau-Zener problem, J. Phys. A 48, 245303 (2015).
  * [58] N. A. Sinitsyn, E. A. Yuzbashyan, V. Y. Chernyak, A. Patra, C. Sun, Integrable time-dependent quantum Hamiltonians, Phys. Rev. Lett. 120, 190402 (2018).
  * [59] Since the BCS interaction is all to all, we can order $\varepsilon_{j}$ on a chain (or arbitrary lattice) so that the $n$ levels $\varepsilon_{k}$ involved are consecutive. Then, the relative length $\frac{n}{N}$ of the $\varepsilon_{k}$ segment contracts to a point in the limit $N\to\infty$ justifying the term “local”. This definition admits $n\sim N^{a}$ with $0\leq a<1$. In all our examples of local operators $a=0$ and we therefore often assume $\frac{n}{N}\sim\frac{1}{N}$ without loss of generality. Let us also emphasize that our notion of locality is unrelated to the property of being local in the real space.
  * [60] P. W. Anderson, Theory of dirty superconductors, J. Phys. Chem. Solids 11, 26 (1959).
  * [61] V. Gurarie and L. Radzihovsky, Resonantly paired fermionic superfluids, Ann. Phys. 322, 2 (2007).
  * [62] E. A. Yuzbashyan and O. Tsyplyatyev, Dynamics of emergent Cooper pairing at finite temperatures, Phys. Rev. B 79, 132504 (2009).
  * [63] J. M. Zhang and R. X. Dong, Exact diagonalization: the Bose-Hubbard model as an example, Eur. J. Phys. 31 591 (2010).
  * [64] R. W. Richardson and N. Sherman, Exact eigenstates of the pairing-force Hamiltonian, Nucl. Phys. 52, 221 (1964).
  * [65] R. W. Richardson, Exact Eigenstates of the Pairing-Force Hamiltonian. II, J. Math. Phys. 6, 1034 (1965).
  * [66] M. Gaudin, Diagonalisation d’une classe d’hamiltoniens de spin, J. Phys. France 37, 1087 (1976).
  * [67] M. Gaudin, The Bethe Wavefunction (Cambridge University Press, 2014).
  * [68] E. K. Sklyanin, Separation of variables in the Gaudin model, J. Sov. Math. 47, 2473 (1989).
  * [69] G. Ortiz, R. Somma, J. Dukelsky, S. Rombouts, Exactly-solvable models derived from a generalized Gaudin algebra, Nucl. Phys. B 707, 421 (2005).
  * [70] E. A. Yuzbashyan, A. A. Baytin, and B. L. Altshuler, Strong-coupling expansion for the pairing Hamiltonian for small superconducting metallic grains, Phys. Rev. B 68, 214509 (2003).
  * [71] B. F. Bayman, A derivation of the pairing-correlation method, Nucl. Phys. 15, 33 (1960).
  * [72] K. Dietrich, H. J. Mang, and J. H. Pradal, Conservation of Particle Number in the Nuclear Pairing Model, Phys. Rev. 135, B22 (1964).
  * [73] J. H. Spencer, L. Florescu, Asymptopia (American Mathematical Society, Providence, Rhode Island, 2014).
  * [74] J. von Delft and D. C. Ralph, Spectroscopy of discrete energy levels in ultrasmall metallic grains, Phys. Rep. 345, 61 (2001).
  * [75] B. Sciolla and G. Biroli, Quantum quenches, dynamical transitions, and off-equilibrium quantum criticality, Phys. Rev. B 88, 201110(R) (2013).
  * [76] P. de Buyl, G. De Ninno, D. Fanelli, C. Nardini, A. Patelli, F. Piazza, and Y. Y. Yamaguchi, Absence of thermalization for systems with long-range interactions coupled to a thermal bath, Phys. Rev. E 87, 042110 (2013).
  * [77] M. P. Müler, E. Adlam, L. Masanes, and N. Wiebe, Thermalization and Canonical Typicality in Translation-Invariant Quantum Lattice Systems, Commun. Math. Phys. 340, 499 (2015).
  * [78] There is also a global phase proportional to $\ln t$ (see note [83] below) and more generally to $\int^{t}dt^{\prime}g(t^{\prime})$, which cancels in $\left|C_{\\{s^{z}\\}}\right|$.
  * [79] A. Polkovnikov, K. Sengupta, A. Silva, M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * [80] E. A. Yuzbashyan, Generalized microcanonical and Gibbs ensembles in classical and quantum integrable dynamics, Ann. Phys. 367, 288 (2016).
  * [81] More rigorously, we verified that the term neglected in making the adiabatic approximation (see p. 330 of Ref. 82) is of higher order in $t$ at small $t$.
  * [82] J. Sakurai and J. Napolitano, Modern Quantum Mechanics, 3rd ed. (Cambridge University Press, 2020).
  * [83] Let us mention for completeness a nonessential overall factor $\exp[-i\eta\frac{N_{\uparrow}}{N}\ln t]$ in $\Psi_{\infty}(N_{\uparrow})$ that we omitted in Eq. (42) for simplicity. We obtain this factor using Eq. (37), which includes the $t^{-1}$ correction, as the saddle point value of $\lambda_{\alpha}$ instead of $\lambda_{\alpha}=\varepsilon_{\alpha}$. An alternative route is the standard time-dependent perturbation theory in the limit $t\to+\infty$ treating $\hat{H}_{\mathrm{int}}$ as a perturbation to $\hat{H}_{0}$. The same term is present in the classical solution. To incorporate it into our answers for the late-time quantum and classical (mean-field) asymptotic states, simply replace $\varphi_{k}\to\varphi_{k}+\frac{\eta}{N}\ln t$ in the expressions for $S_{k}^{\pm}$, $u_{k}$, $v_{k}$, and $\langle\hat{s}_{k}^{-}\rangle_{\mathrm{mf}}$.
  * [84] We showed this for the quantum BCS time evolution with a definite particle number in Sec. IX. Similar arguments apply to the classical BCS dynamics considered in this section.
  * [85] H. Casini and M. Huerta, Entanglement entropy in free quantum field theory, J. Phys. A: Math. Theor. 42 504007
  * [86] E. Witten, A mini-introduction to information theory, Riv. Nuovo Cim. 43, 187 (2020).
  * [87] A. Lerose and S. Pappalardi, Origin of the slow growth of entanglement entropy in long-range interacting spin systems, Phys. Rev. Research 2, 012041(R).
  * [88] A. Georges and G. Kotliar, Hubbard model in infinite dimensions, Phys. Rev. B 45, 6479 (1992).
  * [89] Y. M. Bunkov and G. E. Volovik, Magnon Bose-Einstein condensation and spin superfluidity, J. Phys.: Condens. Matter 22 164210 (2010).
  * [90] M. R. Schweizer, A. J.E. Kreil, G. von Freymann, B. Hillebrands, A. A. Serga, Confinement of Bose-Einstein magnon condensates in adjustable complex magnetization landscapes, arXiv:2208.13507
  * [91] R. M. Gray, Toeplitz and Circulant Matrices: A Review, Found. Trends Commun. Inf. Theory 2(3), 155-239