# Quantum quench in a driven Ising chain

Neil J. Robinson<EMAIL_ADDRESS>Institute for Theoretical
Physics, University of Amsterdam, Postbus 94485, 1090 GL Amsterdam, The
Netherlands  Present address: UKRI EPSRC, Polaris House, North Star Avenue,
Swindon SN2 1ET, United Kingdom  Isaac Pérez Castillo<EMAIL_ADDRESS>Departamento de Física, Universidad Autónoma Metropolitana-Iztapalapa, San
Rafael Atlixco 186, Ciudad de México 09340, Mexico  Edgar Guzmán-González
<EMAIL_ADDRESS>London Mathematical Laboratory, 8 Margravine Gardens,
London W6 8RH, United Kingdom

###### Abstract

We consider the Ising chain driven by oscillatory transverse magnetic fields.
For certain parameter regimes, we reveal a hidden integrable structure in the
problem, which allows access to the exact time-evolution in this driven
quantum system. We compute time-evolved one- and two-point functions following
a quench that activates the driving. It is shown that this model does not heat
up to infinite temperature, despite the absence of energy conservation, and we
further discuss the generalization to a family of driven Hamiltonians that do
not suffer heating to infinite temperature, despite the absence of
integrability and disorder. The particular model studied in detail also
presents a route for realising exotic physics (such as the E8 perturbed
conformal field theory) via driving in quantum chains that could otherwise
never realise such behaviour. In particular we numerically confirm that the
ratio of the meson excitations masses is given by the golden ratio.

Introduction.—Over the last decade, the nonequilibrium dynamics of quantum
systems has attracted a great deal of attention Gogolin and Eisert (2016);
D’Alessio  _et al._ (2016); Calabrese and Cardy (2016); Cazalilla and Chung
(2016); Caux (2016); Essler and Fagotti (2016); Bernard and Doyon (2016);
Ilievski  _et al._ (2016); Vasseur and Moore (2016); De Luca and Mussardo
(2016); Langen  _et al._ (2016); Vidmar and Rigol (2016), motivated by the
desire to address fundamental questions: When and how do quantum systems relax
to equilibrium? How does one describe this equilibrium? What influences the
dynamics and equilibration? Understanding these issues is important when
developing descriptions of a growing number of experiments that examine
nonequilibrium dynamics, both in cold atomic gases Polkovnikov  _et al._
(2011); Langen  _et al._ (2015) and the solid state Bovensiepen and Kirchmann
(2012); Beye  _et al._ (2013). The insights gained may play an important role
in the development of quantum computing resources, especially when considering
how to protect quantum information from the scrambling associated with
thermalization.

Recently attention has turned to understanding driven quantum systems,
partially due to the realization that such systems can host interesting
topological phases (see, e.g., Refs. von Keyserlingk and Sondhi (2016a, b);
Else and Nayak (2016); Potter  _et al._ (2016); Roy and Harper (2016, 2017a,
2017b)) and other exotic behaviors (such as time crystal phases Khemani  _et
al._ (2016); von Keyserlingk  _et al._ (2016); Zhang  _et al._ (2017); Choi
_et al._ (2017); Khemani  _et al._ (2017); Yao  _et al._ (2017); Khemani  _et
al._ (2019)). These studies have generated much discussion of how to extend
and apply the concepts of equilibrium statistical mechanics in the presence of
driving. A particular issue is that, generically, driven quantum systems do
not conserve energy. As a result, in the long time limit entropy maximization
leads them to heat up to infinite temperature, leading to trivial ergodic
behavior. As a result, quantum information is completely scrambled D’Alessio
and Rigol (2014); Lazarides  _et al._ (2014); Ponte  _et al._ (2015a). Routes
to avoid this behavior include introducing disorder to induce a many body
localization transition (see, e.g., Refs. Ponte  _et al._ (2015b); Lazarides
_et al._ (2015); Abanin  _et al._ (2016); Khemani  _et al._ (2016)), or to
consider models that are, in some sense, integrable Gritsev and Polkovnikov
(2017).

In this Letter, we consider a driven model that at each point in time is
nonintegrable but nonetheless possesses the dynamics which is governed by a
hidden integrability. Using this, we compute the nonequilibrium dynamics of
equal-time correlation functions following a quench in which the driving is
initiated. The method for attacking this problem can be generalized to a
(infinite) family of Hamiltonians, opening the door for future
nonperturbative, exact studies. We will see that this whole family of driven
quantum systems, each of which is generically nonintegrable, does not undergo
heating to infinite temperature. We will also see that breaking the special
structure of this family leads to thermalization to infinite temperature.

The driven Ising chain.—We consider a one-dimensional spin-1/2 Ising magnet,
driven by oscillatory transverse fields. The Hamiltonian reads

$\begin{split}H(t)=&-J\sum_{l=1}^{L}\sigma^{z}_{l}\sigma^{z}_{l+1}+h^{z}\sum_{l=1}^{L}\sigma^{z}_{l}\\\
&-g\sum_{l=1}^{L}\Big{(}e^{-i\Omega t}\sigma^{+}_{l}+e^{i\Omega
t}\sigma^{-}_{l}\Big{)},\end{split}$ (1)

with $J>0$ the Ising exchange parameter, $h^{z}$ a static longitudinal field,
$g$ the strength of the transverse fields, which oscillate at frequency
$\Omega$, and $L$ the system size. The spin operators $\sigma^{\alpha}_{l}$
act at the $l$th site of the lattice, $\sigma^{\pm}_{l}=(\sigma^{x}_{l}\pm
i\sigma^{y}_{l})/2$, and we impose periodic boundary conditions
$\sigma^{\alpha}_{L+1}=\sigma^{\alpha}_{1}$. The Hamiltonian (1) is periodic
in time $H(t)=H(t+T)$ with period $T=2\pi/\Omega$, and could be realized in
the quasi-1D ferromagnet CoNb2O6 Coldea  _et al._ (2010); Robinson  _et al._
(2014) by application of oscillating transverse fields.

At a generic time, the Hamiltonian consists of an Ising interaction term and
fields in all ($x,y,z$) directions. Thus instantaneously the Hamiltonian
$H(t)$ is nonintegrable, and the exact computation of quantities seems
unlikely. In the following we will see that this is in fact not the case –
there exists a hidden integrable line within this model where exact results
can be obtained. Furthermore, away from this integrability we will draw
general insights.

Time evolution of observables.—We will now consider how a state
$|\Psi_{0}\rangle$ evolves under the Hamiltonian (1) at times $t>0$. The time-
evolved state $|\Psi(t)\rangle$ will be a solution of the time-dependent
Schrödinger equation

$\Big{[}i\hbar\partial_{t}-H(t)\Big{]}|\Psi(t)\rangle=0,$ (2)

subject to the initial condition $|\Psi(t=0)\rangle=|\Psi_{0}\rangle$. Herein
we set $\hbar=1$, which defines our units. The formal solution of Eq. (2) is
well-known:

$|\Psi(t)\rangle=\mathbb{T}\exp\left(-i\int_{0}^{t}\mathrm{d}t^{\prime}\,H(t^{\prime})\right)|\Psi_{0}\rangle,$
(3)

however, using this to compute time-evolution is a challenge due to the
explicit time-ordering ($\mathbb{T}$) of the exponential. To make some headway
on this problem we apply a time-dependent unitary transformation $U(t)$ 111As
discussed in Refs. Kolodrubetz  _et al._ (2013, 2017); Weinberg  _et al._
(2017), there is nice geometric interpretation of unitary transformations that
depend on a continuous parameter (e.g., time $t$) in terms of gauge
potentials., multiplying both sides of Eq. (2) from the left by $U(t)$ and
inserting a factor of $\mathbb{1}=U(t)^{\dagger}U(t)$ between the wave
function and the operators:

$U(t)\Big{[}i\partial_{t}-H(t)\Big{]}U^{\dagger}(t)U(t)|\Psi(t)\rangle=0.$ (4)

The problem can become much simpler if there is a choice of $U(t)$ such that
this reduces to an effective time-independent Schrödinger equation. Choosing
Takayoshi  _et al._ (2014)

$U(t)=\exp\Big{(}\frac{i\Omega t}{2}\sum_{l}\sigma^{z}_{l}\Big{)}\equiv
e^{\frac{i\Omega t}{2}\sigma^{z}_{\rm tot}},$ (5)

we map Eq. (2) to a time-independent Schrödinger equation,
$(i\partial_{t}-H_{\rm st})|\Phi(t)\rangle=0$, with an effective static
Hamiltonian

$H_{\rm
st}=\sum_{l=1}^{L}\left[-J\sigma^{z}_{l}\sigma^{z}_{l+1}+\left(h^{z}-\frac{\Omega}{2}\right)\sigma^{z}_{l}-g\sigma^{x}_{l}\right].$
(6)

The wave function transforms as $|\Phi(t)\rangle=U(t)|\Psi(t)\rangle$. This
reduction to a static problem is not evident in the Magnus expansion 222See
the Supplemental Material, which also contains the references Blanes  _et al._
(2009); Jordan and Wigner (1928); Fisher and Hartwig (2007); Au-Yang and McCoy
(1974); Basor and Tracy (1991); Forrester and Frankel (2004); Albrecht
Böttcher and Harold Widom (2006); Karlovich (2007); Deift  _et al._ (2011),
for: (i) A discussion of the Magnus Expansion for the driven problem
considered; (ii) details of how two-point correlation functions transform
under action of the time-dependent unitary transformation; (iii) a discussion
of the full-time evolution of observables for a special quench; (iv) detailed
derivations of the required “sudden quench” correlation functions; (v) details
of the numerical algorithm for computing time-evolution of the driven system
and a numerical check of the absence of heating to infinite temperature
outside the integrable line..

Diagonalizing (6) to obtain eigenstates $|E_{n}\rangle$ with energies $E_{n}$,
the time-evolved state can be written as

$|\Psi(t)\rangle=\sum_{n}\exp\Big{[}-i\left(E_{n}+\frac{\Omega}{2}\sigma^{z}_{\text{tot}}\right)t\Big{]}|E_{n}\rangle\langle
E_{n}|\Psi_{0}\rangle\,.$ (7)

The states $|E_{n}\rangle$ are not eigenstates of $\sigma^{z}_{\text{tot}}$
and thus each term in Eq. (7) undergoes nontrivial dynamics. While (7) is
highly nontrivial, there is no need to despair. Our problem reduces to a
tractable one if we focus on equal-time correlation functions, as one can use
that the operator $U(t)$ acts in a simple manner on the spin operators:

$\begin{split}U(t)\sigma^{x,y}_{l}U(t)^{\dagger}&=\cos(\Omega
t)\sigma^{x,y}_{l}\mp\sin(\Omega t)\sigma^{y,x}_{l},\\\
U(t)\sigma^{z}_{l}U(t)^{\dagger}&=\sigma^{z}_{l}.\end{split}$ (8)

Mapping to a “sudden quench”.—Let us now consider the time-evolution of one-
point functions
$s^{\alpha}(t)=\langle\Psi(t)|\sigma^{\alpha}_{l}|\Psi(t)\rangle$, where the
result is independent of $l$ by translational invariance. Using (8) these
become

$\begin{split}s^{z}(t)&=s^{z}_{\text{st}}(t),\\\ s^{x}(t)&=\cos(\Omega
t)s^{x}_{\text{st}}(t)-\sin(\Omega t)s^{y}_{\text{st}}(t),\\\
s^{y}(t)&=\cos(\Omega t)s^{y}_{\text{st}}(t)+\sin(\Omega
t)s^{x}_{\text{st}}(t).\\\ \end{split}$ (9)

Here each time-dependent expectation value on the right hand side describes
time-evolution induced by a sudden quench to the static Hamiltonian (6) when
starting from the initial state $|\Psi_{0}\rangle$:

$s^{\alpha}_{\text{st}}(t)=\sum_{n,m}e^{i(E_{n}-E_{m})t}\langle\Psi_{0}|E_{n}\rangle\langle
E_{n}|\sigma^{\alpha}_{l}|E_{m}\rangle\langle E_{m}|\Psi_{0}\rangle.$ (10)

Equations similar to (9) can be written for the two-point functions,
$s^{\alpha\beta}(\ell;t)=\langle\Psi(t)|\sigma^{\alpha}_{j}\sigma^{\beta}_{j+\ell}|\Psi(t)\rangle$.
These are tractable, but a little unwieldy, so are given in Note (2). All
time-evolved correlation functions are reduced to oscillatory factors
multiplying “sudden quench” correlation functions. Thus for this driven
problem, we can apply the techniques developed for sudden quantum quenches to
compute the time-evolution of observables.

Having reduced the problem from one with driving to an effective sudden
quench, let us return to the static Hamiltonian (6). This describes a quantum
Ising chain with both transverse $g$ and longitudinal $h=h^{z}-\Omega/2$
fields. The two fields can be independently controlled via the amplitude $g$
and frequency $\Omega$ of the driving, see Eq. (1). Two interesting cases are
immediately apparent. Firstly, if the frequency of the driving is tuned to a
$\Omega=2h^{z}$, the longitudinal field is removed from the static
Hamiltonian, which then describes the integrable quantum Ising chain Pfeuty
(1970). Secondly, one can consider tuning both the amplitude and the frequency
such that $g=J$ and $|h^{z}-\Omega/2|\ll g$, where one realizes the lattice
limit of the exotic $E_{8}$ perturbed Ising conformal field theory
Zamolodchikov (1989) (which has recently received renewed attention thanks to
its nonthermal properties Rakovszky  _et al._ (2016); Hódsági  _et al._
(2018); James  _et al._ (2019); Robinson  _et al._ (2019), despite an absence
of integrability). In this work, we will focus on the first scenario and
describe the full time-evolution of one- and two-point functions in this
driven problem. We will touch upon the second case towards the end.

When $\Omega=2h^{z}$, the static Hamiltonian reads:

$H^{0}_{\text{st}}=-J\sum_{l=1}^{L}\sigma^{z}_{l}\sigma^{z}_{l+1}-g\sum_{l=1}^{L}\sigma^{x}_{l}\,.$
(11)

This is the quantum Ising chain, which can be mapped to free fermions and so
is exactly solvable Pfeuty (1970). This reveals that, along the line
$\Omega=2h^{z}$, there is a hidden integrability in the problem (despite,
instantaneously, the Hamiltonian $H(t)$ being nonintegrable). Sudden quenches
in the transverse field Ising model have been extensively studied, with many
exact results being known, see in particular the works of Calabrese, Essler
and Fagotti Calabrese  _et al._ (2011, 2012a, 2012b). We will exploit some of
these results, alongside some new ones, to analytically compute the dynamics
of observables starting from an initial state $|\Psi_{0}\rangle$ that is then
time-evolved with the driven Hamiltonian (1). The derivation of these results
is rather technical, so we provide the details in the Supplemental Material
Note (2).

Time-evolution in the driven model.— Let us now present the time-evolution of
correlation functions in the driven model (1) governed by the effective static
Hamiltonian (11). We compare our analytical results to numerical results
obtained on small finite lattices (our numerical algorithm is explained in
Note (2)).

(a)
---
(b)
Figure 1: (a) Time-evolution of one-point functions starting from the ground
state of the quantum Ising chain, $H(t=0)$ with $J=1,h^{z}=0,g=2$ and time-
evolved with the driven Hamiltonian $H(t)$ (1) with
$J=1,\Omega=1,h^{z}=0.5,g=1.5$ for a system with $L=16$ sites. The behaviour
remains the same for bigger system sizes and larger times. Lines represent
analytical results, while points show numerically exact time-evolution. (b)
Time evolution of the two-point function $s^{zz}(1,t)$ of two adjacent sites
for the same quench, for different system sizes $L$. We see that $s^{zz}(1,t)$
converges to a non-zero value. The revival of the fluctuations is a finite
size effect, as can be seen by increasing the system size.

In Fig. 1 we present the results for one- and two-point functions for a
particular quench. We see that the one-point functions synchronize to the
driving frequency $\Omega$ and no heating to infinite temperature occurs, not
even if we restrict the study of the system to stroboscopic times. For the
two-point functions, we see that $s^{zz}(1,t)$ converges to a non-zero
stationary value, confirming the absence of infinite heating. Although not
included in the figure, we mention that the remaining two-point functions
synchronize to the period $\Omega$ like the one-point functions, except for
$s^{yz}(1,t)$ and $s^{xz}(1,t)$ that converge to zero.

A particularly simple, and solvable in closed form, scenario is realized when
$H_{\rm st}$ coincides with the initial Hamiltonian $H(t<0)$. In this case the
“sudden quench” correlation functions in expressions such as (9) reduce to
equilibrium correlation functions, known since the seminal works of Barouch
_et al._ (1970) and Barouch and McCoy (1971a, b) in the 1970s. Detailed
results in this case are presented in Note (2) and are, to our knowledge, some
of the few closed form exact results known for correlation functions in models
with driving.

Absence of heating to infinite temperature.—With observables mapping in a
simple manner to those from a sudden quench, it is immediately clear that the
system cannot undergo heating to infinite temperature, as is usually assumed
to occur in driven systems D’Alessio and Rigol (2014); Lazarides  _et al._
(2014); Ponte  _et al._ (2015a)). This is easily seen for observables that
feature only $\sigma^{z}_{l}$ operators, which map exactly to “sudden quench”
observables (see, e.g., the first line of Eq. (9)). The long-time limit of
observables after a sudden quench will be described via the relevant
statistical ensemble; for the case detailed above this is the generalized
Gibbs ensemble Rigol  _et al._ (2007); Ilievski  _et al._ (2015); Vidmar and
Rigol (2016). Generically, when $H_{\text{st}}$ is nonintegrable, this will be
a finite-temperature Gibbs ensemble Rigol  _et al._ (2008). It is worth noting
that the absence of heating to infinite temperature is not as result of
integrability, but instead is due to the structure of the driving term. In
Fig. S1 of Note (2), we show an explicit example of a non-integrable system
with absence of heating to infinite temperature by working outside the
integrable line $\Omega=2h^{z}$.

We can then ask, what happens if this structure is broken such that we do not
map to an effective sudden quench problem? We then expect that in the long
time limit the system thermalizes to infinite temperature, due to the absence
of both energy conservation and the mapping to a sudden quench problem,
combined with entropy maximization. We can examine this numerically by adding
terms to our Hamiltonian (1), for example:

$H_{X}(t)=H(t)+J_{X}\sum_{l=1}^{L}\sigma^{x}_{l}\sigma^{x}_{l+1}.$ (12)

The added term breaks $\sigma^{z}$ conservation, and thus evolves non-
trivially under the transformation $U(t)$. This breaks the mapping to a static
Hamiltonian, and hence we expect heating to infinite temperature. It is worth
noting that the thermalization time scale in Floquet systems can be very
large, see e.g. Refs. Abanin  _et al._ (2015, 2017); Else  _et al._ (2017);
Machado  _et al._ (2017). (The Floquet model studied in Ref. Machado  _et al._
(2017) bears some similarity to (12).)

Figure 2: Numerically exact time-evolution of one-point functions with the
driven Hamiltonian $H_{X}(t)$, Eq. (12). This shows that, at the level of one-
point functions, breaking the structure of the drive leads to thermalization
to infinite temperature, as all the expectation values converge to zero. The
parameters considered were the ones of Fig. (1) with $J_{x}=0.5$, and a
Chebyshev expansion of order 64 with a time step $\Delta t=$0.001 (see Note
(2) for the details of the numerical algorithm employed).

In Fig. (2) we present the time-evolution of one-point functions in the driven
model (12). With the addition of the $J_{X}$ term, we see that the system
evolves towards a state with
$\lim_{t\to\infty}\langle\sigma^{\alpha}_{j}(t)\rangle=0$, corresponding to
infinite temperature, at least at the level of one-point subsystems.

Realizing a perturbed critical model.—Let us finish with an illustration of
the second interesting case discussed above. We consider tuning the driving
such that the static Hamiltonian describes the perturbed critical Ising chain:

$H^{1}_{\rm
st}=-J\sum_{l=1}^{L}\sigma^{z}_{l}\sigma^{z}_{l+1}-J\sum_{l=1}^{L}\sigma^{x}_{l}+h\sum_{l=1}^{L}\sigma^{z}_{l}\,.$
(13)

When $h=0$, $H^{1}_{\rm st}$ realizes the critical point of the Ising chain.
For $h\neq 0$, $H^{1}_{\rm st}$ is no longer integrable, but its low-energy
physics is well-understood thanks to Zamolodchikov Zamolodchikov (1989). Pairs
of fermions (corresponding to domain walls in the ordered phase) are confined
by the presence of the longitudinal field $h$ and form “meson” excitations. In
the scaling limit, the algebraic structure of the theory allows the prediction
of these meson masses, including the beautiful result that the ratio of masses
of the first and second meson states realizes the golden ratio,
$m_{2}/m_{1}=\varphi$. With integrability absent, we are limit to performing
small system numerics, such as in Fig. 3.

Figure 3: Dynamical correlation function $s^{xx}(k=0,\omega_{1},\omega_{2})$
(13) for the perturbed critical model with $J=1$ ,$\Omega=1$,
$h_{z}=\Omega/2+0.1$, $g=1.5$ and $L=16$. The data was normalized so that the
maximum value of the plot is one. The four dominant peaks at
$(\pm\Omega,\pm\Omega)$ and $(\mp\Omega,\pm\Omega)$ come from the driving
frequency $\Omega$, while the next three at $\bm{p}_{i}=(m_{i}-\Omega,\Omega-
m_{i})$ come from the masses of the mesons excitations $m_{i}$. Although
$m_{2}/m_{1}$ is not equal to the golden ratio $\varphi$, we verify at the
inset that, as the size of the system is increased, $m_{2}/m_{1}$ gets closer
to $\varphi$.

In Fig. 3 we plot the dynamical correlation function

$\begin{split}s^{xx}(k=0,\omega_{1},\omega_{2})&=\sum_{\ell}\int
dt_{1}dt_{2}e^{i(\omega_{1}t_{1}+\omega_{2}t_{2})}\\\
&\times\langle\Psi_{0}|\sigma^{x}_{j+\ell}(t_{1})\sigma^{x}_{j}(t_{2})|\Psi_{0}\rangle,\end{split}$
(14)

where $\sigma^{x}_{n}(t)$ denotes the time-evolution of $\sigma^{x}_{n}$ in
the Heisenberg picture and, for simplicity, we assume the initial state of the
system was prepared to be the ground state of the static Hamiltonian (13).
Note that, because of the driving, energy is not conserved, and
$\langle\Psi_{0}|\sigma^{x}_{j+\ell}(t_{1})\sigma^{x}_{j}(t_{2})|\Psi_{0}\rangle$
is no longer a function of the time difference $(t_{1}-t_{2})$, therefore we
considered the Fourier transform of both times.

Note that there are four dominant peaks in Fig. 3, these correspond to the
driving frequency $\Omega$ and are located at
$(\omega_{1},\omega_{2})=(\pm\Omega,\pm\Omega),(\mp\Omega,\pm\Omega)$. The
remaining dominant peaks (marked as $\bm{p}_{1}$, $\bm{p}_{2}$ and
$\bm{p}_{3}$ in the figure) correspond to the first excitations or masses of
the static system, and their coordinates are $\bm{p}_{i}=(m_{i}-\Omega,\Omega-
m_{i})$, where $m_{i}$ denotes the masses of the meson excitations. Although
these masses do not satisfy the equality $m_{2}/m_{1}=\varphi$, we verify that
this is a finite size effect in the inset of the figure, as $m_{2}/m_{1}$ gets
closer to $\varphi$ when the size of the system is increased. This reiterates
the fact that driven systems, such as (1), can realise exotic physics that is
inaccessible in its undriven state.

Discussion.—In this Letter, we have explored an example of driven system that
is instantaneously nonintegrable but can nonetheless be solved exactly. This
is due to a hidden integrability in the problem, that is not apparent from the
time-dependent Schrödinger equation: the instantaneous Hamiltonian $H(t)$ is
nonintegrable, but dynamics of observables are nonetheless controlled by an
effective static, integrable Hamiltonian. This may provide a route to
protecting quantum information from the scrambling associated with
thermalization through the addition of driving; this is an interesting
direction for future studies.

The methods applied within this Letter can be used to tackle the dynamics of
an infinite family of Hamiltonians (not necessarily integrable). This family
of driven systems does not undergo heating to infinite temperature, even
though they are absent disorder and (generically) integrability. For example,
consider the Hamiltonian $\widetilde{H}$ of any spin-1/2 chain that conserves
total $\sigma^{z}_{\text{tot}}$ magnetization (this need not be
translationally invariant), which is driven as in (1):

$\widetilde{H}(t)=\widetilde{H}-\widetilde{g}\sum_{l}\left(e^{-i\widetilde{\Omega}t}\sigma^{+}_{l}+\text{H.c.}\right).$
(15)

The transformation (5) still maps Eq. (2) to a time-independent Schrödinger
equation with the new effective static Hamiltonian
$\widetilde{H}_{\text{st}}=\widetilde{H}-(\widetilde{\Omega}/2)\sum_{l}\sigma^{z}_{l}-\widetilde{g}\sum_{l}\sigma^{x}_{l}$.
Time-evolution of observables in a driven system has once again been mapped to
a sudden quench problem. It would be interesting to explore this idea further
in interacting models, such as when $\tilde{H}_{\text{st}}$ describes the
Heisenberg or XXZ model, where potentially integrability can be harnessed to
perform exact calculations.

Another scenario worthy of attention is to consider a problem in which the
parameters of the static Hamiltonian describe a different phase to the initial
Hamiltonian. One may then expect to see signatures of dynamical phase
transitions in the nonequilibrium dynamics, such as kinks in the Lochsmidt
echo Heyl  _et al._ (2013); Karrasch and Schuricht (2013). Further exploring
the lattice limit of Zamolodchikov’s perturbed Ising field theory
Zamolodchikov (1989), which features interesting collective excitations
related to an exotic hidden $E_{8}$ algebraic structure, is interesting. Such
studies would require detailed numerical analysis (perhaps in the scaling
limit James  _et al._ (2018)), an avenue left to future works.

###### Acknowledgements.

Acknowledgments.— Our thanks go to Jean-Sébastien Caux, Mario Collura, Andrew
James, Robert Konik, Tamás Pálmai, and Sergio Tapias Arze for useful
discussions. I.P.C. and E.G.-G. thanks the London Mathematical Laboratory for
financial support. N.J.R. was supported by funding from the EU’s Horizon 2020
research and innovation programme, under grant agreement No 745944, and the
European Research Council under ERC Advanced Grant No 743032 (Dynamint).
E.G.-G. acknowledges the hospitality of the Institute of Physics, University
of Amsterdam, during the completion of this work. N.J.R. thanks the Institute
of Physics of the National Autonomous University of Mexico for hospitality
during a visit, where part of this work was undertaken.

## References

  * Gogolin and Eisert (2016) C. Gogolin and J. Eisert, “Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems,” Rep. Prog. Phys. 79, 056001 (2016).
  * D’Alessio  _et al._ (2016) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, “From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics,” Adv. Phys. 65, 239–362 (2016).
  * Calabrese and Cardy (2016) P. Calabrese and J. Cardy, “Quantum quenches in 1+1 dimensional conformal field theories,” J. Stat. Mech. 2016, P064003 (2016).
  * Cazalilla and Chung (2016) M. A. Cazalilla and M.-C. Chung, “Quantum quenches in the Luttinger model and its close relatives,” J. Stat. Mech. 2016, P064004 (2016).
  * Caux (2016) J.-S. Caux, “The Quench Action,” J. Stat. Mech. 2016, P064006 (2016).
  * Essler and Fagotti (2016) F. H. L. Essler and M. Fagotti, “Quench dynamics and relaxation in isolated integrable quantum spin chains,” J. Stat. Mech. 2016, P064002 (2016).
  * Bernard and Doyon (2016) D. Bernard and B. Doyon, “Conformal field theory out of equilibrium: a review,” J. Stat. Mech. 2016, P064005 (2016).
  * Ilievski  _et al._ (2016) E. Ilievski, M. Medenjak, T. Prosen, and L. Zadnik, “Quasilocal charges in integrable lattice systems,” J. Stat. Mech. 2016, P064008 (2016).
  * Vasseur and Moore (2016) R. Vasseur and J. E. Moore, “Nonequilibrium quantum dynamics and transport: from integrability to many-body localization,” J. Stat. Mech. 2016, P064010 (2016).
  * De Luca and Mussardo (2016) A. De Luca and G. Mussardo, “Equilibration properties of classical integrable field theories,” J. Stat. Mech. 2016, P064011 (2016).
  * Langen  _et al._ (2016) T. Langen, T. Gasenzer, and J. Schmiedmayer, “Prethermalization and universal dynamics in near-integrable quantum systems,” J. Stat. Mech. 2016, P064009 (2016).
  * Vidmar and Rigol (2016) L. Vidmar and M. Rigol, “Generalized Gibbs ensemble in integrable lattice models,” J. Stat. Mech. 2016, P064007 (2016).
  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, “Colloquium: Nonequilibrium dynamics of closed interacting quantum systems,” Rev. Mod. Phys. 83, 863–883 (2011).
  * Langen  _et al._ (2015) T. Langen, R. Geiger, and J. Schmiedmayer, “Ultracold Atoms Out of Equilibrium,” Annu. Rev. Cond. Matter Phys. 6, 201–217 (2015).
  * Bovensiepen and Kirchmann (2012) U. Bovensiepen and P. S. Kirchmann, “Elementary relaxation processes investigated by femtosecond photoelectron spectroscopy of two-dimensional materials,” Laser Photon. Rev. 6, 589–606 (2012).
  * Beye  _et al._ (2013) M. Beye, Ph. Wernet, C. Schüßler-Langeheine, and A. Föhlisch, “Time resolved resonant inelastic X-ray scattering: A supreme tool to understand dynamics in solids and molecules,” J. Electron. Spectrosc. Relat. Phenom. 188, 172 – 182 (2013).
  * von Keyserlingk and Sondhi (2016a) C. W. von Keyserlingk and S. L. Sondhi, “Phase structure of one-dimensional interacting Floquet systems. I. Abelian symmetry-protected topological phases,” Phys. Rev. B 93, 245145 (2016a).
  * von Keyserlingk and Sondhi (2016b) C. W. von Keyserlingk and S. L. Sondhi, “Phase structure of one-dimensional interacting Floquet systems. II. Symmetry-broken phases,” Phys. Rev. B 93, 245146 (2016b).
  * Else and Nayak (2016) D. V. Else and C. Nayak, “Classification of topological phases in periodically driven interacting systems,” Phys. Rev. B 93, 201103 (2016).
  * Potter  _et al._ (2016) A. C. Potter, T. Morimoto, and A. Vishwanath, “Classification of Interacting Topological Floquet Phases in One Dimension,” Phys. Rev. X 6, 041001 (2016).
  * Roy and Harper (2016) R. Roy and F. Harper, “Abelian Floquet symmetry-protected topological phases in one dimension,” Phys. Rev. B 94, 125105 (2016).
  * Roy and Harper (2017a) R. Roy and F. Harper, “Periodic table for Floquet topological insulators,” Phys. Rev. B 96, 155118 (2017a).
  * Roy and Harper (2017b) Rahul Roy and Fenner Harper, “Floquet topological phases with symmetry in all dimensions,” Phys. Rev. B 95, 195128 (2017b).
  * Khemani  _et al._ (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, “Phase Structure of Driven Quantum Systems,” Phys. Rev. Lett. 116, 250401 (2016).
  * von Keyserlingk  _et al._ (2016) C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, “Absolute stability and spatiotemporal long-range order in Floquet systems,” Phys. Rev. B 94, 085112 (2016).
  * Zhang  _et al._ (2017) J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, “Observation of a discrete time crystal,” Nature (London) 543, 217–220 (2017).
  * Choi  _et al._ (2017) S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, “Observation of discrete time-crystalline order in a disordered dipolar many-body system,” Nature (London) 543, 221–225 (2017).
  * Khemani  _et al._ (2017) Vedika Khemani, C. W. von Keyserlingk, and S. L. Sondhi, “Defining time crystals via representation theory,” Phys. Rev. B 96, 115127 (2017).
  * Yao  _et al._ (2017) N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, “Discrete Time Crystals: Rigidity, Criticality, and Realizations,” Phys. Rev. Lett. 118, 030401 (2017).
  * Khemani  _et al._ (2019) Vedika Khemani, Roderich Moessner, and S. L. Sondhi, “A Brief History of Time Crystals,” arXiv e-prints , arXiv:1910.10745 (2019), arXiv:1910.10745 [cond-mat.str-el] .
  * D’Alessio and Rigol (2014) L. D’Alessio and M. Rigol, “Long-time Behavior of Isolated Periodically Driven Interacting Lattice Systems,” Phys. Rev. X 4, 041048 (2014).
  * Lazarides  _et al._ (2014) A. Lazarides, A. Das, and R. Moessner, “Equilibrium states of generic quantum systems subject to periodic driving,” Phys. Rev. E 90, 012110 (2014).
  * Ponte  _et al._ (2015a) P. Ponte, A. Chandran, Z. Papić, and A. A. Abanin, “Periodically driven ergodic and many-body localized quantum systems,” Ann. Phys. (N.Y.) 353, 196 – 204 (2015a).
  * Ponte  _et al._ (2015b) P. Ponte, Z. Papić, F. Huveneers, and D. A. Abanin, “Many-Body Localization in Periodically Driven Systems,” Phys. Rev. Lett. 114, 140401 (2015b).
  * Lazarides  _et al._ (2015) A. Lazarides, A. Das, and R. Moessner, “Fate of Many-Body Localization Under Periodic Driving,” Phys. Rev. Lett. 115, 030402 (2015).
  * Abanin  _et al._ (2016) D. A. Abanin, W. De Roeck, and F. Huveneers, “Theory of many-body localization in periodically driven systems,” Ann. Phys. (N.Y.) 372, 1–11 (2016).
  * Gritsev and Polkovnikov (2017) V. Gritsev and A. Polkovnikov, “Integrable Floquet dynamics,” SciPost Phys. 2, 021 (2017).
  * Coldea  _et al._ (2010) R. Coldea, D. A. Tennant, E. M. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl, and K. Kiefer, “Quantum Criticality in an Ising Chain: Experimental Evidence for Emergent E8 Symmetry,” Science 327, 177 (2010).
  * Robinson  _et al._ (2014) Neil J. Robinson, Fabian H. L. Essler, Ivelisse Cabrera, and Radu Coldea, “Quasiparticle breakdown in the quasi-one-dimensional Ising ferromagnet ${\mathrm{CoNb}}_{2}{\mathrm{O}}_{6}$,” Phys. Rev. B 90, 174406 (2014).
  * Note (1) As discussed in Refs. Kolodrubetz  _et al._ (2013, 2017); Weinberg  _et al._ (2017), there is nice geometric interpretation of unitary transformations that depend on a continuous parameter (e.g., time $t$) in terms of gauge potentials.
  * Takayoshi  _et al._ (2014) S. Takayoshi, H. Aoki, and T. Oka, “Magnetization and phase transition induced by circularly polarized laser in quantum magnets,” Phys. Rev. B 90, 085150 (2014).
  * Note (2) See the Supplemental Material, which also contains the references Blanes  _et al._ (2009); Jordan and Wigner (1928); Fisher and Hartwig (2007); Au-Yang and McCoy (1974); Basor and Tracy (1991); Forrester and Frankel (2004); Albrecht Böttcher and Harold Widom (2006); Karlovich (2007); Deift  _et al._ (2011), for: (i) A discussion of the Magnus Expansion for the driven problem considered; (ii) details of how two-point correlation functions transform under action of the time-dependent unitary transformation; (iii) a discussion of the full-time evolution of observables for a special quench; (iv) detailed derivations of the required “sudden quench” correlation functions; (v) details of the numerical algorithm for computing time-evolution of the driven system and a numerical check of the absence of heating to infinite temperature outside the integrable line.
  * Pfeuty (1970) Pierre Pfeuty, “The one-dimensional Ising model with a transverse field,” Ann. Phys. (N.Y.) 57, 79 – 90 (1970).
  * Zamolodchikov (1989) A. B. Zamolodchikov, “Integrals of motion and S-matrix of the (scaled) T = Tc Ising model with magnetic field,” Int. J. Mod. Phys. A 04, 4235–4248 (1989).
  * Rakovszky  _et al._ (2016) T. Rakovszky, M. Mestyán, M. Collura, M. Kormos, and G. Takács, “Hamiltonian truncation approach to quenches in the Ising field theory,” Nucl. Phys. B 911, 805 – 845 (2016).
  * Hódsági  _et al._ (2018) Kristóf Hódsági, Márton Kormos, and Gábor Takács, “Quench dynamics of the Ising field theory in a magnetic field,” SciPost Phys. 5, 27 (2018).
  * James  _et al._ (2019) Andrew J. A. James, Robert M. Konik, and Neil J. Robinson, “Nonthermal States Arising from Confinement in One and Two Dimensions,” Phys. Rev. Lett. 122, 130603 (2019).
  * Robinson  _et al._ (2019) Neil J. Robinson, Andrew J. A. James, and Robert M. Konik, “Signatures of rare states and thermalization in a theory with confinement,” Phys. Rev. B 99, 195108 (2019).
  * Calabrese  _et al._ (2011) P. Calabrese, F. H. L. Essler, and M. Fagotti, “Quantum Quench in the Transverse-Field Ising Chain,” Phys. Rev. Lett. 106, 227203 (2011).
  * Calabrese  _et al._ (2012a) P. Calabrese, F. H. L. Essler, and M. Fagotti, “Quantum quench in the transverse field Ising chain: I. Time evolution of order parameter correlators,” J. Stat. Mech. 2012, P07016 (2012a).
  * Calabrese  _et al._ (2012b) P. Calabrese, F. H. L. Essler, and M. Fagotti, “Quantum quenches in the transverse field Ising chain: II. Stationary state properties,” J. Stat. Mech. 2012, P07022 (2012b).
  * Barouch  _et al._ (1970) E. Barouch, B. M. McCoy, and M. Dresden, “Statistical Mechanics of the XY Model. I,” Phys. Rev. A 2, 1075–1092 (1970).
  * Barouch and McCoy (1971a) E. Barouch and B. M. McCoy, “Statistical Mechanics of the XY Model. II. Spin-Correlation Functions,” Phys. Rev. A 3, 786–804 (1971a).
  * Barouch and McCoy (1971b) E. Barouch and B. M. McCoy, “Statistical Mechanics of the XY Model. III,” Phys. Rev. A 3, 2137–2140 (1971b).
  * Rigol  _et al._ (2007) Marcos Rigol, Vanja Dunjko, Vladimir Yurovsky, and Maxim 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).
  * Ilievski  _et al._ (2015) E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler, and T. Prosen, “Complete Generalized Gibbs Ensembles in an Interacting Theory,” Phys. Rev. Lett. 115, 157201 (2015).
  * Rigol  _et al._ (2008) Marcos Rigol, Vanja Dunjko, and Maxim Olshanii, “Thermalization and its mechanism for generic isolated quantum systems,” Nature (London) 452, 854 (2008).
  * Abanin  _et al._ (2015) Dmitry A. Abanin, Wojciech De Roeck, and Fran çois Huveneers, “Exponentially Slow Heating in Periodically Driven Many-Body Systems,” Phys. Rev. Lett. 115, 256803 (2015).
  * Abanin  _et al._ (2017) Dmitry Abanin, Wojciech De Roeck, Wen Wei Ho, and François Huveneers, “A Rigorous Theory of Many-Body Prethermalization for Periodically Driven and Closed Quantum Systems,” Commun. Math. Phys. 354, 809–827 (2017).
  * Else  _et al._ (2017) Dominic V. Else, Bela Bauer, and Chetan Nayak, “Prethermal Phases of Matter Protected by Time-Translation Symmetry,” Phys. Rev. X 7, 011026 (2017).
  * Machado  _et al._ (2017) Francisco Machado, Gregory D. Meyer, Dominic V. Else, Chetan Nayak, and Norman Y. Yao, “Exponentially Slow Heating in Short and Long-range Interacting Floquet Systems,” arXiv e-prints , arXiv:1708.01620 (2017), arXiv:1708.01620 [quant-ph] .
  * Heyl  _et al._ (2013) M. Heyl, A. Polkovnikov, and S. Kehrein, “Dynamical Quantum Phase Transitions in the Transverse-Field Ising Model,” Phys. Rev. Lett. 110, 135704 (2013).
  * Karrasch and Schuricht (2013) C. Karrasch and D. Schuricht, “Dynamical phase transitions after quenches in nonintegrable models,” Phys. Rev. B 87, 195104 (2013).
  * James  _et al._ (2018) Andrew J. A. James, Robert M. Konik, Philippe Lecheminant, Neil J. Robinson, and Alexei M. Tsvelik, “Non-perturbative methodologies for low-dimensional strongly-correlated systems: From non-Abelian bosonization to truncated spectrum methods,” Rep. Prog. Phys. 81, 046002 (2018).
  * Kolodrubetz  _et al._ (2013) M. Kolodrubetz, V. Gritsev, and A. Polkovnikov, “Classifying and measuring geometry of a quantum ground state manifold,” Phys. Rev. B 88, 064304 (2013).
  * Kolodrubetz  _et al._ (2017) M. Kolodrubetz, D. Sels, P. Mehta, and A. Polkovnikov, “Geometry and non-adiabatic response in quantum and classical systems,” Phys. Rep. 697, 1–87 (2017), arXiv:1602.01062 [cond-mat.quant-gas] .
  * Weinberg  _et al._ (2017) P. Weinberg, M. Bukov, L. D’Alessio, A. Polkovnikov, S. Vajna, and M. Kolodrubetz, “Adiabatic perturbation theory and geometry of periodically-driven systems,” Phys. Rep. 688, 1–35 (2017).
  * Blanes  _et al._ (2009) S. Blanes, F. Casas, J.A. Oteo, and J. Ros, “The Magnus expansion and some of its applications,” Phys. Rep. 470, 151 – 238 (2009).
  * Jordan and Wigner (1928) P. Jordan and E. Wigner, “Über das Paulische Äquivalenzverbot,” Z. Phys. 47, 631–651 (1928).
  * Fisher and Hartwig (2007) Michael E. Fisher and Robert E. Hartwig, “Toeplitz Determinants: Some Applications, Theorems, and Conjectures,” in  _Advances in Chemical Physics_ (John Wiley & Sons, Ltd, 2007) pp. 333–353.
  * Au-Yang and McCoy (1974) Helen Au-Yang and Barry McCoy, “Theory of layered Ising models. II. Spin correlation functions parallel to the layering,” Phys. Rev. B 10, 3885–3905 (1974).
  * Basor and Tracy (1991) Estelle L. Basor and Craig A. Tracy, “The Fisher-Hartwig conjecture and generalizations,” Physica A 177, 167 – 173 (1991).
  * Forrester and Frankel (2004) P. J. Forrester and N. E. Frankel, “Applications and generalizations of Fisher–Hartwig asymptotics,” J. Math. Phys. 45, 2003–2028 (2004).
  * Albrecht Böttcher and Harold Widom (2006) Albrecht Böttcher and Harold Widom, “Szegö via Jacobi,” Lin. Alg. Appl. 419, 656 – 667 (2006).
  * Karlovich (2007) Alexei Yu. Karlovich, “Asymptotics of block Toeplitz determinants generated by factorable matrix functions with equal partial indices,” Math. Nachr. 280, 1118–1127 (2007).
  * Deift  _et al._ (2011) P. Deift, A. Its, and I. Krasovsky, “Asymptotics of Toeplitz, Hankel, and Toeplitz+Hankel determinants with Fisher-Hartwig singularities,” Ann. Math. 174, 1243–1299 (2011).

Supplemental Material for: “Quantum quench in a driven Ising chain”
Neil J. Robinson1, Isaac Pérez Castillo,2 and Edgar Guzmán-González3

1Institute for Theoretical Physics, University of Amsterdam,
Postbus 94485, 1090 GL Amsterdam, The Netherlands
2Departamento de Física, Universidad Autónoma Metropolitana-Iztapalapa,
San Rafael Atlixco 186, Ciudad de México 09340, Mexico
3London Mathematical Laboratory, 18 Margravine Gardens, London W6 8RH, United
Kingdom

In this Supplemental Material, we present:

  1. 1.

The Magnus expansion for the problem in the main text.

  2. 2.

The exact mapping between time-dependent two-point functions of the driven
systems and a sudden quench problem.

  3. 3.

Exact full-time evolution of one- and two-point functions for a special case
of the driven system quench.

  4. 4.

Details of the computation of “sudden quench” correlation functions,
summarizing both previously known results and presenting new ones.

  5. 5.

Details of our numerical algorithm for simulating time-evolution in the driven
system.

  6. 6.

Details of the procedure used to compare finite-size numerical results with
analytical expressions.

##  S1 Magnus Expansion

An alternative approach to compute the time-evolution of the driven system is
given by the Magnus expansion. For a time-dependent Hamiltonian $H(t)$, the
time-evolution operator ${\cal U}(t)$ is defined via the differential equation

$i\frac{{\rm d}}{{\rm d}t}{\cal U}(t)=H(t){\cal U}(t),\qquad{\cal
U}(t=0)=\mathbb{1}.$ (S1)

This is formally solved in terms of the well-known time-ordered exponential

${\cal U}(t)=\mathbb{T}\,\exp\left[-i\int_{0}^{t}{\rm
d}t^{\prime}H(t^{\prime})\right].$ (S2)

However, the time-ordering here (denoted by the operator $\mathbb{T}$) is
difficult to treat and so it may be beneficial to take an alternative
approach. One such approximate approach is the Magnus expansion [68], where
the time-evolution operator is instead approximated by

${\cal
U}(t)=\exp\left[\Xi(t)\right],\qquad\Xi(t)=\sum_{k=1}^{\infty}\Xi_{k}(t),$
(S3)

where terms in the series are obtained from the differential equation

$i\frac{{\rm d}}{{\rm
d}t}\Xi(t)=\frac{\text{ad}_{\Xi}}{\exp(\text{ad}_{\Xi})-1}\,H(t).$ (S4)

Here $\text{ad}_{\Xi}$ is the adjoint action of $\Xi$:
$\text{ad}_{\Xi}(O)=[\Xi,O]$. The terms of the series $\Xi_{k}(t)$ have well-
known forms

$\begin{split}\Xi_{1}(t)=&-i\int_{0}^{t}{\rm
d}t_{1}H(t_{1}),\qquad\Xi_{2}(t)=-\frac{1}{2}\int_{0}^{t}{\rm
d}t_{1}\int_{0}^{t_{1}}{\rm d}t_{2}\left[H(t_{1}),H(t_{2})\right],\\\
\Xi_{3}(t)=&\frac{i}{6}\int_{0}^{t}{\rm d}t_{1}\int_{0}^{t_{1}}{\rm
d}t_{2}\int_{0}^{t_{2}}{\rm
d}t_{3}\Big{(}\big{[}H(t_{1}),[H(t_{2}),H(t_{3})]\big{]}+\big{[}H(t_{3}),[H(t_{2}),H(t_{1})]\big{]}\Big{)}.\end{split}$
(S5)

Here the first term simply corresponds to taking the time-averaged
Hamiltonian, as might well be expected.

Let us now turn to the problem at hand in the main body of the manuscript, and
examine the Magnus expansion for the time-evolution operator, working at small
$T$, i.e. high frequency for the driving. We will see that the integrability
that we harnessed in the main text to perform exact calculations is completely
hidden in the Magnus expansion. Indeed, the presence of a simple static limit
(independent of integrability of that static Hamiltonian) is not at all
apparent. To illustrate this, we consider the time-evolution operator over a
single period $T=2\pi/\Omega$ of the driving, ${\cal U}(T)$. Treating the
Magnus expansion as a power series approximation in $T$, the first order term
in the expansion is just the average of the Hamiltonian over a single period,
so the driving vanishes:

$i\Xi_{1}(T)=-iT\left(-J\sum_{l=1}^{L}\sigma^{z}_{l}\sigma^{z}_{l+1}+h^{z}\sum_{l=1}^{L}\sigma^{z}_{l}\right).$
(S6)

For the second order term in the expansion, we require the commutator

$\begin{split}[H(t_{1}),H(t_{2})]=&\
2ig\sum_{l=1}^{L}\Bigg{\\{}g\sin\Big{(}\Omega(t_{2}-t_{1})\Big{)}\sigma^{z}_{l}\\\
&-h^{z}\Big{[}(c_{2}-c_{1})\sigma^{y}_{l}-(s_{2}-s_{1})\sigma^{x}_{l}\Big{]}\\\
&+J(c_{2}-c_{1})\Big{(}\sigma^{z}_{l}\sigma^{y}_{l+1}+\sigma^{y}_{l}\sigma^{z}_{l+1}\Big{)}\\\
&-J(s_{2}-s_{1})\Big{(}\sigma^{z}_{l}\sigma^{x}_{l+1}+\sigma^{x}_{l}\sigma^{z}_{l+1}\Big{)}\Bigg{\\}}.\end{split}$

where for conciseness we define the short hand notations $c_{a}=\cos(\Omega
t_{a})$, $s_{a}=\sin(\Omega t_{a})$. Performing the double integral, we obtain
$\Xi_{2}(T)$:

$i\Xi_{2}(T)=\frac{gT^{2}}{\pi}\sum_{l=1}^{L}\bigg{[}h^{z}\sigma^{x}_{l}-\frac{g}{2}\sigma^{z}_{l}-J\Big{(}\sigma^{z}_{l}\sigma^{x}_{l+1}+\sigma^{x}_{l}\sigma^{z}_{l+1}\Big{)}\bigg{]}.$
(S7)

In this second order term (in $T$), we see already that the Magnus expansion
for the time-evolution over one period is rather messy. Continuing to third
order, we expect to generate terms of the form
${\sigma^{x}_{l}\sigma^{y}_{l+1},\,\sigma^{y}_{l}\sigma^{z}_{l+1}}$ too, so
higher orders are unlikely to simplify things further. We see that the
existence of a simple static Hamiltonian, as given in the main body of the
manuscript, is not obvious from the Magnus expansion.

##  S2 Mapping between driven and “sudden quench” two-point functions

Let us consider equal-time two-point functions in the nonequilibrium time-
evolved wave function [see Eq. (7) of the main text],

$s^{\alpha\beta}(\ell;t)=\langle\Psi(t)|\sigma^{\alpha}_{j}\sigma^{\beta}_{j+\ell}|\Psi(t)\rangle.$
(S8)

Using that $U(t)$ is an element of SU(2) leads to the transformations
described in Eq. (8) of the main text, and we arrive at the following mapping
between driven and “sudden quench” correlation functions:

$\begin{split}s^{zz}(\ell;t)=&\ s^{zz}_{\text{st}}(\ell;t),\\\
s^{xx}(\ell;t)=&\ \cos^{2}(\Omega t)s^{xx}_{\text{st}}(\ell;t)+\sin^{2}(\Omega
t)s^{yy}_{\text{st}}(\ell;t)\\\ &-\frac{1}{2}\sin(2\Omega
t)\big{[}s^{xy}_{\text{st}}(\ell;t)+s^{yx}_{\text{st}}(\ell;t)\big{]},\\\
s^{yy}(\ell;t)=&\ \cos^{2}(\Omega t)s^{yy}_{\text{st}}(\ell;t)+\sin^{2}(\Omega
t)s^{xx}_{\text{st}}(\ell;t)\\\ &+\frac{1}{2}\sin(2\Omega
t)\big{[}s^{xy}_{\text{st}}(\ell;t)+s^{yx}_{\text{st}}(\ell;t)\big{]},\\\
s^{xy}(\ell;t)=&\ \cos^{2}(\Omega t)s^{xy}_{\text{st}}(\ell;t)-\sin^{2}(\Omega
t)s^{yx}_{\text{st}}(\ell;t)\\\ &+\frac{1}{2}\sin(2\Omega
t)\big{[}s^{xx}_{\text{st}}(\ell;t)-s^{yy}_{\text{st}}(\ell;t)\big{]},\\\
s^{xz}(\ell;t)=&\ \cos(\Omega t)s^{xz}_{\text{st}}(\ell;t)-\sin(\Omega
t)s^{yz}_{\text{st}}(\ell;t),\\\ s^{yz}(\ell;t)=&\ \cos(\Omega
t)s^{yz}_{\text{st}}(\ell;t)+\sin(\Omega
t)s^{xz}_{\text{st}}(\ell;t).\end{split}$ (S9)

Here, as in the main body of the manuscript, we define the “sudden quench”
correlation functions

$\displaystyle s^{\alpha\beta}_{\text{st}}(\ell;t)=\sum_{n,m}$ $\displaystyle
e^{i(E_{n}-E_{m})t}\langle
E_{n}|\sigma^{\alpha}_{j}\sigma^{\beta}_{j+\ell}|E_{m}\rangle$
$\displaystyle\times\langle\Psi_{0}|E_{n}\rangle\langle
E_{m}|\Psi_{0}\rangle,$ (S10)

where $|E_{n}\rangle$ are eigenstates of the static Hamiltonian
$H_{\text{st}}$ with energy $E_{n}$ and $|\Psi_{0}\rangle=|\Psi(t=0)\rangle$
is the initial state for the time-evolution. We note that Eq. (S10) is the
same as time-evolving the initial state with $H_{\text{st}}$ alone.

##  S3 Full time-evolution of a special case

In this section, we consider a slightly surprising scenario from our mapping.
We consider starting in the ground state of the Hamiltonian

$H_{\text{init}}=-J\sum_{l=1}^{L}\sigma^{z}_{l}\sigma^{z}_{l+1}-g\sum_{l}\sigma^{x}_{l}.$
(S11)

Let us work with $J>0$ and $g>0$. We perform a quench where we start to rotate
the transverse field and apply an additional longitudinal field

$\displaystyle H(t>0)=$
$\displaystyle-J\sum_{l=1}^{L}\sigma^{z}_{l}\sigma^{z}_{l+1}+h^{z}\sum_{l=1}^{N}\sigma^{z}_{l}$
$\displaystyle-\sum_{l=1}^{L}g\Big{(}e^{-i\Omega t}\sigma^{+}_{l}+e^{i\Omega
t}\sigma_{l}^{-}\Big{)}.$ (S12)

with frequency $\Omega$. We choose $h^{z}=\Omega/2$, in which case the
effective static Hamiltonian is

$H_{\text{st}}=-J\sum_{l=1}^{L}\sigma^{z}_{l}\sigma^{z}_{l+1}-g\sum_{l=1}^{L}\sigma^{x}_{l}.$
(S13)

That is, the effective state Hamiltonian is identical to the initial
Hamiltonian! In this case, the “sudden quench” correlation functions reduce to
the equilibrium ones. These can be computed with standard free fermion
techniques, see Ref. [43], and read:

$\displaystyle s^{z}_{\text{st}}(\ell;t)$
$\displaystyle=\left(1-\lambda^{-2}\right)^{1/8}\Theta(\lambda-1),$ (S14)
$\displaystyle s^{x}_{\text{st}}(\ell;t)$ $\displaystyle=G(0),\quad
s^{y}_{\text{st}}(\ell;t)=0,$ (S15)

for the one-point functions, while the diagonal two-point functions read:

$\begin{split}s^{zz}_{\text{st}}(\ell;t)&=\text{det}\left[\begin{array}[]{cccc}G(-1)&G(-2)&\ldots&G(-\ell)\\\
G(0)&G(-1)&\ldots&G(-\ell+1)\\\ \vdots&\vdots&&\vdots\\\
G(\ell-2)&G(\ell-1)&\ldots&G(-1)\end{array}\right]\\\
s^{xx}_{\text{st}}(\ell;t)&=-G(\ell)G(-\ell)+G(0)^{2},\\\
s^{yy}_{\text{st}}(\ell;t)&=\text{det}\left[\begin{array}[]{cccc}G(1)&G(0)&\ldots&G(2-\ell)\\\
G(2)&G(1)&\ldots&G(3-\ell)\\\ \vdots&\vdots&&\vdots\\\
G(\ell)&G(\ell-1)&\ldots&G(1)\end{array}\right]\end{split}$ (S16)

where we use notations, $\lambda=J/g$, $\Theta(x\leq 0)=0$, $\Theta(x>0)=1$,
and we define the functions

$G(n)=L(n)+\lambda L(n+1),$ (S17)

with

$L(n)=\frac{1}{\pi}\int_{0}^{\pi}{\rm
d}k\,\frac{\cos(kn)}{\sqrt{1+\lambda^{2}+2\lambda\cos(k)}}.$ (S18)

The off-diagonal (in spin indices) two-point functions require some further
work, not usually being considered in equilibrium treatments of the quantum
Ising chain. By symmetry arguments
$s^{\alpha\beta}(\ell,t)=s^{\beta\alpha}(\ell,t)$. Using the Heisenberg
equation of motion,
$s^{zy}_{\text{st}}(\ell,t)\propto\dot{s}^{zz}_{\text{st}}(\ell,t)=0$, as we
are working in an static case. The terms $s^{xz}_{\text{st}}(\ell,t)$ and
$s^{xy}_{\text{st}}(\ell,t)$ can not be computed in general using the standard
techniques, as they mix the even sector with the odd sector in the fermionic
theory (see the next section for more details). However, in the paramagnetic
phase, where $\lambda<1$, using the symmetry of the Hamiltonian under $\pi$
rotations of the quantization axis about the $x$-direction, after some tedious
algebra, we find

$\displaystyle s^{zy}_{\text{st}}(\ell;t)$ $\displaystyle=0,$ (S19)
$\displaystyle\left.s^{xz}_{\text{st}}(\ell;t)\right|_{\lambda<1}$
$\displaystyle=0,$ (S20)
$\displaystyle\left.s^{xy}_{\text{st}}(\ell;t)\right|_{\lambda<1}$
$\displaystyle=0,$ (S21)

Putting these known equilibrium results together with the mapping between
driven correlation functions and “sudden quench” ones, we obtain the full
time-dependent correlation functions for the $H_{\text{init}}\to H(t>0)$
quench

$\begin{split}s^{z}(\ell;t)&=s^{z}_{\text{st}}(\ell;t)\\\
s^{x}(\ell;t)&=\cos(\Omega t)s^{x}_{\text{st}}(\ell;t),\\\
s^{y}(\ell;t)&=\sin(\Omega t)s^{x}_{\text{st}}(\ell;t),\\\
s^{zz}(\ell;t)&=s^{zz}_{\text{st}}(\ell;t),\\\
s^{xx}(\ell;t)|_{\lambda<1}&=\cos^{2}(\Omega
t)s^{xx}_{\text{st}}(\ell;t)+\sin^{2}(\Omega t)s^{yy}_{\text{st}}(\ell;t),\\\
s^{yy}(\ell;t)|_{\lambda<1}&=\cos^{2}(\Omega
t)s^{yy}_{\text{st}}(\ell;t)+\sin^{2}(\Omega t)s^{xx}_{\text{st}}(\ell;t),\\\
s^{xy}(\ell;t)|_{\lambda<1}&=\sin(\Omega t)\cos(\Omega
t)\Big{[}s^{xx}_{\text{st}}(\ell;t)-s^{yy}_{\text{st}}(\ell;t)\Big{]},\\\
s^{xz}(\ell;t)|_{\lambda<1}&=0\\\
s^{yz}(\ell;t)|_{\lambda<1}&=0\,.\end{split}$ (S22)

Furthermore, in the limit $\ell\gg 1$ one can simplify some of these results
using an asymptotic analysis via Szegö’s theorem, which leads to (some) closed
form analytical results [52–54]

$\begin{split}s^{zz}_{\text{st}}(\ell\gg
1;t)\bigg{|}_{\lambda<1}\approx&\frac{1}{\sqrt{\pi}}\left(1-\lambda^{2}\right)^{-1/4}\frac{\lambda^{\ell}}{\sqrt{\ell}}\left[1-\frac{1}{8\ell}\frac{1+\lambda^{2}}{1-\lambda^{2}}+O\left(\frac{1}{\ell^{2}}\right)\right],\\\
s^{zz}_{\text{st}}(\ell\gg
1;t)\bigg{|}_{\lambda>1}\approx&\left(1-\lambda^{-2}\right)^{1/4}+O\left(\frac{1}{\ell}\right),\\\
s^{xx}_{\text{st}}(\ell\gg
1;t)\bigg{|}_{\lambda<1}\approx&\left(s^{x}_{\text{st}}(\ell;t)\right)^{2}-\frac{\lambda^{2\ell}}{2\ell^{2}\pi}\left(1-\frac{1}{\ell}\frac{3\lambda^{-2}-1}{\lambda^{-2}+1}\right)\\\
s^{xx}_{\text{st}}(\ell\gg
1;t)\bigg{|}_{\lambda>1}\approx&\left(s^{x}_{\text{st}}(\ell;t)\right)^{2}-\frac{1}{2\lambda^{2\ell+2}\pi}\left(1-\frac{1}{\ell}\frac{\lambda^{-2}-3}{\lambda^{-2}-1}\right)\\\
s^{yy}_{\text{st}}(\ell\gg
1;t)\bigg{|}_{\lambda<1}\approx&-\frac{1}{2\sqrt{\pi}}(1-\lambda^{2})^{3/4}\frac{\lambda^{\ell}}{\ell^{3/2}}\left[1+\frac{3}{2\ell}\left(\frac{1}{2(1-\lambda^{2})}-\frac{1}{4}\right)+O\left(\frac{1}{\ell^{2}}\right)\right]\\\
s^{yy}_{\text{st}}(\ell\gg
1;t)\bigg{|}_{\lambda>1}\approx&-\frac{1}{\pi}\left(1-\lambda^{-2}\right)^{-3/4}\frac{\lambda^{-2\ell}}{\ell^{3}}\left[\frac{1}{2}-\frac{3}{8\ell}\frac{1+\lambda^{-2}}{1-\lambda^{-2}}+O\left(\frac{1}{\ell^{2}}\right)\right].\\\
\end{split}$ (S23)

Putting together these results, we obtain exact, closed form results for one-
point functions and the asymptotics of two-point correlation functions in a
nonequilibrium driven problem. The non-equilibrium one-point functions read:

$\displaystyle s^{z}(t)$
$\displaystyle=\left(1-\lambda^{-2}\right)^{1/8}\Theta(\lambda-1),$ (S24)
$\displaystyle s^{x}(t)\Big{|}_{\lambda\neq 1}$
$\displaystyle=\frac{\cos(\Omega
t)}{\pi|1-\lambda|}\Bigg{[}\left(\lambda-1\right)^{2}{\cal
E}\left(-\frac{4\lambda}{(1-\lambda)^{2}}\right)-\left(\lambda^{2}-1\right){\cal
K}\left(-\frac{4\lambda}{(1-\lambda)^{2}}\right)\Bigg{]},$ (S25)
$\displaystyle s^{y}(t)\Big{|}_{\lambda\neq 1}$
$\displaystyle=\frac{\sin(\Omega
t)}{\pi|1-\lambda|}\Bigg{[}\left(\lambda-1\right)^{2}{\cal
E}\left(-\frac{4\lambda}{(1-\lambda)^{2}}\right)-\left(\lambda^{2}-1\right){\cal
K}\left(-\frac{4\lambda}{(1-\lambda)^{2}}\right)\Bigg{]}$ (S26)

where ${\cal
E}(m)=\int_{0}^{\pi/2}\mathrm{d}\theta\,[1-m\sin^{2}(\theta)]^{1/2}$ and
${\cal K}(m)=\int_{0}^{\pi/2}\mathrm{d}\theta\,[1-m\sin^{2}(\theta)]^{-1/2}$
are elliptic functions. The asymptotics of non-equilibrium two-point functions
read:

$\displaystyle s^{zz}(\ell\gg 1;t)\bigg{|}_{\lambda<1}=$
$\displaystyle~{}\frac{1}{\sqrt{\pi}}\left(1-\lambda^{2}\right)^{-1/4}\frac{\lambda^{\ell}}{\sqrt{\ell}}\left[1-\frac{1}{8\ell}\frac{1+\lambda^{2}}{1-\lambda^{2}}+O\left(\frac{1}{\ell^{2}}\right)\right],$
(S27) $\displaystyle s^{zz}(\ell\gg 1;t)\bigg{|}_{\lambda>1}=$
$\displaystyle~{}\left(1-\lambda^{-2}\right)^{1/4}+O\left(\frac{1}{\ell}\right),$
(S28) $\displaystyle s^{xx}(\ell\gg 1;t)\bigg{|}_{\lambda<1}=$
$\displaystyle~{}\cos^{2}(\Omega
t)\Bigg{[}\left(s^{x}_{\text{st}}(\ell;t)\right)^{2}-\frac{\lambda^{2\ell}}{2\ell^{2}\pi}\left(1-\frac{1}{\ell}\frac{3-\lambda^{2}}{1+\lambda^{2}}\right)\Bigg{]}$
$\displaystyle~{}-\sin^{2}(\Omega
t)\Bigg{[}\frac{1}{2\sqrt{\pi}}(1-\lambda^{2})^{3/4}\frac{\lambda^{\ell}}{\ell^{3/2}}\left[1+\frac{3}{2\ell}\left(\frac{1}{2(1-\lambda^{2})}-\frac{1}{4}\right)+O\left(\frac{1}{\ell^{2}}\right)\right]$
(S29) $\displaystyle s^{yy}(\ell;t)\bigg{|}_{\lambda<1}=$
$\displaystyle~{}\sin^{2}(\Omega
t)\Bigg{[}\left(s^{x}_{\text{st}}(\ell;t)\right)^{2}-\frac{\lambda^{2\ell}}{2\ell^{2}\pi}\left(1-\frac{1}{\ell}\frac{3-\lambda^{2}}{1+\lambda^{2}}\right)\Bigg{]}$
$\displaystyle~{}-\cos^{2}(\Omega
t)\Bigg{[}\frac{1}{2\sqrt{\pi}}(1-\lambda^{2})^{3/4}\frac{\lambda^{\ell}}{\ell^{3/2}}\left[1+\frac{3}{2\ell}\left(\frac{1}{2(1-\lambda^{2})}-\frac{1}{4}\right)+O\left(\frac{1}{\ell^{2}}\right)\right]$
(S30) $\displaystyle s^{xy}(\ell;t)\bigg{|}_{\lambda<1}=$
$\displaystyle~{}\sin(\Omega t)\cos(\Omega
t)\Bigg{[}\left(s^{x}_{\text{st}}(\ell;t)\right)^{2}-\frac{\lambda^{2\ell}}{2\ell^{2}\pi}\left(1-\frac{1}{\ell}\frac{3-\lambda^{2}}{1+\lambda^{2}}\right)$
$\displaystyle~{}+\frac{1}{2\sqrt{\pi}}(1-\lambda^{2})^{3/4}\frac{\lambda^{\ell}}{\ell^{3/2}}\left[1+\frac{3}{2\ell}\left(\frac{1}{2(1-\lambda^{2})}-\frac{1}{4}\right)+O\left(\frac{1}{\ell^{2}}\right)\right]\Bigg{]}$
(S31) $\displaystyle s^{xz}(\ell;t)\bigg{|}_{\lambda<1}=$ $\displaystyle~{}0$
(S32) $\displaystyle s^{yz}(\ell;t)\bigg{|}_{\lambda<1}=$
$\displaystyle~{}0\,.$ (S33)

##  S4 The “sudden quench” correlation functions

In this appendix we detail calculations of the “sudden quench” correlation
functions. We follow the approach of Calabrese, Essler and Fagotti [49–51] in
computing one- and two-point correlation functions following a quantum quench
of the transverse field in the quantum Ising chain. Much of the discussion
will follow similar lines, although we will consider a number of observables
that were not computed in their work.

###  S4.1 Fermionization of the model

Let us begin by describing the exact solution of the transverse field Ising
model in terms of free fermions. The transverse field Ising model has the
Hamiltonian

$H=-J\sum_{l=1}^{L}\sigma^{z}_{l}\sigma^{z}_{l+1}-g\sum_{l=1}^{L}\sigma^{x}_{l}\,,$
(S34)

where $\sigma^{\alpha}_{l}$ are the spin-1/2 operators acting on the $l$th
site of an $L$ chain (we take $L$ even), $J$ is the exchange interaction, and
$g$ is the transverse field strength. We consider periodic boundary
conditions, identifying $\sigma^{\alpha}_{L+1}=\sigma^{\alpha}_{1}$. To
express the model in terms of free fermions, we first perform a $3\pi/2$
rotation of the quantization axes about the $y$ axis to be consistent with
standard notations

$\displaystyle H$
$\displaystyle=-J\sum_{l=1}^{L}\tau^{x}_{l}\tau^{x}_{l+1}-g\sum_{l=1}^{L}\tau^{z}_{l}$
$\displaystyle\equiv-\bar{J}\sum_{l=1}^{L}\Bigg{(}\tau^{x}_{l}\tau^{x}_{l+1}+\bar{g}\tau^{z}_{l}\Bigg{)}.$
(S35)

Here $\tau^{\alpha}_{l}$ are the Pauli matrices with the rotated quantization
axes, $\bar{J}=J$ and $\bar{g}=g/J$.

We proceed to fermionize the problem via the Jordan-Wigner transformation [69]

$\begin{split}\tau^{z}_{l}&=1-2a^{\dagger}_{l}a_{l},\\\
\tau^{-}_{l}&=\exp\left(i\pi\sum_{j<l}a^{\dagger}_{j}a_{j}\right)a^{\dagger}_{l},\\\
\tau^{+}_{l}&=a_{l}\exp\left(i\pi\sum_{j<l}a^{\dagger}_{j}a_{j}\right).\end{split}$
(S36)

where we define the raising/lowering operators $\tau^{\pm}=(\tau^{x}\pm
i\tau^{y})/2$ and the fermions have canonical anticommutation relations
$\\{a_{l},a^{\dagger}_{l^{\prime}}\\}=\delta_{l,l^{\prime}}$. Following this
transformation, the Hamiltonian reads

$\begin{split}H=&-\bar{J}\sum_{l=1}^{L-1}\Big{(}a^{\dagger}_{l}-a_{l}\Big{)}\Big{(}a_{l+1}+a_{l+1}^{\dagger}\Big{)}\\\
&+\bar{J}e^{i\pi\hat{N}}\Big{(}a_{L}-a_{L}^{\dagger}\Big{)}\Big{(}a_{1}+a_{1}^{\dagger}\Big{)}\\\
&+\bar{J}\bar{g}\sum_{l=1}^{L}\Big{(}a^{\dagger}_{l}a_{l}-a_{l}a^{\dagger}_{l}\Big{)},\end{split}$
(S37)

where $\hat{N}=\sum_{l=1}^{L}a^{\dagger}_{l}a_{l}$ is the total number
operator for the fermions. It is easy to see that fermion number parity is
conserved by the Hamiltonian, so the exponential in the second term is

$\exp(i\pi\hat{N})=\left\\{\begin{array}[]{lcl}+1&&\text{for even numbers of
fermions},\\\ -1&&\text{for odd numbers of fermions}.\end{array}\right.$ (S38)

Thus the Hilbert space ${\cal H}={\cal H}_{e}\otimes{\cal H}_{o}$ of the model
splits into two sectors, for $\langle\hat{N}\rangle$ even or odd, which
correspond to antiperiodic or periodic boundary conditions on the fermions. We
need to address each of these cases in turn.

####  S4.1.1 Even numbers of fermions

In the even sector, we have $\exp(i\pi\hat{N})=+1$ and the Hamiltonian (S37)
reads

$\begin{split}H_{e}=&-\bar{J}\sum_{l=1}^{L}\Big{(}a^{\dagger}_{l}-a_{l}\Big{)}\Big{(}a_{l+1}+a_{l+1}^{\dagger}\Big{)}\\\
&+\bar{J}\bar{g}\sum_{l=1}^{L}\Big{(}a^{\dagger}_{l}a_{l}-a_{l}a^{\dagger}_{l}\Big{)},\end{split}$
(S39)

where the fermions obey antiperiodic boundary conditions $a_{L+1}=-a_{1}$. The
Hamiltonian can be diagonalized by first Fourier transforming and then
performing a Bogoliubov rotation. The Fourier transform reads

$c_{k_{n}}=\frac{1}{\sqrt{L}}\sum_{l=1}^{L}e^{ik_{n}l}a_{l},$ (S40)

with $k_{n}=\pi(2n+1)/L$ with $n=-L/2,-L/2+1,\ldots,L/2-1$. In terms of the
Fourier modes, the Hamiltonian reads

$\begin{split}H_{e}=-\bar{J}\sum_{n}&\bigg{[}\cos(k_{n})\left(c^{\dagger}_{k_{n}}c_{k_{n}}-c_{-k_{n}}c_{-k_{n}}^{\dagger}\right)\\\
&-i\sin(k_{n})\left(c^{\dagger}_{k_{n}}c^{\dagger}_{-k_{n}}-c_{-k_{n}}c_{k_{n}}\right)\\\
&-\bar{g}\left(c^{\dagger}_{k_{n}}c_{k_{n}}-c_{-k_{n}}c_{-k_{n}}^{\dagger}\right)\bigg{]}.\end{split}$
(S41)

This can be diagonalized via the Bogoliubov rotation

$\left(\begin{array}[]{c}c_{k_{n}}\\\
c^{\dagger}_{-k_{n}}\end{array}\right)=\left(\begin{array}[]{cc}\cos\left(\frac{\theta_{k_{n}}}{2}\right)&i\sin\left(\frac{\theta_{k_{n}}}{2}\right)\\\
i\sin\left(\frac{\theta_{k_{n}}}{2}\right)&\cos\left(\frac{\theta_{k_{n}}}{2}\right)\end{array}\right)\left(\begin{array}[]{c}\alpha_{k_{n}}\\\
\alpha^{\dagger}_{-k_{n}}\end{array}\right),$ (S42)

where the Bogoliubov angle $\theta_{k_{n}}$ satisfies

$\theta_{k_{n}}=-\text{sgn}(k_{n})\arccos\left(\frac{\bar{g}-\cos(k_{n})}{\sqrt{\bar{g}^{2}-2\bar{g}\cos(k_{n})+1}}\right)$
(S43)

to obtain the Hamiltonian

$H_{e}=\sum_{n}\epsilon_{\bar{g}}(k_{n})\Big{(}\alpha^{\dagger}_{k_{n}}\alpha_{k_{n}}-\frac{1}{2}\Big{)},$
(S44)

with dispersion relation

$\epsilon_{\bar{g}}(k_{n})=2\bar{J}\sqrt{\bar{g}^{2}+2\bar{g}\cos(k_{n})+1}.$
(S45)

Thus the Hamiltonian has been diagonalized.

Let us now say a few words about the nature of the ground state in the even
sector. As is expected from a Bogoliubov rotation (and as we well know from
physics of the Ising chain), the dispersion is generically gapped, i.e.
$\epsilon_{\bar{g}}(k_{n})>0$. As a result, the ground state
$|\Omega_{e},\bar{g}\rangle$ is the vacuum for the Bogoliubov fermions:

$\alpha_{k_{n}}|\Omega_{e},\bar{g}\rangle=0,\quad\forall n.$ (S46)

It immediately follows that the ground state energy is

$E_{\Omega_{e},\bar{g}}=-\frac{1}{2}\sum_{n}\epsilon_{\bar{g}}(k_{n}).$ (S47)

It will be useful to understand the ground state $|\Omega_{e},\bar{g}\rangle$
in terms of the Jordan-Wigner fermions $c_{k_{n}}$. Inverting Eq. (S42) to
express $\alpha_{k_{n}}$ in terms $c_{k_{n}}$, one finds that Eq. (S46)
becomes ($\forall n$)

$\Big{[}\cos\left(\frac{\theta_{k_{n}}}{2}\right)c_{k_{n}}-i\sin\left(\frac{\theta_{k_{n}}}{2}\right)c^{\dagger}_{-k_{n}}\Big{]}|\Omega_{e},\bar{g}\rangle=0.$
(S48)

As $|\Omega_{e},\bar{g}\rangle$ is in the even sector with zero momentum, it
must be of the form $\prod_{k_{n}\geq
0}(a+bc^{\dagger}_{-k_{n}}c^{\dagger}_{k_{n}})|0\rangle$, where $|0\rangle$ is
the vacuum for Jordan-Wigner fermions. Working through some algebra, one finds

$|\Omega_{e},\bar{g}\rangle=\prod_{k_{n}\geq
0}\left[\cos\left(\frac{\theta_{k_{n}}}{2}\right)-i\sin\left(\frac{\theta_{k_{n}}}{2}\right)c^{\dagger}_{-k_{n}}c_{k_{n}}^{\dagger}\right]|0\rangle.$
(S49)

This can be written in a “squeezed state” form following some simple
manipulations

$\begin{split}|\Omega_{e},\bar{g}\rangle=&\left[\prod_{k_{n}\geq
0}\cos\left(\frac{\theta_{k_{n}}}{2}\right)\right]\\\
&\times\exp\left[-i\sum_{k_{n}\geq
0}\tan\left(\frac{\theta_{k_{n}}}{2}\right)c^{\dagger}_{-k_{n}}c^{\dagger}_{k_{n}}\right]|0\rangle.\end{split}$
(S50)

####  S4.1.2 Odd number of fermions

The story with the odd fermion number sector is similar to the even one. This
time we have $\exp(i\pi\hat{N})=-1$ and instead of antiperiodic boundary
conditions, we have periodic ones $a_{L+1}=a_{1}$. The Fourier transform of
the Jordan-Wigner fermions then reads

$c_{p_{n}}=\frac{1}{\sqrt{L}}\sum_{l}e^{ip_{n}l}a_{l}$ (S51)

with the momenta $p_{n}=2n\pi/L$ and $n=-L/2,-L/2+1,\ldots...$. We again
perform a Bogoliubov rotation

$\left(\begin{array}[]{c}c_{p_{n}}\\\
c^{\dagger}_{-p_{n}}\end{array}\right)=\left(\begin{array}[]{cc}\cos\left(\frac{\theta_{p_{n}}}{2}\right)&i\sin\left(\frac{\theta_{p_{n}}}{2}\right)\\\
i\sin\left(\frac{\theta_{p_{n}}}{2}\right)&\cos\left(\frac{\theta_{p_{n}}}{2}\right)\end{array}\right)\left(\begin{array}[]{c}\alpha_{p_{n}}\\\
\alpha^{\dagger}_{-p_{n}}\end{array}\right),$ (S52)

with the Bogoliubov angle as defined in the even sector, with the additional
definition of $\theta_{p_{0}}=0$. The Hamiltonian (S37) is then

$H_{o}=\sum_{n\neq
0}\epsilon_{\bar{g}}(p_{n})\Big{(}\alpha^{\dagger}_{p_{n}}\alpha_{p_{n}}-\frac{1}{2}\Big{)}-2\bar{J}(1-\bar{g})\Big{(}\alpha^{\dagger}_{0}\alpha_{0}-\frac{1}{2}\Big{)}.$
(S53)

We see now that the ground state in the odd sector is

$|\Omega_{o}\rangle=\alpha^{\dagger}_{0}|0_{\alpha}\rangle,$ (S54)

where $|0_{\alpha}\rangle$ is the vacuum state for $\alpha_{p_{n}}$. The
ground state $|\Omega_{o}\rangle$ has energy

$E_{\Omega_{o}}=-\bar{J}(1-\bar{g})-\frac{1}{2}\sum_{n\neq
0}\epsilon_{\bar{g}}(p_{n}).$ (S55)

As with the even sector, it will be useful to know how to express the odd
ground state, Eq. (S54), in terms of the Jordan-Wigner fermions. Following the
same sequence of manipulations as the even sector, one can express the state
in the form

$\begin{split}|\Omega_{o},\bar{g}\rangle=&\left(\prod_{p_{n}\geq
0}\cos\left(\frac{\theta_{p_{n}}}{2}\right)\right)\\\
&\times\exp\left(-i\sum_{p_{n}\geq
0}\tan\left(\frac{\theta_{p_{n}}}{2}\right)c^{\dagger}_{-p_{n}}c^{\dagger}_{p_{n}}\right)c^{\dagger}_{0}|0\rangle.\end{split}$
(S56)

As expected, this state has odd fermion parity.

###  S4.2 Time-evolution following a quantum quench

Let us now turn to the problem of describing time-evolution following a
quantum quench. In particular, we consider the scenario where for time $t<0$
the system is in an eigenstate of the Hamiltonian $H_{i}$ with $\bar{g}=g_{i}$
and at time $t=0$ this is suddenly changed to $\bar{g}=g_{f}$ and the state
evolves under this new Hamiltonian $H_{f}$ for all times $t>0$.

At times $t<0$ we can diagonalize the Hamiltonian in terms of fermions $\beta$
by performing the Bogoliubov rotation with angle $\phi$. For times $t>0$ we
can similarly diagonalize the Hamiltonian in terms of fermions $\alpha$ via
the Bogoliubov rotation with angle $\theta$. Equating the expressions for the
Jordan-Wigner fermions $a_{k}$ in terms of the Bogoliubov fermions, we have

$\displaystyle\cos\left(\frac{\theta_{k_{n}}}{2}\right)\alpha_{k_{n}}+i\sin\left(\frac{\theta_{k_{n}}}{2}\right)\alpha_{-k_{n}}^{\dagger}$
$\displaystyle=\cos\left(\frac{\phi_{k_{n}}}{2}\right)\beta_{k_{n}}+i\sin\left(\frac{\phi_{k_{n}}}{2}\right)\beta^{\dagger}_{-k_{n}},$
(S57) $\displaystyle
i\sin\left(\frac{\theta_{k_{n}}}{2}\right)\alpha_{k_{n}}+\cos\left(\frac{\theta_{k_{n}}}{2}\right)\alpha^{\dagger}_{-k_{n}}$
$\displaystyle=i\sin\left(\frac{\phi_{k_{n}}}{2}\right)\beta_{k_{n}}+\cos\left(\frac{\phi_{k_{n}}}{2}\right)\beta^{\dagger}_{-k_{n}}.$
(S58)

The initial state, an eigenstate of $H_{i}$ will be expressed
straightforwardly in terms of the $\beta$ fermions. Time-evolution, however,
will be easy to compute in terms of the $\alpha$ fermions via

$\alpha_{k_{n}}(t)=\exp\left(-i\epsilon_{g_{f}}(k_{n})t\right)\alpha_{k_{n}}.$
(S59)

Thus we write the $\beta$ fermions in terms of $\alpha$ ones and time-evolve
them

$\beta_{k_{n}}(t)=\cos\left(\frac{\theta_{k_{n}}-\phi_{k_{n}}}{2}\right)e^{-i\epsilon_{g_{f}}(k_{n})t}\alpha_{k_{n}}+i\sin\left(\frac{\theta_{k_{n}}-\phi_{k_{n}}}{2}\right)e^{i\epsilon_{g_{f}}(k_{n})t}\alpha^{\dagger}_{-k_{n}},$
(S60)

before re-expressing the $\alpha$ fermions in terms of the $\beta$ ones:

$\beta_{k_{n}}(t)=\bigg{[}\cos\Big{(}\epsilon_{g_{f}}(k_{n})t\Big{)}-i\sin\Big{(}\epsilon_{g_{f}}(k_{n})t\Big{)}\cos\Big{(}\theta_{k_{n}}-\phi_{k_{n}}\Big{)}\bigg{]}\beta_{k_{n}}-\sin\Big{(}\epsilon_{g_{f}}(k_{n})t\Big{)}\sin\Big{(}\theta_{k_{n}}-\phi_{k_{n}}\Big{)}\beta^{\dagger}_{-k_{n}}.$
(S61)

We see that the difference in Bogoliubov angles often appears, so we denote
this by the short hand

$\Delta_{k}=\theta_{k}-\phi_{k}.$ (S62)

With this, and with some work, we can compute time-evolution of observables.
In the following we first consider one-point functions, before continuing to
two-point functions.

###  S4.3 Time-evolution of one-point functions

We will consider the quench where we start from the ground state of the
initial Hamiltonian $H_{i}\equiv H(t<0)$ with $\bar{g}=g_{i}$ and time-evolve
according to the Hamiltonian $H_{f}$ where $\bar{g}=g_{f}$. In the
thermodynamic limit, the initial state will be a superposition of the ground
states in the even and odd sectors

$|\Omega,\bar{g}\rangle=\frac{1}{\sqrt{2}}\Big{(}|\Omega_{e},\bar{g}\rangle\\!\rangle+|\Omega_{o},\bar{g}\rangle\\!\rangle\Big{)}.$
(S63)

(Here we work with normalized states such that
$\langle\\!\langle\Omega_{i},\bar{g}|\Omega_{i},\bar{g}\rangle\\!\rangle=1$.)
Thus, generally, there will be three sorts of terms when computing
observables: expectation values within the even sector, those within the odd
sector, and those that mix the two sectors. Depending on the operator being
considered, the mixing terms may be forbidden.

Let us begin with the easiest case, and consider the time-evolution of the
one-point function of $\tau^{z}_{l}$. Under the Jordan-Wigner transformation
(S36) this maps to the local number operator of the Jordan-Wigner fermions.
Thus we want to compute

$t^{z}(t)=\langle\Omega,g_{i}|e^{iH_{f}t}\Big{(}1-2a^{\dagger}_{l}a_{l}\Big{)}e^{-iH_{f}t}|\Omega,g_{i}\rangle.$
(S64)

Clearly $t^{z}(t)$ is independent of $l$ by translational invariance. As
$\tau^{z}$ is number conserving in terms of the fermions, the terms mixing
even and odd sectors from Eq. (S63) will vanish. Thus we have

$\begin{split}t^{z}(t)=&\langle\\!\langle\Omega_{e},g_{i}|\bigg{(}\frac{1}{2}-a^{\dagger}_{l}(t)a_{l}(t)\bigg{)}|\Omega_{e},g_{i}\rangle\\!\rangle\\\
&+\langle\\!\langle\Omega_{o},g_{i}|\bigg{(}\frac{1}{2}-a^{\dagger}_{l}(t)a_{l}(t)\bigg{)}|\Omega_{o},g_{i}\rangle\\!\rangle.\end{split}$
(S65)

Let us take the first term, where the ground states is in the even sector of
the Hilbert space. We express the time-evolved Jordan-Wigner fermions in terms
of the time-evolved Bogoliubov fermions (S61). This yields

$\displaystyle
a_{l}(t)_{e}=\frac{1}{\sqrt{L}}\sum_{n}e^{-ik_{n}l}c_{k_{n}}(t)=\frac{1}{\sqrt{L}}\sum_{n}e^{-ik_{n}l}\left[\cos\left(\frac{\phi_{k_{n}}}{2}\right)\beta_{k_{n}}(t)+i\sin\left(\frac{\phi_{k_{n}}}{2}\right)\beta^{\dagger}_{-k_{n}}(t)\right].$
(S66)

Working through some tedious algebra, one finds

$\displaystyle t^{z}(t)_{e}$
$\displaystyle=\frac{1}{2}-\frac{1}{L}\sum_{n}\left[\sin^{2}\left(\epsilon_{g_{f}}(k_{n})t\right)\sin^{2}\left(\Delta_{k_{n}}+\frac{\phi_{k_{n}}}{2}\right)+\sin^{2}\left(\frac{\phi_{k_{n}}}{2}\right)\cos^{2}\left(\epsilon_{g_{f}}(k_{n})t\right)\right],$
(S67) $\displaystyle=\frac{1}{2}-\frac{1}{2\pi}\int_{-\pi}^{\pi}{\rm
d}k\left[\sin^{2}\left(\epsilon_{g_{f}}(k)t\right)\sin^{2}\left(\Delta_{k}+\frac{\phi_{k}}{2}\right)+\sin^{2}\left(\frac{\phi_{k}}{2}\right)\cos^{2}\left(\epsilon_{g_{f}}(k)t\right)\right].$
(S68)

In the second line we use the thermodynamic limit to replace the sum over
discrete momentum modes by an integral. Performing a similar calculation in
the odd sector, we obtain the same result in the thermodynamic limit, leading
to

$t^{z}(t)=1-\frac{1}{\pi}\int_{-\pi}^{\pi}{\rm
d}k\left[\sin^{2}\left(\epsilon_{g_{f}}(k)t\right)\sin^{2}\left(\Delta_{k}+\frac{\phi_{k}}{2}\right)+\sin^{2}\left(\frac{\phi_{k}}{2}\right)\cos^{2}\left(\epsilon_{g_{f}}(k)t\right)\right].$
(S69)

Let us now move on to considering $\tau^{x}_{l}$ and $\tau^{y}_{l}$. These
operators are similar in that they change the number of Jordan-Wigner fermions
by one, and hence their expectation values consist only of the “cross terms”
between even and odd ground states. This is actually quite problematic for the
free fermion techniques discussed here, as it is not clear how to deal with
such terms. Instead, we use the following trick: we relate the one-point
function to the two-point function in the limit of infinite separation of the
operators

$\lim_{\ell\to\infty}\langle\\!\langle\Omega;g_{i}|\tau^{x}_{j}(t)\tau^{x}_{j+\ell}(t)|\Omega;g_{i}\rangle\\!\rangle\to\langle\\!\langle\Omega;g_{i}|\tau^{x}_{j}(t)|\Omega;g_{i}\rangle\\!\rangle^{2},$
(S70)

and similarly for $\tau^{y}$. As such, this naturally brings us on to the
topic of computing two-point functions.

###  S4.4 Time-evolution of two-point functions

We now turn our attention to computing two-point functions

$t^{\alpha\beta}(\ell;t)=\langle\Omega,g_{i}|e^{iH_{f}t}\tau^{\alpha}_{j}\tau^{\beta}_{j+\ell}e^{-iH_{f}t}|\Omega,g_{i}\rangle.$
(S71)

Two examples of these, $\alpha=\beta=x$ and $\alpha=\beta=z$, were considered
in Refs. [49–51]. We will also consider these here, as well as some other
cases that will be of use in this work. As with the one-point functions, we
are unable to compute two-point functions when the operator within the
expectation value is odd under spin inversion (and thus connects states in the
even and odd sectors of the Hilbert space). Thus we restrict our attention to
operators that are even under spin inversion.

This restriction, in the context of the main text, leads us to focus our
results on one of two cases. Firstly, we can consider quenches starting in
states that are even under spin inversion (and thus the expectation value
vanishes by symmetry). This allows us, for example, to compute the full time
evolution of two-point functions in the driven problem where the “effective
sudden quench” is within the paramagnetic phase of the static Hamiltonian.
Secondly, we can restrict attention to stroboscopic times, where contributions
from these problematic expectation values vanish.

####  S4.4.1 The $t^{zz}(\ell;t)$ two-point function

Let us begin with the easiest case, where the operators in the two-point
function are solely in the transverse field direction. Then, according to the
Jordan-Wigner transformation (S36), we need to compute the two-point function
of the (time-dependent) Jordan-Wigner fermion density

$t^{zz}(\ell;t)=\langle\Omega,g_{i}|e^{iH_{f}t}\Big{(}1-2a^{\dagger}_{j}a_{j}\Big{)}\Big{(}1-2a^{\dagger}_{j+\ell}a_{j+\ell}\Big{)}e^{-iH_{f}t}|\Omega,g_{i}\rangle.$
(S72)

The operator under consideration is even under spin inversion, thus matrix
elements between states in the even and odd sectors vanish, leaving only

$t^{zz}(\ell;t)=\frac{1}{2}\langle\langle\Omega_{e},g_{i}|e^{iH_{f}t}\Big{(}1-2a^{\dagger}_{j}a_{j}\Big{)}\Big{(}1-2a^{\dagger}_{j+\ell}a_{j+\ell}\Big{)}e^{-iH_{f}t}|\Omega_{e},g_{i}\rangle\rangle+(e\leftrightarrow
o).$ (S73)

Using Eq. (S66) (and similar for the odd sector) and taking the thermodynamic
limit, one finds

$\displaystyle t^{zz}(\ell;t)$ $\displaystyle=-\int_{-\pi}^{\pi}\frac{{\rm
d}k}{2\pi}\frac{{\rm
d}k^{\prime}}{2\pi}\Bigg{\\{}e^{i(k-k^{\prime})\ell}\Big{[}\sin(2\epsilon_{g_{f}}(k)t)\sin(\Delta_{k})\sin(2\epsilon_{g_{f}}(k^{\prime})t)\sin(\Delta_{k^{\prime}})+\Big{(}e^{i(k-k^{\prime})\ell}-1\Big{)}$
$\displaystyle\qquad\times
e^{i(\theta_{k}+\theta_{k^{\prime}})}\Big{(}\cos(\Delta_{k})-i\sin(\Delta_{k})\cos(2\epsilon_{g_{f}}(k)t)\Big{)}\Big{(}\cos(\Delta_{k^{\prime}})-i\sin(\Delta_{k^{\prime}})\cos(2\epsilon_{g_{f}}(k^{\prime})t)\Big{)}\Bigg{\\}}.$
(S74)

####  S4.4.2 The $t^{xx}(\ell;t)$ two-point function

Let us now turn our attention towards the two-point function where both
operators are in the Ising direction. This is significantly more complicated
than $t^{zz}(\ell;t)$ as the Jordan-Wigner transformation (S36) introduces
“strings” of density operators between $j$ and $j+\ell$. This can easily be
seen by Taylor expanding the exponential factors in Eq. (S36):

$\displaystyle\exp\left(\pm
i\pi\sum_{j<\ell}a^{\dagger}_{j}a_{j}\right)=\prod_{j<\ell}\exp\left(\pm i\pi
a^{\dagger}_{j}a_{j}\right)=\prod_{j<\ell}\left(1-2a_{j}^{\dagger}a_{j}\right).$
(S75)

Thus, in terms of the Jordan-Wigner fermions, we need to compute

$\displaystyle t^{xx}(\ell;t)$
$\displaystyle=\langle\Omega,g_{i}|e^{iH_{f}t}(\tau^{+}_{j}+\tau^{-}_{j})(\tau^{+}_{j+\ell}+\tau^{-}_{j+\ell})e^{-iH_{f}t}|\Omega,g_{i}\rangle,$
(S76)
$\displaystyle=\langle\Omega,g_{i}|e^{iH_{f}t}\left(a_{j}+a^{\dagger}_{j}\right)\left[\prod_{k=0}^{\ell-1}\left(1-2a^{\dagger}_{j+k}a_{j+k}\right)\right]\left(a_{j+\ell}+a^{\dagger}_{j+\ell}\right)e^{-iH_{f}t}|\Omega,g_{i}\rangle,$
(S77)
$\displaystyle=\langle\Omega,g_{i}|e^{iH_{f}t}\left(a^{\dagger}_{j}-a_{j}\right)\left[\prod_{k=1}^{\ell-1}\left(1-2a^{\dagger}_{j+k}a_{j+k}\right)\right]\left(a_{j+\ell}+a^{\dagger}_{j+\ell}\right)e^{-iH_{f}t}|\Omega,g_{i}\rangle.$
(S78)

The calculation will proceed more easily if we introduce Majorana fermions,
which describe the real and imaginary parts of the Jordan-Wigner fermions,

$\gamma^{x}_{l}=a_{l}+a^{\dagger}_{l},\qquad\gamma^{y}_{l}=-i(a_{l}-a_{l}^{\dagger}),$
(S79)

which obey the anticommutation relations
$\\{\gamma^{a}_{j},\gamma^{b}_{l}\\}=2\delta_{a,b}\delta_{j,l}$. In terms of
these Majorana operators we have

$a^{\dagger}_{j}a_{j}=\frac{1}{2}\left(1+i\gamma_{j}^{x}\gamma^{y}_{j}\right),\quad\exp(i\pi
a^{\dagger}_{j}a_{j})=-i\gamma^{x}_{j}\gamma^{y}_{j},$ (S80)

and hence

$t^{xx}(\ell;t)=\langle\langle\Omega_{e},g_{i}|e^{iH_{f}t}\prod_{k=1}^{\ell}(-i\gamma^{y}_{k}\gamma^{x}_{k+1})e^{-iH_{f}t}|\Omega_{e},g_{i}\rangle\rangle.$
(S81)

This can be written as the Pfaffian of a $2\ell\times 2\ell$ matrix

$\displaystyle t^{xx}(\ell;t)$
$\displaystyle=(-i)^{\ell}\text{Pf}(\Gamma)=\text{Pf}(-i\Gamma),$ (S82)
$\displaystyle\Gamma$
$\displaystyle=\left(\begin{array}[]{cccc}\Gamma^{\prime}_{11}&\Gamma^{\prime}_{12}&\cdots&\Gamma^{\prime}_{1\ell}\\\
\Gamma^{\prime}_{21}&\Gamma^{\prime}_{22}&\cdots&\Gamma^{\prime}_{2\ell}\\\
\vdots&\vdots&\ddots&\vdots\\\ \Gamma^{\prime}_{\ell 1}&\Gamma^{\prime}_{\ell
2}&\cdots&\Gamma^{\prime}_{\ell\ell}\end{array}\right),$ (S87)

where $\Gamma^{\prime}$ is a $2\times 2$ sub-matrix

$\displaystyle\Gamma^{\prime}_{mn}=\left(\begin{array}[]{cc}\langle\langle\gamma^{y}_{m-n+j}\gamma^{y}_{j}\rangle\rangle_{t}-\delta_{m,n}&\langle\langle\gamma^{y}_{m-n-1+j}\gamma^{x}_{j}\rangle\rangle_{t}\\\
\langle\langle\gamma^{x}_{m-n+1+j}\gamma^{y}_{j}\rangle\rangle_{t}&\langle\langle\gamma^{x}_{m-n+j}\gamma^{x}_{j}\rangle\rangle_{t}-\delta_{m,n}\end{array}\right),$
(S90)

where

$\langle\langle
O\rangle\rangle_{t}=\langle\langle\Omega_{e},g_{i}|e^{iH_{f}t}Oe^{-iH_{f}t}|\Omega_{e},g_{i}\rangle\rangle.$
(S91)

Using reflection symmetry and translational invariance, we can rewrite this
sub-matrix as

$\displaystyle\bar{\Gamma}_{\eta}$
$\displaystyle=-i\Gamma^{\prime}_{m+\eta,m},$ (S92)
$\displaystyle=-i\left(\begin{array}[]{cc}\langle\langle\gamma^{y}_{\eta+j}\gamma^{y}_{j}\rangle\rangle_{t}-\delta_{\eta,0}&\langle\langle\gamma^{y}_{\eta-1+j}\gamma^{x}_{j}\rangle\rangle_{t}\\\
-\langle\langle\gamma^{y}_{-\eta-1+j}\gamma^{x}_{j}\rangle\rangle_{t}&-\langle\langle\gamma^{y}_{\eta+j}\gamma^{y}_{j}\rangle\rangle_{t}+\delta_{\eta,0}\end{array}\right)$
(S95) $\displaystyle=\left(\begin{array}[]{cc}-f_{\eta}&g_{\eta}\\\
-g_{-\eta}&f_{\eta}\end{array}\right),$ (S98)

where $f_{\eta}$ and $g_{\eta}$ can be obtained using Eq. (S66), and they
read:

$\displaystyle f_{\eta}$ $\displaystyle=\int_{-\pi}^{\pi}\frac{{\rm
d}k}{2\pi}e^{ik\eta}f(k),$ (S99) $\displaystyle f(k)$
$\displaystyle=i\sin\Big{(}2\epsilon_{g_{f}}(k)t\Big{)}\sin\Big{(}\Delta_{k}\Big{)},$
(S100) $\displaystyle g_{\eta}$ $\displaystyle=\int_{-\pi}^{\pi}\frac{{\rm
d}k}{2\pi}e^{ik\eta}g(k),$ (S101) $\displaystyle g(k)$
$\displaystyle=-e^{i\Delta_{k}}\bigg{[}\cos\Big{(}\Delta_{k}\Big{)}-i\sin\Big{(}\Delta_{k}\Big{)}\cos\Big{(}2\epsilon_{g_{f}}(k)t\Big{)}\bigg{]}.$
(S102)

When the two-point correlation function is written in terms of (S98), it is
the Pfaffian of a block Toeplitz matrix:

$t^{xx}(\ell;t)=\text{Pf}\left(\begin{array}[]{cccc}\bar{\Gamma}_{0}&\bar{\Gamma}_{-1}&\cdots&\bar{\Gamma}_{1-\ell}\\\
\bar{\Gamma}_{1}&\bar{\Gamma}_{0}&\cdots&\bar{\Gamma}_{2-\ell}\\\
\vdots&\vdots&\ddots&\vdots\\\
\bar{\Gamma}_{\ell-1}&\bar{\Gamma}_{\ell-2}&\cdots&\bar{\Gamma}_{0}\end{array}\right).$
(S103)

The matrix here is real and antisymmetric, thus it is an anti-Hermitian
matrix. All eigenvalues are imaginary and appear in conjugate pairs.
Furthermore, for an antisymmetric matrix the Pfaffian is the square root of
the determinant, up to a difficult-to-determine sign. In other words, if an
$2\ell\times 2\ell$ antisymmetric matrix $\Lambda=-\Lambda^{T}$ has
eigenvalues $\pm\lambda_{j}$ ($j=1,\ldots,\ell$), the Pfaffian reads
$\text{Pf}(\Lambda)=\prod_{j=1}^{\ell}i\lambda_{j}$.

One can reduce computing the Pfaffian to computing the determinant of an
$\ell\times\ell$ matrix $Q$

$\displaystyle t^{xx}(\ell;t)$
$\displaystyle=(-1)^{\ell(\ell-1)/2}\text{det}(Q),$ (S104) $\displaystyle
Q_{nm}$ $\displaystyle=if_{n-m}+g_{n+m-\ell-1},$ (S105)

which is convenient to do numerically. One can also make analytical progress
by computing the determinant of the block Toeplitz matrix, as considered in
Ref. [49–51], via applications of Szegö’s theorem and Fisher-Hartwig
conjectures [70–76].

####  S4.4.3 The $t^{yy}(\ell;t)$ two-point function

The computation of this two-point function is very similar to the previous
one. In terms of the Majorana fermions,

$t^{yy}(\ell;t)=\langle\langle\Omega_{e},g_{i}|e^{iH_{f}t}\prod_{k=1}^{\ell}(i\gamma^{x}_{k}\gamma^{y}_{k+1})e^{-iH_{f}t}|\Omega_{e},g_{i}\rangle\rangle.$
(S106)

The previous expression can be written as the Pfaffian of a block Toeplitz
matrix,

$t^{yy}(\ell;t)=\text{Pf}\left(\begin{array}[]{cccc}\tilde{\Gamma}_{0}&\tilde{\Gamma}_{-1}&\cdots&\tilde{\Gamma}_{1-\ell}\\\
\tilde{\Gamma}_{1}&\tilde{\Gamma}_{0}&\cdots&\tilde{\Gamma}_{2-\ell}\\\
\vdots&\vdots&\ddots&\vdots\\\
\tilde{\Gamma}_{\ell-1}&\tilde{\Gamma}_{\ell-2}&\cdots&\tilde{\Gamma}_{0}\end{array}\right).$
(S107)

where $\tilde{\Gamma}_{\eta}$ is a $2\times 2$ sub-matrix

$\tilde{\Gamma}_{\eta}=\left(\begin{array}[]{cc}-f_{\eta}&g_{-\eta+2}\\\
-g_{\eta+2}&f_{\eta}\end{array}\right),$ (S108)

where $f_{\eta}$ and $g_{\eta}$ are defined in (S99)-(S101).

The computation of the Pfaffian is equivalent to the one the determinant of an
$\ell\times\ell$ matrix $\tilde{Q}$

$\displaystyle t^{yy}(\ell;t)$
$\displaystyle=(-1)^{\ell(\ell-1)/2}\text{det}(\tilde{Q}),$ (S109)
$\displaystyle\tilde{Q}_{nm}$ $\displaystyle=if_{n-m}+g_{\ell+3-(n+m)}.$
(S110)

####  S4.4.4 The $t^{xy}(\ell;t)$ two-point function

This two point-function can be obtained from $t^{xx}(\ell;t)$. Indeed, by
considering the Heisenberg equation of motion,

$\displaystyle\frac{{\rm d}}{{\rm d}t}\tau^{x}_{j}(t)\tau^{x}_{j+\ell}(t)$
$\displaystyle=i[H_{f},\tau^{x}_{j}(t)\tau^{x}_{j+\ell}(t)]$
$\displaystyle=2\bar{J}\bar{g}_{f}\left(\tau_{j}^{y}(t)\tau_{j+\ell}^{x}(t)+\tau_{j}^{x}(t)\tau_{j+\ell}^{y}(t)\right)\,.$
(S111)

By computing the expectation value and considering inversion symmetry we
conclude,

$t^{xy}(\ell;t)=\frac{1}{4\bar{J}\bar{g}_{f}}\frac{{\rm d}}{{\rm
d}t}t^{xx}(\ell;t)\,.$ (S112)

##  S5 Details of the numerical algorithm for time-evolution of a driven
system

In this section of the Supplemental Material, we discuss details of the
numerical algorithm used to compute time-evolution of the initial state used
in the main text. As we have already shown in the main text, the results of
the simulations have an excellent agreement with the analytical results for
the integrable case when $\Omega=2h^{z}$. We make the following numerical
studies outside this integrable line.

We consider the time-evolution of a given initial state $|\Psi_{0}\rangle$
induced by a time-dependent Hamiltonian, Eq. (1) of the main text. Formal
solution of the time-dependent Schrödinger equation gives

$|\Psi(t)\rangle=\mathbb{T}\exp\left(-i\int_{0}^{t}{\mathrm{d}}t^{\prime}\,H(t^{\prime})\right)|\Psi_{0}\rangle,$
(S113)

where $\mathbb{T}$ is the time-ordering operator. Our Hamiltonian $H(t)$
depends smoothly on time $t$, which makes dealing with the time-ordering
tricky. To proceed we “Trotterize” the time-evolution into $N$ small steps (of
size $\Delta t$, such that $N\Delta t=t$), each of which has a fixed
Hamiltonian

$\displaystyle|\Psi(t)\rangle$
$\displaystyle\approx\mathbb{T}\prod_{j=0}^{N-1}\exp\Big{(}-iH\big{(}j\Delta
t\big{)}\Delta t\Big{)}|\Psi_{0}\rangle,$ (S114)
$\displaystyle\equiv\mathbb{T}\prod_{j=1}^{N-1}{\cal U}_{j}|\Psi_{0}\rangle.$
(S115)

The problem now becomes how to describe the action of the time-evolution
operator ${\cal U}_{j}$ on a state.

We compute the action of a single step time-evolution operator on a state via
the Chebyshev expansion (see, e.g., Ref. [45]). This avoids the need to
diagonalize or exponentiate the Hamiltonian, and it is computationally
efficient. The Chebyshev expansion of the time-evolution operator reads:

${\cal U}_{j}=J_{0}\big{(}\widetilde{\Delta
t}\big{)}\mathbb{1}+2\sum_{n=1}^{\infty}(-i)^{n}J_{n}\big{(}\widetilde{\Delta
t}\big{)}T_{n}\Big{(}\widetilde{H}(j\Delta t)\Big{)}.$ (S116)

Here we have introduced the rescaled time $\widetilde{\Delta t}$ and
Hamiltonian $\widetilde{H}(j\Delta t)$, as well as the Bessel functions
$J_{n}(x)$ and the Chebyshev matrices $T_{n}(x)$, which satisfy the recursion
relation

$\displaystyle T_{0}(x)$ $\displaystyle=\mathbb{1},\quad T_{1}(x)=x,$ (S117)
$\displaystyle T_{n}(x)$ $\displaystyle=2xT_{n-1}(x)-T_{n-2}(x).$ (S118)

The rescaled time and Hamiltonian satisfy $\widetilde{H}(j\Delta
t)\widetilde{\Delta t}=H(j\Delta t)\Delta t$, and enforce that the eigenvalues
of the $\widetilde{H}(j\Delta t)$ lie in the interval $[-1,1]$.

In practice, the Chebyshev expansion (S116) has to be truncated to finite
order. In the following subsections we present some convergence checks of our
“Trotterization+Chebyshev” procedure, where we examine the effects of Trotter
step size and Chebyshev expansion order. We find that only a low order
Chebyshev expansion is required for the problem at hand, leading to a very
efficient algorithm for simulating the time-evolution of the continuously
(smoothly) driven quantum system. We finish our discussions with a brief
inspection of finite system size effects on the nonequilibrium dynamics.

###  S5.1 Convergence with order of the Chebyshev expansion

(a)
---
(b)
Figure S1: The time-evolution of the one-point function (a)
$s^{x}(t)=\langle\Psi(t)|\sigma^{x}_{1}|\Psi(t)\rangle$; (b)
$s^{zz}(1;t)=\langle\Psi(t)|\sigma^{z}_{j}\sigma^{z}_{j+1}|\Psi(t)\rangle$
starting from the ground state of the Hamiltonian (1) at $t=0$ with $J=1$,
$h^{z}=0$, and $h^{x}=2$ and time-evolved with $J=1$, $h^{z}=\Omega=1$, and
$g=1.5$. The time-evolution is computed via “Trotterization+Chebyshev” with
Trotter step $\Delta t=0.005$ and Chebyshev expansion order $N_{c}$ for ten
periods of the drive, $T=2\pi/\Omega$. We see that results are rapidly
converging with increased $N_{c}$. Note that there is no heating to infinite
temperature, despite the fact we are outside the integrable line
$\Omega=2h^{z}$.

Let us first examine how results of the time-evolution vary with the order of
the Chebyshev expansion (S116). We fix the Trotter step size $\Delta t=0.005$
and then proceed to compute an approximation of the time-evolution:

${\cal U}_{j}\approx J_{0}\big{(}\widetilde{\Delta
t}\big{)}\mathbb{1}+2\sum_{n=1}^{N_{c}}(-i)^{n}J_{n}\big{(}\widetilde{\Delta
t}\big{)}T_{n}\Big{(}\widetilde{H}(j\Delta t)\Big{)},$ (S119)

for given values of $N_{c}$. We present an example in Fig. S1(a) for the time-
evolution of a one-point function following a quench within the disordered
(paramagnetic) phase of the quantum Ising chain (details of parameters are
given in the figure caption). We observed that there is rapid convergence of
the result for increasing order of the Chebyshev expansion: after ten periods
of the drive, the results with $N_{c}=8$ and $N_{c}=32$ agree to seven decimal
places. We also stress the fact that, despite being outside the integrable
line $\Omega=2h^{z}$, there is no heating to infinite temperature.

One may question whether similarly good convergence is observed for more
complicated observables. In Fig. S1(b), we present $s^{zz}(1;t)$ for the same
quench, where it is apparent that two-point functions also converge similarly
well with increasing $N_{c}$. Within the main body of the text, we have
explicitly checked that the obtained results are independent of the expansion
order $N_{c}$.

###  S5.2 Convergence with Trotter step size

Let us now turn our attention to the size of the time discretization step.

###  S5.3 Finite size effects

With convergence established for Trotter step size and order of the Chebyshev
expansion, we finally examine finite size effects in our simulations of the
time-evolution. Fixing $\Delta t=0.005$ and $N_{c}=64$, we examine the
nonequilibrium time-evolution of $\sigma^{x}$ following the same quench as in
the previous subsections.

Figure S2: The time-evolution of the one-point function
$s^{x}(t)=\langle\Psi(t)|\sigma^{x}_{1}|\Psi(t)\rangle$ following the quench
presented in Fig. S1. Results are presented for a number of system sizes,
illustrating finite size effects in numerical data due to small accessible
systems. We see that for $L=14,16$ the time-evolution is well matched for
almost the whole period of the driving.

We first examine a local observable,
$\langle\Psi(t)|\sigma^{x}|\Psi(t)\rangle$, and check how it behaves with
changes in the system size $L$, i.e. what the “finite size effects” are.
Example data is presented in Fig. S2, where we see the expected behaviour: at
short times observables for all system sizes are in agreement. Under time-
evolution, where propagating excitations are generated, results for different
system sizes eventually diverge due to finite size revivals. For a single,
sudden quench the picture for this is simple: a quantum quench generates
excitations that can propagate around the system. By conservation of momentum,
such excitations must be generated in pairs, with momentum $k$ and $-k$. These
can propagate around the system, eventually meeting back at the start and
interfering – a purely finite size effect. The time for which this occurs
increases linearly with system size. Here we are seeing the driven system
analogue of this finite size revival behaviour. Over a single period, we see
that results for the two largest system sizes, $L=14$ and $L=16$, match over
almost the whole drive period.

Figure S3: The return amplitude
$\Big{|}\langle\Psi_{0}|\Psi(t)\rangle\Big{|}^{2}$ as a function of time $t$
over one period following the same quench as in Fig. S1 and Fig. S2. Here data
is presented for a number of system sizes, $L$, revealing significant finite
size effects in this nonlocal property. Each simulation was performed with
$\Delta t=0.005$ and $N_{c}=64$.

Finite size effects can more easily be revealed in nonlocal/global
measurements of the system. We illustrate this in Fig. S3, where we consider
the return amplitude $\Big{|}\langle\Psi_{0}|\Psi(t)\rangle\Big{|}^{2}$
following a quench. This rapidly decays as a function of time towards a close-
to-zero value. This close-to-zero value depends exponentially on the system
size. After one period of driving (the right hand side of the plot) we see
that the state $|\Psi(T)\rangle\approx|\Psi_{0}\rangle$, with the overlap
decreasing with system size.

We see that nonlocal observables experiencing severe finite size effects does
not carry through to local observables, cf. Figs. S2 and S3.

##  S6 Comparison of analytical results with finite volume numerics

The equations (S69),(S74),(S104), and (S109) for the one and two point-
functions are exact in the thermodynamic limit. However, to compare with our
numerical results, it is better to consider analytical expressions for finite
size systems. These can be easily obtained by recalling that the integrals in
the analytical expressions were obtained by taking the continuos limit of sums
over the allowed momenta. Therefore we only have to make the substitution
$\int dk\rightarrow(2\pi/L)\sum_{k_{n}}$ in the corresponding equations.
Concretely, we have, the following results for $t^{z}(t)$ and $t^{zz}(\ell,t)$

$\displaystyle t^{z}(t)$
$\displaystyle=\frac{1}{2}-\frac{1}{L}\sum_{n}\Big{[}\sin^{2}\left(\epsilon_{g_{f}}(k_{n})t\right)\sin^{2}\left(\Delta_{k_{n}}+\frac{\phi_{k_{n}}}{2}\right)+\sin^{2}\left(\frac{\phi_{k_{n}}}{2}\right)\cos^{2}\left(\epsilon_{g_{f}}(k_{n})t\right)\Big{]},$
$\displaystyle t^{zz}(\ell;t)$
$\displaystyle=-\frac{1}{L^{2}}\sum_{n}\sum_{n^{\prime}}\Bigg{\\{}e^{i(k_{n}-k_{n^{\prime}})\ell}\Big{[}\sin(2\epsilon_{g_{f}}(k_{n})t)\sin(\Delta_{k_{n}})\sin(2\epsilon_{g_{f}}(k_{n^{\prime}})t)\sin(\Delta_{k_{n^{\prime}}})+\Big{(}e^{i(k_{n}-k_{n^{\prime}})\ell}-1\Big{)}$
$\displaystyle\qquad\times
e^{i(\theta_{k_{n}}+\theta_{k_{n^{\prime}}})}\Big{(}\cos(\Delta_{k_{n}})-i\sin(\Delta_{k_{n}})\cos(2\epsilon_{g_{f}}(k_{n})t)\Big{)}\Big{(}\cos(\Delta_{k_{n^{\prime}}})-i\sin(\Delta_{k_{n^{\prime}}})\cos(2\epsilon_{g_{f}}(k_{n^{\prime}})t)\Big{)}\Bigg{\\}},$

where the set of allowed momenta is,

$k_{n}=\frac{2\pi n}{L},\quad n=0,\dots,L-1\,,$ (S120)

while the expressions for $t^{xx}(\ell,t)$ and $t^{yy}(\ell,t)$ can be
computed using (S104) and (S109), where $f_{\eta}$ and $g_{\eta}$ are given by
the following equations in a finite size system,

$\displaystyle f_{\eta}$
$\displaystyle=\frac{1}{L}\sum_{n}e^{ik_{n}\eta}f(k_{n}),$ (S121)
$\displaystyle g_{\eta}$
$\displaystyle=\frac{1}{L}\sum_{n}e^{ik_{n}\eta}g(k_{n}).$ (S122)