Emergent conservation laws and nonthermal states in the mixed-field Ising model
EEmergent conservation laws and nonthermal states in the mixed-field Ising model
Jonathan Wurtz ∗ and Anatoli Polkovnikov Department of Physics, Boston University, 590 Commonwealth Ave., Boston, MA 02215, USA
This paper presents a method of computing approximate conservation laws and eigenstates ofintegrability-broken models using the concept of adiabatic continuation. Given some Hamiltonian,eigenstates and conserved operators may be computed by using those of a simple Hamiltonian closeby in parameter space, dressed by some unitary rotation. However, most adiabatic continuationanalyses only use this unitary implicitly. In this work, approximate adiabatic gauge potentials areused to construct a state dressing using variational methods, to compute eigenstates via a rotatedtruncated spectrum approximation. These methods allow construction of both low and high-energyapproximate nonthermal eigenstates, as well as quasi-local almost-conserved operators, in modelswhere integrability may be non-perturbatively broken. These concepts will be demonstrated in themixed-field Ising model.
Adiabatic continuation is a well-used concept [1].Given some Hamiltonian of interest, low energy eigen-states of that system may be computed by finding a “sim-ple” Hamiltonian “nearby” in parameter space, whichcan be easily diagonalized. Then, those simple eigen-states can be mapped to interacting ones using somedressing by a unitary U . This procedure can be used toshow that ground states can be mapped to other groundstates, as long as there is some path in parameter spacewhich does not go through some gap-closing critical point[2–4]. For example, low energies of interacting fermionsmay be described by a Fermi gas or Fermi liquid withdressed quasiparticle excitations [5]. However, most ofthese proofs are non-constructive in the sense that theunitary U is generally never actually computed, or onlydone so perturbatively.This work details explicitly connecting non-interactingtrivial states to nontrivial interacting ones, by variation-ally computing an approximate, local adiabatic gaugepotential (AGP) to construct the unitary transforma-tion. The latter modifies non-interacting particle statesto quasiparticle states which are “dressed” within somelocal span of sites. Importantly, the dressing is non-perturbative and not limited to low-energy states. Thisconstruction leads to long-lasting quasiparticles and sta-ble nonthermal states, even at finite energy densities.As a consequence of computing unitary rotations forapproximate eigenstates, this procedure also allows oneto compute local almost-conserved operators from the ap-proximate eigenstates and dressed non-interacting sym-metries. In the presence of integrability breaking termsin some system, undressed symmetries and conserved op-erators are generally no longer conserved [6–8]: instead,these operators may be “dressed” locally by the unitary U to restore approximate conservation of some quasi-local and long lived operator, even though the full systemmay no longer be integrable [9–11]. Similarly, for particu-lar initial wavefunctions with large overlap with good ap-proximate eigenstates, one may compute long time effec-tive quench dynamics within a small subspace of states,in spirit of Schrieffer-Wolff transformations [12].The presence of such unitary dressings, quasi-local con-servation laws, and good approximate eigenstates in in- teracting models may suggest that not all integrability-broken models should be treated equally. Certain modelsmay be “close to integrable”, in the sense that there existsa good local dressing of particular eigenstates, and strongETH may be violated. Two particular cases of such ETHviolation are many body localization [13] and quantumscars in the PXP model [14; 15]. Conversely, other mod-els may be far from any simple system, in the sense thatthere is no path in parameter space that admits a good lo-cal adiabatic dressing, and the model is quantum chaotic[16]. In fact, this unitary rotation may potentially beseen as the analog of the canonical transformation ofKAM theory which restores integrability in classical mod-els [17]. While this paper focuses on a quantum model,the whole methodology, including variationally computedcanonical transformations generated by the AGP [18], isfully applicable to classical non-integrable systems.The rest of this paper is structured as follows. Firstwill describe the methods of computing approximateeigenstates with the variational adiabatic gauge potentialand the rotated truncated spectrum approach (rTSA).The method is a different perspective of Schrieffer-Wolffblock diagonalization methods, as discussed in Ref. [12].Next will introduce the specific model used in this work,the non-integrable mixed-field Ising spin 1/2 chain, anddemonstrate the performance of computing approximateeigenstates. Finally will be an example of an approxi-mately conserved local operator for the mixed-field Isingchain, the total quasiparticle number. I. COMPUTING LOCAL APPROXIMATEEIGENSTATES
In general, directly computing eigenstates and quasi-particle excitations is hard, as computational difficultyscales exponentially in system size, and normally thereare no well-defined quantum numbers. Instead, let uscontinue in spirit of adiabatic continuation [1]. Supposesome parameterized Hamiltonian H ( µ ), with µ = 1 be-ing the particular system of interest, and µ = 0 being ex-actly solvable, in the sense that the eigenstates are easilycomputable via symmetry, integrability, or other means. a r X i v : . [ c ond - m a t . s t r- e l ] F e b µ describes some choice of path in a (multi)-parameterspace of Hamiltonians between 0 and 1. H (0) = ⇒ H ( µ ) = ⇒ H (1)Path from 0 to 1Exactly Solvable ⇓ System of InterestAGP | E n (0) (cid:105) = ⇒ A ( µ ) = ⇒ | E n (1) (cid:105) At all points along this path, there are parameterizedeigenstates {| E n ( µ ) (cid:105)} with eigenenergies { E n ( µ ) } . Howthese eigenstates change as a function of parameter µ isgiven by the adiabatic gauge potential (AGP) A ( µ ) [18],computed from the instantaneous Hamiltonian H ( µ ) i∂ µ | E n ( µ ) (cid:105) = A ( µ ) | E n ( µ ) (cid:105) . (1)In principle, given the (parameter dependent) AGP A ( µ ), one could use Eq. (1) to evolve the simple eigen-states | E n (0) (cid:105) into exact interacting eigenstates | E n (1) (cid:105) ,or equivalently use a unitary rotation U † = T exp (cid:18) i (cid:90) A ( µ ) dµ (cid:19) (2)such that | E n (1) (cid:105) = U | E n (0) (cid:105) , with T indicating pathordering. In practice, this recipe runs into the same in-feasibility as directly computing eigenstates: computingthe exact AGP is generally as difficult as computing theexact eigensystem. It is nonlocal, exponentially large,and highly parameter dependent [18] for generic inter-acting systems, making computing the unitary a similarlydifficult task.To circumvent the problems of computing the ex-act AGP, one can instead use approximate adiabaticgauge potentials. One might hope that simple eigen-states “dressed” by a local approximation of the AGPwill closely resemble eigenstates of the full system. Asthe complexity of the approximation grows, so too shouldthe approximate AGP approach the exact one, and theapproximate eigenstates become exact, at the expense ofthem being highly entangled and nonlocal.From an adiabatic continuation standpoint, the con-nection with this approximate gauge potential is clear;in fact M. Hastings in [1] has is a particular implemen-tation of an approximate AGP for gapped ground states.A simplified version of Eq. (17) in Ref. [1] is equivalentto Eq. (2) with A ( µ ) ≈ − (cid:90) ∞−∞ SGN[ t ] f ( t )[ ∂ µ H ]( t ) dt, (3) f ( t ) = erfc (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) tτ q (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (4) where ∂ µ H ( t ) is the operator ∂ µ H in the Heisenberg rep-resentation, erfc is the complementary error function forwhich erfc(0) = 1 and erfc( ∞ ) = 0, and SGN( t ) is ± t (see Appendix A for more de-tails). This expression is nothing but an approximationof the gauge potential with a particular choice of regular-izer f ( t ) [19]. For the regularization time τ q → ∞ , thisapproximate AGP becomes exact. For a finite regular-ization time, the AGP is approximate but local withinsome span of sites.Instead of computing an approximate gauge potentialvia some choice of regulator, which is computationallydifficult, one can instead compute it variationally [20].Here, some ansatz for the variational gauge potential ischosen based on the system at hand, then variationallyoptimized to best approximate the exact AGP. In thiscase the variational AGP is chosen to be the sum of someset of (local) operators { B i } for variational parameters { α i } A ( { α } ) = (cid:88) i α i B i . (5)The particular choice of operators { B i } depends on theproblem at hand; if the set is expanded to include all op-erators in the space modulo symmetries (a potentially ex-ponential number), the exact AGP can be reconstructedfrom the set, and the variational version becomes exact.Thus, one should expect improvement as the size of theansatz is expanded. It has been found [19] that a localvariational AGP can accurately reproduce the exact one,at least for matrix elements which are far separated inenergy or protected by approximate symmetries. Thisis shown by using a Taylor series expansion of Eq. (3).In fact, this argument is mirrored by an equivalent jus-tification from an adiabatic continuation perspective forgapped ground states: As long as there is a gap, theground state will only entangle within some light cone[2; 3; 21] and thus a local approximation for the rotationgenerator is always possible.The optimization of the AGP is computed by findingthe minimum of [18; 20] S ( α ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:2) H, ∂ µ H + i [ A ( α ) , H ] (cid:3)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (6)where parameter dependence of α and H on µ is implicit,and || Q || = Tr (cid:2) Q (cid:3) / D is the Hilbert-Schmidt norm. Eq.(6) has the property that S = 0 for the exact AGP, and issolvable even for very large system sizes. This is becausethe function S ( α ) can be computed using trace identities,and is quadratic in α , making optimization simple. Thisminimization leads to a quasi-local operator A ( µ ) whichapproximates the exact AGP.One can then use this approximate AGP to computeapproximate eigenstates, in the same manner of the exactAGP computing exact eigenstates. First, choose someset of D p states {| q (cid:105)} of the exactly solvable Hamiltonian H (0). This choice depends on the system at hand. Onechoice could be, for example, all eigenstates below someenergy cutoff. Another choice would be all states withinsome particular symmetry sector such as fixed particlenumber, which form a subset of eigenstates not neces-sarily sorted by energy. The set of states forms someprojective subspace P with projector P = (cid:80) q | q (cid:105)(cid:104) q | .Approximate eigenstates of the interacting model H (1) can be computed by “dressing” each state viaSchr¨odinger evolution of Eq. (1) to implement the uni-tary of Eq. (2). This gives some set of “dressed” states {| q (1) (cid:105)} ≡ { U | q (cid:105)} .As a comment on implementation, care must be takenin the direction of evolution: the AGP is parameter de-pendent so generally A (1) (cid:54) = A (0). For perturbative cou-plings this parameter dependence is very weak so that theAGP approximately commutes with itself for all µ andthe directionality doesn’t matter; however strong cou-pling may lead to nonsensical answers. One may startby accidentally acting on a non-interacting wavefunctionwith the AGP from the interacting point, which may bemuch different than the correct non-interacting AGP.An effective Hamiltonian within that subspace can becomputed via matrix elements in the dressed subspace: (cid:0) H eff (cid:1) pq = (cid:104) p (1) | H (1) | q (1) (cid:105) = (cid:104) p | U † H (1) U | q (cid:105) . (7)Note that this is equivalent to computing the effectiveSchrieffer-Wolff Hamiltonian [12], where one rotates theoperator ˜ H = U † H (1) U instead of the states. If theAGP is exact, this effective Hamiltonian will be exactlydiagonal. However, if the AGP is not exact or the initialsubspace was a degenerate symmetry sector, the effectiveHamiltonian will not be diagonal. One can then computethe eigensystem of the D p ×D p matrix via standard linearalgebra techniques to find eigenenergies E i and eigenvec-tors V i such that (cid:0) H eff (cid:1) nm V mi = E i V ni . (8)The approximate eigenvectors of the system are thengiven by | E i (cid:105) = (cid:88) q V qi | q (1) (cid:105) (9)with eigenvalues E i . This final step is functionally equiv-alent to the truncated spectrum Approach (TSA) [22],except instead of using a non-interacting subspace P ,the subspace is first rotated by the approximate AGPto obtain some improved subspace ˜ P . This basis bet-ter resembles eigenvectors of the interacting system, andmay be exponentially orthogonal from the original basisdue to the finite rotation. This corresponds to a rotatedtruncated spectrum approach (rTSA) The abbreviated method is as follows:1. Define some Hamiltonian H ( µ ), with H (0) being ex-actly solvable and H (1) being a system of interest, withsome path in parameter space linking the two.2. Given some ansatz, compute a variational adiabaticgauge potential A ( µ ) along the points µ ∈ [0 , H (0), either withinsome energy window or within some symmetry sector(s)such as particle number.4. Evolve the set of states via the Schr¨odinger equationfrom µ = 0 to µ = 1 with the variational AGP.5. Compute the effective Hamiltonian and its eigensystemto find approximate eigenstates and eigenvalues. Because the variational AGP is local by construction,various tricks can be employed to go to large or eventhermodynamic system sizes. As commented above, com-puting the variational AGP is not a problem for a largenumber of sites, as complexity scales linearly with systemsize. It is reasonable to compute an AGP for hundredsof sites for all operators spanning less than 5-6 sites on amodern desktop. This locality can also be used for step 4when evolving the basis states: due to the finite evolution“time”, states are only entangled within some finite re-gion (using intuition of Lieb-Robinson bounds [3]). Thissuppresses the finite size effects of evolving some small(typically 15-20) number of sites exactly, and enables ten-sor methods such as matrix product states (MPS) andvariational evolution [23]. This work employs the formermethod of exact evolution on small systems [24].Choosing a larger subspace should also be expected toimprove the computation of the approximate eigensys-tem. In the limit where the projective subspace is thefull Hilbert space or within one of the symmetries of thefull Hamiltonian, the effective Hamiltonian is the exactone and likewise the eigenstates are exact, independentof the rotation. Choosing a subspace within some largerenergy window should also be expected to improve thevariational dressing: the AGP fails to suppress excita-tions close together in energy, but those can then be re-captured within the subspace P via off-diagonal elementsof the effective Hamiltonian.Because both the variational AGP ansatz and the pro-jective subspace can be systematically expanded, thismethod gives a controllable approximation to computeeigenstates: as the complexity increases, the eigenstateswill asymptotically approach the exact ones.The eigenstates computed in this manner are approxi-mate, in that they are not exact eigenstates of the Hamil-tonian H (1). The simplest indicator of the closeness toan exact eigenstate is the energy variance of the state∆ n ≡ (cid:104) E n | H | E n (cid:105) − (cid:12)(cid:12) (cid:104) E n | H | E n (cid:105) (cid:12)(cid:12) . (10)Exact eigenstates have zero energy variance, and so ap-proximate eigenstates should have minimal energy vari-ance ∆ n ≈
0. The average energy variance of theseeigenstates within the subblock corresponds to the av-erage block-off-diagonal matrix elements in the Hamil-tonian and thus indicates the performance of the blockdiagonalization procedure (see Appendix B).
II. MODEL
As a concrete example, suppose the following system,the mixed field Ising model H = N (cid:88) i Jσ iz σ i +1 z + h x σ ix + h z σ iz . (11)For h z = 0 the model is integrable via a Jordan-Wignertransformation to free fermions [25–27] with a criticalpoint at h x = 1 and small h z being an integrable E(8)field theory [28]. For h x = 0 the model is a purely classi-cal Ising model. For h x + h z → ∞ the model is an exactlysolvable collection of single spins with an onsite field. At h z = 2 , h x = 0 there is a first-order multicritical point[29] and for small h x , the low energy effective Hamilto-nian is the PXP model [14; 30]. Elsewhere, the modelhas no apparent conservation laws or symmetries beyondgeometric ones and is generally quantum chaotic [31; 32].However, this does not prevent approximate conservationlaws or nonthermal states, as will indeed be seen.The variational ansatz is chosen to be that of Jordan-Wigner strings, i.e. strings of Pauli operators which mapto fermion bilinear operators, plus all operators localwithin a span of n sites { B } = { σ y , σ x σ y , . . . , σ x σ ny , . . . , σ y σ z σ x , . . . , (12) σ z σ y , σ z σ x σ y , σ z σ x σ x σ y , σ z σ x σ x σ x σ y . . . } . Additional symmetries and properties reduce the sizeof the ansatz: the AGP has all of the symmetries of thefull Hamiltonian [33]. By gauge choice the AGP canbe completely imaginary for real Hamiltonians [18], con-straining { B } to only include terms with an odd numberof σ y . Because the Hamiltonian is translation and reflec-tion invariant, the ansatz can be chosen to be as well.The inclusion of Jordan-Wigner strings is motivated bythis ansatz being exact for the transverse Ising model[18], due to its extra symmetries and mapping to freefermions.While a Hamiltonian of interest is given by particularchoice of parameters h x , h z , there is relative freedomfor choice of the path h x ( µ ), h z ( µ ) in the 2d parameterspace H ( µ ), and especially choice of simple Hamiltonian H (0). This is because there are many “simple” points inthe ( h x , h z ) parameter space which might be considered“close” to the Hamiltonian of interest. The ( h x ,
0) line isthe transverse Ising model; the (0 , h z ) line is the classical Ising model; and the ( h x , h z ) → ∞ line are independentspins with onsite fields.What starting points, and which path in parameterspace, is optimal for computing approximate eigenstates,given ansatz { B } , Hamiltonian H (1), and subspace P?This is a question of a path-dependent Schrieffer Wolfftransformation, as the performance of computing approx-imate eigenstates, or equivalently block diagonalization,may depend on these choices. This work chooses from alimited set of parameterized Hamiltonians with particu-lar starting and ending points. Two additional parame-terizations are discussed in Appendix C. H ( µ ) = N (cid:88) i σ iz σ i +1 z + µ (cid:0) h x σ ix + h z σ iz (cid:1) , (13) H ( µ ) = N (cid:88) i − σ iz σ i +1 z + µ (cid:0) h x σ ix + h z σ iz (cid:1) , (14) H ( µ ) = (cid:80) Ni µh x σ ix + h z σ iz , µ ∈ [0 , . (cid:80) Ni (2 µ − σ iz σ i +1 z µ ∈ [0 . , (cid:0) h x σ ix + h z σ iz (cid:1) . (15)The first and second parameterizations start from the σ z σ z point, whose eigenstates are Z polarized spins. De-pending on the sign, the ground state could be an anti-ferromagnetic (AFM) Ne´el (13) or a polarized ferromag-netic (FM) state (14). Low energy particle excitationsare boundary walls of spin flips (see Fig. 1) [34].The third parameterization (15) is split into two parts.The first leg is simply rotating the on-site field and thusthe AGP is exact A ( µ ) ∼ (cid:80) i σ y , and is an example of theLandau-Zener problem, rotating the spin in the XZ plane.The second leg has no such local exact representation.The ground state is a product state of spins pointing inZ. Low energy particle excitations are spin flips (see Fig.1) to the opposite direction.In all cases, H ∗ (0) is degenerate, with a natural choiceof projective subspace being fixed particle number on topof the ground state. Thus, P is 0 and 2 boundary wallson top of an AFM ground state; P is 0 and 2 boundarywalls on top of a FM state; and P is the 0, 1 and 2particle spin flips on top of a polarized state. P = (cid:26) | ↑↓ . . . ↑↓(cid:105) , (cid:0) σ ix σ i +1 x . . . σ i + nx (cid:1) | ↑↓ . . . ↑↓(cid:105) (cid:27) ,P = (cid:26) | ↑↑ . . . ↑↑(cid:105) , (cid:0) σ ix σ i +1 x . . . σ i + nx (cid:1) | ↑↑ . . . ↑↑(cid:105) (cid:27) ,P = (cid:26)(cid:12)(cid:12) ↓↓ . . . ↓↓(cid:105) , (cid:0) σ ix (cid:1) | ↓↓ . . . ↓↓(cid:105) , (cid:0) σ ix σ jx (cid:1) | ↓↓ . . . ↓↓(cid:105) (cid:27) . (16) FIG. 1. Example basis states of two particles separated by6 sites. Top and middle are states with two boundary walls(red dashes), which are low energy eigenstates of H (0) and H (0), with excitation energy 4 J . Bottom is a state with twospin flip particles, which is a low-energy state of H (0) withexcitation energy 4 h z . Because the system is translation invariant, the zero-momentum sector is chosen as a numerical simplification.Under these constraints, each subspace P has N +1 stateseach out of total Hilbert space dimension ≈ N /N .These basis states are each dressed by the variationalAGP to create dressed boundary wall states: the hardboundary is softened by the dressing procedure to betterdescribe the interacting quasiparticle excitations. III. APPROXIMATE EIGENSTATES ANDSPECTRUM
As an explicit example, let us choose the parameters h x = 0 . h z , and coupling J = ±
1. These parame-ters are non-perturbative, in the sense 0 . O (1) awayfrom any simple point. For J = −
1, the ground stateis ferromagnetic, and the h z term acts as a constant at-tractive force between two boundary wall particles. Thisleads to “meson” bound states of the two boundary walls[35]. For J = +1, the h z term does not change the AFMground state energy [36; 37]; however it will affect en-ergies of spin flips on up/down Ne´el sublattices, which,from a band theory context, leads to two free particlespecies.One can then go through the process of computing ap-proximate eigenstates, as outlined in section I, for theseparticular choices of subspaces. Here, Hamiltonians (13),(14) are chosen starting from the FM and AFM sub-spaces, with an ansatz of all operators local to 3 sitesplus Jordan Wigner strings, with 18 total sites in the 0momentum sector. The form of the AGP is shown inAppendix E.Results for these parameters are shown in Fig. 2. Topplots the effective unrotated and rotated Hamiltonian,or equivalently the Hamiltonian in the projective androtated projective subspaces, for AFM (left) and FM(right) excitations. It can be clearly seen that the rotatedeffective Hamiltonian becomes slightly more nonlocal: adressed boundary wall of width 3 may hop to become width 5, for example. These effects are especially pro-nounced when the two boundary walls are close together,which is an indicator of a 2-particle interaction. Whenthe two particles are far apart the Hamiltonian becomesindependent of distance.The spectrum is shown on the bottom plots of Fig. 2.Clearly, there is remarkable improvement over na¨ıve TSA(red) with the unrotated basis, and the rotated version(blue) is almost identical to the exactly computed eigen-spectrum (black). The error bars are the energy varianceof the approximate eigenstates, as computed from Eq.(10). Note that the exact eigenvectors are matched withapproximate ones by choosing those which have maxi-mum fidelity |(cid:104) E n | E exact m (cid:105)| ; normally this value is > . . J ) has a higher energy then a single FM boundary wallof width 6 (with excitation energy 8 . J ). This means itis energetically possible for the width-6 boundary wallto decay into two width-1 boundary walls, for exam-ple. However, the small energy variance of these dressedstates indicates that such a process is almost completelysuppressed.Because the dressing is local, it is possible to take acontinuum or large system size limit. Numerically, thisis done by duplicating the dressed Hamiltonians of Fig.2 over thousands of sites: The 19 ×
19 matrix is extendedto a N × N matrix, where the middle elements are theduplicated middle elements of the smaller matrix. Then,the eigensystem of that Hamiltonian is computed. Re-sults for the continuum dispersion relation of excitationson top of the AFM ground state are shown in Fig. 3.There are two particle species which have mass 1 .
11 and2 .
56. Like the meson case, the states are not necessarilythe lowest-energy states: for example, 2 heavy particlescan have equal energy to 4 light particles, and may poten-tially decay as such. Note that two light particles couldnot decay into 1 heavy particle, as that is disallowed bythe particles being domain walls.These approximate eigenstates can be compared withthe general bulk eigenstates, as is shown in Fig. 4Top. Here, all eigenstates in the 18-site FM model, 0-momentum sector (14,602 total) are computed exactly,and their half-cut entanglement entropy is found. Ther-mal states are extensively entangled states at finite en-ergy density, while nonthermal states are weakly entan-gled and generally have zero or low energy density [30].Red circles are the states which have maximal overlapwith the approximate eigenstates: low energies are FMstates while high energies are AFM states.One particular approximate eigenstate merits morestudy: the dressed all-up eigenstate, indicated by thegreen circle and arrow in Fig. 4 Top. This is a groundstate of σ z σ z , but the most excited state of σ z ; it hasfinite energy density given roughly by 2 h z . But, in par-ticular, it is an explicit example of a highly non-thermal Unrotated
Rotated
Unrotated
Rotated E i g e n e n e r g y Exact SpectrumUnrotated SpectrumRotated Spectrum 0 5 10 15Eigenindex05101520 E i g e n e n e r g y Exact SpectrumUnrotated SpectrumRotated Spectrum
FIG. 2. Eigenspectrum of the zero momentum mixed-field Ising model of Eq. (11). Left is for 18-site AFM states ( J = +1)while right is for 18-site FM states ( J = − /2 0 /2Momentum1.251.501.752.002.252.502.75 P a r t i c l e E n e r g y FIG. 3. Dispersion relation of the two particle species on topof an AFM ground state. The 2-heavy particle energy liesabove the 4-light particle continuum. state far from the edges of the spectrum [13; 30; 38; 39].It is locally entangled with a half-chain entanglement en-tropy of ≈ .
25 bits. It has very high fidelity of 0 .
995 with an exact eigenstate. Note that in the thermody-namic limit the dressed-all-up state is exponentially or-thogonal to the original all-up state due to the finite localrotation.
A. Quasiparticle Lifetimes
The energy variance of these approximate eigenstatestakes special meaning when they can be interpreted asdressed particles. In this case, the energy variance givesa lower bound on the quasiparticle lifetime. For an ap-proximate particle eigenstate | E n (cid:105) , the time-dependentstate overlap under second order perturbation theory is (cid:12)(cid:12) (cid:104) E n | E n ( t ) (cid:105) (cid:12)(cid:12) ≈ − t τ + O ( t ) , (17)where τ − = ∆ = (cid:104) H (cid:105) − (cid:104) H (cid:105) is the energy vari-ance of the state. In other words, the characteristictime for an (eigen)state of some particles to decay intosome other particle state is given by the energy vari-ance. This timescale is very crude as it assumes all otherstates have the same energy: a more refined timescale | ( t ) | ( ) | RotatedUnrotated
FIG. 4.
Top:
Comparison of half-cut eigenstate entangle-ment entropy for the 18-site FM chain and 0 momentum.Red circles are the states with maximum overlap with theapproximate eigenstates; left are dressed ferromagnetic stateswhile right are dressed anti-ferromagnetic states. Note thathigh-energy states of the FM model are not the same as thelow-energy states of the AFM model. Green circle and arrowindicates the dressed all-up state, which is a nonthermal state.
Bottom:
Fidelity of the dressed-all-up state with the initialstate |(cid:104) E ( t ) | E (0) (cid:105)| which indicates the rotated all-up stateis very close to an eigenstate, while the unrotated version isnot, and is well preserved in time. can be computed using the Fermi golden rule [25] for thedressed states, but is not generally possible without a pri-ori knowledge of the energy of the other states. As such,the energy variance serves as a lowest bound on (inverse)quasiparticle lifetime.As an explicit example of these timescales, a dressedsingle flipped down spin on a FM ground state, whichcorresponds to the lowest energy meson excitation, hasa characteristic lifetime of τ = 110, far longer than anylocal timescale. Excitations on top of an AFM groundstate have lifetimes in excess of τ >
60. The dressedall-up state has a lifetime of τ = 53. These lifetimes arelonger when adding more parameters to the variationalAGP. Explicit time dynamics of Eq. (17) for this dressedall-up state is shown in Fig. 4 Bottom. Clearly, it is muchcloser to an eigenstate than expected, as it is close to 1 atall times. This further indicates the genuineness of this nonthermal eigenstate, especially when compared to theundressed version of the same. Note that the undressedup state is exponentially orthogonal in system size fromthe dressed up state, due to the finite rotation.One application of such a dressed all-up state is for in-formation protection in quantum systems. For a classicalHamiltonian, the all-down [ground] state may be labeledas a logical 0, while the all-up state is labeled as a logi-cal 1. An X-field will generally change these two states,destroying the encoded bit. If this bit is instead encodedin the dressed nonthermal states (the all down groundstate, and the all up nonthermal state), they are muchmore stable, encoding the information for a much longertime by suppressing transitions. B. Quasiparticle parameter dependence
A general measure of the quasiparticle lifetime within aparticular subspace is given by their normalized averageinverse lifetime, or equivalently average energy varianceΓ = 1 (cid:112) D p (cid:112) h x + h z (cid:115)(cid:88) n ∆ n . (18)Γ equivalently is the block-off-diagonal weight of therotated subspace ˜ P (See Appendix B). A small value in-dicates a good block-diagonalization procedure with well-defined quasiparticles within the subspace. A large valueindicates a failure to block diagonalize the Hamiltonian.This error can be computed for various values of h x and h z by evolving with parameterized Hamiltonians (13) and(15), computing approximate eigenvalues, then comput-ing their normalized average energy variance Γ( h x , h z ).Results are shown in Fig. 5, for a 3-site ansatz and 14sites. States are dressed from one of two directions. Oneis H ( µ ), dressing 2-particle AFM boundary wall statesout from the σ z σ z only point, indicated by the radial ar-rows in the bottom left. The other is H ( µ ), dressing1 and 2-particle spin-flip states from the h x σ x + h z σ z only point(s), indicated by the arrows pointing radiallyinwards.In the region where h x , h z is small, the error fromdressing boundary walls is enormously low. With nodressing, the error grows linearly in | h | , while with dress-ing, the error grows sub-linearly, which indicates that thedressing is exact asymptotically. In fact, this dressingaccumulates very small errors even for non-perturbativevalues of | h | , as shown in Fig. 5 Bottom, which is dressingalong the h x = h z line. This indicates that in the whiteareas, there is a good effective quasiparticle descriptionof the low energies of this otherwise quantum chaoticmodel, described by dressed boundary wall particles.Although it is not generally so, the error accumulatesmonotonically with increasing | h | . This means that atsome critical value, the dressing going outwards from h From Z point
Null AnsatzSize-2 AnsatzSize-3 AnsatzSize-4 Ansatz0.0 0.5 1.0 1.5 h I m p r o v e m e n t From ZZ point h A v e r a g e E n e r g y V a r i a n c e FIG. 5. Results for a 3-site VGP ansatz for the 14-site TLIModel.
Top:
Average inverse 2-particle lifetime Γ or equiv-alently average energy variance. Blue line indicates transi-tion between different directions. Star is the (0 . , .
4) pointstudied more in-depth. Arrows indicate direction of dressing.White and yellow indicate areas with a good quasiparticledescription.
Bottom:
Average energy variance in the di-rection π/ h ≈ .
6, the error is vanishinglysmall.
Middle:
Energy variance improvement compared tothe Null ansatz Γ / Γ along the h x = h z line. As the ansatzsize increases, so too does the improvement, as expected. the σ z σ z point will have a larger error then the dress-ing going inwards from the ( h x , h z ) → ∞ point. At thisboundary, the best description of quasiparticles changesfrom dressed pairs of boundary walls , to dressed spin flips .This does not mean that there is no effective descriptionof certain states in terms of quasiparticles: there couldbe some other subspace (say, of doubly-flipped spins) andother path through parameter space which gives a bet-ter quasiparticle picture in that the energy variance issmaller. This cross-over point may be an indicator of an inter-acting phase transition. Around h z = 0 , h x = 1, which isthe transverse Ising phase transition, it has been foundthat local variational adiabatic dressing begins to fail[18]. This finding is now extended to the interacting case:the cross-over gives a rough region where the interactingcritical point may occur, as a local AGP fails to reproducethe long-range entanglement of a critical ground state [3].With increasing ansatz size, this point decreases in totalerror, and shifts in critical parameter (see Fig. 5 Bot-tom) which may eventually converge to some particularvalue, indicating the interacting critical point. This ideais backed up by the convergence in the non-interactinglimit: For h z small, the crossover is around the h x ≈ h x small, thecrossover is around h z ≈
2, which is the first-order phasetransition in h z to change the ground state from AFM topolarized [29]. IV. LOCAL ALMOST-CONSERVEDOPERATORS
Given some set of approximate local eigenvectors {| E n (cid:105)} generated by this adiabatic dressing scheme, it isa relatively simple procedure to construct approximatelyconserved local operators. An operator is conserved if itcommutes with the Hamiltonian, or equivalently if it isconstructed from eigenstates of the Hamiltonian. Givenapproximate eigenstates, one should then be able to com-pute approximately conserved quantities. Conserved op-erators are defined as O = (cid:88) n O n | E n (cid:105)(cid:104) E n | (19)where O n are the eigenstates of the operator, and {| E n (cid:105)} are exact eigenstates. An operator of this form has theproperty that the symmetric time correlation function isconserved (for all initial states) (cid:104){O ( t ) , O (0) }(cid:105) . (20)There are two ways to construct approximate versionsof these operators. The first way is to explicitly use Eq.(19) using only the subspace of dressed eigenstates whichwere directly computed. In this case, the sum is of size D P , the subspace size, as opposed to D , the total Hilbertspace size. Due to the global projective structure of par-ticle excitations on top of a ground state, the resultingoperator is not necessary local. However, this may be im-plemented with some local operator plus post-selection ofstates.In the case of such an operator directly constructedfrom approximate eigenstates, the symmetric correla-tion function is not conserved in time. Under pertur-bation theory, the characteristic timescale is given by a P a r t i c l e nu m b e r c o rr e l a t o r Antiferromagnetic
Infinite Temperature
Ferromagnetic
UnrotatedRotated
FIG. 6. Symmetric time correlation function of dressed (blue) and undressed (red) particle number of Eq. (25) and Eq. (24)for a 2 particle AFM subspace (left), 2 particle FM subspace (right) and all states (middle), and 14 sites. Dashed are theinfinite temperature long-time values. weighted-average energy variance (see Appendix D forderivation) (cid:104){ O ( t ) , O (0) }(cid:105) (cid:104) O (cid:105) (cid:18) − t τ O (cid:19) , τ O ≡ (cid:80) n ρ n O n ∆ n (cid:80) n ρ n O n . (21)Here, ρ n = (cid:104) E n | ρ | E n (cid:105) is the density matrix for theexpectation value, assumed to be diagonal, and ∆ n isthe energy variance of the n th approximate eigenstate.For good eigenstates with low energy variance, the decaytimescale can be very long.The second way to construct approximately conservedoperators is to dress conservation laws of the simple sys-tem H (0) with the unitary. Conservation laws, suchas particle number and particle current, are constructedfrom simple eigenstates in the form of Eq. (19), withparticular choice of O n , and are generally local [26; 40].Then, the dressed operator is constructed from statesbetter resembling eigenstates O = U (cid:18) (cid:88) n O n (cid:12)(cid:12) E n (0) (cid:11)(cid:10) E n (0) (cid:12)(cid:12)(cid:19) U † (22) (cid:109)O = U O U † . (23)Importantly, the eigenstates of operator O are not nec-essarily the same as the constructed approximate eigen-states, as it is missing the re-diagonalization step of Eq.(9). The resulting operator is quasi-local, and approxi-mately conserved [9; 41]. As the ansatz span is increased,the AGP approaches the exact one, resulting in better ap-proximate eigenstates and a better-conserved operator,at the expense of it becoming more and more nonlocal.One such conserved operator for the mixed-fieldIsing model is the dressed total particle number N = U (cid:0) (cid:80) i σ iz σ i +1 z (cid:1) U † which (up to a constant) counts thenumber of boundary walls in the system. For H and H , this is also the initial Hamiltonian; thus one wouldexpect that the dressed operator should also approximatedressed versions of particle number eigenstates.For h x = h z = 0 . J = +1, the dressed particlenumber operator becomes N = (cid:88) i σ iz σ i +1 z (24) ⇓ N = (cid:88) i . σ iz ˆ σ i +1 z + (25)0 . σ ix + 0 . σ iz ˆ σ i +1 x ˆ σ i +2 z +0 . σ iz ˆ σ i +1 x ˆ σ i +2 x ˆ σ i +3 z + − . σ iz ˆ σ i +1 x + ˆ σ ix ˆ σ i +1 z )+ − . σ iy ˆ σ i +1 y +0 . σ iz ˆ σ i +1 x ˆ σ i +2 x ˆ σ i +3 x ˆ σ i +4 z + − . σ ix ˆ σ i +1 x ˆ σ i +2 z + ˆ σ iz ˆ σ i +1 x ˆ σ i +2 x )+0 . σ iz ++ . . . where ellipsis represent the more and more nonlocalterms of the operator. As can be seen, this operator is ap-proximately local, with dominant terms coming from 1,2, and 3-spin terms. One can then compute the symmet-ric correlation function in the initial undressed subspaceof two particles to see its conservationTr (cid:2) P { N ( t ) , N (0) } (cid:3) = C ( t ) . (26)Results are shown in Fig. 6 for AFM and FM states, aswell as infinite temperature typical states [42]. For com-parison, the undressed operator is also shown in red. Theconserved operator for AFM states is almost stationary0in time, while the undressed version is not. The infinitetemperature timescale can be computed analytically asTr (cid:2)(cid:0) [ N, H ] (cid:1) (cid:3) (cid:2) N (cid:3) = τ − N = 50 . − , (27)Tr (cid:2)(cid:0) [ N , H ] (cid:1) (cid:3) (cid:2) N (cid:3) =( τ N ) − = 1 . − . (28)Even for an infinite temperature state, this quasilo-cal dressed operator gets a factor of 40 improvement inthe characteristic decay timescale. This indicates thatdressed quasiparticle excitations may persist in this in-teracting model even at infinite temperature, and thatthis model is “closer to integrable” then one might ex-pect. V. CONCLUSION
The existence of good approximate dressings havesome curious implications. Even if a model system isnot necessarily integrable or exactly solvable, that doesnot mean that there are no local long-lived symmetriesand conservation laws. Indeed, if such a model is closeby to an integrable point, a conservation law of the in-tegrable model can be “dressed” by a unitary generatedby the approximate local adiabatic gauge potential to re-store the symmetry approximately in a now quasi-localoperator. Approximate eigenstates may be computed ina similar manner: simple particle excitations of the inte-grable point can be dressed by the approximate AGP toconstruct long-lived quasiparticle excitations of the inter-acting point. These new dressed states need not be lowenergy states and in fact may be used to construct finiteenergy density low-entanglement nonthermal states, asdemonstrated in the dressed-all-up state of the mixed-field Ising model.Similar studies have been done to compute low en-ergy phenomenology of the Meson case, most predom-inantly in recent work by [43; 44] using a truncated spec-trum approach (TSA). These numerical diagonalizationprocedures are functionally equivalent except that herethe projective subspace is first rotated by the variationalAGP, leading to a subspace closer to the exact eigen-states. While this work uses discrete lattices, general-izations to continuous theories is an interesting futuredirection.The restoration of approximate symmetries and con-struction of quasiparticle excitations in interacting mod-els puts a new perspective on integrability breaking. In-stead of reevaluating a Hamiltonian for every new pointin parameter space, one can instead compute propertiesand approximate symmetries based on nearby Hamilto-nians with a potentially simpler structure. This “close-ness” is defined in the sense of being able to computea good approximate AGP along some path between the simple Hamiltonian and interacting one, not in the senseof perturbative parameter changes. Certain perturba-tions away from integrability may rapidly destroy anylocal conservation laws, if there exists no good local ap-proximate AGP. Other perturbations, while still break-ing integrability, may still admit quasi-local conservationlaws, nonthermal states, and quasiparticles, if there doesexist a good local approximate AGP.These unitary rotations restoring approximate integra-bility are similar in spirit to canonical transformations inKAM theory [17]: integrability may be approximately re-stored for particular subsets of initial conditions of par-ticular integrability-broken systems via the unitary ro-tation (eg canonical transformation) of conserved quan-tities. Whether this approach to stability of quantumintegrable systems can be made more concrete remainsto be seen but these variational local dressings may be astep towards a general theory in that direction.
ACKNOWLEDGMENTS
We would like to thank Artem Rakcheev, Sho Sugiuraand Pieter Claeys for useful discussions. J.W. and A.P.were supported by NSF Grants No. DMR-1813499 andNo. AFOSR FA9550-16- 1-0334.1
Appendix A: Hastings 2005 Eq. 17 to Eq. (3)
This section serves as a derivation of Eq. (3) from [1]Eq. 17. Hastings defines the rotation ˜ V ( s ) (analogously U † ) as˜ V ( s ) = S (cid:48) exp (cid:26) − (cid:90) s ds (cid:48) (cid:90) ∞ dτ e − ( τ/τ q ) / [˜ u + s (cid:48) ( iτ ) − H.c.] (cid:27) . (A1) Here, S (cid:48) is parameter ordering (analogously T ) for pa-rameter s (cid:48) (analogously µ ). The object ˜ u + s (cid:48) ( iτ ) is definedas ˜ u ± ( ± iτ ) ≡ π (cid:90) ∞−∞ dt [ ∂ s H s ]( t ) e − ( t/τ q ) / ± it + τ (A2)with [ ∗ ]( t ) denoting time evolution with respect to in-stantaneous Hamiltonian H s (cid:48) (analogously H ( µ )). Sub-stituting this into the inner integrand and simplifying byintegrating over τ yields (cid:90) ∞ dτ e − ( τ/τ q ) / [˜ u + s (cid:48) ( iτ ) − H.c.] = 12 π (cid:90) ∞−∞ dt [ ∂ s H s ]( t ) e − ( t/τ q ) / (cid:90) ∞ dτ e − ( τ/τ q ) / (cid:20) itt + τ (cid:21) = 12 π (cid:90) ∞−∞ dt [ ∂ s H s ]( t ) e − ( t/τ q ) / (cid:20) − iπe ( t/τ q ) erfc (cid:16)(cid:12)(cid:12)(cid:12) t √ τ q (cid:12)(cid:12)(cid:12)(cid:17) SGN( t ) (cid:21) = − i (cid:90) ∞−∞ dt [ ∂ s H s ]( t )erfc (cid:16)(cid:12)(cid:12)(cid:12) t √ τ q (cid:12)(cid:12)(cid:12)(cid:17) SGN( t ) . Up to a trivial factor of √ τ q , this is the expression in Eq.(3). Appendix B: Off-diagonal matrix elements andenergy variance
The average energy variance of approximate eigen-states analogously gives the average block-off-diagonalelements in the Hamiltonian, as claimed in section III B.This Appendix serves to elaborate on this point.When computing an effective Hamiltonian within a ro-tated subspace, the procedure is an analogous one to aSchrieffer-Wolff transformation: a unitary rotation blockdiagonalizes some Hamiltonian into a subspace P andcomplement Q . A measure of the quality of this diago-nalization is the average strength of the off-diagonal ele-ments: zero strength means exact block diagonalization,while nonzero strength means approximate diagonaliza-tion. The average energy variance is defined asΓ = 1 D P p (cid:88) n (cid:104) E n | H | E n (cid:105) − (cid:16) (cid:104) E n | H | E n (cid:105) (cid:17) , (B1)Γ = 1 D P p (cid:88) n (cid:104) E n | H (cid:16) | q (cid:105)(cid:104) q | + | E p (cid:48) (cid:105)(cid:104) E p (cid:48) | (cid:17) H | E p (cid:105)− (cid:16) (cid:104) E n | H | E n (cid:105) (cid:17) , (B2)Γ = 1 D P (cid:88) nq (cid:12)(cid:12)(cid:12) (cid:104) E n | H | q (cid:105) (cid:12)(cid:12)(cid:12) . (B3) A v e r a g e E n e r g y V a r i a n c e FIG. 7. Three paths to get to the (0.4,0.4) point: AlongX then along Z (blue); along Z then along X (orange) anddiagonally (green). Energy variance is along each point onthe path; the paths turn at parameter value 0 . Step 2 inserts the identity, for complete set of states | q (cid:105) ∈ Q , and complete set of states | E p (cid:105) ∈ P , while step3 simplifies using the fact that | E n (cid:105) are eigenstates of theeffective Hamiltonian within subspace P . Because | E n (cid:105) is a complete set of states in P and similarly for Q , thesum is then over all off-block-diagonal matrix elements,giving an average off diagonal strength. Appendix C: Path Dependence
While this work chooses one particular path to the(0.4,0.4) point, there exist many other options, whichmay give better or worse performance. Here two addi-2tional paths are compared to the diagonal path. Thesecond goes up in h z , for which the gauge potential iszero, and then right in h x . The third goes right in h x asthe transverse Ising model for which the gauge potentialis exact, and then up in h z (See inset in Fig. 7). Thisdirectionality is shown in Fig. 7 for 14 sites, and a size-3ansatz. The diagonal path ends with an average energyvariance of ≈ .
02, while the other paths have varianceof ≈ .
04, a factor of two improvement.
Appendix D: Symmetric correlators and almostconserved quantities
This appendix serves as a derivation of Eq. (21). Sup-pose some set of D P approximate eigenstates | E n (cid:105) di-agonalized within some subspace P , and operator O = (cid:80) n O n | E n (cid:105)(cid:104) E n | for Hamiltonian H . For simplicity, letus choose a subspace [ ρ, O ] = 0 or equivalently ρ = (cid:80) n ρ n | E n (cid:105)(cid:104) E n | . In this case, the symmetric correlationfunction may be equivalently written asTr (cid:2) ρ {O ( t ) , O (0) } (cid:3) (cid:2) ρ O ( t ) O (cid:3) . (D1)Next, expand the operator to second order in a BCHseries O ( t ) ≈ O + it [ H, O ] − t (cid:2) H, [ H, O ] (cid:3) + . . . (D2)and substitute back in. The first order term is zero viatrace identities. What remains isTr (cid:2) ρ {O ( t ) , O (0) } (cid:3) ≈ (cid:104)O (cid:105) − t (cid:2) ρ (cid:2) H, [ H, O ] (cid:3) O (cid:3) + . . . . (D3)Next, inspecting the second term and using the cyclic-ity of the trace and definition of O , findTr (cid:2) ρ (cid:2) H, [ H, O ] (cid:3) O (cid:3) (D4)= Tr (cid:2) H O ρ + H O ρ O − H O H O ρ (cid:3) = 2 ρ n O n (cid:18) O n (cid:104) E n | H | E n (cid:105) − (cid:104) E n | H | E m (cid:105) O m (cid:104) E m | H | E n (cid:105) (cid:19) . (D5)The D P states | E m (cid:105) are constructed such that they arediagonal in H within rotated subspace P , by the TSAprocedure. This means that that (cid:104) E m | H | E n (cid:105) = E n δ mn .However, the operator H will generally be different, as it will include states outside of the subspace (see AppendixB for details). Generally, H may be computed efficientlyanalytically without needing to use a full Hilbert space.Thus this simplifies toTr (cid:2) ρ (cid:2) H, [ H, O ] (cid:3) O (cid:3) = 2 ρ n O n (cid:18) (cid:104) E n | H | E n (cid:105) − (cid:0) (cid:104) E n | H | E n (cid:105) (cid:1) (cid:19) . (D6)The term in parenthesis is the energy variance of eigen-state | E n (cid:105) as defined in Eq. (10), and so the decaytimescale is related to the energy variance as1 τ ≡ (cid:80) n ρ n O n ∆ n (cid:80) n ρ n O n . (D7) Appendix E: AGP Parameters
For completeness, the coefficients for the variationaladiabatic gauge potential are shown in Fig. 8. The valuesare independent of system size and are here computed for20 sites. Some care needs to be taken with the normal-ization of the AGP, as it is derived from the differentialon some path. Note that some terms diverge close to h x = 0, but are cut off numerically. h x A G P C o e ff i c i e n t ZY+YZYXZ+ZXYYX+XYZYZ + YZYZ - Y
FIG. 8. Some of the larger terms in the adiabatic gauge po-tential along the straight-line path from ( h x , h z ) = (0 , ⇒ (0 . , . /h , but is cutoffnumerically when computing the AGP. Notation YX means atranslationally invariant sum of Pauli operators on adjacentsites, eg (cid:80) i σ iy σ i +1 x . ∗ Corresponding author: [email protected] M. B. Hastings and X.-G. Wen, Phys. Rev. B , 045141(2005). B. Nachtergaele and R. Sims, Commun. Math. Phys. ,119 (2006). M. B. Hastings, J. Stat. Mech. Theory Exp. , P08024(2007). X. Chen, Z.-C. Gu, and X.-G. Wen, Phys. Rev. B ,155138 (2010). R. Shankar, Rev. Mod. Phys. , 129 (1994). T. Langen, T. Gasenzer, and J. Schmiedmayer, J. Stat.Mech. Theory Exp. , 064009 (2016). B. Bertini, F. H. L. Essler, S. Groha, and N. J. Robinson,Phys. Rev. B , 245117 (2016). Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol,S. Gopalakrishnan, and B. L. Lev, Phys. Rev. X , 021030(2018). M. Serbyn, Z. Papic, and D. A. Abanin, Phys. Rev. Lett. , 127201 (2013). J. Z. Imbrie, Phys. Rev. Lett. , 027201 (2016). D. A. Abanin, W. De Roeck, W. W. Ho, and F. Huveneers,Phys. Rev. B , 014112 (2017). J. Wurtz, P. W. Claeys, and A. Polkovnikov, Phys. Rev.B , 014302 (2020). D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019). H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Om-ran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres,M. Greiner, V. Vuletic, and M. D. Lukin, Nature ,579 (2017). C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn,and Z. Papic, Nat. Phys. , 745 (2018). M. Srednicki, Phys. Rev. E , 888 (1994). G. P. Brandino, J.-S. Caux, and R. M. Konik, Phys. Rev.X , 041043 (2015). M. Kolodrubetz, D. Sels, P. Mehta, and A. Polkovnikov,Phys. Rep. , 1 (2017). P. W. Claeys, M. Pandey, D. Sels, and A. Polkovnikov,Phys. Rev. Lett. , 090602 (2019). D. Sels and A. Polkovnikov, Proc. Natl. Acad. Sci. U.S.A. , E3909 (2017). E. H. Lieb and D. W. Robinson, Commun. Math. Phys. , 251 (1972). A. J. A. James, R. M. Konik, P. Lecheminant, N. J. Robin-son, and A. M. Tsvelik, Rep. Prog. Phys. , 046002(2017). F. Verstraete, V. Murg, and J. Cirac, Adv. Phys. , 143(2008). P. Weinberg and M. Bukov, SciPost Phys. , 003 (2017). S. Sachdev,
Quantum Phase Transitions , 2nd ed. (Cam-bridge University Press, 2011). P. Calabrese, F. H. L. Essler, and M. Fagotti, J. Stat.Mech. Theory Exp. , P07016 (2012). P. Calabrese, F. H. L. Essler, and M. Fagotti, J. Stat.Mech. Theory Exp. , P07022 (2012). A. B. Zamolodchikov, Int J Mod Phys A , 4235 (1989). J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, andM. Greiner, Nature , 307 (2011). C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn,and Z. Papic, Phys. Rev. B , 155134 (2018). H. Kim, T. N. Ikeda, and D. A. Huse, Phys. Rev. E ,052105 (2014). L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv.Phys. , 239 (2016). This can be seen by the definition of the regularized AGPin Eq. (3). This work exclusively uses an even number of sites andperiodic boundary conditions to avoid any AFM groundstate degeneracy; particles always come in pairs. M. Kormos, M. Collura, G. Tak´acs, and P. Calabrese, Nat.Phys. , 246 (2017). A. A. Ovchinnikov, D. V. Dmitriev, V. Y. Krivnov, andV. O. Cheranovskii, Phys. Rev. B , 214406 (2003). O. F. d. A. Bonfim, B. Boechat, and J. Florencio, Phys.Rev. E , 012122 (2019). C.-J. Lin and O. I. Motrunich, Phys. Rev. Lett. ,173401 (2019). S. Choi, C. J. Turner, H. Pichler, W. W. Ho, A. A. Michai-lidis, Z. Papic, M. Serbyn, M. D. Lukin, and D. A. Abanin,Phys. Rev. Lett. , 220603 (2019). M. Karbach and G. Muller, (1998), arXiv:cond-mat/9809162. M. Serbyn, Z. Papic, and D. A. Abanin, Phys. Rev. B ,174302 (2014). S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zangh`ı,Phys. Rev. Lett. , 050403 (2006). A. J. A. James, R. M. Konik, and N. J. Robinson, Phys.Rev. Lett. , 130603 (2019). N. J. Robinson, A. J. A. James, and R. M. Konik, Phys.Rev. B99