Weak integrability breaking perturbations of integrable models

A quantum integrable system slightly perturbed away from integrability is typically expected to thermalize on timescales of order $\tau\sim \lambda^{-2}$, where $\lambda$ is the perturbation strength. We here study classes of perturbations that violate this scaling, and exhibit much longer thermalization times $\tau\sim \lambda^{-2\ell}$ where $\ell>1$ is an integer. Systems with these"weak integrability breaking"perturbations have an extensive number of quasi-conserved quantities that commute with the perturbed Hamiltonian up to corrections of order $\lambda^\ell$. We demonstrate a systematic construction to obtain families of such weak perturbations of a generic integrable model for arbitrary $\ell$. We then apply the construction to various models, including the Heisenberg, XXZ, and XYZ chains, the Hubbard model, models of spinless free fermions, and the quantum Ising chain. Our analytical framework explains the previously observed evidence of weak integrability breaking in the Heisenberg and XXZ chains under certain perturbations.


I. INTRODUCTION
Understanding how a many-body system reaches thermal equilibrium is one of the fundamental questions in statistical mechanics.A generic (non-integrable) quantum many-body system that evolves under unitary dynamics is typically expected to thermalize [1][2][3][4][5][6]: When thermalization occurs, expectation values of local observables reach stationary values that depend only on few properties of the initial state and can be predicted with usual statistical mechanics ensembles.An exception is represented by integrable models: These models have a large number of extensive local conserved quantities that retain information about the initial state.As a consequence, integrable models do not thermalize in the usual sense, but, in contrast, their time evolution can be described as a relaxation to a stationary ensemble that includes all conserved quantities, called generalized Gibbs ensemble (GGE) [7][8][9][10][11][12].
In the presence of a perturbation that breaks integrability, usual thermalization is again expected to take place.However, if the perturbation is small, the process may require a long time.On a finite time scale, the dynamics is approximately described by the evolution under the integrable unperturbed Hamiltonian.The system initially relaxes to a stationary state of the unperturbed Hamiltonian (prethermalization), while genuine thermalization only occurs at later times [13][14][15].This later thermalization is typically modelled using Fermi's golden rule, which prescribes a thermalization time τ ∼ λ −2 , where λ is the perturbation strength [16,17].However, for specific Hamiltonians this timescale can be much longer.For example, Abanin et.al. [18] proved that for unperturbed Hamiltonians with special structure producing equally-spaced sectors such as the Hubbard model in the limit of zero hopping, the thermalization time has a lower bound τ ∼ e c/λ .The result can be essentially proven by using local unitary (Schrieffer-Wolff) transformations that, order by order, eliminate the perturbations in the rotated frame, stopping at a "sweet" order where the combinatorial growth of the number of perturbation terms becomes overpowering factor.(The rigorous theory of prethermalization was also proved for unperturbed Hamiltonians with energy sectors determined by more than one frequency [19,20].)With this approach, the conserved quantities of the unperturbed Hamiltonian that label the equally-spaced sectors (such as the doublon number in the Hubbard model) are "dressed" by the perturbation, and are conserved up to times that are exponentially large in 1/λ.There is no analog of this rigorous theory of prethermalization for general unperturbed Hamiltonians, and it is not possible, in general, to find such unitary Schrieffer-Wolff transformations in the thermodynamic limit.Nevertheless, the emergence of approximate conserved quantities was observed in certain other models, that do not belong to this special class.
More specifically, Kurlov et.al. [21] recently showed that the Heisenberg chain perturbed with a next-nearest neighbour SU(2) symmetric interaction has several approximately conserved quantities that commute with the perturbed Hamiltonian up to corrections of order λ 2 .A similar case was observed for the Heisenberg and XXZ chains perturbed with isotropic next-nearest neighbor interaction by Jung et. al. in an earlier Ref. [22], where the authors found an anomalously large heat conductivity that they explained from the presence of an approximate conservation law.In such cases, the presence of approximate conservation laws can lead to longer thermalization times, of the order of τ ∼ λ −4 .These approximately conserved quantities were discovered through the use of numerical methods or heuristic procedures, and are typically limited to a small number of operators.It remained unclear why approximate conservation laws appear in these systems.
In this work we use a systematic analytical approach to find models with approximately conserved quantities, and to compute such quantities.This approach is based on recently studied so-called long-range deformations of integrable spin chains.These deformations were introduced in [23,24] in the context of AdS/CFT.It was later shown that some of them can be seen as generalizations of T T deformations of 1 + 1-dimensional integrable quantum field theories [25][26][27].One can view these deformations as produced by continuous unitary transformations generated by special operators which are not necessarily local but which produce local terms order-by-order in the expansion in the continuous parameter.These special operators (that include so-called boosted operators or bi-local operators) are constructed from the original integrals of motion.
In our work, we utilize the idea that truncations of these deformations at finite order can be viewed as special perturbations that break integrability but more weakly than generic perturbations [28,29]; in what follows, we will often refer to such weak integrability breaking perturbations as "weak perturbations" in short (an alternative name which we will sometimes use is "nearlyintegrable" [30]).Since the deformations are generated by unitary flows, the same transformations apply to all integrals of motion of the original integrable model, and the corresponding truncations then produce approximate conserved quantities [which we will often refer to as quasiconserved quantities or quasi-integrals of motion (quasi-IoMs)] of the perturbed model.We show that the Heisenberg and XXZ chains perturbed by the second-neighbor Heisenberg term can be cast as an example of this approach.This explains the previous findings of few quasiintegrals of motion in these chains and also shows that all integrals of motion of the original integrable chains give rise to quasi-integrals of motion of the perturbed chain, allowing us to derive all of them.
The power of this approach is that it allows construction of large families of weak perturbations to a given integrable model.Besides the Heisenberg and XXZ models, we show additional simple examples for a number of wellknown integrable models including the Hubbard model, free-fermion models, and quantum Ising models.Some of these examples -like density-assisted hopping perturbation for spinless fermions or a particular self-dual deformation of the quantum Ising chain -have in fact appeared in the literature in various contexts without appreciating that they are weak integrability breaking perturbations.
As another application of our systematic approach, we provide an explicit demonstration how to construct weak perturbations beyond leading order: Thus, we show that the Heisenberg chain with second-neighbor and third-neighbor interactions (with specific coupling proportional to the square of the second-neighbor coupling) is a second-order weak integrability breaking perturbation.
The existence of approximate conserved quantities implies that the relaxation to thermal equilibrium upon the introduction of a weak perturbation is slow: We provide a rigorous bound on the thermalization rate that is linear in the effective strength of the remaining generic integrability-breaking perturbation.We also review and verify a non-rigorous "Fermi-Golden-Rule" type estimate which is believed to be quadratic in the remaining per-turbation strength.These types of estimates apply to all weak integrability breaking perturbations considered here.
The paper is organized as follows.In Section II we review the families of long-range deformations of integrable models introduced in Refs.[23,24], and the classes of operators (extensive, boosted, bilocal, or a combination of the three) that generate them.These deformations depend smoothly on a parameter λ.In Section III we consider truncations of the deformations to a finite order in λ, and thus derive families of Hamiltonians with quasiconserved quantities.In Section IV we examine several examples of weak perturbations generated by boosted operators to first order in λ: We apply the construction to the Heisenberg, XYZ, XXZ chains and to the Hubbard model.In Section V we consider additional examples of first-order weak perturbations, that are generated by bilocal operators: We focus on free spinless fermions, the quantum Ising chain, and the Heisenberg chain.In Section VI we demonstrate how to apply the procedure to obtain weak perturbations beyond first order, and we consider for concreteness the case of the Heisenberg chain.In Section VII we discuss how weak perturbations imply parametrically longer thermalization times.Finally, in Section VIII we summarize our conclusions and suggest future outlooks.

II. DEFORMATION OF LOCAL CONSERVED CHARGES
We consider a one-dimensional quantum system whose Hamiltonian H commutes with a set of extensive local operators (charges) {Q α } of the form where j labels the sites of a one-dimensional lattice and the charge density q α,j is an operator with finite range (i.e., acting on a finite number of sites around j).We also assume that [Q α , Q β ] = 0 for every pair of charges Q α , Q β , and the set of conserved charges also includes the Hamiltonian H.This setting applies to integrable models, as well as models with only a finite number of conserved quantities.
We will now follow Refs.[23,24,28] and consider deformations of the charges that depend smoothly on a parameter λ.We define a set of conserved charges {Q α (λ)} (and their charge densities {q α,j (λ)}), that for λ = 0 coincides with the original set of (undeformed) charges.The deformed charges are generated by an operator X(λ): With this definition, {Q α (λ)} form a set of mutually commuting charges for any λ, as can be proven by noting that and the initial condition There are various types of operator X(λ) that lead to quasi-local deformations.This is clearly the case, for example, when the operator X(λ) is a local or quasilocal extensive operator.However, more unconventional classes of operators can also generate quasi-local deformations.We now consider three classes of operators that generate quasilocal transformations for a generic set of conserved charges.

A. Local extensive operators
We can consider arbitrary translationally-invariant operators of the form where s labels different such operators R s that we may want to consider (e.g., different types of Pauli strings up to some range), and e s (λ) can be arbitrary functions of λ.

B. Boosted operators
One example of non-trivial generators is the class of boosted operators: we consider an operator X(λ) of the form where f β (λ) is a real function and B[Q β (λ)] is the boosted operator of the charge Q β (λ), defined as While [28]: because of the commutation relations, we can define the generalized currents from which we get The deformed charges and the generalized currents are expected to be quasi-local in a finite range of λ close to λ = 0.
It is important to note that the charge densities (and hence the current operators) are not uniquely defined.A charge density qα,j = q α,j + o α,j+1 − o α,j corresponds to the same extensive operator Q α .The corresponding boosted operator then equals the original boosted operator shifted by an extensive local operator: Thus, the non-uniqueness of the definition of the boost of charge operators when defining X bo (λ) is equivalent to allowing adding an arbitrary extensive local operator to X(λ).In the following, we will make a specific choice when defining the boosted operators and will also make the above allowance where needed.

C. Bilocal operators
Another example of non-trivial generators of quasilocal deformations is the class of bilocal operators, with X(λ) defined as with Using Eq. ( 7) we get We note that similarly to the boosted operators case, the non-uniqueness in the definition of the charge density of the form qβ,j = q β,j + o β,j+1 − o β,j leads to the shift of X bi by an extensive local operator.On the other hand, shifting qβ,j by a constant corresponds essentially to adding a boosted operator to the bilocal operator; equivalently, using a trivially conserved quantity -the identity operator -as one of the operators in the bilocal construction gives In what follows, we are assuming making specific choices for densities and currents when defining the bilocal operators, while the non-uniqueness is taken care of since we are separately including local extensive and boosted operators as generators.

D. Generic deformation
In general we can define Given the deformed charges {Q α (λ)}, we can construct a family of Hamiltonians that commute with them, of the form The coefficients c α (λ) are chosen such that H(λ = 0) coincides with the original Hamiltonian H, but otherwise they can be arbitrary functions of λ.

III. FINITE ORDER
The construction above allows to define deformations of local Hamiltonians that maintain the presence of a set of conservation laws.However, the deformed charges and Hamiltonians obtained with this procedure are not local, but rather quasi-local: they contain arbitrary longrange contributions, with amplitudes that decrease exponentially with the range.In this Section we show how to construct families of Hamiltonians that have strictly finite-range deformed charges that are quasi-conserved.More precisely, for any integer ℓ > 1, the quasi-conserved charges Q <ℓ α (λ) and Hamiltonian H <ℓ (λ) satisfy To define the quasi-conserved charges, we expand the deformed charges Q α (λ) defined by Eq. ( 2) as power series in the small parameter λ: Similarly, we expand the functions e s (λ), f β (λ), g βγ (λ), and c α (λ), yielding Substituting Eqs. ( 4), (5), and (10) in Eq. ( 2) and equating the two sides order by order in λ we get for k = 0, 1, 2, . . . .The above equation shows how to construct the deformed charges order by order in λ, given the real numbers e βγ that parameterize the deformation.The quasi-conserved charges Q <ℓ α (λ) can then be defined as the truncation of the quasi-local charges to a finite order: Since these differ from Q α (λ) by O(λ ℓ ), they indeed satisfy Eq. (15).Similarly, the Hamiltonian H <ℓ (λ) satisfying Eq. ( 16) can be defined as a truncation of H(λ) in Eq. ( 14): with The existence of a set of quasi-conserved charges has practical consequences for the dynamics generated by the Hamiltonian H <ℓ (λ), which we discuss in Sec.VII.

A. First order
We now focus on the case ℓ = 2 and show how to construct the quasi-IoMs α and the perturbed Hamiltonian H <2 (λ) = H (0) + λH (1) .The leading correction to the charge is particularly simple: and involves only "data" of the unperturbed IoMs Q β .
In many examples we consider, the unperturbed Hamiltonian is commonly used to define the charge 2 δ α,2 with some fixed number c (0) 2 , we have to leading order From Eqs. (24) and (26), we see that we can construct a space of possible "weak-perturbations" to first order in λ as a linear space spanned by the unperturbed charges 2 ] (where R s is a generic extensive operator), the boosted-generated deformations 2 ] = J β,2;tot ≡ J β;tot (i.e., the summed current associated to the charge Q (0) β ), and the bilocal- 2 ] [computed as in Eq. ( 12)].The fact that current operators are "nearly integrable" perturbations was recently noted by Durnin et.al. [30] by studying the non-equilibrium dynamics of charges.Our approach extends this class beyond current operators.Note, however, that we do not know if this exhausts the full space of all weak-integrability-breaking perturbations to first order-this is an interesting open question.

IV. EXAMPLES: DEFORMATIONS USING BOOSTED OPERATORS
A. Heisenberg chain As a first example, we here discuss a particular quasiintegrable deformation of the spin-1/2 Heisenberg chain.The undeformed Hamiltonian has the form where ⃗ σ j = (σ x j , σ y j , σ z j ) is the vector of Pauli operators on site j.The model is integrable, and the densities of the first few (i.e., with smallest range) conserved charges have the form The Hamiltonian is included as 2 .A systematic way of obtaining the higher conserved charges in the Heisenberg chain is by applying the commutator with the boosted operator B[Q see Refs.[28,31] and footnote [32].Appendix A 1 shows densities of two more conserved charges obtained this way.
We now want to construct an example of a simple deformation of the Heisenberg Hamiltonian that is generated by a boosted operator and is quasi-integrable to order ℓ = 2.We therefore need to define the parameters of the deformation f (0) β : note that these are the only parameters that define the deformation of the charges to first order in λ, Eq. ( 24), since we are focusing on deformations generated by boosted operators only.
We consider and for the first two charges we obtain We can now define the deformed Hamiltonian using Eqs.( 22) and (23).Since The Hamiltonian H <2 (λ) = H (0) + λH (1) is a quasiintegrable deformation of the Heisenberg chain to order ℓ = 2 for any choice of coefficients c (1)

in which case
Note that we obtained a range-3 term[33] by cancelling the range-4 part in q (1) 2 by a similar part in q 4 ; this is a general property of the deformation generated by Eq. ( 32) when the unperturbed integrable model contains only nearest-neighbor interactions and its IoMs are obtained using Eq. ( 31), see App.A 1 for details.All quantities α are quasi-IoMs of the deformed model in the sense of Eq. ( 16).
The fact that the Heisenberg chain perturbed by the second-neighbor Heisenberg interactions is an example of weak integrability breaking was first noticed in Ref. [22] in their calculations of the thermal conductivity of the perturbed Heisenberg chain, and they pointed presence of a quasi-conserved quantity with density proportional to q3,j = q where we use our convention for writing the leading term.In App.A 1 we show that this quasi-IoM is directly related to q <2 3 = q 3 by adding a combination of the unperturbed conserved densities multiplied by λ, namely q3,j = q <2 3,j + λ 6q which indeed maintains the quasi-conservation property Eq. ( 16).Later work Ref. [21] used brute-force search for quasiconserved quantities in this model and found the same quasi-IoM, Eq. ( 37), as well as several other longerranged quasi-IoMs, and in App.A 1 we show that the next quasi-IoMs are also reproduced by our approach.We thus suggest that the specific unitary transformation explains all these results; our approach proves that there are in fact as many quasi-IoMs as in the original integrable model and provides a straightforward recipe for obtaining all of them, without the need for brute-force searches.

B. XYZ and XXZ Chain
A generalization of the Heisenberg chain is represented by the XYZ chain, with Hamiltonian Similarly to the case of the Heisenberg chain, the set of conserved charges can be obtained by defining 2 and using Eq. ( 31).The first two charge densities have the form Expressions for further conserved charges in the XYZ chain can be found, e.g., in Ref. [34] (with proper translation from their boost definition to ours).Explicit expressions were also calculated in Ref. [35], and in Ref. [36] for the case of the XXZ chain.We again consider a deformation generated by a boost operator, with f α to first order using Eq.(32); we list the expressions for Q with any linear combination of the unperturbed charges is then a generic quasi-integrable perturbation to order ℓ = 2.
A possible choice of coefficients c α of the linear combination that gives a simple (e.g., without three-body terms) perturbation is c 4 , we get where we have defined ⃗ t 2 ≡ α∈{x,y,z} t 2 α .The reason the specific combination cancels range-4 terms so that only range-3 terms remain is the same as in the Heisenberg case, see App. C.
We can readily write the corresponding quasiconserved quantity q <2 3,j = q 3,j .We have verified that we can eliminate range 5 terms by combining with q (0) 5,j as in the Heisenberg case, Eq. (38).Specifically, we find Note that unlike the Heisenberg case, combining with q (0) 3,j does not fully eliminate the first term involving three consecutive sites, so we do not show such combinations.
A particularly relevant case is the XXZ chain, obtained when t x = t y .To simplify the notation, we set t x = t y = 1, so the unperturbed XXZ Hamiltonian reads In this case, we can obtain a simpler perturbation with the choice of coefficients c (1) , which yields The nearest neighbor term σ z j σ z j+1 simply gives a correction O(λ) to the parameter t z , while the next-nearest neighbour Heisenberg interaction can be regarded as the actual weak integrability-breaking perturbation.
To argue this precisely, we need to be careful that the IoMs of the XXZ chain depend on the anisotropy parameter t z , so we need to properly take into account the above shift of the t z by O(λ) when finding the quasi-IoMs for the pure second-neighbor Heisenberg interaction perturbation To accomplish this, we start with the expected commutation for the found quasi-IoMs for the perturbation H (1) : where we have explicitly indicated that these operators have the t z parameter in them.Next, we write , where we introduced the parameter and in the end expanded to the exhibited order in λ.
Repeating the same for the quasi-IoM parts and plugging into the above commutator, we obtain At this point we can drop the prime on the parameter t ′ z and conclude that indeed H(1) is a weak integrability breaking perturbation for the XXZ chain at any value of the anisotropy parameter and also read off the corresponding quasi-IoMs in terms of the ones obtained for H (1) above.
We note that this isotropic next-nearest-neighbor perturbation of the XXZ chain was considered in [22,37], where slow relaxation was signalled by a heat conductivity of order ∼ λ −4 (much larger than the generically expected scaling ∼ λ −2 ), and explained as the consequence of the existence of a quasi-conserved quantity.
Here we found an explicit form of the quasi-IoM derived from Q (0) 3 .Furthermore, our approach shows that there is, in fact, an extensive number of such quasi-conserved quantities.
The above demonstration of a near-integrability of H(1) [Eq.(46)] starting from the near-integrability of H (1) [Eq.(45), derived from the truncated deformations], in fact applies quite generally for integrable models with a continuosuly varying parameter.Specifically, suppose the unperturbed Hamiltonian has the form with a parameter u appearing as a coefficient of some part of the Hamiltonian, where we assume that both H II have no u dependence in them (they may depend on some other parameters, which however are kept fixed throughout).Suppose that H (0) (u) has IoMs Q (0) α (u), which in general depend on u (and in more complicated ways than H (0) ).Then H (1) = H II can be considered as a weak (nearly-integrable) perturbation of H (0) α (u+λ).Of course, by continuing the Taylor expansion we can write quasi-IoMs for the specific perturbation to arbitrary order ℓ, but our main interest is ℓ = 2 where we can combine the above H (1) with general nearly-integrable perturbations since these form a linear space: If H (1)′ and H (1)′′ are nearlyintegrable perturbations with the quasi-IoM corrections Q (1)′ α and Q (1)′′ α respectively, then a ′ H (1)′ + a ′′ H (1)′′ is also a nearly-integrable perturbation with the quasi-IoM correction a ′ Q α , for arbitrary a ′ and a ′′ .As an immediate application, the unperturbed XYZ model has parameters t x , t y , and t z , and hence we can add independent nearest neighbor terms σ x j σ x j+1 , σ y j σ y j+1 , and σ z j σ z j+1 while preserving ℓ = 2 near-integrability.In this way, starting with Eq. ( 42) we conclude that the second-neighbor Heisenberg interaction is a nearly integrable perturbation for the XYZ chain with arbitrary anisotropies (with correspondingly recalculated quasi-IoMs).

C. Hubbard model
The Hamiltonian of the Hubbard model reads We use notation from Ref. [34] for easy referencing to their expressions for the IoMs.We define 2 = H (0) .In contrast to the Heisenberg, XYZ, and XXZ chains, where the conserved charges can be obtained using B[Q 2 ] as a "ladder" operator [Eq.( 31)], no similar systematic construction can be used to find the conserved charges in the Hubbard model, and brute force methods have been used instead [34,38].The first IoM is: where "−s" denotes the opposite spin to s, i.e., −s =↓, ↑ for s =↑, ↓.This IoM is symmetric (even) under the physical spin SU(2) and pseudo-spin (a.k.a.η-pairing) SU(2) symmetries of the Hubbard chain but is odd under the time reversal and inversion symmetries.In what follows, we define the densities q (0) 2,j and q 3,j such that they are respectively even and odd under the inversion in the bond center between j and j + 1, see Eqs. (D1) and (D2) in App.D.
Since B[Q 2 ] does not generate IoMs, we can consider deformations generated by the boost operator B[Q 2 ], as they will correspond to proper perturbations of the model, instead of integrals of motion.We can therefore define Some examples of such first-order weak integrability breaking perturbations H (1) that can be generated in this way are or These terms preserve the spin SU(2) and pseudo-spin SU(2) symmetries but break the time reversal and inversion symmetry.In particular, we see that the simplest hopping modification of the Hubbard model that preserves both SU(2) symmetries-namely, the secondneighbor pure imaginary hopping-is in fact a weak integrability breaking perturbation.Note that to obtain a weak perturbation 2 ] that respects both inversion and time reversal symmetry, we need X to be invariant under inversion but odd under time reversal.Therefore, if X is an extensive local operator, it cannot be a fermion bilinear and, as a consequence, V will contain terms with more than four fermionic operators.Perturbations with only two and four fermionic operators that respect both the time reversal and inversion symmetries can instead be generated using boosted operators, for example, B[Q given by Eq. (D2).We obtain the corresponding operator, Eq. ( 8), This perturbation contains fermion hopping up to range 4 and four-fermion interactions up to range 3. Note also that we can absorb the nearest-neighbor hopping and the on-site Hubbard terms into an O(λ) shift of the parameter U , and hence the remaining parts of J (0) 3,2;tot can also be viewed as weak integrability breaking perturbations.Furthermore, we can remove the range-4 term by combining with Q (0) 4 listed in App.D, Eq. (D3), at the expense of introducing additional range-3 terms, including also six-fermion terms.However, by combining with another weak integrability breaking term generated using A. Free spinless fermions As a simple example to illustrate the deformations induced by bilocal operators, we here consider a model of free spinless fermions hopping on a chain: For simplicity, we consider the case of real nearestneighbor hoppings, but our discussion can be extended to generic complex hoppings of arbitrary ranges.We define the following set of conserved quantities with densities: 1,j = n j , q We are interested in the family of deformations of the Hamiltonian that are weak integrability breaking to first order in λ.In particular, we focus on the contributions in Q  24) and (12) we see that the linear space of deformations Q ′(1) 2 generated by bilocal operators contains in general four-fermion operators.Some examples are: As an illustration, in App.E we show the corresponding quasi-conserved quantities in each case obtained by applying the same deformations to Q 3 .Note that there are total of eight linearly-independent translationally invariant and U(1)-charge conserving Her-mitian four-fermion terms of range up to 3: Thus we see that three of these eight directions in such a space of range-3 perturbations are in fact weak integrability breaking perturbations.In fact, we know one more weak integrability breaking perturbation of range 3 obtained by using an extensive local operator X = j n j n j+1 as a generator: We can further organize these perturbations by their transformation properties under the lattice inversion and time reversal (defined as complex conjugation in the number basis).Thus, out of the eight terms, there are four terms that are invariant under both the inversion and time reversal, namely the two density-density terms, the combination U = V ∈ R, and the term W ∈ R. The U = V ∈ R term is in fact the weak integrability breaking perturbation (c), while the W ∈ R term is a combination of the perturbation (b) and the nearest-neighbor densitydensity term.Furthermore, adding just the nearestneighbor density-density interaction in fact leads to another integrable model equivalent to the XXZ chain, and by argument similar to Sec.IV B any linear combination of this term and the terms (b) and (c) will also be a weak integrability breaking of the free-fermion chain.Hence we conclude that among the inversion and time-reversal invariant range-3 perturbations only the second-neighbor density-density interaction truly breaks the integrability in the leading order.

B. Quantum Ising chain
The transverse field quantum Ising chain has the following Hamiltonian: The model can be solved by a Jordan Wigner transformation, which maps the Hamiltonian to a free fermionic model.We here list the first few conserved quantities, with densities q (0) α,j [39]: Similarly to the case of free spinless fermions, bilocal deformations can be used to generate weak integrability breaking perturbations of the quantum Ising chain that have the form of fermion interactions.Boosted operators, on the other hand, can only generate free-fermion terms.
As a first example of a perturbation induced by a bilocal generator, we consider the operator Dropping the unimportant scale factor, we conclude that ) is a weak integrability breaking perturbation.This perturbation is self-dual under the Kramers-Wannier transformation.Under the Jordan-Wigner transformation, it gives four-fermion interactions consisting of products of Majoranas on four consecutive sites.Interestingly, at the Ising critical point h = 1, the perturbed Hamiltonian Ĥ0 + λ V maps precisely to the interacting Majorana chain studied in Ref. [40].It would be interesting to revisit their study in light of our conclusion that this interaction breaks the integrability only in the next order.For example, one may wonder if this may explain the very large coupling needed to reach the tricritical point, but this requires detailed considerations which we leave for future work.
Another example of a weak perturbation of range 3 induced by a bilocal deformation is Under the Ising duality [41], the first Pauli product (including the sign) is interchanged with the fourth one and the second is interchanged with the third.Thus, this perturbation is self-dual only when h = 1, which is different from the previous perturbation Eq. (69).Furthermore, this perturbation is not invariant under an anti-unitary symmetry of the original Ising model defined as a complex conjugation in the σ x basis, so it is a bit less natural perturbation to consider.

C. Heisenberg chain
In Secs.IV A and IV B we considered weak perturbations of the Heisenberg, XYZ, and XXZ chains generated by boosted operators.Bilocal operators can be used as generators to obtain more examples of weak perturbations of these models.Once again, we focus on the perturbations that have the same symmetries as the original Hamiltonian.For example, with the choice of generator X(λ) = [Q 2 (λ)|Q 3 (λ)], the deformed charges preserve their parity under inversion and time reversal.For concreteness, we consider the deformed charges of the Heisenberg chain obtained with this generator.The first order deformation This perturbation is manifestly invariant under inversion and time reversal and preserves the SU(2) symmetry of the Heisenberg Hamiltonian.We conjecture that it is not possible to generate this weak-integrability-breaking perturbation using the boosted IoMs as generators, i.e., it genuinely requires the bilocal operators as generators.
To give some perspective on the above result, we note that the space of range-4 terms that are symmetric under spin-SU(2), lattice translation and inversion, and time reversal has dimension six.Two directions in this space (Q 4 ) are integrable, while the above four-spin term and the second-neighbor Heisenberg interaction considered in Sec.IV A are independent nearlyintegrable directions.Thus, only two out of six directions truly break the integrability at the leading order, so even including range-4 perturbations the integrability of the Heisenberg chain is more robust than one would naively expect.

VI. EXAMPLES: DEFORMATIONS TO HIGHER ORDERS
In the examples discussed so far we have examined various types of weak perturbations of integrable Hamiltonians of the form such that a set of quasi-conserved quantities Q <2 α can be defined with the property that ).It is possible to apply the general procedure in Sec.III to obtain perturbations of the Hamiltonian and of the quasi-conserved charges, such that they commute up to terms of order λ ℓ with ℓ > 2. Specifically, for ℓ = 3 we consider where the H (k) are defined as in Eq. ( 23): Here we specialized to c 2 δ α,2 intending to work around the unperturbed Hamiltonian, which in the examples we consider is commonly used to define the charge

A. Heisenberg chain
In Sec.IV A we showed that a weak first-order perturbation of the Heisenberg chain [H (0) in Eq. ( 27), with convention c 3 and c (1) α = 8δ α,2 +δ α,4 [cf.Eq. ( 36)].We now want to show how to construct H (2) such that H <3 has a set of quasi-conserved quantities up to order λ 3 .From Eq. ( 76) we get where using Eq. ( 20) with f α ] and With this definition, H <3 satisfies the desired property for any choice of coefficients f (1) β and c α .We can use this freedom to look for perturbations that have small range and involve a small number of spins.An example of a particularly simple deformation of the Heisenberg chain that is quasi-integrable to order ℓ = 3 is We refer the readers to App.A 2 for details on the specific choices of coefficients f (2) α and some intermediate steps.The degree of simplification that we managed to achieve is quite surprising given the complicated intermediate expressions, and we are wondering if there may be some reason for this.It would be interesting to check if one can achieve comparable simplification at the next order.

VII. THERMALIZATION TIME
We now discuss the implications of the existence of quasi-conserved quantities for the thermalization time.In a generic quench (with an arbitrary initial state), the local quantities Q <ℓ α (λ) are conserved by the Hamiltonian H <ℓ (λ) for at least a time t ∼ O(λ −ℓ ) (see App. F).We remark that this lower bound on the thermalization time is a completely rigorous bound using only the locality of the Hamiltonian and of the quasi-conserved observable.
On the other hand, if we use a perturbative calculation for the time-averaged rate of decay of the conserved quantity in the spirit of the Fermi's Golden Rule (see App. F for precise meaning), we would get a (non-rigorous) estimate for the thermalization time as O(λ −2ℓ ).This estimate can be intuitively understood by noting that in our construction the Hamiltonian is H <ℓ (λ) = H(λ) + O(λ ℓ ), with H(λ) integrable, and therefore the effective integrability-breaking perturbation strength is λ ℓ .On a time scale t ≪ λ −ℓ the dynamics is determined by H(λ), and the system prethermalizes to a generalized Gibbs ensemble with conserved charges Q α (λ).Genuine thermalization is triggered by the effective perturbation O(λ ℓ ), whose effect on the dynamics becomes non-negligible at longer times.Using standard estimates of the rate, but noting that the perturbation strength is λ ℓ , we then expect thermalization after time τ ∼ λ −2ℓ .
We note that this argument agrees with the direct perturbative estimates of the rate in the important case ℓ = 2: For such special weak perturbations the formal Fermi's Golden Rule rate O(λ 2 ) vanishes after a time O(1) (see App. F).The vanishing of the λ 2 order of the rate indicates that after a decay at small times, the expectation value of an observable Q (0) α (i.e., one of the original charges) reaches a plateau.Using this observation, it was argued in [30], that these perturbations do not lead to thermalization in the Boltzmann regime (i.e., in the limit λ → 0, t → ∞ with λ 2 t = const.),but lead to hydrodynamic diffusion.
The observable Q α reaching a plateau after O(1) time is consistent with a picture where α that does not thermalize until a much later time.In fact, it is much easier to see the vanishing O(λ 2 ) rate by thinking directly about the quasi-conserved is identically zero at any time (and not just after O(1) time).
In App.F we also consider formal O(λ 3 ) term in the rate of change of Q <2 α and show that it vanishes after O(1) time.Hence the leading non-vanishing rate after O(1) time is actually O(λ 4 ), in agreement with our picture and the intuition based on the effective perturbation of the integrable model H(λ).The above statements about the rates are for initial states or ensembles generated by the unperturbed integrable model (i.e., defined by {Q (0) α }).In App.F we also show initial ensembles where the rate of change of the quasi-conserved Q <2 α starts at O(λ 4 ) for all times, which can be viewed as formalizing the above intuition by also finding appropriate "quasi-stationary" initial states.Thus we have established a good understanding of the connection between direct perturbation theory calculations in the case of such nearly-integrable perturbations and thinking in terms of the quasi-conserved quantities, observing both the physical intuition and analytical power of the latter framework.

VIII. CONCLUSIONS
We have here demonstrated a general construction for producing weak (i.e., nearly integrable) perturbations of integrable lattice models to arbitary order λ ℓ .These weakly-perturbed models have an extensive number of (extensive local) approximate conserved quantities that commute with the Hamiltonian up to corrections O(λ ℓ ).
We have applied the construction to several well-known spin chains and fermionic models, focusing on finding particularly simple such examples and also putting some previously known instances and their physics into a unified framework (e.g., the relation between the vanishing of the Fermi's Golden Rule rate and the presence of the quasi-IoM).
In the case of truncation to linear order, i.e., achieving commutation up to corrections O(λ 2 ), this approach produces a linear space of weak integrability breaking perturbations ("diffusive subspace" [30], as opposed to generic "thermalizing" perturbations).While this subspace is of course measure zero in the full space of possible perturbations, it may nevertheless play an important role in realistic physics applications.For example, for the spin-1/2 Heisenberg chain, it turns out that all SU(2)-spin symmetric perturbations that are translationally invariant and have range up to 3 (i.e., involving up to three consecutive spins) are in fact either integrable or weak integrability breaking perturbations.As discussed at the end of Sec.V C, for range-4 perturbations, if we also require lattice inversion and time reversal symmetry, three out of five perturbations (not counting 2 ) of the Heisenberg chain are either integrable or weak integrability breaking.Such unexpected paucity of natural true integrability breaking perturbations of the Heisenberg chain may help explain the robustness of the superdiffusion signatures in numerical studies [42] as well as in physical world experiments [43,44].
As another example in the same spirit, for the spinless fermion chain, in the four-dimensional space of inversion and time-reversal symmetric four-fermion interactions up to range 3, three directions are in fact weak integrability breaking perturbations.In general, when we study effects of a given perturbation, we would want to understand/remove appropriate "components" onto the weak integrability breaking ones, since only the remaining part represents generic integrability breaking perturbation that dominates the thermalization rate of the sys-tem.Furthermore, in some situations such a removal may not be possible, e.g., under sufficiently restrictive range and symmetry conditions, and it is important to understand the thermalization rates in such cases as well.This shows importance of systematic constructions of weak integrability breaking perturbations, and we hope that our work will stimulate further such studies.
A very interesting Ref. [45] proposed to analyze perturbations to integrable and non-integrable systems by constructing so-called Adiabatic Gauge Potential (AGP) [46] which is an analog of the generator X(λ) in our formalism.While a regularized AGP always exists, for true integrability breaking perturbations it is expected to be highly non-local, and they proposed that a particular norm of the AGP scales exponentially with the system size and provided strong evidences for this.On the other hand, for perturbations along exact integrabilitypreserving directions like varying the anisotropy in the XXZ chain, they found that the AGP norm scales polynomially with the system size.They also noted that perturbations like the ones here generated by X e or X bo , while breaking the exact integrability and resulting eventually in the chaotic behavior, nevertheless also have the AGP norm scaling polynomially with system size, and they suggested to exclude such special perturbations when checking for quantum chaos [47].
Interestingly, in our formalism both the integrabilitypreserving perturbations and the integrability-breaking perturbations generated by X e , X bo , and X bi are formally ℓ = 2 weak perturbations, i.e., can be viewed as falling into the same group; likewise, they share a similar polynomial scaling of the AGP norm.Of course, the latter perturbations will eventually lead to thermalization, and while they do not behave like other "truly generic" perturbations, it is also interesting to explore how thermalization happens under such "less generic" weak perturbations: We discuss this in Sec.F, where we argue that the relaxation times have different parametric dependence on the perturbation strength λ.Furthermore, our formalism provides recipes to tabulate families of such weak integrability breaking perturbations beforehand, and such tables can then be used to systematically exclude them when needed for studies of more generic thermalization phenomena.While we do not provide complete tabulations (which is important future work), we showcase many examples; interestingly, already we find that the number of such special perturbations can be significant enough for them to be of practical importance, as discussed earlier.
Our results suggest several intriguing directions for further studies.One interesting possibility is the study of transport properties of integrable models with weak perturbations [48,49].Recent developments in the study of transport in integrable systems, in particular within the framework of generalized hydrodynamics, revealed different regimes, including ballistic, diffusive, and also anomalous superdiffusive [50][51][52].Small integrability breaking perturbations induce scattering of quasi-particles, affect-ing the transport properties of the system.It is an interesting question to understand how weak perturbations compare to ordinary perturbations in this context [42,53].For example, Ref. [42] showed that the superdiffusive behavior of the Heisenberg model persists up to all numerically accessible times when a small to sizable next-nearest neighbor Heisenberg interaction (a weak integrability breaking perturbation) is included.However, it is not clear whether the "weak property" of the perturbation plays a role, since similar persistence is observed for other SU(2) preserving perturbations that do not have this property.
One of the possibilities offered by our approach is the ability to systematically compute all the approximately conserved quantities.The construction could be applied also to the quasi-local charges, that were first discovered in the XXZ chain [54].It was shown that, in many cases, taking into account these quasi-local charges is crucial to obtain the correct results, for example in the computation of the Drude weight using Mazur bounds, and of the expectation values of the generalized Gibbs ensembles [55,56].We expect that the approximatelyconserved quasi-local operators that can be computed with our construction are similarly important in determining the properties of weakly-perturbed integrable models.
Another interesting question regards systems with finite size.We argued that weak perturbations correspond to longer thermalization times in the non-equilibrium dynamics in the thermodynamic limit.However, the signatures of thermalization are present also in the opposite order of limits [57]: The transition from an integrable to a thermalizing regime can be captured from the energy spectrum at finite size, for example, by studying the crossover in the level spacing distribution from a Poisson to a Wigner-Dyson statistics.A recent work [29] proposed that a weak perturbation (such as the ones generated by boosted operators) corresponds to a different scaling of the position of the crossover with system size and provided numerical evidence in a perturbed XXZ chain.The crossover to a Wigner-Dyson statistics was also studied in Ref. [58] for the Heisenberg spin chain perturbed with the next-nearest-neighbor interactions, which found that the crossover occurs at larger coupling than for generic perturbations.It would be interesting to study this for different types of generators X e , X bo , and X bi , for different models and boundary conditions, and also see if there is relation to the AGP norms.
We note that, while we focused on translationally invariant systems, weak integrability breaking perturbations can be also inhomogeneous.The simplest example is when we turn an extensive local generator X e into an inhomogeneous one or even a strictly local operator, which produces an inhomogeneous or a strictly local weak integrability breaking perturbation.There are also variants starting from boosted generators [59], but so far we have not been able to achieve this starting from bilocal generators.We leave a systematic study of possible in-homogeneous weak integrability breaking perturbations for future work.
While we focus here on integrable Hamiltonians, some of the most relevant experimental realizations of integrable models are Floquet integrable models realized with circuits of gates [60][61][62][63].It would be interesting to extend our construction to such models: While the procedure can be immediately generalized for deformations induced by extensive local operators, it is unclear to us how to generate weak perturbations that preserve the simple gate structure using boosted or bilocal operators.
We remark that the construction that is the object of our study applies not only to integrable models, but it can also be used to generate weak perturbations of onedimensional models with a finite number of conserved quantities.A possible question is whether a similar procedure can generate perturbations that are weak only for a specific subspace, as such perturbations would correspond to models with prominent non-exact quantum many-body scars.In fact, the tower of scars of the PXP model is an example of such non-exact quantum manybody scars [64][65][66][67][68][69][70][71], and the mechanism that protects their subspace is still unclear.We also note that numerical results on some exact quantum-many body scars [72] show that they might be robust to lowest order in the perturbation strength [73,74], suggesting that a notion of weak perturbation may be formulated also for single eigenstates.Moreover, a recent experiment has found that some eigenstates of a Floquet integrable model [62] decay very slowly in the presence of an integrability-breaking perturbation [63].It is not yet clear if these states are robust to first order in the perturbation strength.
Finally, while we focused on the case of onedimensional integrable models, one can think of constructing similar weak perturbations in higher dimensions.In this case, extensive local and boosted operators can be used as generators [75], but we are not aware of other classes of operators -analogous to bilocal operators in one dimension -that would produce physical (i.e., local and extensive) perturbations. of course can be added with arbitrary coefficients.Reference [21] found only oneparameter family given by the above equation with the parameter a, probably because their search required the coefficient of ⃗ σ j • ⃗ σ j+1 in the very final expression to be zero.While the appearance and count of such free parameters is somewhat mysterious in their formalism, it is not from our perspective where all operators of the form Eq. (A7) are formally quasi-IoMs to order ℓ = 2 in the sense of Eq. (15).Equation (A7) clearly shows that there can be infinitely many quasi-IoMs associated with each original IoM Q α .While formally for fixed α they are linearly independent for linearly independent d (1) α,β , it is not to say that physically they are equally important (and they are not independent when varying α).The perturbative setup keeping λ small works well when λ multiplies objects that do not get large themselves, and which specific quasi-IoM is best to use can depend on the context (e.g., Ref. [21] used conditions of fixed range and minimization of the Frobenius norm of the commutator with the perturbed Hamiltonian, which is reasonable for certain types of quenches [76]).

Details of constructing ℓ = 3 nearly-integrable perturbation to the Heisenberg chain
Here we provide some details for constructing particularly simple higher-order nearly-integrable perturbation to the Heisenberg chain.To evaluate Eq. ( 78), we first calculate 2 ] and i[B[Q 2 ] are in general not extensive local operators, since [Q the Hamiltonian is an extensive local operator: j jh j,j+1 , H 0 = − j g j,j+1,j+2 . (C4) The underlying reason for this is the energy conservation, and the R.H.S. is (proportional to) a sum of local energy currents.The above is valid for any Hamiltonian.
From now on we will consider integrable Hamiltonians whose IoMs can be generated using the boosted Hamiltonian as as in the Heisenberg and XYZ spin chains in the main text.Following Ref. [78], let us consider when thus defined Q 3 can commute with H 0 .We have where we used the expression for g's in terms of h's to rewrite [h j,j+1 , g j−2,j−1,j ] = −[h j−2,j−1 , g j−1,j,j+1 ] (since h j−2,j−1 and h j,j+2 commute).To evaluate [H 0 , Q 3 ] we need to sum Eq. (C5) over j.The first two terms have a "telescoping structure," hence their contributions will cancel upon the summation.No such structure is present for the last two terms for general Hamiltonians; however, we will have the telescoping structure if so-called Reshetikhin condition [78] is satisfied: There exist two-site operators R j,j+1 such that [h j,j+1 +h j+1,j+2 , g j,j+1,j+2 ] = R j,j+1 −R j+1,j+2 .(C6) This condition is satisfied for the Heisenberg and XYZ chains, and in general it implies [H 0 , Q 3 ] = 0.
If the Reshetikhin condition is satisfied, then it is also clear that Q 4 as defined above is an extensive local operator.A straightforward calculation using Eqs.(C5) and (C6) gives j jh j,j+1 , In the main text we considered perturbed models generated using X ∼ B[Q 3 ], which led to special perturbations ∼ [B[Q 3 ], H 0 ].Reusing some calculations behind Eq. (C5), we have Summation over j using the telescoping structure in the first two terms and the Reshetikhin condition for the second two terms then gives j h j,j+1 , We see that in this total current associated with the conserved quantity Q 3 , the range-4 terms have exactly the same structure as in Q 4 in Eq. (C7).Hence we can obtain a simpler nearly-integrable perturbation by combining Eqs.(C9) and (C7): Connecting with the notation in the main text used in the Heisenberg and XYZ cases, Secs.IV A and IV B, we have: = − i 4 j g j,j+1,j+2 , (C12) 2 ≡ −i = − 1 8 In the preceding discussions and expressions, everything is still general other than the convention Q used to obtain more simple nearly-integrable perturbations in the Heisenberg and XYZ models in the main text, Eqs.(36) and (42) respectively, in fact works for all such integrable models, since it is precisely the combination in Eq. (C10) that cancels the four-site terms: [h j+1,j+2 , g j,j+1,j+2 ]+ 3 2 R j,j+1 .

(C15)
At time t we have where in the intermediate steps L is the system size and we have used the translational invariance of H, M , R, and |Ψ ini ⟩.We have | ⟨Ψ ini | e iHt ′ r j e −iHt ′ |Ψ ini ⟩ | ≤ ∥r j ∥ ≤ c r , and hence Since we expect that the "thermalized" value of the local observable at long time, lim t→∞ ⟨Ψ(t)| m j |Ψ(t)⟩, differs from the initial value, ⟨Ψ ini | m j |Ψ ini ⟩, by a non-zero O(1) amount ∆m th , we conclude that it will take at least time to thermalize, which is parametrically large for small λ.
Note that this is a completely rigorous lower bound on the thermalization time valid for any L (hence also in the thermodynamic limit L → ∞) that makes no assumptions about the thermalization physics of H and M other than their extensive local character and the nearcommutation Eq. (F2).Note also that no assumption is made about the initial state |Ψ ini ⟩ with respect to H and M ; the bound on the rate of change of the observable m j is always valid, and the "thermalization" assumption is only made to convert the rate to a finite time assuming an O(1) change in the observable under the dynamics.
In the case of eigenvalue degeneracies, it is natural to consider an initial ensemble where the states {|ϕ (0) k ⟩} appear with probabilities that depend only on their H (0)  and M (0) eigenvalues, {p k = f (ϵ k )}.For such an ensemble, using Eq.(F12) we can easily see that hence the initial rate of change of the observable M is 0.
Note that in the above we have only used the quasicommutation condition Eq. (F2), specialized to Eq. (F8) in the present case.Suppose we further know that the quasi-commutation derives from a truncated unitary rotation, i.e., there exists X such that A straightforward calculation utilizing the commutation of H (0) and M (0) and the common eigenstate condition on |ϕ k ⟩ then gives ⟨ϕ k ⟩ = 0 without the nondegeneracy assumption.Indeed, we simply write out and observe that the above terms can be grouped in pairs that cancel each other when evaluated on any common eigenstate |ϕ k ⟩, e.g., the first and the last term, etc.
b. Assumptions about the dynamical correlation functions Let us be more specific about the precise conjecture used to argue the vanishing of FGR rate and various assumptions motivating it.
First, because of the locality of H (0) (which we always assume), the Lieb-Robinson bound implies that at any fixed t the dynamical correlation functions like Eq. (F20) are well defined in the thermodynamic limit, and that with such correlation functions the L.H.S. of Eq. (F21) is a convergent sum over the real-space separations |j ′ − j|; this means that we have a well-defined formal O(λ 2 ) FGR rate of change of the local observable m (0) j at any t in the thermodynamic limit.Our main conjecture is that the sum in Eq. (F21) -and hence this FGR rate -vanishes for sufficiently large t.
Next, we discuss the underlying assumptions about the dynamical correlation functions.While the "factorization" in Eq. (F20) (i.e., vanishing of the connected correlation functions) is intuitively reasonable, we do not know of a proof in general.Several recent papers [81][82][83] obtained some rigorous results about such "factorization" for general Hamiltonians.However, these results are either in different regimes (e.g., correlations of local observables in the t → ∞ limit first before the thermodynamic limit), or are too weak for us to use (e.g., bounds on correlations of extensive local observables).Nevertheless, our situation is better in that we actually need the specific difference of correlation functions which corresponds to the commutator in Eq. (F20).For fixed t and large enough |j ′ − j| > v LR t, i.e., outside the Lieb-Robinson cone where v LR is the corresponding velocity, we expect this to decay exponentially with |j ′ − j|, so the sum in Eq. (F21) is convergent at any t.Whether this sum indeed approaches zero at large t depends on the behavior of the correlation functions within the Lieb-Robinson cone |j ′ −j| ≲ v LR t.Unfortunately, we do not know sufficiently strong general results for such correlations.Still, if we can establish that, e.g., correlations at the same location j ′ = j decay sufficiently quickly with time, it is reasonable to assume that the total contribution from all |j ′ − j| ≲ v LR t also decays to zero with time.
Furthermore, as already mentioned, these are solely questions about the unperturbed H (0) and may be directly answerable for solvable H (0) like all the integrable models considered in this paper.Thus, for free-fermion unperturbed Hamiltonians like the ones in Sec.V, either for the ground state or thermal initial ensembles, one can employ Wick's theorem and carry out such calculations exactly.While we have not done the full calculations for the corresponding near-integrable perturbations H (1)  and M (1) as observables, toy calculations with simpler local observables (e.g., a j = c † j c † j+1 , b j ′ = c j ′ c j ′ +1 , for the spinless fermion hopping problem) suggest that the above conjecture indeed holds, and we leave full calculations for future work.
For interacting integrable models, Ref. [22] studied thermal conductivity in the nearest-neighbor Heisenberg chain perturbed by the second-neighbor Heisenberg terms; using numerical high-temperature series and exact diagonalization estimates they found that the corresponding FGR-like O(λ 2 ) contribution indeed vanishes in this case.Also, Ref. [30] appealed to the hydrodynamic projection principle calculations of the dynamical correlation functions in integrable models to argue that the FGR-like rate of change of an unperturbed IoM indeed vanishes after some initial time, for special perturbations that are equivalent to the ones in our work obtained using generators X that are boosted operators.More precisely, one can start with Eq. (F18) and use relation ) t −e iH (0) t M (1) e −iH (0) t ) |Ψ 0 ⟩, which can be checked, e.g., by differentiating both sides with respect to t and utilizing Eq. (F8).The latter form matches expression analyzed in Ref. [30], which then appealed to the hydrodynamic projection principle to suggest that the corresponding positive time and negative time correlations effectively cancel each other at large t.While these calculations are specific for integrable models, we think that our conjecture does not require integrability and holds for generic local Hamiltonians H (0) , perhaps under mild additional conditions, but leave more tests (e.g., by numerical methods) to future work.
Finally, we emphasize that while fully justifying the vanishing of the O(λ 2 ) rate for the M (0) required additional arguments like the ones in this subsection, no such arguments were needed when considering the O(λ 2 ) rate for the properly corrected (i.e., quasi-conserved) observable M .Returning to the dynamics of the full M and Eq.(F17), in the case H (1) = i[X, H (0) ] and using that ⟨Ψ 0 |[i[X, H (0) ], •]|Ψ 0 ⟩ = ⟨Ψ 0 |[X, i[H (0) , •]]|Ψ 0 ⟩, we have ⟨Ψ 0 | [H (1) , e iH (0) s Re −iH (0) s ] |Ψ 0 ⟩ (F22) = ⟨Ψ 0 | [X, i[H (0) , e iH (0) s Re −iH Furthermore, if we assume "factorization" of the dynamical correlation functions at long times, similar to Eq. (F20) and conjecture Eq. (F21), we then have ⟨Ψ 0 | [X, e iH (0) t Re −iH (0) t ] |Ψ 0 ⟩ ≈ 0 (F25) for large t.Hence after some initial time, the formal O(λ 3 ) term in the rate of change of M vanishes.The non-vanishing rate is then O(λ 4 ), in agreement with the expectations of FRG reasoning applied to an effective strength of true integrability breaking being O(λ 2 ).Some remarks are in order.Note that we also needed such arguments appealing to "factorization" of the dynamical correlations to show the vanishing of the O(λ 2 ) rate of change of M (0) after some initial time, while we did not need such arguments for the O(λ 2 ) rate of change of the "corrected" (quasi-IoM) M -the formal O(λ 2 ) rate of change of M is identically zero at any time.We think that we need such arguments for the O(λ 3 ) rate of change of the corrected M because the initial state |Ψ 0 ⟩ is not "corrected" to reflect that the true integrability breaking perturbation has effective strength O(λ 2 ).
Indeed, let us return to Eq. (F9) valid for any initial state, and instead of |Ψ 0 ⟩ [a common eigenstate of H (0)  and M (0) as before], let us start with where the O(λ 2 ) and O(λ 3 ) terms vanish by Eqs.(F10) and (F24).Thus, for such initial states or ensembles, the leading contribution to the rate of change of M is O(λ 4 ) at any time t.One may worry that such |Ψ ini ⟩ or ρ ini are not readily preparable (e.g., experimentally preparable) from |Ψ 0 ⟩ or ρ 0 .However, if X is a sum of on-site terms (as in the example in Sec.VIII of the on-site Hubbard term used as a generator starting with free fermions in any dimension), then e iλX is a product of on-site unitaries and therefore is a depth-1 unitary circuit, hence such initial states or ensembles can be viewed as preparable from the unperturbed ones.For a general extensive local X, we can appropriately Trotterize e iλX and obtain a finite depth unitary circuit U ′ (λ) that reproduces e iλX with O(λ 2 ) accuracy (e.g., e iλX ≈ e iλX even bonds e iλX odd bonds ≡ U ′ (λ) familiar for a 1d chain with only nearest-neighbor bond terms in X); in this case we can initialize with the corresponding |Ψ ′ ini ⟩ = U ′ (λ) |Ψ 0 ⟩ or ρ ′ ini = U ′ (λ)ρ 0 U ′ (λ) † differing from the above |Ψ ini ⟩ or ρ ini by O(λ 2 ), and observe that with such preparable initial states or ensembles we also obtain vanishing formal O(λ 2 ) and O(λ 3 ) terms in the rate of change of M .[Note that these |Ψ ′ ini ⟩ or ρ ′ ini are not the naive truncated series approximations |Ψ 0 ⟩ + iλX |Ψ 0 ⟩ or ρ 0 + iλ[X, ρ 0 ] that would be problematic since X is an extensive operator.]One should worry more about the case where the generator X is a boosted or bilocal operator, hence is not a regular local extensive operator.While we do not have a full resolution of this concern, we observe that the commutator of such X with H (0) and M (0) gives regular extensive local operators; hence, the non-locality of X is mitigated when one considers, e.g., ρ 0 that is an equilibrium ensemble under H (0) and M (0) : −γM (0)  , Z 0 ≡ Tr e −βH (0) −γM (0)  ,

(F32)
Hence, for such ρ ′ ini we expect the formal rate of change of M is O(λ 4 ).Note that this is not an equilibrium ensemble for H = H (0) + H (1) , since H does not commute exactly with M = M (0) + M (1) .Intuitively, we expect such ensembles to describe prethermalized states in our system with quasi-conserved M , while here we use ρ ′ ini to illustrate initial ensembles that will "know" from the outset that the true integrability breaking strength is O(λ 2 ) and where any reference to X has dropped out.Indeed, we can verify this by a direct calculation:

2 ≡ H 0
/2, i.e., applicable for any integrable model satisfying the Reshetikhin condition.The specific combination of Q c. Vanishing of the O(λ 3 ) term in the rate of change of M at long time