A flow equation approach to periodically driven quantum systems

We present a theoretical method to generate a highly accurate {\em time-independent} Hamiltonian governing the finite-time behavior of a time-periodic system. The method exploits infinitesimal unitary transformation steps, from which renormalization group-like flow equations are derived to produce the effective Hamiltonian. Our tractable method has a range of validity reaching into frequency regimes that are usually inaccessible via high frequency $\omega$ expansions in the parameter $h/\omega$, where $h$ is the upper limit for the strength of local interactions. We demonstrate our approach on both interacting and non-interacting many-body Hamiltonians where it offers an improvement over the more well-known Magnus expansion and other high frequency expansions. For the interacting models, we compare our approximate results to those found via exact diagonalization. While the approximation generally performs better globally than other high frequency approximations, the improvement is especially pronounced in the regime of lower frequencies and strong external driving. This regime is of special interest because of its proximity to the resonant regime where the effect of a periodic drive is the most dramatic. Our results open a new route towards identifying novel non-equilibrium regimes and behaviors in driven quantum many-particle systems.


I. INTRODUCTION
Recent years have seen rapid progress in our understanding of dynamics and non-equilibrium phenomena in quantum systems [1,2]. This has been a result of experimental advances in the ability to control cold atom [2][3][4] and condensed matter systems [5][6][7], by developments in time-resolved laser techniques [8,9], and by the fact that stepping into the time domain opens up new ways of ultrafast control of material properties [5,10,11] and access to different phases of matter. These include photoinduced superconductivity [12,13], hidden orders [14], and metastable states [15], but also entirely novel phases, such as time crystals [16,17] and non-equilibrium topological phases [18,19].
In particular, there has been growing interest in periodically driven (or Floquet) [20,21] many-body systems, which can bear a close resemblance to equilibrium systems [22]. The Floquet systems come in three established thermodynamic classes: integrable [23][24][25], manybody localized (MBL) [19,26,27], and generic interacting ones [28]. The first two classes can avoid thermalization, allowing for a notion of a Floquet phase of matter at long stroboscopic times t = nT , where T is the period of the Hamiltonian, H(t + T ) = H(t), and n is an integer. The physics of these phases is captured by an effective, time-independent Floquet Hamiltonian H F , given via the time evolution operator over one period * These authors contributed equally to this work.
The existence of a prethermal regime is important because realistic systems usually contain integrabilitybreaking perturbations that support it, and because the thermalization (or more specifically, the energy absorption) time τ can correspond to experimentally accessible time scales. The existence of such a regime also implies that there is interesting physics to be found at intermediate times 0 < t < τ [52,59], where one may use timedependent perturbations to drive dynamical phase transitions [60][61][62][63], control interactions [64,65], or engineer phase transitions and topological phases [66][67][68][69][70].
To understand the properties of a system in the prethermal regime, it is convenient to use a description in terms of the effective Hamiltonian, H F . It is, however, notoriously difficult to calculate H F or the exact time-evolution operator U (t) for interacting systems, so arXiv:1808.01697v4 [cond-mat.stat-mech] 19 Mar 2019 generally one uses an expansion technique to find an approximate, effective Hamiltonian in the high-frequency limit. These include the Magnus expansion [71][72][73], rotating frames [53], and many more [20,67,[74][75][76][77][78][79][80]. Unfortunately, these methods do not produce a cleanly convergent expansion series for general systems. Instead, they are asymptotic expansions, subject to an optimal cut-off order which prevents them (in principle) from reaching into the lower frequency regimes [53,71]. By this statement we do not mean to imply that methods such as ours may not be subject to their own cut-offs but that these cut-offs may differ [81]. Whether this is the case for the exact version of the flow equations in this paper is a matter that still has to be determined.
One of the more controlled descriptions of a system occurs in the quasiequilibrium regime, W ω ∆, where W is the bandwidth of the system, ω the driving frequency ( Planck's constant), and ∆ is the gap to the continuum of higher energy states. While this separation of energy scales is quite feasible in cold atom systems, it is harder to reach in solid state systems. Mott insulators are the most promising class of systems in this regard, but even there the range of frequencies is limited since we typically have W ∼ 1eV and ∆ ∼ 1eV, which are of the same order of magnitude. In addition, lower frequency regimes are required for certain topological phases [18], and are of interest in cold atom systems [82,83], and in the study of thermalization [84]. Hence, techniques to handle lower frequencies are needed.
In this paper we improve on the limitations of previous methods, and provide better access to lower frequencies and higher driving strengths. To achieve this we introduce a formalism to remove the time-dependent part of a Hamiltonian using infinitesimal unitary transformations. This results in flow equations for different couplings, reminiscent of renormalization group calculations [85] and Wegner's flow-equation approach to diagonalizing Hamiltonians [86,87]. There has also recently been progress in using the Wegner flow to describe the time-evolution of a many-body localized system [88], which however still requires the solution of flow equations for each time-step -a problem we avoid in our construction.
We note that while a flow equation method for finding effective Floquet Hamiltonians exists in the literature, it uses an approximate version of the Wegner generator (keeping only terms proportional to ω in the generator) [89] in Sambe space [90], where the approximation brings up a question as to the range of validity. Our method differs in that we do not need to introduce Sambe space, and our generator is obtained in a constructive manner and differs completely from the Wegner generator. For our method, we describe both the exact flow equations, and ways to approximate them. We apply our method to the Schwinger-Rabi model of a single spin in a magnetic field, and also to four different spin chain Hamiltonians: (i) an integrable XY model with antisymmetric exchange, (ii, iii) two integrability breaking extensions of a J 1 -J 2 -type XXZ model [91], and (iv) the transverse field Ising model. The extended XY model is driven by a transverse magnetic field, the first J 1 -J 2 -type XXZ model is driven locally by a magnetic field in the x-direction, and the second by a nearest neighbor Ising exchange interaction, making for a time-dependent J 1 -J 2 model [92,93]. For the transverse field Ising model we consider (i) a harmonic driving case, and (ii) a case where the time evolution operator factorizes into two matrix exponentials, which allows us to find a family of different resummations of the Baker-Campbell-Hausdorff (BCH) identity. This observation leaves open the question of how to construct the optimal effective Hamiltonian for a given time evolution operator (reverse of the usual situation in which one seeks the optimal time evolution operator approximation for a given Hamiltonian).
In this paper, we study the time evolution of the exact models, and their effective models obtained in our approach. We compare our results with those obtained by the Magnus expansion. The integrability breaking models are studied numerically using full exact diagonalization, which provides an unbiased test of the validity of our approach. We find that our flow method generally outperforms the Magnus expansion, with significantly greater accuracy as the resonant regime is approached, as well as in the case when the time-dependent term in the Hamiltonian is large. Both of these cases are of direct physical relevance and interest. Our method thus opens new possibilities in the analytical and numerical simulation of time-dependent quantum many-particle systems, and will facilitate the search for novel prethermal, and non-equilibrium regimes.
Our paper is organized as follows. In Sec.II we develop the general flow-equation formalism, and discuss its structure and approximations. In Sec.III we relate the general results obtained from the flow equation approach to various high-frequency expansions used in the literature. In Sec. IV we test the flow equations on an exactly solvable two level system and discuss in detail the properties of the fixed points of the flow equations and their stability. This discussion is continued in Sec. V for a many body-system studied via a truncated ansatz where we show it outperforms a high frequency Magnus expansion and rotating wave approximation. In Sec.VI we introduce four different one-dimensional spin chain Hamiltonians we will use to assess the performance of the approximate method described in Sec.III. In Sec.VII we summarize our results for the different models. In Sec.VIII we compare our results to a resummation of the Baker-Campbell-Hausdorff identity that was of recent interest [81]. We also show what advantages our approach has over a standard rotating frame approximation-namely that it can be truncated when a rotating frame transformation is not practically possible and that it still performs well under these circumstances. In Sec.IX we present our main conclusions. Various technical details and formulas appear in the appendices.

II. GENERAL FORMALISM
We take the Schrödinger equation of a periodically driven many-particle system as our starting point. Following Ref. [53], the Hamiltonian H(t) is split into a constant part H 0 = 1 T T 0 dtH(t), and a time-periodic term ) that averages to zero over one period, 1 T T 0 dtV (t) = 0. Thus, the time-dependent Schrödinger equation takes the form, where we have set Planck's reduced constant = 1.
We introduce a unitary transformation, U = e δΩ(t) , generated by an as yet undetermined quantity δΩ that will be chosen to reduce the time-dependent term V (t). The δ in front of the Ω indicates we keep the generator infinitesimal, which ensures that the exponential can be safely expanded to lowest order.
Let us now introduce a new wavefunction |φ δs = U † |ψ 0 = [1 − δΩ(t)] |ψ 0 and act with U (t) † = 1 − δΩ(t) (to leading order in δΩ) from the left on the Schrödinger equation. This new wavefunction now fulfills the modified Schrödinger equation (keeping lowest order in δΩ only), One may read off a new Hamiltonian, which, since δΩ is infinitesimal, can be written as Up to this point, this treatment coincides with the use of time-dependent generators [94]. We now, however, choose δΩ very different from the Wegner generator. We choose it such that it reduces the time dependent part of the Hamiltonian V (t) → (1 − δs)V (t) by some infinitesimal value δs, where the generator in Eq.(4) also has the nice property that it vanishes at stroboscopic times T . Therefore, at stroboscopic times, expectation values Ô of operatorsÔ can be calculated without a change of basis. The behavior at other times can be found by applying the unitary transformation to the operatorÔ. One could now repeat the procedure of splitting the Hamiltonian into a constant and a time average zero part and then apply this infinitesimal unitary transformation to find the Floquet Hamiltonian after an infinite amount of steps (or an approximation to it by stopping after a finite amount of steps). To simplify the process, we recognize that one can track the progress of the unitary transformations by a single flow parameter, s. To do so we extend the functional dependencies of the Hamiltonian to include this parameter, replacing H(t) → H(s, t) andH (t) → H(s+δs, t). Note that H(s, t) represents a family of effective Hamiltonians interpolating between a starting Hamiltonian H(0, t), and a Hamiltonian H(∞, t). H(∞, t) is the Floquet Hamiltonian H F if V (∞, t) = 0. It seems plausible that V (∞, t) = 0 and we find this to be true in an explicit example and some limiting cases but it remains to be shown rigorously. We set appropriate boundary conditions by enforcing that s = 0 corresponds to the initial, non-transformed Hamiltonian.
With this notation, Eq. (3) takes the form with V (s, t) = 1 ). One may note that this leaves a residual time-dependence of δs [V (t), H(t)] in Eq. (3), which is small in magnitude if δs is small.
Taylor expanding the left hand side since δs is infinitesimal we find, which is a central result of this work. We refer to Eq.(6) as the exact flow equation. This equation is similar in spirit to the infinitesimal unitary transforms that Wegner [86] employs to diagonalize an interacting Hamiltonian in the equilibrium case. One can readily see that Eq.(6) has a fixed point with the desired property V (s, t) = 0. This fixed point is guaranteed to be stable for sufficiently large ω because in this case the commutator term can be neglected. Under these circumstances the time independent parts of H(s, t) remain unchanged. More precisely equation (6) then reduces to which trivially has the stable fixed point V (s, t) = 0 since V (s, t) for this case can be treated like a scalar. But what about smaller ω? Because an analytic understanding is difficult to achieve, we will discuss this in the context of an explicit example in Sec.IV. While our discussion gives a mechanism by which the fixed points can be stable in general it does not give a rigorous proof.
How should one interpret the flow of s in Eq.(6)? Note that H(s, t) is a Hamiltonian and therefore a linear sum of the various energy contributions, and can be expressed as a sum of linear operators with coefficients c i (s, t), H(s, t) = i c i (s, t)Ô i (similar in spirit to a Landau-Ginzburg energy functional). TheÔ i operators are nothing other than kinetic and potential energy terms appearing in a Hamiltonian, such as a hopping term c † i c j in a lattice model, an interaction term n i↑ n i↓ on a lattice, or a multiple-spin term ( S i · S j )( S k · S k ) in a spin model, among many other possibilities. The coefficients c i (s, t) describe the coupling constants (strength) of these terms. This mathematical structure of H(s, t) = Here g i has a functional dependence on the c j (s, t ) with t ∈ [0, T ], because V (s, t) itself depends on the c j (s, t) and it appears under an integral.
One may therefore write Eq. (6) as which is just a flow equation for the coupling parameters c i (s, t) at different times. Note that the set of operatorŝ O i may include both the original operators, and ones generated from the kinetic and potential energy terms of the original Hamiltonian, Eq.(1), as the Hamiltonian flows. In general, new terms are generated such as hopping and interaction terms that involve more and more sites of a lattice as the order of the transformation increases. These new terms can in principle change the balance of kinetic and potential energy in the effective time-independent Hamiltonian and therefore may lead to new physical regimes for a periodically driven manyparticle quantum system. The reason we write the flow equations in this form is to emphasize that Eq.(6) actually describes couplings that flow as we reduce out the time dependence and to show how this operator equation corresponds to a numerically tractable scheme to determine couplings.

III. APPROXIMATIONS TO THE FLOW EQUATIONS
It is important to note that Eq.(6) offers a convenient starting point to approximate the Floquet Hamiltonian. In particular, it allows us to improve on the various high frequency expansions of the Floquet Hamiltonian that have appeared in the literature. As an example, we can find an analytically tractable equation if we set s = 0 only for the terms V (s, t). This corresponds to removing the original time dependent part V (t) from the Hamiltonian via the rotating frame transformation [95] e −i t 0 dtV (t) , while generating other new time dependences. (This approximation is for convenience. Indeed in the following section we will present an example in which we exactly solve Eq. (6) without taking s = 0 in V (s, t).) To ensure that this approximation actually corresponds to the aforementioned unitary transformation we also need to restrict the range of s to [0, 1], rather than the previous [0, ∞).
Let us justify this approximation slightly more carefully by using an analogy. One may notice that Eq. (6) is very similar in structure to the classical problem of a first order differential equation, where g(t, f (t)) would correspond to −V (s, t) in our case and all the couplings in H and V correspond to f (t).
A standard method of solving this class of problems [101] is plugging in the initial condition f (t) → f 0 = f (t = 0) on the right hand side. Integrating both sides of the equation one finds a first approximation to f (t), which we call f 1 (t). One may then repeat the procedure and plug successive approximations f n (t) into the right-hand side. This procedure is called Picard iteration. In our case, it is the same as replacing A variant of Picard iteration that quite often works better is to only set f (t) = f n (t) in some places of g(t, f (t)) but keep it as f (t) in others. This is a particularly helpful improvement when this is done in such a way that some symmetries are explicitly kept that would otherwise be destroyed [71]. For our case, if we only replace the first two V (s, t) → V (0, t) but keep H(s, t) then we find approximate flow equations, where s is set to run from zero to one only. As required above, this still implements a unitary transformation, which can be seen explicitly by reconstructing Eq. (10) from a unitary transformation. Introducing the dependence on the flow parameter s, Eq. (3) reads (11) One may plug in the manifestly anti-hermitian generator δΩ(s, t) ≡ δΩ(0, t) = −iδs t 0 dtV (0, t), corresponding to a unitary transformation U . The result is Eq.(10). Therefore, making such an approximation is a particularly convenient improvement on a Picard iteration.
One may ask why s should run from zero to one as claimed above. One reason for this is that in the lowest order improved Picard iteration we neglect terms that are proportional to s. Neglecting such terms is only justified if s ≤ 1. Therefore, we let the flow parameter run from zero to 1. If we reach a fixed point in this range of values or come close to it, then it is a good approximation. Letting s run to higher values would not be justified and may yield a bad result. Another reason we apply this approach is that we know that for infinite frequencies one reaches a stable fixed point for s = 1. This can be seen easily because Eq.(10) is then approximately given as This procedure also works well in other cases because often at s = 1 one may be close to an unstable fixed point (see, for example Fig.1). We should also mention that the multitude of different possible fixed points (all V (s, t) = 0) and their corresponding s value makes it difficult to estimate the size of the error from letting s only run from zero to one. After all, often s = 1 is close to a fixed point but there may be more fixed points further out (for larger values of s). We will see this explicitly in the next section where we work with the exact flow equations. Now let us return to discussing Eq.(10). One finds that this also can be rewritten in terms of coupling constants as, where one can write g i (t, c j (s, t)) = j γ ij (t)c j (s, t) as a linear combination of couplings c j (s, t). We are therefore left with a first order linear differential equation that doesn't couple coefficients c j at different times. Gone is the more complicated structure of a functional in the c i . The effective time-independent Hamiltonian is then given by where we have taken an average over one period, which is physically meaningful if one is only looking at stroboscopic times. If one is interested in micromotions, one could in principle retain the time-dependence of c i (1, t)-the important "flow" having been taken into account in the parameter s, which has now been set to unity.
The approximation in Eq.(10), setting s = 0 in V (s, t), does not make any implicit assumptions, such as V (t) is small. By contrast, many other high frequency approximations do make the assumption of smallness. As a result, our approach like the rotating frame approximation, works especially well in the limit of strong V (t). We will demonstrate this explicitly in later sections of this work.
It is important to pause for a moment and stress the advantages our approximate method, Eq. (10), offers over a rotating frame approximation, if the latter is carried out exactly. Firstly, if the driving is complicated it is often not possible to calculate the matrix exponential needed for a rotating frame approximation or the rotation induced on operators by such a matrix exponential because it will generate infinitely many components of the operator algebra. This is indeed the case with one of our example models namely the square-wave driven Ising model we discuss later. In this case our method allows one to keep all orders in 1/ω with a truncated ansatz for the Hamiltonian. That this method performs well can be seen in the plot shown in Sec.VIII.
It is also important to recognize that, even if a rotating frame approximation can be done exactly, usually most terms in the Hamiltonian become time-dependent. In most cases this makes a second rotating frame approximation not possible. Our method allows one to avoid this issue by truncating the ansatz Hamiltonian. Lastly, in some cases one would like to prevent the generation of any new terms and see what happens to the coupling constants of a restricted set of terms. Thus, our method provides a convenient starting point for many different approximation schemes.
We would also like to stress that Eq. (10) implements a unitary transformation exactly. Its solution therefore still retains the full information of the original Hamiltonian. In this paper we will be content with discussing results from the first order iteration only. Again, the formalism we present here lays the groundwork for further development of approximation schemes.
Let us explicitly relate the first order iteration to the more common high frequency approximations. For the moment, neglect , which assumes that all couplings in the Hamiltonian are negligible compared to the driving frequency. This is an approximation common to many of the high frequency approximations. We then find that Inserting this back into Eq. (10) and taking a time average we find which is the lowest order of many common high frequency approximations. Hence, our approximation agrees with other approximations in the high frequency limit.
One should also note that there are other ways to approximately solve the exact flow equation by directly working with Eq.(6) and a truncated ansatz rather than solving Eq.(10). We we will do this in one example in Sec.V and will find that it indeed offers an improvement to the methods above (rotating frame and high frequency expansion) and opens the door to many semi-analytical schemes.
Next, we turn to an application of our method to a number of different Hamiltonians and compare our results with other approaches. We find the method nearly always provides more accurate evolution than other approximations, and in many cases our method works substantially better, particularly as the strong coupling resonant regime is approached. This is also true if we solve Eq.(10) with a truncated ansatz like on one of the cases in Sec.VIII.

IV. FIXED POINT STABILITY AND THE PROPERTIES OF THE EXACT FLOW EQUATIONS
Because it is difficult to discuss the stability of the flow equations in Eq.(10) analytically in full generality, we consider a simple example model where the exact flow equations can be written down explicitly. This will allow us to identify a mechanism that makes the fixed point stable. It is conjectured, but we stress not rigorously proven, that this mechanism will persist even for more complicated systems. In Sec.V this conjecture will be further supported. The current section serves as a means to gain some insight into how the flow equations work. We consider the Schwinger-Rabi model of a spin in a rotating magnetic field, For this model the Floquet Hamiltonian, can be found for all frequencies (see for instance [102]). Let us discuss how the flow equations apply to this model. After repeatedly inserting the form of the original Hamiltonian in our exact flow equations in Eq.(6) (always including newly generated terms) we find that the Hamiltonian H(s, t) takes the form, and the flow equations for the couplings {Z 0 , X 0 , Y S , X C , Z C } are given as, (where the denotes the derivative with respect to s) with initial conditions, As expected from Eq.(6), we find that the fixed point is Y S = X C = Z C = 0, with arbitrary Z 0 and X 0 . This is the only fixed point. For this fixed point we may carry out a stability analysis. That is, we expand Eq. (20) around the fixed point to find linearized equations C (s) = JC(s), where C = {Z 0 , X 0 , Y S , X C , Z C } is a vector of the couplings. The eigenvalues of the corresponding Jacobian J are given as,  17). This corresponds to a low-frequency regime. Note that while the couplings exhibit a non-trivial dependence on s until sufficiently large s, the unitary evolution remains stable down to small frequencies, as seen in the red curve (exact flow) in Fig. 3. The couplings after the range of the plot do not change within the limits of the line thicknesses.
down. If the form of the Hamiltonian at the fixed point reproduces that of Eq. (18) this could indeed be the case, since there Z 0 and X 0 would be finite for arbitrarily small ω. One might thus expect that flow equations would be unable to reach a stable fixed point for low enough frequencies. However, this outcome is avoided. To see how this works, recall that the Floquet Hamiltonian H F is determined only up to some phases by  Fig. 1.
In Fig.1 we see that the couplings made a few approaches to a fixed point V (s, t) = 0, but it wasn't stable. However, the couplings Z 0 and X 0 kept shrinking until a stable fixed point was reached. The matrix logarithm, log(U (T )), has branches with relatively small H F and the couplings continued flowing until a branch with sufficiently small couplings to have a stable fixed point was reached. In the language of the exact flow equations, Eq.(6), there existed a branch of the matrix logarithm log(U (T )) such that H(s, t) became sufficiently small that the commutator  Fig. 1) for B p = 3, B z = 1. Note that in spite of the rapid oscillations for small ω, the resultant unitary evolution remains stable, as seen in the red curve (exact flow) in Fig. 3.
be neglected when compared to V (s, t) and therefore a stable fixed point was reached. We were able to observe this effect in all cases we studied and it is plausible that this could be a general mechanism that leads to stable fixed points in our flow equations. This is illustrated in Fig. 2.
From Figs. 1-2, one may suspect numerical issues. However, this is not the case. Rather, the oscillations stem from the fact that the flow equations do not consistently stay on one branch of the matrix logarithm for H F . Flowing to a stable fixed point means choosing the branch of the matrix logarithm that corresponds to a stable fixed point.
Indeed, if we take the time independent couplings in Fig. 2 to calculate the time evolution operator at stroboscopic times and compare it to the time evolution operator calculated via the standard method of a Trotter expansion we find them to be identical. More specifically we calculate the l 2 distance between two unitary operators, that was normed such that it takes values between zero and one (D dim is the dimension of the Hilbert space), where one corresponds to the maximum distance between two unitary operators and zero to agreement between the two operators. A comparison is shown in Fig. 3. Details of the rotating frame approximation and Magnus expansion are given in appendix A). We find that the exact flow equations-despite the couplings rapidly changingfully agree with the Trotter expansion as they should. The wildly jumping couplings are therefore not a numerical artifact.

V. EXACT FLOW EQUATIONS WITH A TRUNCATED ANSATZ
In this section we discuss how the results from the previous section seem to be quite generic by considering a many-body system. We limit ourselves to a specific We choose this Hamiltonian because: (i) it has a relatively strong external drive, (ii) a time-dependent term that does not commute with itself at different times, and (iii) because the time-dependence is convenient for studying the flow equations. One may find flow equations by making the truncated ansatz, with where a ∈ {0, +, −} and the s dependence of the coupling constants C a was dropped for notational simplicity. We do not discuss the specific form of the flow equations here because they are rather complicated and not insightful. Let us rather first have a look at how some of the couplings behave for this system. Specifically let us first look at one representative coupling as a function of flow parameter s.  Fig.2 for the two level system we solved exactly in the previous section. In particular, the coupling constant nearly approaches zero for the fixed point multiple times before eventually a stable fixed point is reached. This strengthens our interpretation that our method might be kept stable by the mechanism we provided in Sec.IV. To get further evidence of this we plot in Fig.5 one of the couplings as a function of ω and find it again to be consistent with the mechanism we proposed in Sec.IV and illustrated in Fig.2. We stress that this is not a rigorous proof of our understanding of how the flow equations manage to converge, but it is does provide good evidence for the general structure of the convergence.
Let us now discuss these results further. One finds numerically that letting s → ∞ only certain terms survive. Namely, as expected from the fixed point C ± i → 0, one is left with The couplings in the range ω ∈ [8,40] are well approximated by Such fitted couplings allow for a semi-analytic understanding in some cases. One should note that for smaller ω the expressions become much more complicated because of the non-analytic behavior of the couplings as seen in Fig.5.
Let us show below how well our approximation [also using results for smaller ω and not just the expression in Eq. (29)] does when compared to the rotating frame approximation and the Magnus expansion. We do not explicitly give the expressions for the couplings in the Magnus expansion and the rotating frame approximation because they are cumbersome and do not provide much physical insight. Instead, we refer the interested reader Ref. [103]. From Fig.6, one finds that the flow equations (with a truncated ansatz) perform better than both the Magnus expansion and the rotating frame approximation. To stress that the comparison to the rotating frame approximation is a fair one, we note that the operators in Eq.(28) are the same as those appearing within the rotating frame approximation. From this example, one sees that the exact flow equations allowed one to find better coefficients than those afforded by the rotating frame approximation.

VI. EXAMPLE MODELS
To demonstrate the power and validity of the flow equation approach for a wider range of many-body systems we will next consider a selection of quantum spin chain (S = 1 2 ) models. Recall that the spin operators S x,y,z n fulfill the commutation relations, (j, k, l ∈ {x, y, z} and m, n label lattice sites) with the special condition for S = 1 2 that where 1 H is the unit operator in the many-body Hilbert space. Here jkl is the fully antisymmetric tensor and δ mn is the Kronecker delta function. In this section we introduce four different spin models that exhibit different functional dependences of the time-dependent term. The first model (XY spin chain) is integrable, and in particular one-particle reducible. The next two models are integrability-breaking modifications of the XXZ spin chain, and the final model is a transverse field Ising model which will be discussed independently in Sec.VIII. These models possess a range of different symmetries and form of the driving term. They will illustrate the generality and mathematical structure of the flow equation approach.
A. XY spin chain with antisymmetric exchange in a driven magnetic field As a first example model we choose an XY spin-chain with an antisymmetric Dzyaloshinskii-Moriya exchange interaction and a time-periodic magnetic field that both point along the z-axis, where and Here, J x/y is the strength of the exchange interaction in the x/y−direction, D the strength of the antisymmetric exchange, h 0 the static magnetic field strength, and h the strength of the magnetic field driving. This model has the advantage that its instantaneous Hamiltonian can be diagonalized by applying a Jordan-Wigner transformation, followed by a Bogoliubov transformation [91]. Furthermore, it has multiple coefficients, which can be varied to check the validity of our approximation based on the flow equations in a variety of cases. Note that the driving term does not generally commute with the static part of the Hamiltonian.
B. J1-J2-model with a driven magnetic field in the isotropic plane In order to find out if a new approximation scheme is valuable for more realistic interacting systems, it is important to go beyond non-interacting models. To this end, we study the J 1 -J 2 -model [92,93], where and with a time periodic magnetic field in the x-direction h(t) = B· 1; 2nπ < ωt < 2nπ + π −1; 2nπ + π < ωt < 2(n + 1)π ; n ∈ Z, (38) where the time dependence was chosen to simplify the numerical treatment done by exact diagonalization. None of our physical conclusions-nor our flow equation method-rely on this piecewise constant form of the time-dependence. It should be noted that J n is the strength of the n-th neighbor exchange interaction in the isotropic plane and J z n is the exchange interaction in z-direction. For a more compact notation we defined . We chose this model because the external magnetic field breaks magnetization conservation and it therefore also allows us to see if the flow equation approach works under circumstances where the driving breaks a symmetry of the static part of the Hamiltonian. We also applied the flow equation approach to a model in which one of the spin-spin interaction terms is timedependent. The model we consider is another J 1 -J 2 model given by, where and where mod denotes modulo. In this model, the timedependence is in an interaction term. In Sec.VIII, we will consider one further spin model (Ising model) separately because the structure of the fixed point Hamiltonian is different than the three models introduced in this section. Together, these four spin models and the example given in Sec.IV should provide a compelling picture for the generality and power of our method.

VII. RESULTS
In this section we study how well our flow equation approach performs compared to common high frequency approximations. We compare the approximate time evolution operators obtained through various approximations to the exact time evolution operator (obtained by exact diagonalization) at stroboscopic times.
We adhere to the following procedure: We first make use of the translational invariance of our models and calculate the exact time evolution operator U k ex (T ) and the approximate time evolution operator U k approx (T ) at different points in k-space (momentum space). Then, we calculate the mismatch of the approximate time evolution operator and the exact time evolution operator via, which is a quantity that takes values on the interval [0, 1], with zero meaning perfect agreement and one meaning the largest possible disagreement. Here, D dim is the dimensionality of the Hilbert space for any given k-point, N is the number of k-points that the sum runs over, and A Frob := √ trAA † is the Frobenius norm. Let us motivate this quantity: For a given point in kspace this is just the l 2 distance, Eq. (24), between two unitary operators at this point in k-space divided by the maximum l 2 distance of two unitary operators. We average this quantity over all points of k-space. The Frobenius norm provides us with a basis-independent measure of how accurate unitary evolution of a quantum system will be with various time-independent approximations to the full time-dependent Hamiltonian. Similar formulas are used in the context of quantum information science.
where a labels the approximation scheme, with different coupling constants for different approximation schemes. The details of the derivation are given in Appendix B.
There are newly generated terms in Eq. (43) compared to Eq. (33). We note that a suitably chosen rotation in spin space gives back the original undriven Hamiltonian with ∆J = J x − J y modified.
The coupling constants for the leading order Magnus expansion are, and the results for the flow equation approach are, It should be emphasized that both approximations agree in the limit of ω → ∞ -a general result mentioned previously at the end of Sec.II. We also stress that, in the case that h is much larger than all other coefficients, the flow equation approximation works well even when expanded around 2h/ω 1, which is not what one would normally expect from a high frequency expansion. The flow equation approach does not make the assumption at any point that V (t) is small, and therefore it handles this regime more accurately.
We are particularly interested in quantum manyparticle systems with a large number of degrees of freedom. We therefore compute the mismatch E, Eq.(42), of the time evolution operators for a long spin chain. We plot the relative error E as a function of the number of k-points to find out how many k-points are needed for a stable result. (The details on how the time evolution operator was calculated are given in Appendix C.) The plot for the Magnus expansion and for the flow equation approach are given in Fig.7. From Fig. 7 one can see Note that the flow equation error is much smaller than the Magnus expansion error, particularly at the lowest frequencies. In both approximations, the error decreases as the frequency increases. We consider the case of D = 0.1, J x = 1, J y = 1.1, h 0 = 1 and h = 1 that at 256 k-points the value of the relative error E has stabilized. Therefore, for this model all further plots will be done sampling 256 k-points.
To study the accuracy of the different approximations as a function of frequency, we choose a set of coefficients D = 0.1, J x = 1, J y = 1.1, h 0 = 1, and h = 1, where J x was fixed at unity because one may divide the Hamiltonian by J x to make it dimensionless. The strength of D was chosen to be small since often the anti-symmetric exchange is small when compared to the exchange inter-actions. The other values were chosen to to be in a similar range. The plot of the relative error E as a function of frequency ω is given in Fig. 8. From Fig. 8 one can see that the results from the flow equation approach are valid down to much lower frequencies ω. In fact, one can expect higher order Magnus expansions to become worse at lower frequencies than the first order Magnus expansion we plotted. This is because the optimal cut-off order of the Magnus expansion (and a number of other high frequency expansions) shrinks with decreasing frequencies [53] unless couplings are small enough to suppress this effect. It should also be noted that the stuttering (wiggles) at low frequencies seen in the plot is an effect that happens because the U k matrices are relatively small. For larger matrices this averages out as we will see in interacting models to follow.
In Fig. 9 we show how well the approximation does as a function of various couplings. From the plots it is clear that the results obtained via the flow equation approach are generally more accurate than the results from the Magnus expansion. As expected from general arguments, we find that the approximation does increasingly well for large values of driving h. We now turn to non-integrable models.
where (a) labels the approximation scheme (either flow or Magnus). The details of the calculation are given in Appendix B.
It is important to note that one of the new terms, Γ n , can be removed by a suitable rotation in spin space, which tells us that we went from an XXZ model to a XYZ model followed by a rotation in spin space. The effective coefficients for the Magnus expansion are and for the flow equation approach Calculating higher orders in the Magnus expansion for this model yields extremely complicated effective Hamiltonians. The second order Magnus expansion already gives a Hamiltonian that is a sum of 60 different operators with complicated prefactors. One tractable way to improve on the first order Magnus expansion is via the flow equation approach. The plots in Fig. 10 illustrate the quality of the approximation for different frequencies. These results are obtained numerically using exact diagonalization for finite size systems, as described in Appendix D.
One finds that the flow equations outperform the Magnus expansion for all frequencies. For the plot of strong driving magnetic field h this is especially pronounced. There, the Magnus expansion for a large range of frequencies gives poor results and the flow equations generally give quite precise results.
One may also in this case ask how well the approximation does as a function of all the different coefficients. In Fig. 11, we show a plot for different values of the coefficients. These plots were done only including the sector k = 0 in k-space because this is numerically quicker and because other points in k-space reproduce the same results.
Similar to the previous integrable model, for this non- where S ⊥ i = (S x i , S y i , 0) and (a) labels the approximation scheme.
The last two terms of Eq.(51) are newly generated terms in the Hamiltonian. If S z i has an approximately uniform orientation the terms proportional to D (a) n can be interpreted as different range antisymmetric exchange terms -treating S z i as a mean-field term. By the same token, in a mean field approximation the term proportional to Q (a) n can be interpreted as exchange terms. Beyond the mean-field case it is clear that higher order spin interactions are generated. Such terms can lead to new physics and can drive new phases.
The coupling constants within the flow equation approach [solving Eq.(10) exactly] are given by, and within the Magnus expansion, While the form of the Hamiltonian in Eq.(51) is already complicated (with three and four-spin interactions) it is worth noting that the second order Magnus expansion would become forbiddingly complicated with a sum of over 100 operators, which makes even a numerical implementation impractical. Therefore, the result from the flow equations, while also complicated, is a significant improvement on the first order Magnus expansion.
In Fig.12 we plot the frequency dependence of the approximation. One finds that the flow equation result is much better in the lower frequency regime and outperforms the Magnus approximation significantly when the external drive is relatively strong. The performance of the two approximations as a function of the different couplings is shown in the plots in Fig.13. Consistent with the models previously discussed, the flow equation approximation does substantially better across all parameter regimes. For this case we made use of the QuSpin package [96] to obtain a comparison to the exact result.

VIII. COMPARISON WITH RESUMMATIONS OF THE BAKER-CAMPBELL-HAUSDORFF IDENTITY
In this section, we turn the logic around relative to the conventional Hamiltonian-evolution operator relationship. Up to this point in the manuscript, we have been asking about computing an effective time-independent Hamiltonian for a time-dependent problem, and we have used this effective Hamiltonian to compute the time evolution of the system. Now, we turn our attention to a situation in which the time evolution operator is known (in our case it takes a specific product form) and we wish to determine an optimal Hamiltonian that can be used to produce the desired time evolution. This may be useful in certain quantum computing applications, for example.
A second goal of this section is to show that our method has advantages over the rotating frame approximation in that one can capture most of its features by a truncated ansatz even when an exact rotating frame approximation cannot be calculated because the effective Hamiltonian would include infinitely many long range interacting terms. This highlights an another important dimension to our flow equation approach, beyond the examples illustrating its use in earlier sections of the manuscript.
There has been a recent surge of interest in resummations of the Baker-Campbell-Haussdorff (BCH) identity [81]. An important evolution case where the BCH identity is useful is when the time evolution opera-tor factorizes into a product of matrix exponentials e −iH1t e −iH2t . This structure corresponds to multiple different Schrödinger equations. One possible correspondence is to a delta function time dependence in the Schrödinger equation. For example, the kicked transverse field Ising model that is discussed in Ref. [81] has, and can be put into the form, where to stay close to the notation of Ref. [81] we use Pauli operators σ x,y,z i = 2S x,y,z i rather than the spin operators we used earlier in our manuscript. Here δ(t) is the Dirac delta function.
Another possibility is to rewrite the problem in terms of a Heaviside θ function as, (56) Both choices lead to different flow equations and can therefore be interpreted as leading to different resummations of the BCH identity. Thus, we discuss here these two Hamiltonian choices for a given time evolution operator. As a matter of fact, there are infinitely many ways to make a choice in the time dependence, and likely one is an ideal choice. However, we will not discuss this issue of the optimal choice any further. An important difference between the two formulations is that the flow equations in one case can be solved exactly and in the other case require truncation. This allows us to assess how useful our method is in a case where a rotating frame approximation cannot be calculated exactly. This example helps to illustrate the point that even when the flow equations are not solved exactly, they still give results beyond the Magnus expansion.
One finds that within the lowest order in the BCH expansion, the replica approximation used in Ref. [81] and our flow equation approach lead to an approximate Floquet Hamiltonian of the form, where a labels the approximation scheme. The different approximations only differ in their coefficients (and some coefficients may be zero). The coefficients themselves offer little to illuminate our discussion. Therefore, their derivation is given in Appendix E.
In Fig. 14 we show a comparative plot for the δ-type and the Heaviside-type resummations. The plots are done for spin chains of length L = 14 to get a smooth plot. There are only small numerical differences for longer spin chains.
In the plots one can see that the flow equation approach Eq.(10) does better for small values of coupling strength than the Magnus expansion -in some cases also better than the replica expansion. For large couplings, it outperforms both.
From Fig. 14, one can see that the flow equation approach is the most reliable approximation with the mismatch in some cases plateauing at values of around 0.1. For those values one is still able to capture at least qualitative features of the time-evolution. Thus, the flow equation approach offers a useful numerical strategy for finding a Hamiltonian describing a given time-evolution. This may be of practical importance in a wide variety of applications where it is difficult to determine the underlying Hamiltonian from microscopic considerations, such as may be the case in various types of quantum information scenarios.
We would also like to stress that for the step-wise drive the exact rotating frame transformation was not possible to calculate and therefore a truncated ansatz for the Hamiltonian had to be employed to solve Eq.(10). One can see that this truncated ansatz performs well (red curve). It should be stressed that the truncated ansatz performed similar to the case where an exact rotating frame approximation was possible. Our method therefore allows one to capture properties of a rotating frame approximation even when calculating a rotating frame approximation exactly is not possible.

IX. CONCLUSIONS
In conclusion, we have introduced an accurate "flow equation" approach to compute effective timeindependent Hamiltonians, valid for finite times (which may be exponentially long) for periodically driven quantum many-particle systems. We have demonstrated the power of the flow equation approach by illustrating how one can reach into perturbatively inaccessible frequency regimes, and shown that the approximation generally yields an improvement over the Magnus expansion, and that it can also outperform the rotating frame approximation. Furthermore, in many instances the results from the flow equation approach also yield a practically accessible improvement on the first order Magnus expansion where no other method appears to be available. A straightforward application of the Magnus expansion leads to an explosion in the number of different operators that contribute to the effective Hamiltonian with coefficients that are tedious to evaluate. In our approach, one is able to truncate the number of operators contributing to the flow equations in a controlled way, which allows one to keep fewer terms but find highly accurate coefficients. We have also demonstrated that our method compares favorably to resummations of the Baker-Campbell-Hausdorff identity, illustrating it shows its strength even in niche applications, where more powerful methods are to be expected. Our approach also has a wider range of applicability than standard rotating frame approximations because, even if a rotating frame approximation is impractical or not possible because the matrix exponential or the rotation of operators induced by it cannot be calculated, our method allows for a truncated ansatz that may still capture the important features of the transformation.
In summary, we hope that the demonstration of the validity of our approximate method illustrates its power and potential impact on time-dependent quantum manybody systems. The method is completely general and applicable to any form of time-dependent terms in the Hamiltonian-be it through the potential energy, kinetic energy, or both. With the accurate effective timeindependent Hamiltonians that one obtains, new access is granted to potential prethermal regimes with properties not present in the equilibrium phase diagram of the original Hamiltonian. Our results also open the door to new opportunities for quantum control through Hamiltonian engineering to create desired properties out-ofequilibrium. The effective Hamiltonian can be used to compute any observable over finite times through the standard formulas of statistical mechanics, in addition to accurately governing the evolution of the quantum states themselves. We hope our approach will inspire new studies that exploit its flexibility and expand the range of approximation schemes that can be employed within it. With it, new regimes of cold atom, condensed matter, and other systems will likely be uncovered and manipulated in new ways.

Flow equations for the Heaviside θ-function model
The flow equations for the Heaviside θ-function model, Eq.(56), in our approximation, Eq.(10), are found to generate an infinite amount of terms. This means an exact solution of (10) is impossible in turn this also means that a rotating frame approximation is impossible because matrix exponentials and also the rotation of operators induced by it cannot be calculated. Our method allows to truncate terms and therefore find an approximate rotating frame transformation. The terms that appear in Eq.(57) are generated quickly when using an ansatz that starts with the form of the original Hamiltonian, and subsequently adding the new terms that appear to that ansatz. This motivates one to include as many terms from the Hamiltonian Eq.(57) as possible while still allowing for a compact analytical result. We choose the ansatz Hamiltonian, (E4) The flow equations, Eq.(10), give us the following equations for the coefficients dC F,θ x (s, t) ds = −(4J z C F,θ yz (s, t)f I (t) + h x f (t), ) dC F,θ z (s, t) ds = −h z f (t), dC F,θ yy (s, t) ds = 4h x C F,θ yz (s, t)f I (t), dC F,θ xz (s, t) ds = 2h z C F,θ yz (s, t)f I (t), dC F,θ yz (s, t) ds = 2f I (t)(J z C F,θ x (s, t) − h z C F,θ xz (s, t), + J z C F,θ xzz (s, t) − h x Cyy(s, t) + h x Czz(s, t)), dC F,θ zz (s, t) ds = (J z f (t) − 4h x C F,θ yz (s, t)f I (t)), dC F,θ xzz (s, t) ds = −4J z C F,θ yz (s, t)f I (t), (E5) where f (t) = θ t − 1 2 , f I (t) = t + (1 − 2t)θ t − 1 2 and θ the Heaviside function.

Result for the BCH idenity
For the BCH idenity one finds coefficients (E8)

Result for the replica approximation
The coefficients for the replica case were taken from Ref. [81] as