Replica resummation of the Baker-Campbell-Hausdorff series
Szabolcs Vajna Katja Klobas Tomaz Prosen Anatoli Polkovnikov
cond-mat.str-el cond-mat.quant-gas cond-mat.stat-mech quant-ph



Department of Physics, Faculty of Mathematics and Physics, University of Ljubljana,  SI-1000 Ljubljana, Slovenia
Department of Physics and BME-MTA Exotic  Quantum  Phases Research Group, Budapest University of Technology and Economics, 1521 Budapest, Hungary

Department of Physics, Faculty of Mathematics and Physics, University of Ljubljana,  SI-1000 Ljubljana, Slovenia

Department of Physics, Faculty of Mathematics and Physics, University of Ljubljana,  SI-1000 Ljubljana, Slovenia

Department of Physics, Boston University, 590 Commonwealth Ave., Boston, MA 02215, USA



We developed a novel perturbative expansion based on the replica trick for the Floquet Hamiltonian governing the dynamics of periodically kicked systems  where the kick strength is the small parameter. The expansion is formally equivalent to an infinite resummation of the Baker-Campbell-Hausdorff series in the un-driven (non-perturbed) Hamiltonian, 
while considering terms up to a finite order in the kick strength. As an application of the replica expansion, we 
analyze an Ising spin 1/2 chain periodically kicked with magnetic field of strength h, which has both longitudinal and transverse components.
We demonstrate that even away from the regime of high frequency driving, the heating rate is nonperturbative in the kick strength bounded from above by a stretched exponential: e^- const h^-1/2. This guarantees existence of a very long pre-thermal regime, where the dynamics is governed by the Floquet Hamiltonian obtained from the replica expansion.

Replica resummation of the Baker-Campbell-Hausdorff series
    Anatoli Polkovnikov
    December 30, 2023
==========================================================







Introduction.—
Time-periodic modulation of interactions is a powerful tool to engineer properties of materials in both artificial and condensed matter systems <cit.>. In particular, high frequency driving is the cornerstone of various experiments and proposals inducing interactions such as spin-orbit coupling <cit.>, artificial gauge fields for uncharged particles <cit.>; it has been applied to dynamically tune or suppress hopping amplitude in optical lattices <cit.>, and also to change topological properties of materials <cit.>.

Given a periodic driving protocol, however, determining the Floquet Hamiltonian that governs the stroboscopic evolution is usually a highly non-trivial task. Except for some special integrable cases <cit.>, one is compelled to apply approximate methods, e.g. variants of high frequency expansion (Magnus <cit.>, van Vleck <cit.> or Brillouin-Wigner <cit.> expansions). These provide a local effective Hamiltonian in each order of the expansion, however, until recently, not much had been known about the convergence properties of these series. A conjecture based on the generalization of the eigenstate thermalization hypothesis suggests that generic closed periodically driven systems heat up in the thermodynamic limit, i.e. they approach a completely structureless, infinite temperature steady state <cit.>. The convergence of the expansions of the effective Hamiltonian is intimately related to heating. Recently upper bounds on heating had been reported in the linear response regime <cit.> and for the Magnus expansion <cit.>, with the central result that the heating is at least exponentially suppressed in the driving frequency for periodically driven models characterized by local Hamiltonians with bounded energy spectrum. This theoretical finding implies that one can engineer nontrivial phases of matter, which remain stable for the experimentally relevant timescales. 
In some situations heating seems to be either absent completely or remain well below exponential bounds <cit.>. Another recent theoretical work showed that nontrivial non-equilibrium Floquet phases can be stabilized by weak coupling to environment <cit.>.



One of the most studied driving protocols is a time-periodic sequence of sudden quenches between different Hamiltonians <cit.>,
which is interpreted as kicked dynamics if one of the time intervals on which the Hamiltonians act is much shorter than the other, or equivalently, the strength of one of the Hamiltonians is smaller than the other. Such protocols naturally appear e.g. in the context of digital quantum simulation in trapped ions <cit.> and are frequently realized in other experimental platforms (see e.g. Refs. <cit.>).

In the present work a novel expansion for the effective Floquet Hamiltonian is introduced for periodically kicked systems, which clearly outperforms the traditional high frequency expansions in a wide parameter range, as illustrated using two examples shown in Figure <ref>. Our approach uses the replica trick to calculate the logarithm of the time evolution operator describing a single period. The small parameter is the kick strength and we do not assume high frequency driving. The Magnus expansion is equivalent to the Baker-Campbell-Hausdorff (BCH) series in periodically kicked systems, and the replica expansion can be thought of as an infinite resummation of the BCH formula.  As such, the possible applications of the replica expansion can reach far beyond periodically driven systems, including the theory of differential equations <cit.>, Lie group theory <cit.>, analysis of NMR experiments <cit.> or estimation of Trotterization errors in various numerical integration schemes. As an application of the replica expansion, we establish a conjecture for a non-perturbative – stretched exponential exp(- const h^-1/2) in the kick strength h – upper bound for the heating rate in the kicked tilted field Ising model. Our conjecture supplements the bounds introduced in Refs. <cit.>, which address only the regime of high frequency driving. 

Kicked dynamics.—
The effective Floquet Hamiltonian evolving the kicked system is given by the logarithm of the stroboscopic time evolution operator over one period U=U_0 U_1=e^-i J H_0e^-i h H_1 as 

    H_F=ilog(e^-i J H_0e^-i h H_1)  ,

where we incorporated the time intervals of the Hamiltonians H_0,1 to the coupling constants J,h. The BCH formula provides a series expansion for H_F assuming both J and h are small as 

    H_F=   J H_0+h H_1-i J h 1/2[H_0,H_1]-
       Jh1/12(J[H_0,[H_0,H_1]]+h[H_1,[H_1,H_0]]) +…

In the usual setup of kicked systems, where J≫ h, that is, at intermediate frequencies and weak kick strengths, the terms with low order in H_1 but high order in H_0 are not negligible. A series expansion in the small parameter of h could be obtained formally by a resummation of the BCH series in H_0. The first order resummation is well known <cit.>,

    H_F=J H_0+-i J _H_0e^-i J_H_0/e^-i J_H_0-1h H_1+(h^2)

where _X(Y)=[X,Y] is the Lie derivative. However, to the best knowledge of the authors, closed form expressions for the infinite resummation in higher orders of H_1 have not yet been reported in the literature.

Replica expansion.—
We tackle this problem by constructing a series expansion in h in Eq. (<ref>). Because of the noncommutativity of H_0 and H_1, the higher order derivatives of the logarithm of the time evolution operator cannot be obtained easily.  To circumvent this obstacle, we apply the replica trick to express the logarithm,   

    log U=lim_ρ→ 01/ρ(U^ρ-1).
 
This idea has been proven to be uniquely useful in various fields of science, such as in the statistical physics of spin glasses <cit.>, machine learning <cit.>, and also in calculation of the entanglement entropy <cit.>. 
Assuming that the replica limit (∙)≡lim_ρ→ 01/ρ(∙) commutes with the differentiation, the series expansion of the Floquet Hamiltonian in Eq. (<ref>) reads 

    H_F^(n)=∑_r=0^n h^r Γ_r ,

with 
Γ_0=J H_0 and 
Γ_r= 1/r!∂^r_h U^ρ.
The derivatives of the powers of the time evolution are easy to calculate at integer values of the replica index ρ, and the replica limit is taken following an analytical continuation to arbitrary real values. After evaluating the derivatives and neglecting the prefactor U_0^ρ, which approaches the identity operator in the replica limit, the r^th order correction reads 

    Γ_r=(-i)^r-1/r!∑_0≤ m_1≤… m_r<ρm_rm_r-1…m_1 c_m_1 … m_r

where m=U_0^-m H_1 U_0^m and c_m_1 … m_r=r!/n_0!n_1!… is the multinomial coefficient with n_s being the number of indices taking value s.
Our method provides a remarkably simple derivation of the known first order term in Eq. (<ref>) giving a certain degree of confidence in the replica expansion (see <cit.>).

It is not clear directly from Eq. (<ref>) if higher order terms can be expressed as nested commutators (multiple Lie derivatives), which would be expected from the resummation of the BCH series. One can show that the corrections to the Floquet Hamiltonian can be represented as sums of terms containing the commutators and extra terms, which vanish in the zero replica limit. For example at second order 

    Γ_2=-i/2{∑_0 ≤ m_1 ≤  m_2<ρ [m_2,m_1]+ (∑_0 ≤ m_1 <ρm_1)^2} ,

where the second term is proportional to ρ^2, assuming that the replica trick works at first order, that is, ∑m∼ρ. Similar transformations that produce the nested commutator expression are given in the supplemental material up to the 5^th order.

The resulting expansion of the replica Floquet Hamiltonian expressed in terms of nested commutators reads


    Γ_0   =J H_0
    Γ_1   =∑_0≤ m < ρm
    Γ_r   =(-i)^r-1/r!∑_0≤ m_1≤… m_r<ρ [m_r,… [m_2,m_1]]] c_m_2 … m_r



where c_m_2 … m_r=(r-1)!/n_0!n_1!…. The expansion can be constructed similarly for different initial phases of the driving, U'=e^-i J H_0 (1-φ)e^-i h H_1 e^-i J H_0 φ=e^-i H_F', which results in the same equations as Eq. (<ref>) except for a simple substitution m_i→m_i+φ.
Having established the first main result of this Letter, we now demonstrate its performance in the example of the kicked Ising model in a tilted field <cit.>.

Kicked Ising model.— 
The time evolution is characterized by time-periodic quenches between the Hamiltonians H_0,1,


    H_0   =∑_i σ_i^z σ_i+1^z
    
     H_1   =∑_i σ_i^x+ σ_i^z ,



where  and  are shorthand notations for cosθ and sinθ. The purpose of introducing the tilt angle is to break the integrability of the model  at  θ= 0 and π/2 <cit.>. Figure <ref> shows the performance of the replica expansion for the kicked Ising model in two different limits: when the kick parameter is the magnetic field h, or the Ising interaction J. 
The spectral norm of the difference between the approximate and exact time evolution operators, Δ_n=U-U^(n) shown in Figure <ref>, bounds the accuracy of the expansion for the dynamics of any observable A, as |A(t)-A_n(t)|≤ 2 t AΔ_n +(Δ_n^2), where A(t) (A_n(t)) is the expectation value of the observable following t periods with respect to the exact (approximate) time evolution, starting from an arbitrary initial state [Although Δ_n scales linearly with the system size L, this accounts for expectation values of nonlocal operators. The accuracy of the dynamics of local observables is independent of system size <cit.>.]. 
As the two cases, kicking with H_0 or H_1, show qualitatively similar behavior, we discuss here only the replica expansion for kicking magnetic field.


The time evolution of kicking Hamiltonian H_1 with respect to unperturbed dynamics H_0 for m periods reads explicitly

    m= ∑_i    σ_i^z-(1/2sin(4mJ)(σ_i-1^zσ_i^y+σ_i^yσ_i+1^z)+ 
       sin^2(2mJ)σ_i-1^zσ_i^xσ_i+1^z -cos^2(2mJ)σ_i^x)  ,

and is the main building block of the replica expansion. The computation of the nested commutators of these objects is trivial, following which one has to deal with the multiple sums and the replica limit. It is convenient to separate the operator part of the expansion from the replica coefficients by expressing H̃_m in Fourier harmonics as m=_0+e^i 4 J m_1+e^-i 4 J m_-1, where


    _0   =∑_i σ_i^z+/2(σ_i^x-σ_i^z σ_i+1^xσ_i+2^z)
    _± 1   =∑_i /4[σ_i^x+σ_i^z σ_i+1^xσ_i+2^z± i (σ_i^zσ_i+1^y+σ_i^yσ_i+1^z)],


which brings us to the simplest formulation of the replica expansion,

    Γ_1   =∑_x_1_x_1_x_1
    Γ_r   =(-i)^r-1/r!∑_x_1,x_2,…,x_r_x_1 x_2 … x_r [_x_r,… [_x_2,_x_1]]

where x_i∈{0,± 1} and we introduced the replica sum

    _x_1 x_2 … x_r   =∑_0≤ m_1 ≤… m_r < ρ  e^i 4J m_1 x_1e^i 4J m_2 x_2… e^i 4J m_r x_r c_m_2 … m_r .

These sums are evaluated gradually (with attention to the combinatorial factors) as 

    ∑_m_j=0^m_j+1-1 m_j^y e^i  m_j   =(∂/i ∂_)^ye^i  m_j+1-1/e^i -1 ,

with m_r+1=ρ, and J̃ is an integer multiple of 4J. The prefactor m^y, 0≤ y∈ℕ may arise from the previous sum with respect to m_j-1, e.g. from the sum of constant terms. This way of evaluating the sums already defines the analytical continuation to arbitrary real values of ρ, allowing one to take the replica limit . 
This analytical continuation leads to e^i ρ-1=log e^i =i, which tries to enforce a Floquet Hamiltonian continuous in J at a price of breaking the periodicity H_F(J)=H_F(J+2π). Alternatively, one can choose a different branch of the logarithm, e.g. which folds J into the interval (-π,π] by applying a different analytical continuation <cit.>. This ambiguity in choosing the branch of the logarithm can be potentially used to further improve the expansion, restore the periodicity in J and eliminate divergences which are discussed below.

The sum is especially simple in the first order correction: _0=0, _± 1=2J( 2J ∓ i), yielding 
 
    Γ_1=∑_i     a_+σ_i^x+ a_-σ_i-1^zσ_i^x σ_i+1^z+σ_i^z+
        J (σ_i-1^zσ_i^y+σ_i^yσ_i+1^z)

with a_±=(J 2J ± 1/2). 

The second order correction is written in a compact form by noticing that  _± 1^†=_∓ 1 and _x_1,x_2^*=_-x_1,-x_2,

    Γ_2=-i/2{(_10-_01)[_0,_1]+_1-1[_-1,_1]}+h.c.

The replica coefficients are evaluated as

    _10-_01   =(1-2J  2J)(1+i  2J)
    _1-1   =1/2-isin 4J-4J/4sin^2 2J ,

which finally yields

    Γ_2=   ∑_i  a_- (σ_i^y σ_i+1^x σ_i+2^z+σ_i^z σ_i+1^x σ_i+2^y) - 
        a_- [σ_i^y+σ_i^z σ_i+1^y σ_i+2^z+2J(σ_i^x σ_i+1^z +σ_i^z σ_i+1^x )]+
    
       (b+c)σ_i^z σ_i+1^xσ_i+2^x σ_i+1^z -b σ_i^y σ_i+1^y-c σ_i^z σ_i+1^z .

The coefficients are b=^2/84Jcos 4J-sin 4J/sin^2 4J and c=^2/44J-sin 4J/sin^2 4J. The higher order corrections can be  calculated similarly <cit.>.


Notice that the first order correction diverges near J_k,1=kπ/2, which was identified as a signal of a heating (or nonergodicity–ergodicity) transition in a different spin model <cit.>, similar to the divergence of high-temperature expansion signaling phase transition in statistical physics. 
The source of this divergence is easily identified as the zero of the denominator in Eq. (<ref>). Similar to the high frequency expansion, the higher order corrections become less and less local due to the increasing number of commutators.  
The degree of divergence at J_k,1 also increases with the order, as the denominators from the consecutive sums become multiplied, and it can also increase because of the derivative in Eq. (<ref>), leading to a divergence ∼ |J-k π/2|^-r at the r^th order of expansion. These divergences restrict the convergence radius of the expansion. 
A finite convergence radius would imply no heating in the domain of convergence and would suggest the existence of a heating transition in the parameter J. However, at order r, in addition to the r^th order divergence at J=kπ/2, additional lower order divergencies may appear at J_k,m=kπ/2m, m=1… r. Consider e.g. the replica sum _11=2J ( 4J-i) appearing in the second order expansion, which diverges at J_k,2=k π/4. Many of these possible divergences do not enter the expansion because of the vanishing commutators in the operator part or due to cancellations, e.g. [_1,_1]=0, [_1,[_1,_0]]=0, etc. For instance, the divergence at k π/4 only appears at the 5^th order of the expansion, see Figure <ref>. In spite of the cancellations we conjecture that new divergences keep appearing in increasing orders similar to the dual case with interaction kicks (Figure <ref>(b)), and the expansion blows up near every rational fraction of π/2 (similar to KAM series).

In spite of the divergences at high orders of the expansion, it can provide a very accurate estimate of the Floquet Hamiltonian at low orders. Moreover, the divergences determine the order at which to stop the expansion. That is, given a fixed J, similar to the method introduced in Ref. <cit.>, one can introduce an optimal order of expansion n^*, up to which the corrections increase the accuracy of the approximation of the Floquet Hamiltonian. 

It is natural to assume that the width of the resonances is proportional to the small parameter h, which is further supported by the analysis of the magnitude of the corrections Γ_r
[We give a more rigorous discussion and estimate about the magnitudes of Γ_r in a subsequent publication.]. As an illustration, we give the scaling of the Hilbert-Schmidt norm Γ_r_ HS = √( Γ^†_r Γ_r)∼ f^r_J(J-π/4) near the resonance π/4,

    f^r_π/4(δ J)   =c_π/4(r)/δ J^r-4+𝒪(δ J^-(r-5)) .

The r dependence of the prefactor is illustrated in <cit.>. Up to the highest order we had access to, we found c_J(r) to decrease with r. For our purposes it is enough to assume that it grows at most exponentially ≲α^r, and we expect that at high orders this exponential growth indeed appears as the asymptote of c_J(r). Then the series ∑_r Γ_rh^r diverges for δ J< α h, which gives the width of the resonances.
The optimal order of the expansion is hence estimated by the maximal order at which the closest resonance is located further than ∼α h.
As the resonances appear at the rational fractions of π/2, J=kπ/2m, where m=1,…, n at the n^th order of the expansion, the question is how far one can get in the expansion without having a resonance approaching a fixed J.

Rational approximation of irrational numbers has been thoroughly studied in the mathematical literature <cit.>, and is the cornerstone of the KAM theorem in classical dynamical systems, where the stability of the 
(quasi)periodic motion to integrability-breaking perturbations depends on the irrationality of the corresponding frequencies.
The irrationality of a number is defined by how difficult it is to approximate by rational numbers. The number x is of type (K,ν) if it satisfies
|x-p/q|>K q^-ν
for all integer pairs (p,q) <cit.>.  For example, the most irrational number in this sense is the golden ratio, which is of type (1/√(5),2). Such badly approximable numbers are generic in the sense that for any ν>2, almost all irrational numbers x are of type (K,ν) for some K <cit.>. In the following we choose a J for which 2J/π is of type (K,ν), such that

    |J-kπ/2 m|=π/2|2J/π-k/m|>π/2K/m^ν .

Hence J is not affected by any resonances as long as n<n^*

    n^*=(π K/2α h)^1/ν ,

which we set as the optimal order of expansion. By the construction of the expansion,

    U-e^-i H_F^(n)∼ h^n+1 ,

which gives

    U-e^-i H_F^(n^*)∼ h^n^*+1∼ h^C/h^1/ν≲ e^-C'/h^1/2-ϵ

at the optimal order with some constants C, C' and arbitrary ϵ>0, by choosing ν close enough to 2. Consequently, the Floquet Hamiltonian in the optimal order is conserved for stretched exponentially long times in the inverse kick strength, and, if the steady state is the infinite temperature ensemble, it is approached at least stretched exponentially slowly.
In the above analysis we gave an estimate for the accuracy of the replica expansion. We leave a more rigorous mathematical analysis, similar to the ones in Refs. <cit.>, to future work.


Conclusion.—
We have developed a novel expansion applicable to periodically driven systems where the driving consists of sudden quenches between different Hamiltonians. The expansion takes into account all orders in one of the Hamiltonians and is perturbative in the other. As such, it is an infinite resummation of the BCH formula, whose coefficients can be reproduced by taking the derivatives of the terms in the replica expansion <cit.>. We demonstrated that, similar to the high frequency expansions, the replica expansion is  asymptotic for systems with unbounded Hamilton operators, that is, it may not converge, but performs very well when evaluated at an optimal order. The expansion suffers from resonances near rational frequencies, whose avoidance determines the optimal order of expansion. It is an interesting question whether these resonances have a physical meaning or   
are just an artifact of the expansion, and whether one could remove the resonances by a proper choice of analytical continuation in the replica trick.


This research has been supported by the grants P1-0044 and N1-0025 of Slovenian Research Agency (ARRS), Hungarian-Slovenian (OTKA/ARRS) bilateral  grant N1-0055, ERC AdG grant OMNES, and grants by the Hungarian National Research, Development and Innovation Office - NKFIH K119442, SNN118028. A.P. was supported by NSF DMR-1506340, ARO W911NF1410540 and AFOSR FA9550-16-1-0334.	



apsrev


 


§ SUPPLEMENTARY MATERIAL FOR "REPLICA RESUMMATION OF THE BAKER-CAMPBELL-HAUSDORFF SERIES"








§ EVALUATION OF THE REPLICA LIMIT IN THE FIRST ORDER CORRECTION

The first order correction in the replica trick from Eq. (<ref>) is simply written as

    Γ_1=∑_0≤ m <ρm=∑_0≤ m <ρ (U_0^ad)^m H_1=lim_ρ→ 01/ρ(U_0^ad)^ρ-1/U_0^ad-1=log(U_0^ad)/U_0^ad-1

where we expressed m by the adjoint action of U_0 as m=U_0^-mH_1 U_0^m=(U_0^ad)^m H_1 and the substitution of U_0^ad=e^i Jad_H_0 brings us to Eq. (<ref>) of the main text. 


§ NESTED COMMUTATOR EXPRESSION OF REPLICA EXPANSION

A direct evaluation of the series expansion of the integer powers of the time evolution operator leads to Eq. (<ref>) of the main text. However, similar to the BCH expansion, it is expected that the corrections can be written in terms of nested commutators, which produce local operators if the original Hamiltonians are local. We denote 
the sums in Eq. (<ref>) by η_r before evaluating the replica limit, that is, Γ_r=(-i)^r-1/r!η_r, and

    η_r=∑_0≤ m_1≤… m_r<ρm_rm_r-1…m_1 c_m_1 … m_r .

The corresponding nested commutator expressions from Eq. (<ref>) are denoted by η̃_r,

    η̃_r=∑_0≤ m_1≤… m_r<ρ [m_r,… [m_2,m_1]]] c_m_2 … m_r .

We show that η_r can be expressed as a sum of η̃_̃r̃ and terms which vanish in the replica limit, that is, η_r=η̃_r. 
A necessary condition for the replica trick to work up to the r^th order is to have an analytical continuation for which η_r=(ρ). This means that polynomials of at least second order of {η_1,…,η_r-1} vanish in the replica limit. The second order correction was discussed in the main text, in this notation 

    η̃_2   =∑_0 ≤ m_1 ≤  m_2<ρ [m_2,m_1]
    η_1^2    =∑_0 ≤ m_1 ≤  m_2<ρm_2m_1+∑_0 ≤ m_1 <  m_2<ρm_1m_2
    η_2    =∑_0 ≤ m_1 ≤  m_2<ρm_2m_1(2-δ_m_1,m_2)= η̃_2 +η_1^2

Below we list the combinations of the operators η_s which produce the nested commutators in Eq. (<ref>) of the main text, up to the 5^th order (note the non-commutativity of the different η_ss). 

    η̃_1=   η_1
    η̃_2=   η_2-η_1^2
    η̃_3=   η_3-2 η_2 η_1-η_1 η_2+2η_1^3
    η̃_4=   η_4-3η_3 η_1-η_1 η_3-3 η_2^2+6η_2 η_1^2+3 η_1 η_2 η_1+3 η_1^2η_2-6η_1^4
    η̃_5=   η_5-4η_4 η_1-η_1 η_4-6η_3 η_2-4η_2 η_3+12η_3 η_1^2 +4η_1 η_3 η_1+4η_1^2 η_3+12η_2^2 η_1+12η_2 η_1 η_2 
        +6 η_1 η_2^2 -24η_2 η_1^3-12η_1 η_2 η_1^2 -12 η_1^2 η_2 η_1-12η_1^3 η_2+24η_1^5

We note that if the operators η_s commute then these equations reduce to the cumulant expansion.



§ EXTRACTION OF THE BCH COEFFICIENTS FROM THE REPLICA EXPANSION

It is easy to check that series expansion, or equivalently, the partial derivatives of the corrections Γ_r at J=0 reproduce the commutators from the BCH formula containing r instances of H_1, e.g. 

    ∂_J Γ_1|_J=0   = ∑_j σ_j^yσ_j+1^z+σ_j^zσ_j+1^y=-1/2i [H_0,H_1]
    1/2∂_J^2 Γ_1|_J=0   = -2/3∑_j σ_j^x+σ_j^zσ_j+1^xσ_j+2^z=-1/12[H_0,[H_0,H_1]]
    ∂_J Γ_2|_J=0   = 2/3∑_j  (σ_j^x σ_j+1^z+σ_j^z σ_j+1^x)+   2   (σ_j^y σ_j+1^y-σ_j^z σ_j+1^z)       = -1/12[H_1,[H_1,H_0]]  
    1/2∂_J^2 Γ_2|_J=0   = 2/3∑_j  (σ_j^y +σ_j^zσ_j+1^y σ_j+2^z) -     (σ_j^z σ_j+1^x σ_j+2^y+σ_j^y σ_j+1^x σ_j+2^z)       = -i1/24[H_1,[H_0,[H_0,H_1]]]




§ EXPECTATION VALUE OF LOCAL OPERATORS

An important measure of the accuracy of the replica expansion is how well it reproduces the dynamics of local observables. This can be quantified by the operator norm of the difference between the time evolved observable with respect to the exact and approximate time evolution: Δ̃_n(A)=U^adA-U^ad_(n)A, for a local operator A and notation U^adA ≡ U^-1A U. This measure is independent of the system size up to an exponentially small correction due to Lieb-Robinson bounds. In principle it could be possible that the resonances seen in Figure <ref> correspond to some non-local observables and they do not affect the dynamics of local operators. In contrast, we find that the accuracy measured by the time evolution of single spin operators σ_i^μ, μ∈{x,y,z} follows closely the estimate given by Δ_n, as shown in Figure <ref> for σ_i^x. The other spin components behave in qualitatively the same way.  
As discussed in the main text, the spectral norm of the difference of the exact and approximate time evolution operators Δ_n=U-U^(n) bounds the accuracy of the expectation value of any dynamical observable A, including non-local ones, as Δ̃_n(A)≤ (2 Δ_n +Δ_n^2)A. However, to account for the possibly non-local observables, Δ_n scales linearly with system size, and underestimates the accuracy of the replica expansion for the time evolution of local operators. 






§ HILBERT-SCHMIDT NORM OF THE CORRECTIONS

The Hilbert-Schmidt norm of corrections Γ_r are illustrated in Figure <ref> for both magnetic field and interaction kicks. Note that in the latter case the expansion parameter is J rather than h.

The asymptotic behavior of the Hilbert-Schmidt norm near divergencies can be further analyzed, as in Eq. (<ref>) of the main text,

    f^r_J^*(δ J)   =c_J^*(r)/δ J^r-r_J^*+𝒪(δ J^-(r-r_J^*)+1)
    
    f^r_h^*(δ h)   =c_h^*(r)/δ h^r-r_h^*+𝒪(δ h^-(r-r_h^*)+1) ,

where r_J^*+1 (r_h^*+1) determines the order of expansion at which the divergence at coupling J=J^* (h=h^*) first appears. The coefficient c_π/4(r) is plotted in Figure <ref> (a) as a function of  for kicking magnetic field. It is more informative to analyze the model with interaction kicks, because we have much more resonances at hand in this situation. The coefficients c_h^*(r) are shown for that model in Figure <ref> (b). The new resonances appearing at higher orders are characterized by smaller coefficients than the previous ones, at least up to the 6^th order.





§ PERIODICITY PRESERVING ANALYTICAL CONTINUATION

In the manuscript we used the naive analytical continuation, that is, we evaluated the replica sums in Eq. (<ref>) with the prescription Eq. (<ref>), then we continued the resulting expression to arbitrary real values of the replica index ρ. This results in replica coefficients continuous in J except for the resonances. The denominators in Eq. (<ref>) do not break the 2π periodicity of the time evolution operator U in variable J. The linear in J factors originate in the replica limit of the expressions e^i 4 J z ρ-1,

    e^i 4 J z ρ-1= lim_ρ→ 0e^i 4 J z ρ-1/ρ=log e^i 4 J z ρ=i m_z(4 J z)

where z is an integer. In the naive analytical continuation, for all values of z, we chose m_z as the identity function. However, we can choose a different branch cut of the logarithm, which is equivalent to dividing our analytical continued function by e^i 2πζρ, which is unity at integer values of ρ, and ζ is an integer defining the new analytical continuation. In principle we are allowed to define a different ζ for every single term in the expansion, here we consider it as a function of z. The periodicity H_F(J)=H_F(J+2π) can be enforced by substituting J by J  mod  2π in the replica coefficients. This corresponds to the analytical continuation defined by m_z(4 J z)=4 z (J  mod  2π).

From the point of view of observables, the periodicity in J is shorter, since U only obtains a real phase factor (± 1) following a shift of π:  U(J+π)=(-1)^b U(J), where b is the number of nearest neighbor bonds in the system. In open boundary conditions b=L-1 and in periodic boundary conditions (PBC) b=L. The period is further decreased in PBC, because then U(J+π/2)=i^L U(J) since

    e^i π/2 σ_i^zσ_i+1^z=i σ_i^zσ_i+1^z

and ∏σ_i^zσ_i+1^z is the identity operator if PBC is applied. Thus, at appropriate system sizes, the period of U in J is halved or even quartered, which can be similarly established by a corresponding analytical continuation, that is, by the proper choice of the function m_z.  


