Quantum quench in a driven Ising chain
Neil J. Robinson, Isaac Pérez Castillo, Edgar Guzmán-González
QQuantum quench in a driven Ising chain
Neil J. Robinson, ∗ Isaac Pérez Castillo, † and Edgar Guzmán-González ‡ Institute for Theoretical Physics, University of Amsterdam,Postbus 94485, 1090 GL Amsterdam, The Netherlands § Departamento de Física, Universidad Autónoma Metropolitana-Iztapalapa,San Rafael Atlixco 186, Ciudad de México 09340, Mexico London Mathematical Laboratory, 8 Margravine Gardens, London W6 8RH, United Kingdom (Dated: December 1, 2020)We consider the Ising chain driven by oscillatory transverse magnetic fields. For certain parameterregimes, 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 functionsfollowing a quench that activates the driving. It is shown that this model does not heat up to infinitetemperature, despite the absence of energy conservation, and we further discuss the generalizationto a family of driven Hamiltonians that do not suffer heating to infinite temperature, despite theabsence of integrability and disorder. The particular model studied in detail also presents a routefor realising exotic physics (such as the E8 perturbed conformal field theory) via driving in quantumchains that could otherwise never realise such behaviour. In particular we numerically confirm thatthe ratio of the meson excitations masses is given by the golden ratio.
Introduction. —Over the last decade, the nonequilib-rium dynamics of quantum systems has attracted a greatdeal of attention [1–12], motivated by the desire to ad-dress fundamental questions: When and how do quantumsystems relax to equilibrium? How does one describe thisequilibrium? What influences the dynamics and equili-bration? Understanding these issues is important whendeveloping descriptions of a growing number of experi-ments that examine nonequilibrium dynamics, both incold atomic gases [13, 14] and the solid state [15, 16].The insights gained may play an important role in thedevelopment of quantum computing resources, especiallywhen considering how to protect quantum informationfrom the scrambling associated with thermalization.Recently attention has turned to understanding drivenquantum systems, partially due to the realization thatsuch systems can host interesting topological phases (see,e.g., Refs. [17–23]) and other exotic behaviors (such astime crystal phases [24–30]). These studies have gener-ated much discussion of how to extend and apply the con-cepts of equilibrium statistical mechanics in the presenceof driving. A particular issue is that, generically, drivenquantum systems do not conserve energy. As a result,in the long time limit entropy maximization leads themto heat up to infinite temperature, leading to trivial er-godic behavior. As a result, quantum information is com-pletely scrambled [31–33]. Routes to avoid this behaviorinclude introducing disorder to induce a many body lo-calization transition (see, e.g., Refs. [24, 34–36]), or toconsider models that are, in some sense, integrable [37].In this Letter, we consider a driven model that at eachpoint in time is nonintegrable but nonetheless possessesthe dynamics which is governed by a hidden integrabil-ity. Using this, we compute the nonequilibrium dynam-ics of equal-time correlation functions following a quenchin 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 non-perturbative, exact studies. We will see that this wholefamily of driven quantum systems, each of which is gener-ically nonintegrable, does not undergo heating to infinitetemperature. We will also see that breaking the specialstructure of this family leads to thermalization to infinitetemperature.
The driven Ising chain. —We consider a one-dimensional spin-1/2 Ising magnet, driven by oscillatorytransverse fields. The Hamiltonian reads H ( t ) = − J L (cid:88) l =1 σ zl σ zl +1 + h z L (cid:88) l =1 σ zl − g L (cid:88) l =1 (cid:16) e − i Ω t σ + l + e i Ω t σ − l (cid:17) , (1)with J > the Ising exchange parameter, h z a staticlongitudinal field, g the strength of the transverse fields,which oscillate at frequency Ω , and L the system size.The spin operators σ αl act at the l th site of the lattice, σ ± l = ( σ xl ± iσ yl ) / , and we impose periodic boundaryconditions σ αL +1 = σ α . The Hamiltonian (1) is periodic intime H ( t ) = H ( t + T ) with period T = 2 π/ Ω , and couldbe realized in the quasi-1D ferromagnet CoNb O [38, 39]by application of oscillating transverse fields.At a generic time, the Hamiltonian consists of an Isinginteraction term and fields in all ( x, y, z ) directions. Thusinstantaneously 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 modelwhere exact results can be obtained. Furthermore, awayfrom this integrability we will draw general insights. Time evolution of observables. —We will now considerhow a state | Ψ (cid:105) evolves under the Hamiltonian (1) at a r X i v : . [ c ond - m a t . o t h e r] N ov times t > . The time-evolved state | Ψ( t ) (cid:105) will be asolution of the time-dependent Schrödinger equation (cid:104) i (cid:126) ∂ t − H ( t ) (cid:105) | Ψ( t ) (cid:105) = 0 , (2)subject to the initial condition | Ψ( t = 0) (cid:105) = | Ψ (cid:105) . Hereinwe set (cid:126) = 1 , which defines our units. The formal solutionof Eq. (2) is well-known: | Ψ( t ) (cid:105) = T exp (cid:18) − i (cid:90) t d t (cid:48) H ( t (cid:48) ) (cid:19) | Ψ (cid:105) , (3)however, using this to compute time-evolution is a chal-lenge due to the explicit time-ordering ( T ) of the expo-nential. To make some headway on this problem we ap-ply 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 = U ( t ) † U ( t ) between the wavefunction and the operators: U ( t ) (cid:104) i∂ t − H ( t ) (cid:105) U † ( t ) U ( t ) | Ψ( t ) (cid:105) = 0 . (4)The problem can become much simpler if there is achoice of U ( t ) such that this reduces to an effective time-independent Schrödinger equation. Choosing [41] U ( t ) = exp (cid:16) i Ω t (cid:88) l σ zl (cid:17) ≡ e i Ω t σ z tot , (5)we map Eq. (2) to a time-independent Schrödinger equa-tion, ( i∂ t − H st ) | Φ( t ) (cid:105) = 0 , with an effective static Hamil-tonian H st = L (cid:88) l =1 (cid:20) − Jσ zl σ zl +1 + (cid:18) h z − Ω2 (cid:19) σ zl − gσ xl (cid:21) . (6)The wave function transforms as | Φ( t ) (cid:105) = U ( t ) | Ψ( t ) (cid:105) .This reduction to a static problem is not evident in theMagnus expansion [42].Diagonalizing (6) to obtain eigenstates | E n (cid:105) with en-ergies E n , the time-evolved state can be written as | Ψ( t ) (cid:105) = (cid:88) n exp (cid:104) − i (cid:18) E n + Ω2 σ z tot (cid:19) t (cid:105) | E n (cid:105)(cid:104) E n | Ψ (cid:105) . (7)The states | E n (cid:105) are not eigenstates of σ z tot and thus eachterm in Eq. (7) undergoes nontrivial dynamics. While (7)is highly nontrivial, there is no need to despair. Our prob-lem reduces to a tractable one if we focus on equal-timecorrelation functions, as one can use that the operator U ( t ) acts in a simple manner on the spin operators: U ( t ) σ x,yl U ( t ) † = cos(Ω t ) σ x,yl ∓ sin(Ω t ) σ y,xl ,U ( t ) σ zl U ( t ) † = σ zl . (8) Mapping to a “sudden quench”.—
Let us now con-sider the time-evolution of one-point functions s α ( t ) = (cid:104) Ψ( t ) | σ αl | Ψ( t ) (cid:105) , where the result is independent of l bytranslational invariance. Using (8) these become s z ( t ) = s z st ( t ) ,s x ( t ) = cos(Ω t ) s x st ( t ) − sin(Ω t ) s y st ( t ) ,s y ( t ) = cos(Ω t ) s y st ( t ) + sin(Ω t ) s x st ( t ) . (9)Here each time-dependent expectation value on the righthand side describes time-evolution induced by a suddenquench to the static Hamiltonian (6) when starting fromthe initial state | Ψ (cid:105) : s α st ( t ) = (cid:88) n,m e i ( E n − E m ) t (cid:104) Ψ | E n (cid:105)(cid:104) E n | σ αl | E m (cid:105)(cid:104) E m | Ψ (cid:105) . (10)Equations similar to (9) can be written for the two-point functions, s αβ ( (cid:96) ; t ) = (cid:104) Ψ( t ) | σ αj σ βj + (cid:96) | Ψ( t ) (cid:105) . Theseare tractable, but a little unwieldy, so are given in [42].All time-evolved correlation functions are reduced to os-cillatory factors multiplying “sudden quench” correlationfunctions. Thus for this driven problem, we can applythe techniques developed for sudden quantum quenchesto compute the time-evolution of observables.Having reduced the problem from one with driving toan effective sudden quench, let us return to the staticHamiltonian (6). This describes a quantum Ising chainwith both transverse g and longitudinal h = h z − Ω / fields. The two fields can be independently controlledvia the amplitude g and frequency Ω of the driving, seeEq. (1). Two interesting cases are immediately appar-ent. Firstly, if the frequency of the driving is tuned to a Ω = 2 h z , the longitudinal field is removed from the staticHamiltonian, which then describes the integrable quan-tum Ising chain [43]. Secondly, one can consider tuningboth the amplitude and the frequency such that g = J and | h z − Ω / | (cid:28) g , where one realizes the lattice limit ofthe exotic E perturbed Ising conformal field theory [44](which has recently received renewed attention thanks toits nonthermal properties [45–48], despite an absence ofintegrability). In this work, we will focus on the firstscenario and describe the full time-evolution of one- andtwo-point functions in this driven problem. We will touchupon the second case towards the end.When Ω = 2 h z , the static Hamiltonian reads: H st = − J L (cid:88) l =1 σ zl σ zl +1 − g L (cid:88) l =1 σ xl . (11)This is the quantum Ising chain, which can be mapped tofree fermions and so is exactly solvable [43]. This revealsthat, along the line Ω = 2 h z , there is a hidden integrabil-ity in the problem (despite, instantaneously, the Hamilto-nian H ( t ) being nonintegrable). Sudden quenches in thetransverse field Ising model have been extensively stud-ied, with many exact results being known, see in partic-ular the works of Calabrese, Essler and Fagotti [49–51].We will exploit some of these results, alongside some newones, to analytically compute the dynamics of observ-ables starting from an initial state | Ψ (cid:105) that is then time-evolved with the driven Hamiltonian (1). The derivationof these results is rather technical, so we provide the de-tails in the Supplemental Material [42]. Time-evolution in the driven model.—
Let us nowpresent the time-evolution of correlation functions in thedriven model (1) governed by the effective static Hamilto-nian (11). We compare our analytical results to numeri-cal results obtained on small finite lattices (our numericalalgorithm is explained in [42]).In Fig. 1 we present the results for one- and two-pointfunctions for a particular quench. We see that the one-point functions synchronize to the driving frequency Ω and no heating to infinite temperature occurs, not evenif we restrict the study of the system to stroboscopictimes. For the two-point functions, we see that s zz (1 , t ) converges to a non-zero stationary value, confirming theabsence of infinite heating. Although not included in thefigure, we mention that the remaining two-point func-tions synchronize to the period Ω like the one-point func-tions, except for s yz (1 , t ) and s xz (1 , t ) that converge tozero.A particularly simple, and solvable in closed form,scenario is realized when H st coincides with the initialHamiltonian H ( t < . In this case the “sudden quench”correlation functions in expressions such as (9) reduce to equilibrium correlation functions , known since the semi-nal works of Barouch et al. [52] and Barouch and McCoy[53, 54] in the 1970s. Detailed results in this case are pre-sented in [42] and are, to our knowledge, some of the fewclosed form exact results known for correlation functionsin models with driving. Absence of heating to infinite temperature.—
With ob-servables mapping in a simple manner to those from asudden quench, it is immediately clear that the systemcannot undergo heating to infinite temperature, as is usu-ally assumed to occur in driven systems [31–33]). This iseasily seen for observables that feature only σ zl operators,which map exactly to “sudden quench” observables (see,e.g., the first line of Eq. (9)). The long-time limit of ob-servables after a sudden quench will be described via therelevant statistical ensemble; for the case detailed abovethis is the generalized Gibbs ensemble [12, 55, 56]. Gener-ically, when H st is nonintegrable, this will be a finite-temperature Gibbs ensemble [57]. It is worth noting thatthe absence of heating to infinite temperature is not asresult of integrability, but instead is due to the structureof the driving term. In Fig. S1 of [42], we show an ex-plicit example of a non-integrable system with absenceof heating to infinite temperature by working outside theintegrable line Ω = 2 h z .We can then ask, what happens if this structure isbroken such that we do not map to an effective suddenquench problem? We then expect that in the long time (a) − − . .
51 0 2 4 6 8 10Time, t/Ts x ( t ) s y ( t ) s z ( t ) (b) . . . . . . . . .
42 0 2 4 6 8 10 s zz ( , t ) Time, t/TL = 16 L = 150 FIG. 1. (a) Time-evolution of one-point functions startingfrom the ground state of the quantum Ising chain, H ( t = 0) with J = 1 , h z = 0 , g = 2 and time-evolved with the drivenHamiltonian H ( t ) (1) with J = 1 , Ω = 1 , h z = 0 . , g = 1 . for a system with L = 16 sites. The behaviour remains thesame for bigger system sizes and larger times. Lines repre-sent analytical results, while points show numerically exacttime-evolution. (b) Time evolution of the two-point function s zz (1 , t ) of two adjacent sites for the same quench, for dif-ferent system sizes L . We see that s zz (1 , t ) converges to anon-zero value. The revival of the fluctuations is a finite sizeeffect, as can be seen by increasing the system size. limit the system thermalizes to infinite temperature, dueto the absence of both energy conservation and the map-ping to a sudden quench problem, combined with en-tropy maximization. We can examine this numericallyby adding terms to our Hamiltonian (1), for example: H X ( t ) = H ( t ) + J X L (cid:88) l =1 σ xl σ xl +1 . (12)The added term breaks σ z conservation, and thus evolvesnon-trivially under the transformation U ( t ) . This breaks − − . .
51 0 5 10 15 20 25 30 35 40 45Time, t/Ts x ( t ) s y ( t ) s z ( t ) FIG. 2. Numerically exact time-evolution of one-point func-tions with the driven Hamiltonian H X ( t ) , Eq. (12). Thisshows that, at the level of one-point functions, breaking thestructure of the drive leads to thermalization to infinite tem-perature, as all the expectation values converge to zero. Theparameters considered were the ones of Fig. (1) with J x = 0 . ,and a Chebyshev expansion of order 64 with a time step ∆ t = the mapping to a static Hamiltonian, and hence we ex-pect heating to infinite temperature. It is worth notingthat the thermalization time scale in Floquet systems canbe very large, see e.g. Refs. [58–61]. (The Floquet modelstudied in Ref. [61] bears some similarity to (12).)In Fig. (2) we present the time-evolution of one-pointfunctions in the driven model (12). With the addition ofthe J X term, we see that the system evolves towards astate with lim t →∞ (cid:104) σ αj ( t ) (cid:105) = 0 , corresponding to infinitetemperature, at least at the level of one-point subsys-tems. Realizing a perturbed critical model.—
Let us finish withan illustration of the second interesting case discussedabove. We consider tuning the driving such that thestatic Hamiltonian describes the perturbed critical Isingchain: H = − J L (cid:88) l =1 σ zl σ zl +1 − J L (cid:88) l =1 σ xl + h L (cid:88) l =1 σ zl . (13)When h = 0 , H realizes the critical point of the Isingchain. For h (cid:54) = 0 , H is no longer integrable, but its low-energy physics is well-understood thanks to Zamolod-chikov [44]. Pairs of fermions (corresponding to domainwalls in the ordered phase) are confined by the presenceof the longitudinal field h and form “meson” excitations.In the scaling limit, the algebraic structure of the the-ory allows the prediction of these meson masses, includ-ing the beautiful result that the ratio of masses of thefirst and second meson states realizes the golden ratio, FIG. 3. Dynamical correlation function s xx ( k = 0 , ω , ω ) (13) for the perturbed critical model with J = 1 , Ω = 1 , h z = Ω / . , g = 1 . and L = 16 . The data was normalizedso that the maximum value of the plot is one. The four domi-nant peaks at ( ± Ω , ± Ω) and ( ∓ Ω , ± Ω) come from the drivingfrequency Ω , while the next three at p i = ( m i − Ω , Ω − m i ) come from the masses of the mesons excitations m i . Although m /m is not equal to the golden ratio ϕ , we verify at the in-set that, as the size of the system is increased, m /m getscloser to ϕ . m /m = ϕ . With integrability absent, we are limit toperforming small system numerics, such as in Fig. 3.In Fig. 3 we plot the dynamical correlation function s xx ( k = 0 , ω , ω ) = (cid:88) (cid:96) (cid:90) dt dt e i ( ω t + ω t ) × (cid:104) Ψ | σ xj + (cid:96) ( t ) σ xj ( t ) | Ψ (cid:105) , (14)where σ xn ( t ) denotes the time-evolution of σ xn in theHeisenberg picture and, for simplicity, we assume theinitial state of the system was prepared to be theground state of the static Hamiltonian (13). Note that,because of the driving, energy is not conserved, and (cid:104) Ψ | σ xj + (cid:96) ( t ) σ xj ( t ) | Ψ (cid:105) is no longer a function of the timedifference ( t − t ) , therefore we considered the Fouriertransform of both times.Note that there are four dominant peaks in Fig. 3,these correspond to the driving frequency Ω and arelocated at ( ω , ω ) = ( ± Ω , ± Ω) , ( ∓ Ω , ± Ω) . The re-maining dominant peaks (marked as p , p and p inthe figure) correspond to the first excitations or massesof the static system, and their coordinates are p i =( m i − Ω , Ω − m i ) , where m i denotes the masses of themeson excitations. Although these masses do not satisfythe equality m /m = ϕ , we verify that this is a finitesize effect in the inset of the figure, as m /m gets closerto ϕ when the size of the system is increased. This reiter-ates the fact that driven systems, such as (1), can realiseexotic physics that is inaccessible in its undriven state. Discussion.—
In this Letter, we have explored an ex-ample of driven system that is instantaneously nonin-tegrable but can nonetheless be solved exactly. This isdue to a hidden integrability in the problem, that is notapparent from the time-dependent Schrödinger equation:the instantaneous Hamiltonian H ( t ) is nonintegrable, butdynamics of observables are nonetheless controlled by aneffective static, integrable Hamiltonian. This may pro-vide a route to protecting quantum information from thescrambling associated with thermalization through theaddition of driving; this is an interesting direction forfuture studies.The methods applied within this Letter can be used totackle the dynamics of an infinite family of Hamiltonians(not necessarily integrable). This family of driven sys-tems does not undergo heating to infinite temperature,even though they are absent disorder and (generically)integrability. For example, consider the Hamiltonian (cid:101) H of any spin-1/2 chain that conserves total σ z tot magneti-zation (this need not be translationally invariant), whichis driven as in (1): (cid:101) H ( t ) = (cid:101) H − (cid:101) g (cid:88) l (cid:16) e − i (cid:101) Ω t σ + l + H.c. (cid:17) . (15)The transformation (5) still maps Eq. (2) to a time-independent Schrödinger equation with the new effectivestatic Hamiltonian (cid:101) H st = (cid:101) H − ( (cid:101) Ω / (cid:80) l σ zl − (cid:101) g (cid:80) l σ xl .Time-evolution of observables in a driven system hasonce again been mapped to a sudden quench problem.It would be interesting to explore this idea further in in-teracting models, such as when ˜ H st describes the Heisen-berg or XXZ model, where potentially integrability canbe harnessed to perform exact calculations.Another scenario worthy of attention is to consider aproblem in which the parameters of the static Hamil-tonian describe a different phase to the initial Hamilto-nian. One may then expect to see signatures of dynamicalphase transitions in the nonequilibrium dynamics, suchas kinks in the Lochsmidt echo [62, 63]. Further exploringthe lattice limit of Zamolodchikov’s perturbed Ising fieldtheory [44], which features interesting collective excita-tions related to an exotic hidden E algebraic structure,is interesting. Such studies would require detailed nu-merical analysis (perhaps in the scaling limit [64]), anavenue left to future works. Acknowledgments.—
Our thanks go to Jean-SébastienCaux, Mario Collura, Andrew James, Robert Konik,Tamás Pálmai, and Sergio Tapias Arze for useful dis-cussions. I.P.C. and E.G.-G. thanks the London Math-ematical Laboratory for financial support. N.J.R. wassupported by funding from the EU’s Horizon 2020 re-search and innovation programme, under grant agree-ment No 745944, and the European Research Councilunder ERC Advanced Grant No 743032 (
Dynamint ).E.G.-G. acknowledges the hospitality of the Institute of Physics, University of Amsterdam, during the completionof this work. N.J.R. thanks the Institute of Physics of theNational Autonomous University of Mexico for hospital-ity during a visit, where part of this work was undertaken. ∗ [email protected] † [email protected] ‡ [email protected] § Present address: UKRI EPSRC, Polaris House, NorthStar Avenue, Swindon SN2 1ET, United Kingdom[1] C. Gogolin and J. Eisert, “Equilibration, thermalisa-tion, and the emergence of statistical mechanics in closedquantum systems,” Rep. Prog. Phys. , 056001 (2016).[2] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol,“From quantum chaos and eigenstate thermalization tostatistical mechanics and thermodynamics,” Adv. Phys. , 239–362 (2016).[3] P. Calabrese and J. Cardy, “Quantum quenches in 1+1dimensional conformal field theories,” J. Stat. Mech. , P064003 (2016).[4] M. A. Cazalilla and M.-C. Chung, “Quantum quenchesin the Luttinger model and its close relatives,” J. Stat.Mech. , P064004 (2016).[5] J.-S. Caux, “The Quench Action,” J. Stat. Mech. ,P064006 (2016).[6] F. H. L. Essler and M. Fagotti, “Quench dynamics andrelaxation in isolated integrable quantum spin chains,” J.Stat. Mech. , P064002 (2016).[7] D. Bernard and B. Doyon, “Conformal field theory outof equilibrium: a review,” J. Stat. Mech. , P064005(2016).[8] E. Ilievski, M. Medenjak, T. Prosen, and L. Zadnik,“Quasilocal charges in integrable lattice systems,” J. Stat.Mech. , P064008 (2016).[9] R. Vasseur and J. E. Moore, “Nonequilibrium quantumdynamics and transport: from integrability to many-body localization,” J. Stat. Mech. , P064010 (2016).[10] A. De Luca and G. Mussardo, “Equilibration propertiesof classical integrable field theories,” J. Stat. Mech. ,P064011 (2016).[11] T. Langen, T. Gasenzer, and J. Schmiedmayer, “Prether-malization and universal dynamics in near-integrablequantum systems,” J. Stat. Mech. , P064009 (2016).[12] L. Vidmar and M. Rigol, “Generalized Gibbs ensemble inintegrable lattice models,” J. Stat. Mech. , P064007(2016).[13] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat-tore, “ Colloquium:
Nonequilibrium dynamics of closedinteracting quantum systems,” Rev. Mod. Phys. , 863–883 (2011).[14] T. Langen, R. Geiger, and J. Schmiedmayer, “UltracoldAtoms Out of Equilibrium,” Annu. Rev. Cond. MatterPhys. , 201–217 (2015).[15] U. Bovensiepen and P. S. Kirchmann, “Elementary relax-ation processes investigated by femtosecond photoelec-tron spectroscopy of two-dimensional materials,” LaserPhoton. Rev. , 589–606 (2012).[16] M. Beye, Ph. Wernet, C. Schüßler-Langeheine, andA. Föhlisch, “Time resolved resonant inelastic X-ray scat-tering: A supreme tool to understand dynamics in solids and molecules,” J. Electron. Spectrosc. Relat. Phenom. , 172 – 182 (2013).[17] C. W. von Keyserlingk and S. L. Sondhi, “Phase struc-ture of one-dimensional interacting Floquet systems. I.Abelian symmetry-protected topological phases,” Phys.Rev. B , 245145 (2016).[18] C. W. von Keyserlingk and S. L. Sondhi, “Phase struc-ture of one-dimensional interacting Floquet systems. II.Symmetry-broken phases,” Phys. Rev. B , 245146(2016).[19] D. V. Else and C. Nayak, “Classification of topologicalphases in periodically driven interacting systems,” Phys.Rev. B , 201103 (2016).[20] A. C. Potter, T. Morimoto, and A. Vishwanath, “Classi-fication of Interacting Topological Floquet Phases in OneDimension,” Phys. Rev. X , 041001 (2016).[21] R. Roy and F. Harper, “Abelian Floquet symmetry-protected topological phases in one dimension,” Phys.Rev. B , 125105 (2016).[22] R. Roy and F. Harper, “Periodic table for Floquet topo-logical insulators,” Phys. Rev. B , 155118 (2017).[23] Rahul Roy and Fenner Harper, “Floquet topologicalphases with symmetry in all dimensions,” Phys. Rev. B , 195128 (2017).[24] V. Khemani, A. Lazarides, R. Moessner, and S. L.Sondhi, “Phase Structure of Driven Quantum Systems,”Phys. Rev. Lett. , 250401 (2016).[25] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi,“Absolute stability and spatiotemporal long-range orderin Floquet systems,” Phys. Rev. B , 085112 (2016).[26] 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, “Observa-tion of a discrete time crystal,” Nature (London) ,217–220 (2017).[27] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya,F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. vonKeyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin,“Observation of discrete time-crystalline order in a disor-dered dipolar many-body system,” Nature (London) ,221–225 (2017).[28] Vedika Khemani, C. W. von Keyserlingk, and S. L.Sondhi, “Defining time crystals via representation the-ory,” Phys. Rev. B , 115127 (2017).[29] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vish-wanath, “Discrete Time Crystals: Rigidity, Criticality,and Realizations,” Phys. Rev. Lett. , 030401 (2017).[30] Vedika Khemani, Roderich Moessner, and S. L.Sondhi, “A Brief History of Time Crystals,” arXive-prints , arXiv:1910.10745 (2019), arXiv:1910.10745[cond-mat.str-el].[31] L. D’Alessio and M. Rigol, “Long-time Behavior of Iso-lated Periodically Driven Interacting Lattice Systems,”Phys. Rev. X , 041048 (2014).[32] A. Lazarides, A. Das, and R. Moessner, “Equilibriumstates of generic quantum systems subject to periodicdriving,” Phys. Rev. E , 012110 (2014).[33] P. Ponte, A. Chandran, Z. Papić, and A. A. Abanin, “Pe-riodically driven ergodic and many-body localized quan-tum systems,” Ann. Phys. (N.Y.) , 196 – 204 (2015).[34] P. Ponte, Z. Papić, F. Huveneers, and D. A. Abanin,“Many-Body Localization in Periodically Driven Sys-tems,” Phys. Rev. Lett. , 140401 (2015).[35] A. Lazarides, A. Das, and R. Moessner, “Fate of Many- Body Localization Under Periodic Driving,” Phys. Rev.Lett. , 030402 (2015).[36] D. A. Abanin, W. De Roeck, and F. Huveneers, “The-ory of many-body localization in periodically driven sys-tems,” Ann. Phys. (N.Y.) , 1–11 (2016).[37] V. Gritsev and A. Polkovnikov, “Integrable Floquet dy-namics,” SciPost Phys. , 021 (2017).[38] R. Coldea, D. A. Tennant, E. M. Wheeler, E. Wawrzyn-ska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl,and K. Kiefer, “Quantum Criticality in an Ising Chain:Experimental Evidence for Emergent E Symmetry,” Sci-ence , 177 (2010).[39] Neil J. Robinson, Fabian H. L. Essler, Ivelisse Cabr-era, and Radu Coldea, “Quasiparticle breakdown inthe quasi-one-dimensional Ising ferromagnet
CoNb O ,”Phys. Rev. B , 174406 (2014).[40] As discussed in Refs. [65–67], there is nice geometric in-terpretation of unitary transformations that depend ona continuous parameter (e.g., time t ) in terms of gaugepotentials.[41] S. Takayoshi, H. Aoki, and T. Oka, “Magnetization andphase transition induced by circularly polarized laser inquantum magnets,” Phys. Rev. B , 085150 (2014).[42] See the Supplemental Material, which also contains thereferences [68–76], for: (i) A discussion of the MagnusExpansion for the driven problem considered; (ii) de-tails of how two-point correlation functions transform un-der action of the time-dependent unitary transformation;(iii) a discussion of the full-time evolution of observablesfor a special quench; (iv) detailed derivations of the re-quired “sudden quench” correlation functions; (v) detailsof the numerical algorithm for computing time-evolutionof the driven system and a numerical check of the absenceof heating to infinite temperature outside the integrableline.[43] Pierre Pfeuty, “The one-dimensional Ising model with atransverse field,” Ann. Phys. (N.Y.) , 79 – 90 (1970).[44] A. B. Zamolodchikov, “Integrals of motion and S-matrixof the (scaled) T = T c Ising model with magnetic field,”Int. J. Mod. Phys. A , 4235–4248 (1989).[45] T. Rakovszky, M. Mestyán, M. Collura, M. Kormos,and G. Takács, “Hamiltonian truncation approach toquenches in the Ising field theory,” Nucl. Phys. B ,805 – 845 (2016).[46] Kristóf Hódsági, Márton Kormos, and Gábor Takács,“Quench dynamics of the Ising field theory in a magneticfield,” SciPost Phys. , 27 (2018).[47] Andrew J. A. James, Robert M. Konik, and Neil J.Robinson, “Nonthermal States Arising from Confinementin One and Two Dimensions,” Phys. Rev. Lett. ,130603 (2019).[48] Neil J. Robinson, Andrew J. A. James, and Robert M.Konik, “Signatures of rare states and thermalization ina theory with confinement,” Phys. Rev. B , 195108(2019).[49] P. Calabrese, F. H. L. Essler, and M. Fagotti, “QuantumQuench in the Transverse-Field Ising Chain,” Phys. Rev.Lett. , 227203 (2011).[50] P. Calabrese, F. H. L. Essler, and M. Fagotti, “Quan-tum quench in the transverse field Ising chain: I. Timeevolution of order parameter correlators,” J. Stat. Mech. , P07016 (2012).[51] P. Calabrese, F. H. L. Essler, and M. Fagotti, “Quantumquenches in the transverse field Ising chain: II. Stationary state properties,” J. Stat. Mech. , P07022 (2012).[52] E. Barouch, B. M. McCoy, and M. Dresden, “StatisticalMechanics of the XY Model. I,” Phys. Rev. A , 1075–1092 (1970).[53] E. Barouch and B. M. McCoy, “Statistical Mechanicsof the XY Model. II. Spin-Correlation Functions,” Phys.Rev. A , 786–804 (1971).[54] E. Barouch and B. M. McCoy, “Statistical Mechanics ofthe XY Model. III,” Phys. Rev. A , 2137–2140 (1971).[55] Marcos Rigol, Vanja Dunjko, Vladimir Yurovsky, andMaxim Olshanii, “Relaxation in a Completely IntegrableMany-Body Quantum System: An Ab Initio Study ofthe Dynamics of the Highly Excited States of 1D LatticeHard-Core Bosons,” Phys. Rev. Lett. , 050405 (2007).[56] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L.Essler, and T. Prosen, “Complete Generalized Gibbs En-sembles in an Interacting Theory,” Phys. Rev. Lett. ,157201 (2015).[57] Marcos Rigol, Vanja Dunjko, and Maxim Olshanii,“Thermalization and its mechanism for generic isolatedquantum systems,” Nature (London) , 854 (2008).[58] Dmitry A. Abanin, Wojciech De Roeck, and Fran çoisHuveneers, “Exponentially Slow Heating in PeriodicallyDriven Many-Body Systems,” Phys. Rev. Lett. ,256803 (2015).[59] Dmitry Abanin, Wojciech De Roeck, Wen Wei Ho, andFrançois Huveneers, “A Rigorous Theory of Many-BodyPrethermalization for Periodically Driven and ClosedQuantum Systems,” Commun. Math. Phys. , 809–827(2017).[60] Dominic V. Else, Bela Bauer, and Chetan Nayak,“Prethermal Phases of Matter Protected by Time-Translation Symmetry,” Phys. Rev. X , 011026 (2017).[61] Francisco Machado, Gregory D. Meyer, Dominic V. Else,Chetan Nayak, and Norman Y. Yao, “ExponentiallySlow Heating in Short and Long-range Interacting Flo-quet Systems,” arXiv e-prints , arXiv:1708.01620 (2017),arXiv:1708.01620 [quant-ph].[62] M. Heyl, A. Polkovnikov, and S. Kehrein, “DynamicalQuantum Phase Transitions in the Transverse-Field IsingModel,” Phys. Rev. Lett. , 135704 (2013).[63] C. Karrasch and D. Schuricht, “Dynamical phase transi-tions after quenches in nonintegrable models,” Phys. Rev.B , 195104 (2013). [64] Andrew J. A. James, Robert M. Konik, Philippe Lechem-inant, Neil J. Robinson, and Alexei M. Tsvelik, “Non-perturbative methodologies for low-dimensional strongly-correlated systems: From non-Abelian bosonization totruncated spectrum methods,” Rep. Prog. Phys. ,046002 (2018).[65] M. Kolodrubetz, V. Gritsev, and A. Polkovnikov, “Clas-sifying and measuring geometry of a quantum groundstate manifold,” Phys. Rev. B , 064304 (2013).[66] M. Kolodrubetz, D. Sels, P. Mehta, and A. Polkovnikov,“Geometry and non-adiabatic response in quantumand classical systems,” Phys. Rep. , 1–87 (2017),arXiv:1602.01062 [cond-mat.quant-gas].[67] P. Weinberg, M. Bukov, L. D’Alessio, A. Polkovnikov,S. Vajna, and M. Kolodrubetz, “Adiabatic perturba-tion theory and geometry of periodically-driven systems,”Phys. Rep. , 1–35 (2017).[68] S. Blanes, F. Casas, J.A. Oteo, and J. Ros, “The Magnusexpansion and some of its applications,” Phys. Rep. ,151 – 238 (2009).[69] P. Jordan and E. Wigner, “Über das Paulische Äquiv-alenzverbot,” Z. Phys. , 631–651 (1928).[70] Michael E. Fisher and Robert E. Hartwig, “Toeplitz De-terminants: Some Applications, Theorems, and Conjec-tures,” in Advances in Chemical Physics (John Wiley &Sons, Ltd, 2007) pp. 333–353.[71] Helen Au-Yang and Barry McCoy, “Theory of layeredIsing models. II. Spin correlation functions parallel tothe layering,” Phys. Rev. B , 3885–3905 (1974).[72] Estelle L. Basor and Craig A. Tracy, “The Fisher-Hartwigconjecture and generalizations,” Physica A , 167 –173 (1991).[73] P. J. Forrester and N. E. Frankel, “Applications andgeneralizations of Fisher–Hartwig asymptotics,” J. Math.Phys. , 2003–2028 (2004).[74] Albrecht Böttcher and Harold Widom, “Szegö via Ja-cobi,” Lin. Alg. Appl. , 656 – 667 (2006).[75] Alexei Yu. Karlovich, “Asymptotics of block Toeplitzdeterminants generated by factorable matrix functionswith equal partial indices,” Math. Nachr. , 1118–1127(2007).[76] P. Deift, A. Its, and I. Krasovsky, “Asymptotics ofToeplitz, Hankel, and Toeplitz+Hankel determinantswith Fisher-Hartwig singularities,” Ann. Math. ,1243–1299 (2011). Supplemental Material for: “Quantum quench in a driven Ising chain”
Neil J. Robinson , Isaac Pérez Castillo, and Edgar Guzmán-González Institute for Theoretical Physics, University of Amsterdam,Postbus 94485, 1090 GL Amsterdam, The Netherlands Departamento de Física, Universidad Autónoma Metropolitana-Iztapalapa,San Rafael Atlixco 186, Ciudad de México 09340, Mexico London Mathematical Laboratory, 18 Margravine Gardens, London W6 8RH, United Kingdom
In this Supplemental Material, we present:1. The Magnus expansion for the problem in the maintext.2. The exact mapping between time-dependent two-point functions of the driven systems and a suddenquench problem.3. Exact full-time evolution of one- and two-pointfunctions for a special case of the driven systemquench.4. Details of the computation of “sudden quench” cor-relation functions, summarizing both previouslyknown results and presenting new ones.5. Details of our numerical algorithm for simulatingtime-evolution in the driven system.6. Details of the procedure used to compare finite-sizenumerical results with analytical expressions.
MAGNUS EXPANSION
An alternative approach to compute the time-evolutionof the driven system is given by the Magnus expan- sion. For a time-dependent Hamiltonian H ( t ) , the time-evolution operator U ( t ) is defined via the differentialequation i dd t U ( t ) = H ( t ) U ( t ) , U ( t = 0) = . (S1)This is formally solved in terms of the well-known time-ordered exponential U ( t ) = T exp (cid:20) − i (cid:90) t d t (cid:48) H ( t (cid:48) ) (cid:21) . (S2)However, the time-ordering here (denoted by the oper-ator T ) is difficult to treat and so it may be beneficialto take an alternative approach. One such approximateapproach is the Magnus expansion [68], where the time-evolution operator is instead approximated by U ( t ) = exp [Ξ( t )] , Ξ( t ) = ∞ (cid:88) k =1 Ξ k ( t ) , (S3)where terms in the series are obtained from the differen-tial equation i dd t Ξ( t ) = ad Ξ exp( ad Ξ ) − H ( t ) . (S4)Here ad Ξ is the adjoint action of Ξ : ad Ξ ( O ) = [Ξ , O ] .The terms of the series Ξ k ( t ) have well-known forms Ξ ( t ) = − i (cid:90) t d t H ( t ) , Ξ ( t ) = − (cid:90) t d t (cid:90) t d t [ H ( t ) , H ( t )] , Ξ ( t ) = i (cid:90) t d t (cid:90) t d t (cid:90) t d t (cid:16)(cid:2) H ( t ) , [ H ( t ) , H ( t )] (cid:3) + (cid:2) H ( t ) , [ H ( t ) , H ( t )] (cid:3)(cid:17) . (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 mainbody of the manuscript, and examine the Magnus expan-sion for the time-evolution operator, working at small T ,i.e. high frequency for the driving. We will see thatthe integrability that we harnessed in the main text toperform exact calculations is completely hidden in theMagnus expansion. Indeed, the presence of a simple static limit (independent of integrability of that staticHamiltonian) is not at all apparent. To illustrate this,we consider the time-evolution operator over a single pe-riod T = 2 π/ Ω of the driving, U ( T ) . Treating the Mag-nus expansion as a power series approximation in T , thefirst order term in the expansion is just the average of theHamiltonian over a single period, so the driving vanishes: i Ξ ( T ) = − iT (cid:32) − J L (cid:88) l =1 σ zl σ zl +1 + h z L (cid:88) l =1 σ zl (cid:33) . (S6)For the second order term in the expansion, we requirethe commutator [ H ( t ) , H ( t )] = 2 ig L (cid:88) l =1 (cid:40) g sin (cid:16) Ω( t − t ) (cid:17) σ zl − h z (cid:104) ( c − c ) σ yl − ( s − s ) σ xl (cid:105) + J ( c − c ) (cid:16) σ zl σ yl +1 + σ yl σ zl +1 (cid:17) − J ( s − s ) (cid:16) σ zl σ xl +1 + σ xl σ zl +1 (cid:17)(cid:41) . where for conciseness we define the short hand notations c a = cos(Ω t a ) , s a = sin(Ω t a ) . Performing the doubleintegral, we obtain Ξ ( T ) : i Ξ ( T ) = gT π L (cid:88) l =1 (cid:20) h z σ xl − g σ zl − J (cid:16) σ zl σ xl +1 + σ xl σ zl +1 (cid:17)(cid:21) . (S7)In this second order term (in T ), we see already that theMagnus expansion for the time-evolution over one periodis rather messy. Continuing to third order, we expect togenerate terms of the form σ xl σ yl +1 , σ yl σ zl +1 too, so higherorders are unlikely to simplify things further. We see thatthe existence of a simple static Hamiltonian, as given inthe main body of the manuscript, is not obvious from theMagnus expansion. MAPPING BETWEEN DRIVEN AND “SUDDENQUENCH” TWO-POINT FUNCTIONS
Let us consider equal-time two-point functions in thenonequilibrium time-evolved wave function [see Eq. (7)of the main text], s αβ ( (cid:96) ; t ) = (cid:104) Ψ( t ) | σ αj σ βj + (cid:96) | Ψ( t ) (cid:105) . (S8)Using that U ( t ) is an element of SU(2) leads to the trans-formations described in Eq. (8) of the main text, and wearrive at the following mapping between driven and “sud-den quench” correlation functions: s zz ( (cid:96) ; t ) = s zz st ( (cid:96) ; t ) ,s xx ( (cid:96) ; t ) = cos (Ω t ) s xx st ( (cid:96) ; t ) + sin (Ω t ) s yy st ( (cid:96) ; t ) −
12 sin(2Ω t ) (cid:2) s xy st ( (cid:96) ; t ) + s yx st ( (cid:96) ; t ) (cid:3) ,s yy ( (cid:96) ; t ) = cos (Ω t ) s yy st ( (cid:96) ; t ) + sin (Ω t ) s xx st ( (cid:96) ; t )+ 12 sin(2Ω t ) (cid:2) s xy st ( (cid:96) ; t ) + s yx st ( (cid:96) ; t ) (cid:3) ,s xy ( (cid:96) ; t ) = cos (Ω t ) s xy st ( (cid:96) ; t ) − sin (Ω t ) s yx st ( (cid:96) ; t )+ 12 sin(2Ω t ) (cid:2) s xx st ( (cid:96) ; t ) − s yy st ( (cid:96) ; t ) (cid:3) ,s xz ( (cid:96) ; t ) = cos(Ω t ) s xz st ( (cid:96) ; t ) − sin(Ω t ) s yz st ( (cid:96) ; t ) ,s yz ( (cid:96) ; t ) = cos(Ω t ) s yz st ( (cid:96) ; t ) + sin(Ω t ) s xz st ( (cid:96) ; t ) . (S9) Here, as in the main body of the manuscript, we definethe “sudden quench” correlation functions s αβ st ( (cid:96) ; t ) = (cid:88) n,m e i ( E n − E m ) t (cid:104) E n | σ αj σ βj + (cid:96) | E m (cid:105)× (cid:104) Ψ | E n (cid:105)(cid:104) E m | Ψ (cid:105) , (S10)where | E n (cid:105) are eigenstates of the static Hamiltonian H st with energy E n and | Ψ (cid:105) = | Ψ( t = 0) (cid:105) is the initial statefor the time-evolution. We note that Eq. (S10) is thesame 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 sce-nario from our mapping. We consider starting in theground state of the Hamiltonian H init = − J L (cid:88) l =1 σ zl σ zl +1 − g (cid:88) l σ xl . (S11)Let us work with J > and g > . We perform a quenchwhere we start to rotate the transverse field and apply anadditional longitudinal field H ( t >
0) = − J L (cid:88) l =1 σ zl σ zl +1 + h z N (cid:88) l =1 σ zl − L (cid:88) l =1 g (cid:16) e − i Ω t σ + l + e i Ω t σ − l (cid:17) . (S12)with frequency Ω . We choose h z = Ω / , in which casethe effective static Hamiltonian is H st = − J L (cid:88) l =1 σ zl σ zl +1 − g L (cid:88) l =1 σ xl . (S13)That is, the effective state Hamiltonian is identical to theinitial Hamiltonian! In this case, the “sudden quench”correlation functions reduce to the equilibrium ones.These can be computed with standard free fermion tech-niques, see Ref. [43], and read: s z st ( (cid:96) ; t ) = (cid:0) − λ − (cid:1) / Θ( λ − , (S14) s x st ( (cid:96) ; t ) = G (0) , s y st ( (cid:96) ; t ) = 0 , (S15)0for the one-point functions, while the diagonal two-pointfunctions read: s zz st ( (cid:96) ; t ) = det G ( − G ( − . . . G ( − (cid:96) ) G (0) G ( − . . . G ( − (cid:96) + 1) ... ... ... G ( (cid:96) − G ( (cid:96) − . . . G ( − s xx st ( (cid:96) ; t ) = − G ( (cid:96) ) G ( − (cid:96) ) + G (0) ,s yy st ( (cid:96) ; t ) = det G (1) G (0) . . . G (2 − (cid:96) ) G (2) G (1) . . . G (3 − (cid:96) ) ... ... ... G ( (cid:96) ) G ( (cid:96) − . . . G (1) (S16)where we use notations, λ = J/g , Θ( x ≤
0) = 0 , Θ( x >
0) = 1 , and we define the functions G ( n ) = L ( n ) + λL ( n + 1) , (S17)with L ( n ) = 1 π (cid:90) π d k cos( kn ) (cid:112) λ + 2 λ cos( k ) . (S18)The off-diagonal (in spin indices) two-point functionsrequire some further work, not usually being consideredin equilibrium treatments of the quantum Ising chain.By symmetry arguments s αβ ( (cid:96), t ) = s βα ( (cid:96), t ) . Using theHeisenberg equation of motion, s zy st ( (cid:96), t ) ∝ ˙ s zz st ( (cid:96), t ) = 0 ,as we are working in an static case. The terms s xz st ( (cid:96), t ) and s xy st ( (cid:96), t ) can not be computed in general using thestandard techniques, as they mix the even sector with the odd sector in the fermionic theory (see the next sectionfor more details). However, in the paramagnetic phase,where λ < , using the symmetry of the Hamiltonianunder π rotations of the quantization axis about the x -direction, after some tedious algebra, we find s zy st ( (cid:96) ; t ) = 0 , (S19) s xz st ( (cid:96) ; t ) | λ< = 0 , (S20) s xy st ( (cid:96) ; t ) | λ< = 0 , (S21)Putting these known equilibrium results together withthe mapping between driven correlation functions and“sudden quench” ones, we obtain the full time-dependentcorrelation functions for the H init → H ( t > quench s z ( (cid:96) ; t ) = s z st ( (cid:96) ; t ) s x ( (cid:96) ; t ) = cos(Ω t ) s x st ( (cid:96) ; t ) ,s y ( (cid:96) ; t ) = sin(Ω t ) s x st ( (cid:96) ; t ) ,s zz ( (cid:96) ; t ) = s zz st ( (cid:96) ; t ) ,s xx ( (cid:96) ; t ) | λ< = cos (Ω t ) s xx st ( (cid:96) ; t ) + sin (Ω t ) s yy st ( (cid:96) ; t ) ,s yy ( (cid:96) ; t ) | λ< = cos (Ω t ) s yy st ( (cid:96) ; t ) + sin (Ω t ) s xx st ( (cid:96) ; t ) ,s xy ( (cid:96) ; t ) | λ< = sin(Ω t ) cos(Ω t ) (cid:104) s xx st ( (cid:96) ; t ) − s yy st ( (cid:96) ; t ) (cid:105) ,s xz ( (cid:96) ; t ) | λ< = 0 s yz ( (cid:96) ; t ) | λ< = 0 . (S22)Furthermore, in the limit (cid:96) (cid:29) one can simplify someof these results using an asymptotic analysis via Szegö’stheorem, which leads to (some) closed form analyticalresults [52–54] s zz st ( (cid:96) (cid:29) t ) (cid:12)(cid:12)(cid:12)(cid:12) λ< ≈ √ π (cid:0) − λ (cid:1) − / λ (cid:96) √ (cid:96) (cid:20) − (cid:96) λ − λ + O (cid:18) (cid:96) (cid:19)(cid:21) ,s zz st ( (cid:96) (cid:29) t ) (cid:12)(cid:12)(cid:12)(cid:12) λ> ≈ (cid:0) − λ − (cid:1) / + O (cid:18) (cid:96) (cid:19) ,s xx st ( (cid:96) (cid:29) t ) (cid:12)(cid:12)(cid:12)(cid:12) λ< ≈ ( s x st ( (cid:96) ; t )) − λ (cid:96) (cid:96) π (cid:18) − (cid:96) λ − − λ − + 1 (cid:19) s xx st ( (cid:96) (cid:29) t ) (cid:12)(cid:12)(cid:12)(cid:12) λ> ≈ ( s x st ( (cid:96) ; t )) − λ (cid:96) +2 π (cid:18) − (cid:96) λ − − λ − − (cid:19) s yy st ( (cid:96) (cid:29) t ) (cid:12)(cid:12)(cid:12)(cid:12) λ< ≈ − √ π (1 − λ ) / λ (cid:96) (cid:96) / (cid:20) (cid:96) (cid:18) − λ ) − (cid:19) + O (cid:18) (cid:96) (cid:19)(cid:21) s yy st ( (cid:96) (cid:29) t ) (cid:12)(cid:12)(cid:12)(cid:12) λ> ≈ − π (cid:0) − λ − (cid:1) − / λ − (cid:96) (cid:96) (cid:20) − (cid:96) λ − − λ − + O (cid:18) (cid:96) (cid:19)(cid:21) . (S23)Putting together these results, we obtain exact, closed form results for one-point functions and the asymptotics of1two-point correlation functions in a nonequilibrium driven problem. The non-equilibrium one-point functions read: s z ( t ) = (cid:0) − λ − (cid:1) / Θ( λ − , (S24) s x ( t ) (cid:12)(cid:12)(cid:12) λ (cid:54) =1 = cos(Ω t ) π | − λ | (cid:34) ( λ − E (cid:18) − λ (1 − λ ) (cid:19) − (cid:0) λ − (cid:1) K (cid:18) − λ (1 − λ ) (cid:19) (cid:35) , (S25) s y ( t ) (cid:12)(cid:12)(cid:12) λ (cid:54) =1 = sin(Ω t ) π | − λ | (cid:34) ( λ − E (cid:18) − λ (1 − λ ) (cid:19) − (cid:0) λ − (cid:1) K (cid:18) − λ (1 − λ ) (cid:19) (cid:35) (S26)where E ( m ) = (cid:82) π/ d θ [1 − m sin ( θ )] / and K ( m ) = (cid:82) π/ d θ [1 − m sin ( θ )] − / are elliptic functions. The asymptoticsof non-equilibrium two-point functions read: s zz ( (cid:96) (cid:29) t ) (cid:12)(cid:12)(cid:12)(cid:12) λ< = 1 √ π (cid:0) − λ (cid:1) − / λ (cid:96) √ (cid:96) (cid:20) − (cid:96) λ − λ + O (cid:18) (cid:96) (cid:19)(cid:21) , (S27) s zz ( (cid:96) (cid:29) t ) (cid:12)(cid:12)(cid:12)(cid:12) λ> = (cid:0) − λ − (cid:1) / + O (cid:18) (cid:96) (cid:19) , (S28) s xx ( (cid:96) (cid:29) t ) (cid:12)(cid:12)(cid:12)(cid:12) λ< = cos (Ω t ) (cid:34) ( s x st ( (cid:96) ; t )) − λ (cid:96) (cid:96) π (cid:18) − (cid:96) − λ λ (cid:19) (cid:35) − sin (Ω t ) (cid:34) √ π (1 − λ ) / λ (cid:96) (cid:96) / (cid:20) (cid:96) (cid:18) − λ ) − (cid:19) + O (cid:18) (cid:96) (cid:19)(cid:21) (S29) s yy ( (cid:96) ; t ) (cid:12)(cid:12)(cid:12)(cid:12) λ< = sin (Ω t ) (cid:34) ( s x st ( (cid:96) ; t )) − λ (cid:96) (cid:96) π (cid:18) − (cid:96) − λ λ (cid:19) (cid:35) − cos (Ω t ) (cid:34) √ π (1 − λ ) / λ (cid:96) (cid:96) / (cid:20) (cid:96) (cid:18) − λ ) − (cid:19) + O (cid:18) (cid:96) (cid:19)(cid:21) (S30) s xy ( (cid:96) ; t ) (cid:12)(cid:12)(cid:12)(cid:12) λ< = sin(Ω t ) cos(Ω t ) (cid:34) ( s x st ( (cid:96) ; t )) − λ (cid:96) (cid:96) π (cid:18) − (cid:96) − λ λ (cid:19) + 12 √ π (1 − λ ) / λ (cid:96) (cid:96) / (cid:20) (cid:96) (cid:18) − λ ) − (cid:19) + O (cid:18) (cid:96) (cid:19)(cid:21) (cid:35) (S31) s xz ( (cid:96) ; t ) (cid:12)(cid:12)(cid:12)(cid:12) λ< = 0 (S32) s yz ( (cid:96) ; t ) (cid:12)(cid:12)(cid:12)(cid:12) λ< = 0 . (S33) THE “SUDDEN QUENCH” CORRELATIONFUNCTIONS
In this appendix we detail calculations of the “suddenquench” correlation functions. We follow the approach ofCalabrese, Essler and Fagotti [49–51] in computing one-and two-point correlation functions following a quantumquench of the transverse field in the quantum Ising chain.Much of the discussion will follow similar lines, althoughwe will consider a number of observables that were notcomputed in their work.
Fermionization of the model
Let us begin by describing the exact solution of thetransverse field Ising model in terms of free fermions.The transverse field Ising model has the Hamiltonian H = − J L (cid:88) l =1 σ zl σ zl +1 − g L (cid:88) l =1 σ xl , (S34)where σ αl are the spin-1/2 operators acting on the l th siteof an L chain (we take L even), J is the exchange interac-tion, and g is the transverse field strength. We considerperiodic boundary conditions, identifying σ αL +1 = σ α . Toexpress the model in terms of free fermions, we first per-2form a π/ rotation of the quantization axes about the y axis to be consistent with standard notations H = − J L (cid:88) l =1 τ xl τ xl +1 − g L (cid:88) l =1 τ zl ≡ − ¯ J L (cid:88) l =1 (cid:32) τ xl τ xl +1 + ¯ gτ zl (cid:33) . (S35)Here τ αl are the Pauli matrices with the rotated quanti-zation axes, ¯ J = J and ¯ g = g/J .We proceed to fermionize the problem via the Jordan-Wigner transformation [69] τ zl = 1 − a † l a l ,τ − l = exp iπ (cid:88) j In the even sector, we have exp( iπ ˆ N ) = +1 and theHamiltonian (S37) reads H e = − ¯ J L (cid:88) l =1 (cid:16) a † l − a l (cid:17)(cid:16) a l +1 + a † l +1 (cid:17) + ¯ J ¯ g L (cid:88) l =1 (cid:16) a † l a l − a l a † l (cid:17) , (S39)where the fermions obey antiperiodic boundary conditions a L +1 = − a . The Hamiltonian can be diagonalized byfirst Fourier transforming and then performing a Bogoli-ubov rotation. The Fourier transform reads c k n = 1 √ L L (cid:88) l =1 e ik n l a l , (S40)with k n = π (2 n + 1) /L with n = − L/ , − L/ , . . . , L/ − . In terms of the Fourier modes, the Hamil-tonian reads H e = − ¯ J (cid:88) n (cid:20) cos( k n ) (cid:16) c † k n c k n − c − k n c †− k n (cid:17) − i sin( k n ) (cid:16) c † k n c †− k n − c − k n c k n (cid:17) − ¯ g (cid:16) c † k n c k n − c − k n c †− k n (cid:17) (cid:21) . (S41)This can be diagonalized via the Bogoliubov rotation (cid:18) c k n c †− k n (cid:19) = cos (cid:16) θ kn (cid:17) i sin (cid:16) θ kn (cid:17) i sin (cid:16) θ kn (cid:17) cos (cid:16) θ kn (cid:17) (cid:18) α k n α †− k n (cid:19) , (S42)where the Bogoliubov angle θ k n satisfies θ k n = − sgn ( k n ) arccos (cid:32) ¯ g − cos( k n ) (cid:112) ¯ g − g cos( k n ) + 1 (cid:33) (S43)to obtain the Hamiltonian H e = (cid:88) n (cid:15) ¯ g ( k n ) (cid:16) α † k n α k n − (cid:17) , (S44)with dispersion relation (cid:15) ¯ g ( k n ) = 2 ¯ J (cid:112) ¯ g + 2¯ g cos( k n ) + 1 . (S45)Thus the Hamiltonian has been diagonalized.Let us now say a few words about the nature of theground state in the even sector. As is expected from aBogoliubov rotation (and as we well know from physicsof the Ising chain), the dispersion is generically gapped,i.e. (cid:15) ¯ g ( k n ) > . As a result, the ground state | Ω e , ¯ g (cid:105) isthe vacuum for the Bogoliubov fermions: α k n | Ω e , ¯ g (cid:105) = 0 , ∀ n. (S46)3It immediately follows that the ground state energy is E Ω e , ¯ g = − (cid:88) n (cid:15) ¯ g ( k n ) . (S47)It will be useful to understand the ground state | Ω e , ¯ g (cid:105) in terms of the Jordan-Wigner fermions c k n . InvertingEq. (S42) to express α k n in terms c k n , one finds thatEq. (S46) becomes ( ∀ n ) (cid:104) cos (cid:18) θ k n (cid:19) c k n − i sin (cid:18) θ k n (cid:19) c †− k n (cid:105) | Ω e , ¯ g (cid:105) = 0 . (S48)As | Ω e , ¯ g (cid:105) is in the even sector with zero momentum,it must be of the form (cid:81) k n ≥ ( a + bc †− k n c † k n ) | (cid:105) , where | (cid:105) is the vacuum for Jordan-Wigner fermions. Workingthrough some algebra, one finds | Ω e , ¯ g (cid:105) = (cid:89) k n ≥ (cid:20) cos (cid:18) θ k n (cid:19) − i sin (cid:18) θ k n (cid:19) c †− k n c † k n (cid:21) | (cid:105) . (S49)This can be written in a “squeezed state” form followingsome simple manipulations | Ω e , ¯ g (cid:105) = (cid:89) k n ≥ cos (cid:18) θ k n (cid:19) × exp − i (cid:88) k n ≥ tan (cid:18) θ k n (cid:19) c †− k n c † k n | (cid:105) . (S50) Odd number of fermions The story with the odd fermion number sector is sim-ilar to the even one. This time we have exp( iπ ˆ N ) = − and instead of antiperiodic boundary conditions, we haveperiodic ones a L +1 = a . The Fourier transform of theJordan-Wigner fermions then reads c p n = 1 √ L (cid:88) l e ip n l a l (S51)with the momenta p n = 2 nπ/L and n = − L/ , − L/ , . . . ... . We again perform a Bogoliubov rotation (cid:18) c p n c †− p n (cid:19) = cos (cid:16) θ pn (cid:17) i sin (cid:16) θ pn (cid:17) i sin (cid:16) θ pn (cid:17) cos (cid:16) θ pn (cid:17) (cid:18) α p n α †− p n (cid:19) , (S52) with the Bogoliubov angle as defined in the even sector,with the additional definition of θ p = 0 . The Hamilto-nian (S37) is then H o = (cid:88) n (cid:54) =0 (cid:15) ¯ g ( p n ) (cid:16) α † p n α p n − (cid:17) − J (1 − ¯ g ) (cid:16) α † α − (cid:17) . (S53)We see now that the ground state in the odd sector is | Ω o (cid:105) = α † | α (cid:105) , (S54)where | α (cid:105) is the vacuum state for α p n . The ground state | Ω o (cid:105) has energy E Ω o = − ¯ J (1 − ¯ g ) − (cid:88) n (cid:54) =0 (cid:15) ¯ g ( p n ) . (S55)As with the even sector, it will be useful to know how toexpress the odd ground state, Eq. (S54), in terms of theJordan-Wigner fermions. Following the same sequenceof manipulations as the even sector, one can express thestate in the form | Ω o , ¯ g (cid:105) = (cid:89) p n ≥ cos (cid:18) θ p n (cid:19) × exp − i (cid:88) p n ≥ tan (cid:18) θ p n (cid:19) c †− p n c † p n c † | (cid:105) . (S56)As expected, this state has odd fermion parity. Time-evolution following a quantum quench Let us now turn to the problem of describing time-evolution following a quantum quench. In particular, weconsider the scenario where for time t < the system isin an eigenstate of the Hamiltonian H i with ¯ g = g i andat time t = 0 this is suddenly changed to ¯ g = g f andthe state evolves under this new Hamiltonian H f for alltimes t > .At times t < we can diagonalize the Hamiltonian interms of fermions β by performing the Bogoliubov ro-tation with angle φ . For times t > we can similarlydiagonalize the Hamiltonian in terms of fermions α viathe Bogoliubov rotation with angle θ . Equating the ex-pressions for the Jordan-Wigner fermions a k in terms ofthe Bogoliubov fermions, we have4 cos (cid:18) θ k n (cid:19) α k n + i sin (cid:18) θ k n (cid:19) α †− k n = cos (cid:18) φ k n (cid:19) β k n + i sin (cid:18) φ k n (cid:19) β †− k n , (S57) i sin (cid:18) θ k n (cid:19) α k n + cos (cid:18) θ k n (cid:19) α †− k n = i sin (cid:18) φ k n (cid:19) β k n + cos (cid:18) φ k n (cid:19) β †− k n . (S58)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 α k n ( t ) = exp (cid:0) − i(cid:15) g f ( k n ) t (cid:1) α k n . (S59)Thus we write the β fermions in terms of α ones and time-evolve them β k n ( t ) = cos (cid:18) θ k n − φ k n (cid:19) e − i(cid:15) gf ( k n ) t α k n + i sin (cid:18) θ k n − φ k n (cid:19) e i(cid:15) gf ( k n ) t α †− k n , (S60)before re-expressing the α fermions in terms of the β ones: β k n ( t ) = (cid:20) cos (cid:16) (cid:15) g f ( k n ) t (cid:17) − i sin (cid:16) (cid:15) g f ( k n ) t (cid:17) cos (cid:16) θ k n − φ k n (cid:17)(cid:21) β k n − sin (cid:16) (cid:15) g f ( k n ) t (cid:17) sin (cid:16) θ k n − φ k n (cid:17) β †− k n . (S61)We see that the difference in Bogoliubov angles oftenappears, so we denote this by the short hand ∆ k = θ k − φ k . (S62)With this, and with some work, we can compute time-evolution of observables. In the following we first con-sider one-point functions, before continuing to two-pointfunctions. Time-evolution of one-point functions We will consider the quench where we start from theground state of the initial Hamiltonian H i ≡ H ( t < with ¯ g = g i and time-evolve according to the Hamiltonian H f where ¯ g = g f . In the thermodynamic limit, the initialstate will be a superposition of the ground states in theeven and odd sectors | Ω , ¯ g (cid:105) = 1 √ (cid:16) | Ω e , ¯ g (cid:105)(cid:105) + | Ω o , ¯ g (cid:105)(cid:105) (cid:17) . (S63)(Here we work with normalized states such that (cid:104)(cid:104) Ω i , ¯ g | Ω i , ¯ g (cid:105)(cid:105) = 1 .) Thus, generally, there will be threesorts of terms when computing observables: expectation values within the even sector, those within the odd sec-tor, and those that mix the two sectors. Depending onthe operator being considered, the mixing terms may beforbidden.Let us begin with the easiest case, and consider thetime-evolution of the one-point function of τ zl . Underthe Jordan-Wigner transformation (S36) this maps to thelocal number operator of the Jordan-Wigner fermions.Thus we want to compute t z ( t ) = (cid:104) Ω , g i | e iH f t (cid:16) − a † l a l (cid:17) e − iH f t | Ω , g i (cid:105) . (S64)Clearly t z ( t ) is independent of l by translational in-variance. As τ z is number conserving in terms of thefermions, the terms mixing even and odd sectors fromEq. (S63) will vanish. Thus we have t z ( t ) = (cid:104)(cid:104) Ω e , g i | (cid:18) − a † l ( t ) a l ( t ) (cid:19) | Ω e , g i (cid:105)(cid:105) + (cid:104)(cid:104) Ω o , g i | (cid:18) − a † l ( t ) a l ( t ) (cid:19) | Ω o , g i (cid:105)(cid:105) . (S65)Let us take the first term, where the ground states is inthe 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 a l ( t ) e = 1 √ L (cid:88) n e − ik n l c k n ( t ) = 1 √ L (cid:88) n e − ik n l (cid:20) cos (cid:18) φ k n (cid:19) β k n ( t ) + i sin (cid:18) φ k n (cid:19) β †− k n ( t ) (cid:21) . (S66)5Working through some tedious algebra, one finds t z ( t ) e = 12 − L (cid:88) n (cid:20) sin (cid:0) (cid:15) g f ( k n ) t (cid:1) sin (cid:18) ∆ k n + φ k n (cid:19) + sin (cid:18) φ k n (cid:19) cos (cid:0) (cid:15) g f ( k n ) t (cid:1)(cid:21) , (S67) = 12 − π (cid:90) π − π d k (cid:20) sin (cid:0) (cid:15) g f ( k ) t (cid:1) sin (cid:18) ∆ k + φ k (cid:19) + sin (cid:18) φ k (cid:19) cos (cid:0) (cid:15) g f ( k ) t (cid:1)(cid:21) . (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 − π (cid:90) π − π d k (cid:20) sin (cid:0) (cid:15) g f ( k ) t (cid:1) sin (cid:18) ∆ k + φ k (cid:19) + sin (cid:18) φ k (cid:19) cos (cid:0) (cid:15) g f ( k ) t (cid:1)(cid:21) . (S69)Let us now move on to considering τ xl and τ yl . Theseoperators are similar in that they change the number ofJordan-Wigner fermions by one, and hence their expecta-tion values consist only of the “cross terms” between evenand odd ground states. This is actually quite problem-atic for the free fermion techniques discussed here, as itis not clear how to deal with such terms. Instead, we usethe following trick: we relate the one-point function tothe two-point function in the limit of infinite separationof the operators lim (cid:96) →∞ (cid:104)(cid:104) Ω; g i | τ xj ( t ) τ xj + (cid:96) ( t ) | Ω; g i (cid:105)(cid:105) → (cid:104)(cid:104) Ω; g i | τ xj ( t ) | Ω; g i (cid:105)(cid:105) , (S70)and similarly for τ y . As such, this naturally brings us onto the topic of computing two-point functions. Time-evolution of two-point functions We now turn our attention to computing two-pointfunctions t αβ ( (cid:96) ; t ) = (cid:104) Ω , g i | e iH f t τ αj τ βj + (cid:96) e − iH f t | Ω , g i (cid:105) . (S71)Two examples of these, α = β = x and α = β = z ,were considered in Refs. [49–51]. We will also considerthese here, as well as some other cases that will be of use in this work. As with the one-point functions, we areunable to compute two-point functions when the operatorwithin the expectation value is odd under spin inversion(and thus connects states in the even and odd sectorsof the Hilbert space). Thus we restrict our attention tooperators that are even under spin inversion.This restriction, in the context of the main text, leadsus to focus our results on one of two cases. Firstly, wecan consider quenches starting in states that are even un-der spin inversion (and thus the expectation value van-ishes by symmetry). This allows us, for example, to com-pute the full time evolution of two-point functions in thedriven problem where the “effective sudden quench” iswithin the paramagnetic phase of the static Hamiltonian.Secondly, we can restrict attention to stroboscopic times,where contributions from these problematic expectationvalues vanish. The t zz ( (cid:96) ; t ) two-point function Let us begin with the easiest case, where the opera-tors in the two-point function are solely in the transversefield direction. Then, according to the Jordan-Wignertransformation (S36), we need to compute the two-pointfunction of the (time-dependent) Jordan-Wigner fermiondensity t zz ( (cid:96) ; t ) = (cid:104) Ω , g i | e iH f t (cid:16) − a † j a j (cid:17)(cid:16) − a † j + (cid:96) a j + (cid:96) (cid:17) e − iH f t | Ω , g i (cid:105) . (S72)The operator under consideration is even under spin inversion, thus matrix elements between states in the even andodd sectors vanish, leaving only t zz ( (cid:96) ; t ) = 12 (cid:104)(cid:104) Ω e , g i | e iH f t (cid:16) − a † j a j (cid:17)(cid:16) − a † j + (cid:96) a j + (cid:96) (cid:17) e − iH f t | Ω e , g i (cid:105)(cid:105) + ( e ↔ o ) . (S73)Using Eq. (S66) (and similar for the odd sector) and taking the thermodynamic limit, one finds t zz ( (cid:96) ; t ) = − (cid:90) π − π d k π d k (cid:48) π (cid:40) e i ( k − k (cid:48) ) (cid:96) (cid:104) sin(2 (cid:15) g f ( k ) t ) sin(∆ k ) sin(2 (cid:15) g f ( k (cid:48) ) t ) sin(∆ k (cid:48) ) + (cid:16) e i ( k − k (cid:48) ) (cid:96) − (cid:17) × e i ( θ k + θ k (cid:48) ) (cid:16) cos(∆ k ) − i sin(∆ k ) cos(2 (cid:15) g f ( k ) t ) (cid:17)(cid:16) cos(∆ k (cid:48) ) − i sin(∆ k (cid:48) ) cos(2 (cid:15) g f ( k (cid:48) ) t ) (cid:17)(cid:41) . (S74)6 The t xx ( (cid:96) ; t ) two-point function Let us now turn our attention towards the two-pointfunction where both operators are in the Ising direction. This is significantly more complicated than t zz ( (cid:96) ; t ) as theJordan-Wigner transformation (S36) introduces “strings”of density operators between j and j + (cid:96) . This can easilybe seen by Taylor expanding the exponential factors inEq. (S36): exp ± iπ (cid:88) j<(cid:96) a † j a j = (cid:89) j<(cid:96) exp (cid:16) ± iπa † j a j (cid:17) = (cid:89) j<(cid:96) (cid:16) − a † j a j (cid:17) . (S75)Thus, in terms of the Jordan-Wigner fermions, we need to compute t xx ( (cid:96) ; t ) = (cid:104) Ω , g i | e iH f t ( τ + j + τ − j )( τ + j + (cid:96) + τ − j + (cid:96) ) e − iH f t | Ω , g i (cid:105) , (S76) = (cid:104) Ω , g i | e iH f t (cid:16) a j + a † j (cid:17) (cid:34) (cid:96) − (cid:89) k =0 (cid:16) − a † j + k a j + k (cid:17)(cid:35) (cid:16) a j + (cid:96) + a † j + (cid:96) (cid:17) e − iH f t | Ω , g i (cid:105) , (S77) = (cid:104) Ω , g i | e iH f t (cid:16) a † j − a j (cid:17) (cid:34) (cid:96) − (cid:89) k =1 (cid:16) − a † j + k a j + k (cid:17)(cid:35) (cid:16) a j + (cid:96) + a † j + (cid:96) (cid:17) e − iH f t | Ω , g i (cid:105) . (S78)The calculation will proceed more easily if we introduceMajorana fermions, which describe the real and imagi-nary parts of the Jordan-Wigner fermions, γ xl = a l + a † l , γ yl = − i ( a l − a † l ) , (S79)which obey the anticommutation relations { γ aj , γ bl } =2 δ a,b δ j,l . In terms of these Majorana operators we have a † j a j = 12 (cid:0) iγ xj γ yj (cid:1) , exp( iπa † j a j ) = − iγ xj γ yj , (S80)and hence t xx ( (cid:96) ; t ) = (cid:104)(cid:104) Ω e , g i | e iH f t (cid:96) (cid:89) k =1 ( − iγ yk γ xk +1 ) e − iH f t | Ω e , g i (cid:105)(cid:105) . (S81)This can be written as the Pfaffian of a (cid:96) × (cid:96) matrix t xx ( (cid:96) ; t ) = ( − i ) (cid:96) Pf (Γ) = Pf ( − i Γ) , (S82) Γ = Γ (cid:48) Γ (cid:48) · · · Γ (cid:48) (cid:96) Γ (cid:48) Γ (cid:48) · · · Γ (cid:48) (cid:96) ... ... . . . ... Γ (cid:48) (cid:96) Γ (cid:48) (cid:96) · · · Γ (cid:48) (cid:96)(cid:96) , (S83)where Γ (cid:48) is a × sub-matrix Γ (cid:48) mn = (cid:18) (cid:104)(cid:104) γ ym − n + j γ yj (cid:105)(cid:105) t − δ m,n (cid:104)(cid:104) γ ym − n − j γ xj (cid:105)(cid:105) t (cid:104)(cid:104) γ xm − n +1+ j γ yj (cid:105)(cid:105) t (cid:104)(cid:104) γ xm − n + j γ xj (cid:105)(cid:105) t − δ m,n (cid:19) , (S84)where (cid:104)(cid:104) O (cid:105)(cid:105) t = (cid:104)(cid:104) Ω e , g i | e iH f t Oe − iH f t | Ω e , g i (cid:105)(cid:105) . (S85) Using reflection symmetry and translational invariance,we can rewrite this sub-matrix as ¯Γ η = − i Γ (cid:48) m + η,m , (S86) = − i (cid:18) (cid:104)(cid:104) γ yη + j γ yj (cid:105)(cid:105) t − δ η, (cid:104)(cid:104) γ yη − j γ xj (cid:105)(cid:105) t −(cid:104)(cid:104) γ y − η − j γ xj (cid:105)(cid:105) t −(cid:104)(cid:104) γ yη + j γ yj (cid:105)(cid:105) t + δ η, (cid:19) (S87) = (cid:18) − f η g η − g − η f η (cid:19) , (S88)where f η and g η can be obtained using Eq. (S66), andthey read: f η = (cid:90) π − π d k π e ikη f ( k ) , (S89) f ( k ) = i sin (cid:16) (cid:15) g f ( k ) t (cid:17) sin (cid:16) ∆ k (cid:17) , (S90) g η = (cid:90) π − π d k π e ikη g ( k ) , (S91) g ( k ) = − e i ∆ k (cid:20) cos (cid:16) ∆ k (cid:17) − i sin (cid:16) ∆ k (cid:17) cos (cid:16) (cid:15) g f ( k ) t (cid:17)(cid:21) . (S92)When the two-point correlation function is written interms of (S88), it is the Pfaffian of a block Toeplitz ma-trix: t xx ( (cid:96) ; t ) = Pf ¯Γ ¯Γ − · · · ¯Γ − (cid:96) ¯Γ ¯Γ · · · ¯Γ − (cid:96) ... ... . . . ... ¯Γ (cid:96) − ¯Γ (cid:96) − · · · ¯Γ . (S93)The matrix here is real and antisymmetric, thus it isan anti-Hermitian matrix. All eigenvalues are imaginary7and appear in conjugate pairs. Furthermore, for an an-tisymmetric matrix the Pfaffian is the square root of thedeterminant, up to a difficult-to-determine sign. In otherwords, if an (cid:96) × (cid:96) antisymmetric matrix Λ = − Λ T has eigenvalues ± λ j ( j = 1 , . . . , (cid:96) ), the Pfaffian readsPf (Λ) = (cid:81) (cid:96)j =1 iλ j .One can reduce computing the Pfaffian to computingthe determinant of an (cid:96) × (cid:96) matrix Qt xx ( (cid:96) ; t ) = ( − (cid:96) ( (cid:96) − / det ( Q ) , (S94) Q nm = if n − m + g n + m − (cid:96) − , (S95)which is convenient to do numerically. One can also makeanalytical progress by computing the determinant of theblock Toeplitz matrix, as considered in Ref. [49–51], viaapplications of Szegö’s theorem and Fisher-Hartwig con-jectures [70–76]. The t yy ( (cid:96) ; t ) two-point function The computation of this two-point function is verysimilar to the previous one. In terms of the Majoranafermions, t yy ( (cid:96) ; t ) = (cid:104)(cid:104) Ω e , g i | e iH f t (cid:96) (cid:89) k =1 ( iγ xk γ yk +1 ) e − iH f t | Ω e , g i (cid:105)(cid:105) . (S96)The previous expression can be written as the Pfaffian ofa block Toeplitz matrix, t yy ( (cid:96) ; t ) = Pf ˜Γ ˜Γ − · · · ˜Γ − (cid:96) ˜Γ ˜Γ · · · ˜Γ − (cid:96) ... ... . . . ... ˜Γ (cid:96) − ˜Γ (cid:96) − · · · ˜Γ . (S97)where ˜Γ η is a × sub-matrix ˜Γ η = (cid:18) − f η g − η +2 − g η +2 f η (cid:19) , (S98)where f η and g η are defined in (S89)-(S91).The computation of the Pfaffian is equivalent to theone the determinant of an (cid:96) × (cid:96) matrix ˜ Qt yy ( (cid:96) ; t ) = ( − (cid:96) ( (cid:96) − / det ( ˜ Q ) , (S99) ˜ Q nm = if n − m + g (cid:96) +3 − ( n + m ) . (S100) The t xy ( (cid:96) ; t ) two-point function This two point-function can be obtained from t xx ( (cid:96) ; t ) .Indeed, by considering the Heisenberg equation of mo- tion, dd t τ xj ( t ) τ xj + (cid:96) ( t ) = i [ H f , τ xj ( t ) τ xj + (cid:96) ( t )]= 2 ¯ J ¯ g f (cid:16) τ yj ( t ) τ xj + (cid:96) ( t ) + τ xj ( t ) τ yj + (cid:96) ( t ) (cid:17) . (S101)By computing the expectation value and considering in-version symmetry we conclude, t xy ( (cid:96) ; t ) = 14 ¯ J ¯ g f dd t t xx ( (cid:96) ; t ) . (S102) DETAILS OF THE NUMERICAL ALGORITHMFOR TIME-EVOLUTION OF A DRIVEN SYSTEM In this section of the Supplemental Material, we dis-cuss details of the numerical algorithm used to computetime-evolution of the initial state used in the main text.As we have already shown in the main text, the resultsof the simulations have an excellent agreement with theanalytical results for the integrable case when Ω = 2 h z .We make the following numerical studies outside this in-tegrable line.We consider the time-evolution of a given initial state | Ψ (cid:105) induced by a time-dependent Hamiltonian, Eq. (1)of the main text. Formal solution of the time-dependentSchrödinger equation gives | Ψ( t ) (cid:105) = T exp (cid:18) − i (cid:90) t d t (cid:48) H ( t (cid:48) ) (cid:19) | Ψ (cid:105) , (S103)where T is the time-ordering operator. Our Hamiltonian H ( t ) depends smoothly on time t , which makes dealingwith the time-ordering tricky. To proceed we “Trotterize”the time-evolution into N small steps (of size ∆ t , suchthat N ∆ t = t ), each of which has a fixed Hamiltonian | Ψ( t ) (cid:105) ≈ T N − (cid:89) j =0 exp (cid:16) − iH (cid:0) j ∆ t (cid:1) ∆ t (cid:17) | Ψ (cid:105) , (S104) ≡ T N − (cid:89) j =1 U j | Ψ (cid:105) . (S105)The problem now becomes how to describe the action ofthe time-evolution operator U j on a state.We compute the action of a single step time-evolutionoperator on a state via the Chebyshev expansion (see,e.g., Ref. [45]). This avoids the need to diagonalize orexponentiate the Hamiltonian, and it is computationallyefficient. The Chebyshev expansion of the time-evolutionoperator reads: U j = J (cid:0)(cid:102) ∆ t (cid:1) +2 ∞ (cid:88) n =1 ( − i ) n J n (cid:0)(cid:102) ∆ t (cid:1) T n (cid:16) (cid:101) H ( j ∆ t ) (cid:17) . (S106)8Here we have introduced the rescaled time (cid:102) ∆ t and Hamil-tonian (cid:101) H ( j ∆ t ) , as well as the Bessel functions J n ( x ) andthe Chebyshev matrices T n ( x ) , which satisfy the recur-sion relation T ( x ) = , T ( x ) = x, (S107) T n ( x ) = 2 xT n − ( x ) − T n − ( x ) . (S108)The rescaled time and Hamiltonian satisfy (cid:101) H ( j ∆ t ) (cid:102) ∆ t = H ( j ∆ t )∆ t , and enforce that the eigenvalues of the (cid:101) H ( j ∆ t ) lie in the interval [ − , .In practice, the Chebyshev expansion (S106) has to betruncated to finite order. In the following subsectionswe present some convergence checks of our “Trotteriza-tion+Chebyshev” procedure, where we examine the ef-fects of Trotter step size and Chebyshev expansion or-der. We find that only a low order Chebyshev expansionis required for the problem at hand, leading to a very ef-ficient algorithm for simulating the time-evolution of thecontinuously (smoothly) driven quantum system. We fin-ish our discussions with a brief inspection of finite systemsize effects on the nonequilibrium dynamics. Convergence with order of the Chebyshev expansion Let us first examine how results of the time-evolutionvary with the order of the Chebyshev expansion (S106).We fix the Trotter step size ∆ t = 0 . and then proceedto compute an approximation of the time-evolution: U j ≈ J (cid:0)(cid:102) ∆ t (cid:1) +2 N c (cid:88) n =1 ( − i ) n J n (cid:0)(cid:102) ∆ t (cid:1) T n (cid:16) (cid:101) H ( j ∆ t ) (cid:17) , (S109)for given values of N c . We present an example inFig. S1(a) for the time-evolution of a one-point functionfollowing a quench within the disordered (paramagnetic)phase of the quantum Ising chain (details of parametersare given in the figure caption). We observed that thereis rapid convergence of the result for increasing order ofthe Chebyshev expansion: after ten periods of the drive,the results with N c = 8 and N c = 32 agree to seven dec-imal places. We also stress the fact that, despite beingoutside the integrable line Ω = 2 h z , there is no heatingto infinite temperature.One may question whether similarly good conver-gence is observed for more complicated observables. InFig. S1(b), we present s zz (1; t ) for the same quench,where it is apparent that two-point functions also con-verge similarly well with increasing N c . Within the mainbody of the text, we have explicitly checked that the ob-tained results are independent of the expansion order N c . (a) − − . . 51 0 2 4 6 8 10 s x ( t ) Time, t/TN c = 2 N c = 4 N c = 8 N c = 16 N c = 32 (b) . . . . . . . . s zz ( ; t ) Time, t/TN c = 8 N c = 32 FIG. S1. The time-evolution of the one-point function (a) s x ( t ) = (cid:104) Ψ( t ) | σ x | Ψ( t ) (cid:105) ; (b) s zz (1; t ) = (cid:104) Ψ( t ) | σ zj σ zj +1 | Ψ( t ) (cid:105) 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 . . The time-evolution iscomputed via “Trotterization+Chebyshev” with Trotter step ∆ t = 0 . and Chebyshev expansion order N c for ten peri-ods of the drive, T = 2 π/ Ω . We see that results are rapidlyconverging with increased N c . Note that there is no heatingto infinite temperature, despite the fact we are outside theintegrable line Ω = 2 h z . Convergence with Trotter step size Let us now turn our attention to the size of the timediscretization step. Finite size effects With convergence established for Trotter step size andorder of the Chebyshev expansion, we finally examine fi-nite size effects in our simulations of the time-evolution.9 − − . . 51 0 0 . . . . s x ( t ) Time, t/TL = 6 L = 8 L = 10 L = 12 L = 14 L = 16 FIG. S2. The time-evolution of the one-point function s x ( t ) = (cid:104) Ψ( t ) | σ x | Ψ( t ) (cid:105) following the quench presented in Fig. S1. Re-sults are presented for a number of system sizes, illustratingfinite size effects in numerical data due to small accessiblesystems. We see that for L = 14 , the time-evolution is wellmatched for almost the whole period of the driving. Fixing ∆ t = 0 . and N c = 64 , we examine the nonequi-librium time-evolution of σ x following the same quenchas in the previous subsections.We first examine a local observable, (cid:104) Ψ( t ) | σ x | Ψ( t ) (cid:105) ,and check how it behaves with changes in the systemsize L , i.e. what the “finite size effects” are. Exampledata is presented in Fig. S2, where we see the expectedbehaviour: at short times observables for all system sizesare in agreement. Under time-evolution, where propagat-ing excitations are generated, results for different systemsizes eventually diverge due to finite size revivals. Fora single, sudden quench the picture for this is simple: aquantum quench generates excitations that can propa-gate around the system. By conservation of momentum,such excitations must be generated in pairs, with momen-tum k and − k . These can propagate around the system,eventually meeting back at the start and interfering – apurely finite size effect. The time for which this occursincreases linearly with system size. Here we are seeingthe driven system analogue of this finite size revival be-haviour. Over a single period, we see that results for thetwo largest system sizes, L = 14 and L = 16 , match overalmost the whole drive period.Finite size effects can more easily be revealed in non-local/global measurements of the system. We illustrate this in Fig. S3, where we consider the return amplitude (cid:12)(cid:12)(cid:12) (cid:104) Ψ | Ψ( t ) (cid:105) (cid:12)(cid:12)(cid:12) following a quench. This rapidly decays asa function of time towards a close-to-zero value. Thisclose-to-zero value depends exponentially on the systemsize. After one period of driving (the right hand side ofthe plot) we see that the state | Ψ( T ) (cid:105) ≈ | Ψ (cid:105) , with theoverlap decreasing with system size. − − − − − − − − . . . . (cid:12)(cid:12)(cid:12) h Ψ | Ψ ( t ) i (cid:12)(cid:12)(cid:12) Time, t/TL = 8 L = 10 L = 12 L = 14 L = 16 FIG. S3. The return amplitude (cid:12)(cid:12)(cid:12) (cid:104) Ψ | Ψ( t ) (cid:105) (cid:12)(cid:12)(cid:12) as a function oftime t over one period following the same quench as in Fig. S1and Fig. S2. Here data is presented for a number of systemsizes, L , revealing significant finite size effects in this nonlocal property. Each simulation was performed with ∆ t = 0 . and N c = 64 . We see that nonlocal observables experiencing severefinite size effects does not carry through to local observ-ables, cf. Figs. S2 and S3. COMPARISON OF ANALYTICAL RESULTSWITH FINITE VOLUME NUMERICS The equations (S69),(S74),(S94), and (S99) for the oneand two point-functions are exact in the thermodynamiclimit. However, to compare with our numerical results, itis better to consider analytical expressions for finite sizesystems. These can be easily obtained by recalling thatthe integrals in the analytical expressions were obtainedby taking the continuos limit of sums over the allowedmomenta. Therefore we only have to make the substitu-tion (cid:82) dk → (2 π/L ) (cid:80) k n in the corresponding equations.Concretely, we have, the following results for t z ( t ) and t zz ( (cid:96), t ) t z ( t ) = 12 − L (cid:88) n (cid:104) sin (cid:0) (cid:15) g f ( k n ) t (cid:1) sin (cid:18) ∆ k n + φ k n (cid:19) + sin (cid:18) φ k n (cid:19) cos (cid:0) (cid:15) g f ( k n ) t (cid:1) (cid:105) ,t zz ( (cid:96) ; t ) = − L (cid:88) n (cid:88) n (cid:48) (cid:40) e i ( k n − k n (cid:48) ) (cid:96) (cid:104) sin(2 (cid:15) g f ( k n ) t ) sin(∆ k n ) sin(2 (cid:15) g f ( k n (cid:48) ) t ) sin(∆ k n (cid:48) ) + (cid:16) e i ( k n − k n (cid:48) ) (cid:96) − (cid:17) × e i ( θ kn + θ kn (cid:48) ) (cid:16) cos(∆ k n ) − i sin(∆ k n ) cos(2 (cid:15) g f ( k n ) t ) (cid:17)(cid:16) cos(∆ k n (cid:48) ) − i sin(∆ k n (cid:48) ) cos(2 (cid:15) g f ( k n (cid:48) ) t ) (cid:17)(cid:41) , where the set of allowed momenta is, k n = 2 πnL , n = 0 , . . . , L − , (S110)while the expressions for t xx ( (cid:96), t ) and t yy ( (cid:96), t ) can be com-puted using (S94) and (S99), where f η and g η are givenby the following equations in a finite size system, f η = 1 L (cid:88) n e ik n η f ( k n ) , (S111) g η = 1 L (cid:88) n e ik n η g ( k n ) ..