{'fraction_non_alphanumeric': 0.03303040894792031, 'fraction_numerical': 0.008713237129874669, 'mean_word_length': 4.743475767135073, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Emergent conservation laws and nonthermal states in the mixed-field Ising model

This paper presents a method of computing approximate conservation laws and eigenstates of integrability-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 close by in parameter space, dressed by some unitary rotation. However, most adiabatic continuation analyses only use this unitary implicitly. In this work, approximate adiabatic gauge potentials are used to construct a state dressing using variational methods, to compute eigenstates via a rotated truncated spectrum approximation. These methods allow construction of both low and high-energy approximate nonthermal eigenstates, as well as quasi-local almost-conserved operators, in models where integrability may be non-perturbatively broken. These concepts will be demonstrated in the mixed-field Ising model.

Adiabatic continuation is a well-used concept [1]. Given some Hamiltonian of interest, low energy eigenstates of that system may be computed by finding a "simple" Hamiltonian "nearby" in parameter space, which can be easily diagonalized. Then, those simple eigenstates can be mapped to interacting ones using some dressing by a unitary U . This procedure can be used to show that ground states can be mapped to other ground states, as long as there is some path in parameter space which does not go through some gap-closing critical point [2][3][4]. For example, low energies of interacting fermions may be described by a Fermi gas or Fermi liquid with dressed quasiparticle excitations [5]. However, most of these proofs are non-constructive in the sense that the unitary U is generally never actually computed, or only done so perturbatively.
This work details explicitly connecting non-interacting trivial states to nontrivial interacting ones, by variationally computing an approximate, local adiabatic gauge potential (AGP) to construct the unitary transformation. The latter modifies non-interacting particle states to quasiparticle states which are "dressed" within some local span of sites. Importantly, the dressing is nonperturbative and not limited to low-energy states. This construction leads to long-lasting quasiparticles and stable nonthermal states, even at finite energy densities.
As a consequence of computing unitary rotations for approximate eigenstates, this procedure also allows one to compute local almost-conserved operators from the approximate eigenstates and dressed non-interacting symmetries. In the presence of integrability breaking terms in some system, undressed symmetries and conserved operators are generally no longer conserved [6][7][8]: instead, these operators may be "dressed" locally by the unitary U to restore approximate conservation of some quasilocal and long lived operator, even though the full system may no longer be integrable [9][10][11]. Similarly, for particular initial wavefunctions with large overlap with good approximate eigenstates, one may compute long time effective quench dynamics within a small subspace of states, in spirit of Schrieffer-Wolff transformations [12].
The presence of such unitary dressings, quasi-local conservation laws, and good approximate eigenstates in in-teracting models may suggest that not all integrabilitybroken models should be treated equally. Certain models may be "close to integrable", in the sense that there exists a good local dressing of particular eigenstates, and strong ETH may be violated. Two particular cases of such ETH violation are many body localization [13] and quantum scars in the PXP model [14; 15]. Conversely, other models may be far from any simple system, in the sense that there is no path in parameter space that admits a good local adiabatic dressing, and the model is quantum chaotic [16]. In fact, this unitary rotation may potentially be seen as the analog of the canonical transformation of KAM theory which restores integrability in classical models [17]. While this paper focuses on a quantum model, the whole methodology, including variationally computed canonical transformations generated by the AGP [18], is fully applicable to classical non-integrable systems.
The rest of this paper is structured as follows. First will describe the methods of computing approximate eigenstates with the variational adiabatic gauge potential and the rotated truncated spectrum approach (rTSA). The method is a different perspective of Schrieffer-Wolff block 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, and demonstrate the performance of computing approximate eigenstates. Finally will be an example of an approximately conserved local operator for the mixed-field Ising chain, the total quasiparticle number.

I. COMPUTING LOCAL APPROXIMATE EIGENSTATES
In general, directly computing eigenstates and quasiparticle excitations is hard, as computational difficulty scales exponentially in system size, and normally there are no well-defined quantum numbers. Instead, let us continue in spirit of adiabatic continuation [1]. Suppose some parameterized Hamiltonian H(µ), with µ = 1 being the particular system of interest, and µ = 0 being exactly solvable, in the sense that the eigenstates are easily computable via symmetry, integrability, or other means. µ describes some choice of path in a (multi)-parameter space of Hamiltonians between 0 and 1.
At all points along this path, there are parameterized eigenstates {|E n (µ) } with eigenenergies {E n (µ)}. How these eigenstates change as a function of parameter µ is given by the adiabatic gauge potential (AGP) A(µ) [18], computed from the instantaneous Hamiltonian H(µ) (1) In principle, given the (parameter dependent) AGP A(µ), one could use Eq. (1) to evolve the simple eigenstates |E n (0) into exact interacting eigenstates |E n (1) , or equivalently use a unitary rotation such that |E n (1) = U |E n (0) , with T indicating path ordering. In practice, this recipe runs into the same infeasibility as directly computing eigenstates: computing the exact AGP is generally as difficult as computing the exact eigensystem. It is nonlocal, exponentially large, and highly parameter dependent [18] for generic interacting systems, making computing the unitary a similarly difficult task.
To circumvent the problems of computing the exact AGP, one can instead use approximate adiabatic gauge potentials. One might hope that simple eigenstates "dressed" by a local approximation of the AGP will closely resemble eigenstates of the full system. As the complexity of the approximation grows, so too should the approximate AGP approach the exact one, and the approximate eigenstates become exact, at the expense of them being highly entangled and nonlocal.
From an adiabatic continuation standpoint, the connection with this approximate gauge potential is clear; in fact M. Hastings in [1] has is a particular implementation of an approximate AGP for gapped ground states. A simplified version of Eq. (17) where ∂ µ H(t) is the operator ∂ µ H in the Heisenberg representation, erfc is the complementary error function for which erfc(0) = 1 and erfc(∞) = 0, and SGN(t) is ±1 depending on sign of t (see Appendix A for more details). This expression is nothing but an approximation of the gauge potential with a particular choice of regularizer f (t) [19]. For the regularization time τ q → ∞, this approximate AGP becomes exact. For a finite regularization time, the AGP is approximate but local within some span of sites. Instead of computing an approximate gauge potential via some choice of regulator, which is computationally difficult, one can instead compute it variationally [20]. Here, some ansatz for the variational gauge potential is chosen based on the system at hand, then variationally optimized to best approximate the exact AGP. In this case the variational AGP is chosen to be the sum of some set of (local) operators {B i } for variational parameters The particular choice of operators {B i } depends on the problem at hand; if the set is expanded to include all operators in the space modulo symmetries (a potentially exponential number), the exact AGP can be reconstructed from the set, and the variational version becomes exact. Thus, one should expect improvement as the size of the ansatz is expanded. It has been found [19] that a local variational AGP can accurately reproduce the exact one, at least for matrix elements which are far separated in energy or protected by approximate symmetries. This is shown by using a Taylor series expansion of Eq. (3). In fact, this argument is mirrored by an equivalent justification from an adiabatic continuation perspective for gapped ground states: As long as there is a gap, the ground state will only entangle within some light cone [2; 3; 21] and thus a local approximation for the rotation generator is always possible.
The optimization of the AGP is computed by finding the minimum of [18; 20] where parameter dependence of α and H on µ is implicit, and ||Q|| = Tr Q 2 /D is the Hilbert-Schmidt norm. Eq. (6) has the property that S = 0 for the exact AGP, and is solvable even for very large system sizes. This is because the function S(α) can be computed using trace identities, and is quadratic in α, making optimization simple. This minimization leads to a quasi-local operator A(µ) which approximates the exact AGP.
One can then use this approximate AGP to compute approximate eigenstates, in the same manner of the exact AGP computing exact eigenstates. First, choose some set of D p states {|q } of the exactly solvable Hamiltonian H(0). This choice depends on the system at hand. One choice could be, for example, all eigenstates below some energy cutoff. Another choice would be all states within some particular symmetry sector such as fixed particle number, which form a subset of eigenstates not necessarily sorted by energy. The set of states forms some projective subspace P with projector P = q |q q|.
Approximate eigenstates of the interacting model H(1) can be computed by "dressing" each state via Schrödinger evolution of Eq. (1) to implement the unitary of Eq. (2). This gives some set of "dressed" states {|q(1) } ≡ {U |q }.
As a comment on implementation, care must be taken in the direction of evolution: the AGP is parameter dependent so generally A(1) = A(0). For perturbative couplings this parameter dependence is very weak so that the AGP approximately commutes with itself for all µ and the directionality doesn't matter; however strong coupling may lead to nonsensical answers. One may start by accidentally acting on a non-interacting wavefunction with the AGP from the interacting point, which may be much different than the correct non-interacting AGP.
An effective Hamiltonian within that subspace can be computed via matrix elements in the dressed subspace: Note that this is equivalent to computing the effective Schrieffer-Wolff Hamiltonian [12], where one rotates the operatorH = U † H(1)U instead of the states. If the AGP is exact, this effective Hamiltonian will be exactly diagonal. However, if the AGP is not exact or the initial subspace was a degenerate symmetry sector, the effective Hamiltonian will not be diagonal. One can then compute the eigensystem of the D p ×D p matrix via standard linear algebra techniques to find eigenenergies E i and eigenvectors V i such that The approximate eigenvectors of the system are then given by with eigenvalues E i . This final step is functionally equivalent to the truncated spectrum Approach (TSA) [22], except instead of using a non-interacting subspace P , the subspace is first rotated by the approximate AGP to obtain some improved subspaceP . This basis better resembles eigenvectors of the interacting system, and may be exponentially orthogonal from the original basis due to the finite rotation. This corresponds to a rotated truncated spectrum approach (rTSA) The abbreviated method is as follows: 1. Define some Hamiltonian H(µ), with H(0) being exactly solvable and H(1) being a system of interest, with some path in parameter space linking the two.
3. Define some set of eigenstates of H(0), either within some energy window or within some symmetry sector(s) such as particle number.
4. Evolve the set of states via the Schrödinger equation from µ = 0 to µ = 1 with the variational AGP.
5. Compute the effective Hamiltonian and its eigensystem to find approximate eigenstates and eigenvalues.
Because the variational AGP is local by construction, various tricks can be employed to go to large or even thermodynamic system sizes. As commented above, computing the variational AGP is not a problem for a large number of sites, as complexity scales linearly with system size. It is reasonable to compute an AGP for hundreds of sites for all operators spanning less than 5-6 sites on a modern desktop. This locality can also be used for step 4 when evolving the basis states: due to the finite evolution "time", states are only entangled within some finite region (using intuition of Lieb-Robinson bounds [3]). This suppresses the finite size effects of evolving some small (typically 15-20) number of sites exactly, and enables tensor methods such as matrix product states (MPS) and variational evolution [23]. This work employs the former method of exact evolution on small systems [24].
Choosing a larger subspace should also be expected to improve the computation of the approximate eigensystem. In the limit where the projective subspace is the full Hilbert space or within one of the symmetries of the full Hamiltonian, the effective Hamiltonian is the exact one and likewise the eigenstates are exact, independent of the rotation. Choosing a subspace within some larger energy window should also be expected to improve the variational dressing: the AGP fails to suppress excitations close together in energy, but those can then be recaptured within the subspace P via off-diagonal elements of the effective Hamiltonian.
Because both the variational AGP ansatz and the projective subspace can be systematically expanded, this method gives a controllable approximation to compute eigenstates: as the complexity increases, the eigenstates will asymptotically approach the exact ones.
The eigenstates computed in this manner are approximate, in that they are not exact eigenstates of the Hamiltonian H(1). The simplest indicator of the closeness to an exact eigenstate is the energy variance of the state Exact eigenstates have zero energy variance, and so ap-proximate eigenstates should have minimal energy variance ∆ 2 n ≈ 0. The average energy variance of these eigenstates within the subblock corresponds to the average block-off-diagonal matrix elements in the Hamiltonian and thus indicates the performance of the block diagonalization procedure (see Appendix B).

II. MODEL
As a concrete example, suppose the following system, the mixed field Ising model For h z = 0 the model is integrable via a Jordan-Wigner transformation to free fermions [25][26][27] with a critical point 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 classical Ising model. For h 2 x +h 2 z → ∞ the model is an exactly solvable 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 Hamiltonian is the PXP model [14; 30]. Elsewhere, the model has no apparent conservation laws or symmetries beyond geometric ones and is generally quantum chaotic [31; 32]. However, this does not prevent approximate conservation laws 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 map to fermion bilinear operators, plus all operators local within a span of n sites Additional symmetries and properties reduce the size of the ansatz: the AGP has all of the symmetries of the full Hamiltonian [33]. By gauge choice the AGP can be completely imaginary for real Hamiltonians [18], constraining {B} to only include terms with an odd number of σ y . Because the Hamiltonian is translation and reflection invariant, the ansatz can be chosen to be as well. The inclusion of Jordan-Wigner strings is motivated by this ansatz being exact for the transverse Ising model [18], due to its extra symmetries and mapping to free fermions.
While a Hamiltonian of interest is given by particular choice of parameters h x , h z , there is relative freedom for choice of the path h x (µ), h z (µ) in the 2d parameter space H(µ), and especially choice of simple Hamiltonian H(0). This is because there are many "simple" points in the (h x , h z ) parameter space which might be considered "close" to the Hamiltonian of interest. The (h x , 0) line is the transverse Ising model; the (0, h z ) line is the classical Ising model; and the (h x , h z ) → ∞ line are independent spins with onsite fields.
What starting points, and which path in parameter space, is optimal for computing approximate eigenstates, given ansatz {B}, Hamiltonian H(1), and subspace P? This is a question of a path-dependent Schrieffer Wolff transformation, as the performance of computing approximate eigenstates, or equivalently block diagonalization, may depend on these choices. This work chooses from a limited set of parameterized Hamiltonians with particular starting and ending points. Two additional parameterizations are discussed in Appendix C.
The first and second parameterizations start from the σ z σ z point, whose eigenstates are Z polarized spins. Depending on the sign, the ground state could be an antiferromagnetic (AFM) Neél (13) or a polarized ferromagnetic (FM) state (14). Low energy particle excitations are 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 thus the AGP is exact A(µ) ∼ i σ y , and is an example of the Landau-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 in Z. Low energy particle excitations are spin flips (see Fig.  1) to the opposite direction.
In all cases, H * (0) is degenerate, with a natural choice of projective subspace being fixed particle number on top of the ground state. Thus, P 1 is 0 and 2 boundary walls on top of an AFM ground state; P 2 is 0 and 2 boundary walls on top of a FM state; and P 3 is the 0, 1 and 2 particle spin flips on top of a polarized state.
FIG. 1. Example basis states of two particles separated by 6 sites. Top and middle are states with two boundary walls (red dashes), which are low energy eigenstates of H1(0) and H2(0), with excitation energy 4J. Bottom is a state with two spin flip particles, which is a low-energy state of H3(0) with excitation energy 4hz.
Because the system is translation invariant, the zeromomentum sector is chosen as a numerical simplification. Under these constraints, each subspace P has N +1 states each out of total Hilbert space dimension ≈ 2 N /N . These basis states are each dressed by the variational AGP to create dressed boundary wall states: the hard boundary is softened by the dressing procedure to better describe the interacting quasiparticle excitations.

III. APPROXIMATE EIGENSTATES AND SPECTRUM
As an explicit example, let us choose the parameters h x = 0.4 = h z , and coupling J = ±1. These parameters are non-perturbative, in the sense 0.4 is O(1) away from any simple point. For J = −1, the ground state is ferromagnetic, and the h z term acts as a constant attractive force between two boundary wall particles. This leads to "meson" bound states of the two boundary walls [35]. For J = +1, the h z term does not change the AFM ground state energy [36; 37]; however it will affect energies of spin flips on up/down Neél sublattices, which, from a band theory context, leads to two free particle species.
One can then go through the process of computing approximate eigenstates, as outlined in section I, for these particular choices of subspaces. Here, Hamiltonians (13), (14) are chosen starting from the FM and AFM subspaces, with an ansatz of all operators local to 3 sites plus Jordan Wigner strings, with 18 total sites in the 0 momentum sector. The form of the AGP is shown in Appendix E.
Results for these parameters are shown in Fig. 2. Top plots the effective unrotated and rotated Hamiltonian, or equivalently the Hamiltonian in the projective and rotated projective subspaces, for AFM (left) and FM (right) excitations. It can be clearly seen that the rotated effective Hamiltonian becomes slightly more nonlocal: a dressed boundary wall of width 3 may hop to become width 5, for example. These effects are especially pronounced when the two boundary walls are close together, which is an indicator of a 2-particle interaction. When the two particles are far apart the Hamiltonian becomes independent 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 eigenspectrum (black). The error bars are the energy variance of the approximate eigenstates, as computed from Eq. (10). Note that the exact eigenvectors are matched with approximate ones by choosing those which have maximum fidelity | E n |E exact m | 2 ; normally this value is > 0.9. Importantly, the eigenstates are not necessarily all the lowest energy states. For example, two of the lowestenergy FM boundary walls (each with excitation energy 4.4J) has a higher energy then a single FM boundary wall of width 6 (with excitation energy 8.7J). This means it is energetically possible for the width-6 boundary wall to decay into two width-1 boundary walls, for example. However, the small energy variance of these dressed states indicates that such a process is almost completely suppressed.
Because the dressing is local, it is possible to take a continuum or large system size limit. Numerically, this is done by duplicating the dressed Hamiltonians of Fig.  2 over thousands of sites: The 19 × 19 matrix is extended to a N × N matrix, where the middle elements are the duplicated middle elements of the smaller matrix. Then, the eigensystem of that Hamiltonian is computed. Results for the continuum dispersion relation of excitations on top of the AFM ground state are shown in Fig. 3. There are two particle species which have mass 1.11 and 2.56. Like the meson case, the states are not necessarily the lowest-energy states: for example, 2 heavy particles can have equal energy to 4 light particles, and may potentially decay as such. Note that two light particles could not decay into 1 heavy particle, as that is disallowed by the particles being domain walls.
These approximate eigenstates can be compared with the general bulk eigenstates, as is shown in Fig. 4 Top. Here, all eigenstates in the 18-site FM model, 0momentum sector (14,602 total) are computed exactly, and their half-cut entanglement entropy is found. Thermal states are extensively entangled states at finite energy density, while nonthermal states are weakly entangled and generally have zero or low energy density [30]. Red circles are the states which have maximal overlap with the approximate eigenstates: low energies are FM states while high energies are AFM states.
One particular approximate eigenstate merits more study: the dressed all-up eigenstate, indicated by the green circle and arrow in Fig. 4 Top. This is a ground state of σ z σ z , but the most excited state of σ z ; it has finite energy density given roughly by 2h z . But, in particular, it is an explicit example of a highly non-thermal  with an exact eigenstate. Note that in the thermodynamic limit the dressed-all-up state is exponentially orthogonal to the original all-up state due to the finite local rotation.

A. Quasiparticle Lifetimes
The energy variance of these approximate eigenstates takes special meaning when they can be interpreted as dressed particles. In this case, the energy variance gives a lower bound on the quasiparticle lifetime. For an approximate particle eigenstate |E n , the time-dependent state overlap under second order perturbation theory is where τ −2 = ∆ 2 = H 2 − H 2 is the energy variance of the state. In other words, the characteristic time for an (eigen)state of some particles to decay into some other particle state is given by the energy variance. This timescale is very crude as it assumes all other states have the same energy: a more refined timescale can be computed using the Fermi golden rule [25] for the dressed states, but is not generally possible without a priori 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 dressed single flipped down spin on a FM ground state, which corresponds to the lowest energy meson excitation, has a characteristic lifetime of τ = 110, far longer than any local timescale. Excitations on top of an AFM ground state have lifetimes in excess of τ > 60. The dressed all-up state has a lifetime of τ = 53. These lifetimes are longer when adding more parameters to the variational AGP. Explicit time dynamics of Eq. (17) for this dressed all-up state is shown in Fig. 4 Bottom. Clearly, it is much closer to an eigenstate than expected, as it is close to 1 at all times. This further indicates the genuineness of this nonthermal eigenstate, especially when compared to the undressed version of the same. Note that the undressed up state is exponentially orthogonal in system size from the dressed up state, due to the finite rotation.
One application of such a dressed all-up state is for information protection in quantum systems. For a classical Hamiltonian, the all-down [ground] state may be labeled as a logical 0, while the all-up state is labeled as a logical 1. An X-field will generally change these two states, destroying the encoded bit. If this bit is instead encoded in the dressed nonthermal states (the all down ground state, and the all up nonthermal state), they are much more stable, encoding the information for a much longer time by suppressing transitions.

B. Quasiparticle parameter dependence
A general measure of the quasiparticle lifetime within a particular subspace is given by their normalized average inverse lifetime, or equivalently average energy variance Γ equivalently is the block-off-diagonal weight of the rotated subspaceP (See Appendix B). A small value indicates a good block-diagonalization procedure with welldefined quasiparticles within the subspace. A large value indicates 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 computing their normalized average energy variance Γ(h x , h z ). Results are shown in Fig. 5, for a 3-site ansatz and 14 sites. States are dressed from one of two directions. One is H 1 (µ), dressing 2-particle AFM boundary wall states out from the σ z σ z only point, indicated by the radial arrows in the bottom left. The other is H 3 (µ), dressing 1 and 2-particle spin-flip states from the h x σ x + h z σ z only point(s), indicated by the arrows pointing radially inwards.
In the region where h x , h z is small, the error from dressing boundary walls is enormously low. With no dressing, the error grows linearly in |h|, while with dressing, the error grows sub-linearly, which indicates that the dressing is exact asymptotically. In fact, this dressing accumulates very small errors even for non-perturbative values of |h|, as shown in Fig. 5 Bottom, which is dressing along the h x = h z line. This indicates that in the white areas, there is a good effective quasiparticle description of the low energies of this otherwise quantum chaotic model, described by dressed boundary wall particles.
Although it is not generally so, the error accumulates monotonically with increasing |h|. This means that at some critical value, the dressing going outwards from  White and yellow indicate areas with a good quasiparticle description. Bottom: Average energy variance in the direction π/4 from vertical. Middle is improvement from the undressed subspace. Below h ≈ 0.6, the error is vanishingly small. Middle: Energy variance improvement compared to the Null ansatz Γ/Γ0 along the hx = hz line. As the ansatz size increases, so too does the improvement, as expected.
the σ z σ z point will have a larger error then the dressing going inwards from the (h x , h z ) → ∞ point. At this boundary, the best description of quasiparticles changes from dressed pairs of boundary walls, to dressed spin flips. This does not mean that there is no effective description of certain states in terms of quasiparticles: there could be some other subspace (say, of doubly-flipped spins) and other path through parameter space which gives a better quasiparticle picture in that the energy variance is smaller.
This cross-over point may be an indicator of an interacting phase transition. Around h z = 0, h x = 1, which is the transverse Ising phase transition, it has been found that 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 interacting critical point may occur, as a local AGP fails to reproduce the long-range entanglement of a critical ground state [3]. With increasing ansatz size, this point decreases in total error, and shifts in critical parameter (see Fig. 5 Bottom) which may eventually converge to some particular value, indicating the interacting critical point. This idea is backed up by the convergence in the non-interacting limit: For h z small, the crossover is around the h x ≈ 1 transverse Ising critical point. Similarly, for h x small, the crossover is around h z ≈ 2, which is the first-order phase transition in h z to change the ground state from AFM to polarized [29].

IV. LOCAL ALMOST-CONSERVED OPERATORS
Given some set of approximate local eigenvectors {|E n } generated by this adiabatic dressing scheme, it is a relatively simple procedure to construct approximately conserved local operators. An operator is conserved if it commutes with the Hamiltonian, or equivalently if it is constructed from eigenstates of the Hamiltonian. Given approximate eigenstates, one should then be able to compute approximately conserved quantities. Conserved operators are defined as where O n are the eigenstates of the operator, and{|E n } are exact eigenstates. An operator of this form has the property that the symmetric time correlation function is conserved (for all initial states) There are two ways to construct approximate versions of these operators. The first way is to explicitly use Eq. (19) using only the subspace of dressed eigenstates which were directly computed. In this case, the sum is of size D P , the subspace size, as opposed to D, the total Hilbert space size. Due to the global projective structure of particle excitations on top of a ground state, the resulting operator is not necessary local. However, this may be implemented with some local operator plus post-selection of states.
In the case of such an operator directly constructed from approximate eigenstates, the symmetric correlation function is not conserved in time. Under perturbation theory, the characteristic timescale is given by a weighted-average energy variance (see Appendix D for derivation) Here, ρ n = E n |ρ|E n is the density matrix for the expectation value, assumed to be diagonal, and ∆ n is the energy variance of the nth approximate eigenstate. For good eigenstates with low energy variance, the decay timescale can be very long.
The second way to construct approximately conserved operators is to dress conservation laws of the simple system H(0) with the unitary. Conservation laws, such as particle number and particle current, are constructed from simple eigenstates in the form of Eq. (19), with particular choice of O n , and are generally local [26; 40]. Then, the dressed operator is constructed from states better resembling eigenstates Importantly, the eigenstates of operator O are not necessarily the same as the constructed approximate eigenstates, as it is missing the re-diagonalization step of Eq. (9). The resulting operator is quasi-local, and approximately conserved [9; 41]. As the ansatz span is increased, the AGP approaches the exact one, resulting in better approximate eigenstates and a better-conserved operator, at the expense of it becoming more and more nonlocal.
One such conserved operator for the mixed-field Ising model is the dressed total particle number N = z U † which (up to a constant) counts the number of boundary walls in the system. For H 1 and H 2 , this is also the initial Hamiltonian; thus one would expect that the dressed operator should also approximate dressed versions of particle number eigenstates.
For h x = h z = 0.4 and J = +1, the dressed particle number operator becomes where ellipsis represent the more and more nonlocal terms of the operator. As can be seen, this operator is approximately local, with dominant terms coming from 1, 2, and 3-spin terms. One can then compute the symmetric correlation function in the initial undressed subspace of two particles to see its conservation Results are shown in Fig. 6 for AFM and FM states, as well as infinite temperature typical states [42]. For comparison, the undressed operator is also shown in red. The conserved operator for AFM states is almost stationary in time, while the undressed version is not. The infinite temperature timescale can be computed analytically as Even for an infinite temperature state, this quasilocal dressed operator gets a factor of 40 improvement in the characteristic decay timescale. This indicates that dressed quasiparticle excitations may persist in this interacting model even at infinite temperature, and that this model is "closer to integrable" then one might expect.

V. CONCLUSION
The existence of good approximate dressings have some curious implications. Even if a model system is not necessarily integrable or exactly solvable, that does not mean that there are no local long-lived symmetries and conservation laws. Indeed, if such a model is close by to an integrable point, a conservation law of the integrable model can be "dressed" by a unitary generated by the approximate local adiabatic gauge potential to restore the symmetry approximately in a now quasi-local operator. Approximate eigenstates may be computed in a similar manner: simple particle excitations of the integrable point can be dressed by the approximate AGP to construct long-lived quasiparticle excitations of the interacting point. These new dressed states need not be low energy states and in fact may be used to construct finite energy density low-entanglement nonthermal states, as demonstrated in the dressed-all-up state of the mixedfield Ising model. Similar studies have been done to compute low energy phenomenology of the Meson case, most predominantly in recent work by [43; 44] using a truncated spectrum approach (TSA). These numerical diagonalization procedures are functionally equivalent except that here the projective subspace is first rotated by the variational AGP, leading to a subspace closer to the exact eigenstates. While this work uses discrete lattices, generalizations to continuous theories is an interesting future direction.
The restoration of approximate symmetries and construction of quasiparticle excitations in interacting models puts a new perspective on integrability breaking. Instead of reevaluating a Hamiltonian for every new point in parameter space, one can instead compute properties and approximate symmetries based on nearby Hamiltonians with a potentially simpler structure. This "closeness" is defined in the sense of being able to compute a good approximate AGP along some path between the simple Hamiltonian and interacting one, not in the sense of perturbative parameter changes. Certain perturbations away from integrability may rapidly destroy any local conservation laws, if there exists no good local approximate AGP. Other perturbations, while still breaking integrability, may still admit quasi-local conservation laws, nonthermal states, and quasiparticles, if there does exist a good local approximate AGP.
These unitary rotations restoring approximate integrability are similar in spirit to canonical transformations in KAM theory [17]: integrability may be approximately restored for particular subsets of initial conditions of particular integrability-broken systems via the unitary rotation (eg canonical transformation) of conserved quantities. Whether this approach to stability of quantum integrable systems can be made more concrete remains to be seen but these variational local dressings may be a step towards a general theory in that direction.
Up to a trivial factor of √ 2 which can be absorbed by the regularization time τ q , this is the expression in Eq. (3).

Appendix B: Off-diagonal matrix elements and energy variance
The average energy variance of approximate eigenstates analogously gives the average block-off-diagonal elements in the Hamiltonian, as claimed in section III B. This Appendix serves to elaborate on this point.
When computing an effective Hamiltonian within a rotated subspace, the procedure is an analogous one to a Schrieffer-Wolff transformation: a unitary rotation block diagonalizes some Hamiltonian into a subspace P and complement Q. A measure of the quality of this diagonalization is the average strength of the off-diagonal elements: zero strength means exact block diagonalization, while nonzero strength means approximate diagonalization. The average energy variance is defined as Step 2 inserts the identity, for complete set of states |q ∈ Q, and complete set of states |E p ∈ P , while step tional paths are compared to the diagonal path. The second goes up in h z , for which the gauge potential is zero, and then right in h x . The third goes right in h x as the transverse Ising model for which the gauge potential is exact, and then up in h z (See inset in Fig. 7). This directionality is shown in Fig. 7 for 14 sites, and a size-3 ansatz. The diagonal path ends with an average energy variance of ≈ 0.02, while the other paths have variance of ≈ 0.04, a factor of two improvement.
The D P states |E m are constructed such that they are diagonal in H within rotated subspace P , by the TSA procedure. This means that that E m |H|E n = E n δ mn . However, the operator H 2 will generally be different, as it will include states outside of the subspace (see Appendix B for details). Generally, H 2 may be computed efficiently analytically without needing to use a full Hilbert space. Thus this simplifies to (D6) The term in parenthesis is the energy variance of eigenstate |E n as defined in Eq. (10), and so the decay timescale is related to the energy variance as (D7)

Appendix E: AGP Parameters
For completeness, the coefficients for the variational adiabatic gauge potential are shown in Fig. 8. The values are independent of system size and are here computed for 20 sites. Some care needs to be taken with the normalization of the AGP, as it is derived from the differential on some path. Note that some terms diverge close to h x = 0, but are cut off numerically. Note that one term diverges as 1/h, but is cutoff numerically when computing the AGP. Notation YX means a translationally invariant sum of Pauli operators on adjacent sites, eg i σ i y σ i+1 x .