{'fraction_non_alphanumeric': 0.04533906134515891, 'fraction_numerical': 0.012541574279379157, 'mean_word_length': 4.412125, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum quench in a driven Ising chain

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 \textit{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.

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 timeevolution 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 [1][2][3][4][5][6][7][8][9][10][11][12], 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 [13,14] and the solid state [15,16]. 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. [17][18][19][20][21][22][23]) and other exotic behaviors (such as time crystal phases [24][25][26][27][28][29][30]). 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 [31][32][33]. Routes to avoid this behavior include introducing disorder to induce a many body localization transition (see, e.g., Refs. [24,[34][35][36]), or to consider models that are, in some sense, integrable [37].
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 at-tacking 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 onedimensional spin-1/2 Ising magnet, driven by oscillatory transverse fields. The Hamiltonian reads with J > 0 the Ising exchange parameter, h z a static longitudinal field, g the strength of the transverse fields, which oscillate at frequency Ω, and L the system size.
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 |Ψ 0 evolves under the Hamiltonian (1) at times t > 0. The time-evolved state |Ψ(t) will be a solution of the time-dependent Schrödinger equation subject to the initial condition |Ψ(t = 0) = |Ψ 0 . Herein we set = 1, which defines our units. The formal solution of Eq. (2) is well-known: however, using this to compute time-evolution is a challenge due to the explicit time-ordering (T) of the exponential. To make some headway on this problem we apply a time-dependent unitary transformation U (t) [40], multiplying both sides of Eq. (2) from the left by U (t) and inserting a factor of 1 = U (t) † U (t) between the wave function and the operators: The problem can become much simpler if there is a choice of U (t) such that this reduces to an effective timeindependent Schrödinger equation. Choosing [41] we map Eq. (2) to a time-independent Schrödinger equation, (i∂ t −H st )|Φ(t) = 0, with an effective static Hamiltonian The wave function transforms as |Φ(t) = U (t)|Ψ(t) . This reduction to a static problem is not evident in the Magnus expansion [42]. Diagonalizing (6) to obtain eigenstates |E n with energies E n , the time-evolved state can be written as The states |E n are not eigenstates of σ z 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: Mapping to a "sudden quench".-Let us now consider the time-evolution of one-point functions s α (t) = Ψ(t)|σ α l |Ψ(t) , where the result is independent of l by translational invariance. Using (8) these become 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 |Ψ 0 : (10) Equations similar to (9) can be written for the twopoint functions, s αβ ( ; t) = Ψ(t)|σ α j σ β j+ |Ψ(t) . These are tractable, but a little unwieldy, so are given in [42]. 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 − Ω/2 fields. The two fields can be independently controlled via the amplitude g and frequency Ω of the driving, see Eq. (1). Two interesting cases are immediately apparent. Firstly, if the frequency of the driving is tuned to a Ω = 2h z , the longitudinal field is removed from the static Hamiltonian, which then describes the integrable quantum Ising chain [43]. Secondly, one can consider tuning both the amplitude and the frequency such that g = J and |h z −Ω/2| g, where one realizes the lattice limit of the exotic E 8 perturbed Ising conformal field theory [44] (which has recently received renewed attention thanks to its nonthermal properties [45][46][47][48], 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 Ω = 2h z , the static Hamiltonian reads: This is the quantum Ising chain, which can be mapped to free fermions and so is exactly solvable [43]. This reveals that, along the line Ω = 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 [49][50][51].
We will exploit some of these results, alongside some new ones, to analytically compute the dynamics of observables starting from an initial state |Ψ 0 that is then timeevolved with the driven Hamiltonian (1). The derivation of these results is rather technical, so we provide the details in the Supplemental Material [42]. 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 [42]).
In Fig. 1 we present the results for one-and two-point functions for a particular quench. We see that the onepoint functions synchronize to the driving frequency Ω 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 Ω 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 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. [52] and Barouch and McCoy [53,54] in the 1970s. Detailed results in this case are presented in [42] 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 [31][32][33]). This is easily seen for observables that feature only σ 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 [12,55,56]. Generically, when H st is nonintegrable, this will be a finitetemperature Gibbs ensemble [57]. 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 [42], we show an explicit example of a non-integrable system with absence of heating to infinite temperature by working outside the integrable line Ω = 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 with J = 1, h z = 0, g = 2 and time-evolved with the driven Hamiltonian H(t) (1) with J = 1, Ω = 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.
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: The added term breaks σ z conservation, and thus evolves non-trivially under the transformation U (t). This breaks .001 (see [42] for the details of the numerical algorithm employed).
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. [58][59][60][61]. (The Floquet model studied in Ref. [61] bears some similarity to (12).) 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→∞ σ α j (t) = 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: When h = 0, H 1 st realizes the critical point of the Ising chain. For h = 0, H 1 st is no longer integrable, but its lowenergy physics is well-understood thanks to Zamolodchikov [44]. 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, The data was normalized so that the maximum value of the plot is one. The four dominant peaks at (±Ω, ±Ω) and (∓Ω, ±Ω) come from the driving frequency Ω, while the next three at pi = (mi − Ω, Ω − mi) come from the masses of the mesons excitations mi. Although m2/m1 is not equal to the golden ratio ϕ, we verify at the inset that, as the size of the system is increased, m2/m1 gets closer to ϕ. m 2 /m 1 = ϕ. With integrability absent, we are limit to performing small system numerics, such as in Fig. 3.
In Fig. 3 we plot the dynamical correlation function where σ x n (t) denotes the time-evolution of σ 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 , 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 Ω and are located at (ω 1 , ω 2 ) = (±Ω, ±Ω), (∓Ω, ±Ω). The remaining dominant peaks (marked as p 1 , p 2 and p 3 in the figure) correspond to the first excitations or masses of the static system, and their coordinates are where m i denotes the masses of the meson excitations. Although these masses do not satisfy the equality m 2 /m 1 = ϕ, we verify that this is a finite size effect in the inset of the figure, as m 2 /m 1 gets closer to ϕ 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 H of any spin-1/2 chain that conserves total σ z tot magnetization (this need not be translationally invariant), which is driven as in (1): The transformation (5) 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 whenH 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 [62,63]. Further exploring the lattice limit of Zamolodchikov's perturbed Ising field theory [44], 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 [64]), an avenue left to future works.
Acknowledgments In this Supplemental Material, we present: 1. The Magnus expansion for the problem in the main text.
2. The exact mapping between time-dependent twopoint functions of the driven systems and a sudden quench problem.
3. Exact full-time evolution of one-and two-point functions for a special case of the driven system quench.
4. Details of the computation of "sudden quench" correlation functions, summarizing both previously known results and presenting new ones.
5. Details of our numerical algorithm for simulating time-evolution in the driven system.
6. Details of the procedure used to compare finite-size numerical results with analytical expressions.

MAGNUS EXPANSION
An alternative approach to compute the time-evolution of the driven system is given by the Magnus expan-sion. For a time-dependent Hamiltonian H(t), the timeevolution operator U(t) is defined via the differential equation This is formally solved in terms of the well-known timeordered exponential However, the time-ordering here (denoted by the operator 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 timeevolution operator is instead approximated by where terms in the series are obtained from the differential equation Here ad Ξ is the adjoint action of Ξ: The terms of the series Ξ k (t) have well-known forms Here the first term simply corresponds to taking the timeaveraged 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π/Ω of the driving, 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: For the second order term in the expansion, we require the commutator where for conciseness we define the short hand notations c a = cos(Ωt a ), s a = sin(Ωt a ). Performing the double integral, we obtain Ξ 2 (T ): (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 σ x l σ y l+1 , σ y l σ 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.
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: Here, as in the main body of the manuscript, we define the "sudden quench" correlation functions where |E n are eigenstates of the static Hamiltonian H st with energy E n and |Ψ 0 = |Ψ(t = 0) is the initial state for the time-evolution. We note that Eq. (S10) is the same as time-evolving the initial state with H st alone.

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 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 (S12) with frequency Ω. We choose h z = Ω/2, in which case the effective static Hamiltonian is 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: for the one-point functions, while the diagonal two-point functions read: where we use notations, λ = J/g, Θ(x ≤ 0) = 0, Θ(x > 0) = 1, and we define the functions with L(n) = 1 π π 0 dk cos(kn) 1 + λ 2 + 2λ 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 αβ ( , t) = s βα ( , t). Using the Heisenberg equation of motion, s zy st ( , t) ∝ṡ zz st ( , t) = 0, as we are working in an static case. The terms s xz st ( , t) and s xy st ( , 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 λ < 1, using the symmetry of the Hamiltonian under π rotations of the quantization axis about the xdirection, after some tedious algebra, we find 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 init → H(t > 0) quench (S22) Furthermore, in the limit 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][53][54] 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: where E(m) =

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][50][51] in computing oneand 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.
where σ α l are the spin-1/2 operators acting on the lth 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 σ α L+1 = σ α 1 . To express the model in terms of free fermions, we first per-form a 3π/2 rotation of the quantization axes about the y axis to be consistent with standard notations Here τ α l are the Pauli matrices with the rotated quantization axes,J = J andḡ = g/J.
We proceed to fermionize the problem via the Jordan-Wigner transformation [69] where we define the raising/lowering operators τ ± = (τ x ± iτ y )/2 and the fermions have canonical anticommutation relations {a l , a † l } = δ l,l . Following this transformation, the Hamiltonian reads whereN = L l=1 a † 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πN ) = +1 for even numbers of fermions, −1 for odd numbers of fermions.
(S38) Thus the Hilbert space H = H e ⊗ H o of the model splits into two sectors, for N even or odd, which correspond to antiperiodic or periodic boundary conditions on the fermions. We need to address each of these cases in turn.
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. ḡ (k n ) > 0. As a result, the ground state |Ω e ,ḡ is the vacuum for the Bogoliubov fermions: α kn |Ω e ,ḡ = 0, ∀n. (S46) It immediately follows that the ground state energy is It will be useful to understand the ground state |Ω e ,ḡ in terms of the Jordan-Wigner fermions c kn . Inverting Eq. (S42) to express α kn in terms c kn , one finds that Eq. (S46) becomes (∀n) As |Ω e ,ḡ is in the even sector with zero momentum, it must be of the form kn≥0 (a + bc † −kn c † kn )|0 , where |0 is the vacuum for Jordan-Wigner fermions. Working through some algebra, one finds (S49) This can be written in a "squeezed state" form following some simple manipulations (S50)

Odd number of fermions
The story with the odd fermion number sector is similar to the even one. This time we have exp(iπ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 with the Bogoliubov angle as defined in the even sector, with the additional definition of θ p0 = 0. The Hamiltonian (S37) is then (S53) We see now that the ground state in the odd sector is where |0 α is the vacuum state for α pn . The ground state |Ω o has energy 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 As expected, this state has odd fermion parity.

Time-evolution following a quantum quench
Let us now turn to the problem of describing timeevolution 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ḡ = g i and at time t = 0 this is suddenly changed toḡ = 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 β by performing the Bogoliubov rotation with angle φ. For times t > 0 we can similarly diagonalize the Hamiltonian in terms of fermions α via the Bogoliubov rotation with angle θ. Equating the expressions for the Jordan-Wigner fermions a k in terms of the Bogoliubov fermions, we have The initial state, an eigenstate of H i will be expressed straightforwardly in terms of the β fermions. Time-evolution, however, will be easy to compute in terms of the α fermions via Thus we write the β fermions in terms of α ones and time-evolve them before re-expressing the α fermions in terms of the β ones: We see that the difference in Bogoliubov angles often appears, so we denote this by the short hand With this, and with some work, we can compute timeevolution of observables. In the following we first consider one-point functions, before continuing to two-point functions.

Time-evolution of one-point functions
We will consider the quench where we start from the ground state of the initial Hamiltonian H i ≡ H(t < 0) withḡ = g i and time-evolve according to the Hamiltonian H f whereḡ = g f . In the thermodynamic limit, the initial state will be a superposition of the ground states in the even and odd sectors (Here we work with normalized states such that Ω i ,ḡ|Ω i ,ḡ = 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 τ 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 Clearly t z (t) is independent of l by translational invariance. As τ z is number conserving in terms of the fermions, the terms mixing even and odd sectors from Eq. (S63) will vanish. Thus we have Let us take the first term, where the ground states is in the even sector of the Hilbert space. We express the timeevolved Jordan-Wigner fermions in terms of the timeevolved Bogoliubov fermions (S61). This yields Working through some tedious algebra, one finds 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 Let us now move on to considering τ x l and τ 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 (S70) and similarly for τ y . As such, this naturally brings us on to the topic of computing two-point functions.

Time-evolution of two-point functions
We now turn our attention to computing two-point functions Two examples of these, α = β = x and α = β = z, were considered in Refs. [49][50][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.
The operator under consideration is even under spin inversion, thus matrix elements between states in the even and odd sectors vanish, leaving only Using Eq. (S66) (and similar for the odd sector) and taking the thermodynamic limit, one finds The t xx ( ; 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 ( ; t) as the Jordan-Wigner transformation (S36) introduces "strings" of density operators between j and j + . This can easily be seen by Taylor expanding the exponential factors in Eq. (S36): Thus, in terms of the Jordan-Wigner fermions, we need to compute The calculation will proceed more easily if we introduce Majorana fermions, which describe the real and imaginary parts of the Jordan-Wigner fermions, which obey the anticommutation relations {γ a j , γ b l } = 2δ a,b δ j,l . In terms of these Majorana operators we have a † j a j = and hence (S81) This can be written as the Pfaffian of a 2 × 2 matrix where Γ is a 2 × 2 sub-matrix where Using reflection symmetry and translational invariance, we can rewrite this sub-matrix as where f η and g η can be obtained using Eq. (S66), and they read: When the two-point correlation function is written in terms of (S88), it is the Pfaffian of a block Toeplitz matrix: 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 × 2 antisymmetric matrix Λ = −Λ T has eigenvalues ±λ j (j = 1, . . . , ), the Pfaffian reads Pf(Λ) = j=1 iλ j . One can reduce computing the Pfaffian to computing the determinant of an × matrix Q 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][50][51], via applications of Szegö's theorem and Fisher-Hartwig conjectures [70][71][72][73][74][75][76].
The t yy ( ; t) two-point function The computation of this two-point function is very similar to the previous one. In terms of the Majorana fermions, (S96) The previous expression can be written as the Pfaffian of a block Toeplitz matrix, where f η and g η are defined in (S89)-(S91). The computation of the Pfaffian is equivalent to the one the determinant of an × matrixQ t yy ( ; t) = (−1) ( −1)/2 det(Q), (S99) (S100) The t xy ( ; t) two-point function This two point-function can be obtained from t xx ( ; t). Indeed, by considering the Heisenberg equation of mo-tion, (S101) By computing the expectation value and considering inversion symmetry we conclude, (S102)

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 Ω = 2h z . We make the following numerical studies outside this integrable line.
We consider the time-evolution of a given initial state |Ψ 0 induced by a time-dependent Hamiltonian, Eq. (1) of the main text. Formal solution of the time-dependent Schrödinger equation gives where 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 ∆t, such that N ∆t = t), each of which has a fixed Hamiltonian The problem now becomes how to describe the action of the time-evolution operator 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: (−i) n J n ∆t T n H(j∆t) . (S106) Here we have introduced the rescaled time ∆t and Hamiltonian H(j∆t), as well as the Bessel functions J n (x) and the Chebyshev matrices T n (x), which satisfy the recursion relation (S108) The rescaled time and Hamiltonian satisfy H(j∆t) ∆t = H(j∆t)∆t, and enforce that the eigenvalues of the H(j∆t) lie in the interval [−1, 1].
In practice, the Chebyshev expansion (S106) has to be truncated to finite order. In the following subsections we present some convergence checks of our "Trotteriza-tion+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.

Convergence with order of the Chebyshev expansion
Let us first examine how results of the time-evolution vary with the order of the Chebyshev expansion (S106). We fix the Trotter step size ∆t = 0.005 and then proceed to compute an approximation of the time-evolution: (−i) n J n ∆t T n H(j∆t) , (S109) 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 Ω = 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 . FIG. S1. The time-evolution of the one-point function (a) s x (t) = Ψ(t)|σ x 1 |Ψ(t) ; (b) s zz (1; t) = Ψ(t)|σ z j σ z j+1 |Ψ(t) 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 = Ω = 1, and g = 1.5. The time-evolution is computed via "Trotterization+Chebyshev" with Trotter step ∆t = 0.005 and Chebyshev expansion order Nc for ten periods of the drive, T = 2π/Ω. We see that results are rapidly converging with increased Nc. Note that there is no heating to infinite temperature, despite the fact we are outside the integrable line Ω = 2h z .

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

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.  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.
Fixing ∆t = 0.005 and N c = 64, we examine the nonequilibrium time-evolution of σ x following the same quench as in the previous subsections.
We first examine a local observable, Ψ(t)|σ x |Ψ(t) , 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.
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 Ψ 0 |Ψ(t) 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 |Ψ(T ) ≈ |Ψ 0 , with the overlap decreasing with system size. 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 ∆t = 0.005 and Nc = 64.
We see that nonlocal observables experiencing severe finite size effects does not carry through to local observables, cf. Figs. S2 and S3.

COMPARISON OF ANALYTICAL RESULTS
WITH FINITE VOLUME NUMERICS The equations (S69),(S74),(S94), and (S99) 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 dk → (2π/L) kn in the corresponding equations. Concretely, we have, the following results for t z (t) and t zz ( , t) t z (t) = 1 2 − 1 L n sin 2 g f (k n )t sin 2 ∆ kn + φ kn 2 + sin 2 φ kn 2 cos 2 g f (k n )t , t zz ( ; t) = − 1 L 2 n n e i(kn−k n ) sin(2 g f (k n )t) sin(∆ kn ) sin(2 g f (k n )t) sin(∆ k n ) + e i(kn−k n ) − 1 × e i(θ kn +θ k n ) cos(∆ kn ) − i sin(∆ kn ) cos(2 g f (k n )t) cos(∆ k n ) − i sin(∆ k n ) cos(2 g f (k n )t) , where the set of allowed momenta is, while the expressions for t xx ( , t) and t yy ( , t) can be computed using (S94) and (S99), where f η and g η are given by the following equations in a finite size system, g η = 1 L n e iknη g(k n ). (S112)