Estimating heating times in periodically driven quantum many-body systems via avoided crossing spectroscopy

Periodic driving of a quantum (or classical) many-body system can alter the systems properties significantly and therefore has emerged as a promising way to engineer exotic quantum phases, such as topological insulators and discrete time crystals. A major limitation in such setups, is that generally interacting, driven systems will heat up over time and lose the desired properties. Understanding the relevant time scales is thus an important topic in the field and so far, there have only been few approaches to determine heating times for a concrete system quantitatively, and in a computationally efficient way. In this article we propose a new approach, based on building the heating rate from microscopic processes, encoded in avoided level crossings of the Floquet propagator. We develop a method able to resolve individual crossings and show how to construct the heating rate based on these. The method is closely related to the Fermi Golden Rule approach for weak drives, but can go beyond it, since it captures non-perturbative effects by construction. This enables our method to be applicable in scenarios such as the heating time of discrete time crystals or frequency dependent couplings, which are very relevant for Floquet engineering, where previously no efficient methods for estimating heating times were available.

All these developments rely on the insight, that a system subject to a periodic drive can be described by an effective Hamiltonian, which may be related to the system Hamiltonian and simply feature renormalized couplings, but may also be completely different. However, there are strong arguments indicating that for generic (ergodic) interacting QMBS the periodic driving leads the system to heat up to infinite temperature resulting in a featureless state at late times [34,35]. The only exceptions known so far are many body localized systems, which are believed to resist heating to infinite temperature [36][37][38][39] and an O(N ) model in the N → ∞ limit [40]. While these arguments are not disputed in principle, over time * artem.rakcheev@psi.ch a number of numerical studies have observed and reported absence of thermalization to infinite temperature in clean systems, which was either attributed to dynamical localization phenomena [41,42] or threshold behavior [43][44][45], seemingly challenging the heating to infinite temperature paradigm.
Another major challenge is the actual determination of the time it takes for a particular system of interest to heat up to infinite temperature. This question is particularly important for the Floquet engineering of e.g. topological phases, where the prethermal regime governed by the effective Hamiltonian should be long enough to observe the transient stabilization of interesting phases, well before the heating dynamics takes over. Understanding the time scales in this setup has thus gained attention in the last years and there has been corresponding theoretical progress. Most significant perhaps are proofs that the heating time is exponentially large, typically t h ∝ exp(ω/J) with some microscopic energy scale J, in the high-frequency regime ω J ( = 1 throughout the paper). The proofs use different analytical techniques such as an analysis of the errors in linear response theory [46], the Magnus expansion [47] or multiple rotating frame transformations [48,49]. These approaches rely on bounds such as Lieb-Robinson bounds [50,51] and can typically not be used to obtain numerical estimates of the actual heating time in a specific system. It is therefore of significant importance to have an accurate, flexible, reliable and computationally efficient method at disposal to predict the heating times in driven quantum many body systems.
So far only relatively few quantitative analyses of heating times (not necessarily focusing specifically on heating to infinite temperature) in Floquet systems have been performed. Most rely on explicit and computationally expensive real-time simulations of either sufficiently large finitesize systems [44,52] or systems treated within truncated schemes such as Density Matrix Truncation [53], methods based on the Density Matrix Renormalization Group [54,55], non-equilibrium Dynamical Mean Field Theory [56,57] or a Keldysh approach [58]. For effectively weakly driven systems, the Fermi Golden Rule (FGR) approach provides an accurate picture [59]. To our best knowledge, the latter one is the only studied method for generic systems, that is (significantly) computationally less expensive than the real-time simulations and applicable to generic systems but is restricted to effectively weak coupling. Part of the reason for the reduced cost is, that it can obtain accurate predictions from smaller systems than real-time simulations. The problem in working with small systems, as discussed in Refs. [34,42,60] and also in later parts of this article, is that the Hilbert space is too small to support heating at large frequencies. This though does however not mean, that the information about heating time scales is not yet contained within small systems, as we will demonstrate here.
In this work, we analyze the appearance of avoided level crossings in the eigenvalues of the Floquet propagator and will be able to infer and quantitatively predict also very long heating times. The significance of avoided crossings in Floquet (and generally many-body) systems is well-known and documented in [61][62][63], but to our best knowledge there have not yet been efforts to resolve individual (as we will see often very narrow) crossings systematically and to link them to the heating rate.
The article is structured as follows: in Sec. II we introduce some basic notions in Floquet theory necessary to follow the arguments in further sections. In Sec. III we review some of the prior work centered around heating in Floquet systems. In particular we explain how heating rates show in real-time simulations and how they are predicted using the FGR. Our method of predicting heating rates, based on avoided crossing spectroscopy, is described in Sec. IV, where we also make the connection to the FGR for weak drives. In Sec. V we discuss in detail the application of our method to a particular driven spin chain [44,45]. For this model we also identify certain commensurate parameter points, akin to discrete time crystals, where the effective Hamiltonian cannot be obtained by a perturbative expansion and show that our method detects these features and is still applicable and accurate. Finally in Sec. VI, we use the driven spin chain with frequency dependent couplings, an important scenario in Floquet engineering, to illustrate the applicability of our method in this case as well.
A different spin chain model [52,53] with weakly broken spin inversion symmetry is discussed in Appendix F as a further illustration of the power of our method.

II. ELEMENTS OF FLOQUET THEORY
Floquet theory is concerned with the study of time-periodic quantum many body systems with Hamiltonian H(t) with period T = 2π/ω. A standard setup consists of an average Hamiltonian H 0 and a drive Hamiltonian V as with a T -periodic function f (t) = f (t + T ) with zero mean.
The propagator over a single period, formally given by the time-ordered exponential (with time-ordering operator T ) has the eigenvalues λ i = exp(−iθ i ). We call θ i the eigenangles in the following. The Floquet Hamiltonian H Fl is the generator of the propagator over one period with eigenvalues θ i /T . It is not unique since the angles can be chosen modulo 2π. A common choice, that we also make, is to restrict the eigenangles to the first Floquet zone θ i ∈ (−π, π]. In this work, for simplicity, we focus on a square wave drive (also known as a switched or bang-bang protocol) which however is more naturally understood in a discrete sense using the product with the half period τ = T /2 and the Hamiltonians H ± = H 0 ± V . Such setups are often used in theoretical studies, since they can be simpler to analyze analytically and numerically. They also arise naturally in digital quantum simulation, for example via a Trotter decomposition of a time independent Hamiltonian, see e.g. [44,45].

A. Floquet Hamiltonian
The Floquet Hamiltonian, as defined in (3), governs the stroboscopic evolution between periods. However, using another starting point within a period (i.e. a different phase of the square wave), would lead to a different Floquet Hamiltonian. For this reason, notions such as Floquet or effective Hamiltonian are not always used to denote the generator of the evolution operator. Some literature rather reserves these names for a gauge invariant formulation, moving influences such as the initial phase to the so called kick operator [1, 64,65]. We will not make use this formalism, but would like inform that the Floquet Hamiltonian as defined in (3) is not gauge invariant [64].
The appearance of the Floquet Hamiltonian, which can have a non-trivial dependence on the average and drive Hamiltonians, is what makes Floquet systems an interesting research subject. Correspondingly, large efforts have been devoted to obtain approximations to the Floquet Hamiltonian, usually at large driving frequencies ω [1, 48,[64][65][66][67][68][69][70][71]. In the digital setup (5) the Baker-Campbell-Hausdorff (BCH) series [72] can be used at small half-period τ (high frequency ω) where we recognize that the average and the Floquet Hamiltonian correspond to each other to first order in T . Figure 1. Floquet diagram of the spin chain with L = 8 with eigenangles of Usw(τ ) colored by their energy density with respect to H0. The location of the first occurrence of an avoided crossing in this diagram is indicated by a black circle. Due to specifics of the protocol the eigenangles concentrate onto two points (one point) at τ = π (at τ = 2π).

B. Floquet Diagram
A visualization tool, used at times in the literature and providing a lot of insight for our method, is the Floquet diagram wherein the eigenangles θ i of U sw (τ ) for a finite size system are plotted as a function of the half-period τ . For a model spin chain, to be specified in Sec. III A, such a diagram is shown in Fig. 1, where we only show the relevant symmetry sector. The lines denote the eigenangles θ i (τ ). The color of the lines highlights the expectation value of the energy density with respect to the average Hamiltonian in the eigenstates of the Floquet propagator: 0 = θ i (τ )|H 0 /L|θ i (τ ) . Let us walk through some of the features of the Floquet diagram, which will be important in developing our method. At small τ , starting from τ = 0, the lines are almost straight as the Floquet propagator eigenstates and the eigenstates of H 0 basically coincide, and thus their slope in the diagram is proportional to the energy given by H 0 . Incidentally, for the protocol family at hand (5), this quantity is proportional to the derivative of the eigenangles with respect to τ at all values of τ [73]. Due to the 2π periodicity of the eigenangles, at a certain value τ =τ c the continuation of the lowest and highest energy state of H 0 seem to cross as indicated by a circle in Fig. 1. A more refined analy-sis shows that the two states undergo an avoided crossing with a very small minimal angular gap of ∆θ(τ c ) ≈ 10 −7 rad.
In a small τ −θ region around this avoided crossing the many-body system can be sketched as an effective two-level system, whose dynamics can be understood as Rabi oscillations [74] (see Appendix A for an example). If we were to sit right on the crossing at τ c and initialize the system in one of the two eigenstates of H 0 , the off-diagonal matrix element responsible for the minimal gap would drive a (resonant) Rabi oscillation between the two H 0 eigenstates, therefore violating the energy conservation with respect to H 0 in an explicit manner. This exemplary first crossing thus provides an initial seed on a small system for the proliferation of many-body heating processes in larger systems. Our proposed method will build on this important intuition and consists of an automated analysis of all finite size level crossings in a certain window of τ and the H 0 energy transfer at each of them.
Looking again at the Floquet diagram, we notice that for our specific choice of H 0 and V , the Floquet propagator U sw (τ ) shows unusual behavior in the considered τ window at τ = π and τ = 2π, where the eigenangles join at θ = 0, π and θ = 0 respectively. This behavior is closely related to discrete time crystals, where at least for τ = π one observes a recurrent dynamics with period 2τ (i.e. a period doubling) for generic initial states. While we relegate the discussion of the specific properties of the Floquet propagators at these values of τ to Sec. V and Appendix D, an important feature is that the number of avoided crossings as well as the magnitude of angular gaps, is strongly suppressed in the vicinity of those special points leading to a reduced heating rate, which competing methods such as the FGR treatment cannot easily access.

A. Driven Spin Chain
As an illustration for heating in QMBS and for our method in later sections, we focus on a model which was recently studied in the context of digital quantum simulation [44,45] and argued to exhibit a threshold behavior as a function of the half period τ , i.e. the absence of detectable heating below a threshold value of τ . The model is defined by with spin one-half operators s α i and periodic boundary conditions 1 , leading to Any operator in the protocol has a spatial translation symmetry and a spatial reflection symmetry, which allows us to reduce the Hilbert space dimension by working in the zero momentum and even spatial parity sector throughout this article. This restriction is only possible for initial states lying fully within the given sector, for example translation invariant products states, which we use throughout the article. H 0 is an Ising model with transverse and longitudinal field, with parameter values not too far from other instances which have been reported to obey "eigenstates thermalization hypothesis" (ETH) properties [75,76]. The average and the drive Hamiltonians can thus be characterized as generic (non-integrable) QMBS.

B. Phenomenology and Earlier Diagnostics of Heating
Suppose that we evolve a pure state with the Floquet propagator U sw (τ ) and monitor 0 after each cycle. Since the system is a Floquet system with only discrete time translation invariance, the average energy 0 need not be conserved. The initial value is given by the expectation value of the average Hamiltonian in the initial (product) state. If the hypothesis of heating to infinite temperature holds true, then 0 is supposed to approach zero at late times for our Hamiltonian (in the thermodynamic limit).
It has been predicted analytically [46][47][48][49]77] and observed in numerical simulations [52,53,59], that for large parts of the dynamics the energy density changes exponentially slowly 0 ∼ exp(−Γt) with the heating rate Γ, or equivalently the heating time t h ≡ 1/Γ, and that this time increases exponentially with the frequency of the drive. For certain small to intermediate system sizes, the heating rates can be obtained from real-time simulations using numerically exact methods. We perform the time evolution using Krylov subspace methods [78] with partial reorthogonalization [79] and appropriate error bounds [80] and show the results for a product state along the x-axis |x, + in Fig. 2 for four different values of the half-period τ . In the figure we can clearly observe the almost perfect exponential decay of | 0 |, until the curves reach a finite size plateau with fluctuating | 0 |, here of the order of 10 −4 for the system size of L = 28 spins. In the thermodynamic limit the energy plateau would be at 0, but due to the finite system size, it has a finite value decreasing with system size. More precisely, the steady state can be described as a random pure state, which can be inferred from computing the Shannon entropy − i P i ln P i with probabilities in the computational basis. In a Hilbert space with large dimension N , this quantity is given by approximately ln(N ) − 0.4228 [81], which we also observe for the steady states with an error of ≈ 10 −3 .
The heating times can be extracted by exponential fits in an appropriate window, which are indicated by black dashed lines. Even without extracting the rates, one can see clearly that the heating time decreases from τ = 1.1 to τ = 1.35 and τ = 2.2 as expected, but then the heating time increases again at τ = 2.3. In Sec. V we will explain this unusual behavior in more detail.
In complementary previous work the heating to infinite For larger τ this trend is reversed due to the specifics of the model. temperature was often diagnosed not from the actual real-time evolution of the energy, but instead through diagnostics which build on properties of the set of all eigenstates of the Floquet propagator U sw (τ ). A prominent example is to determine the level spacing statistics of the eigenvalues of the propagator, or the (inverse) participation ratio of eigenfunctions [34][35][36]44]. These diagnostics build on the idea that for systems which heat up to infinite temperature the propagator is effectively an instance of a random unitary matrix (in the circular unitary (CUE) or circular orthogonal (COE) ensemble ) [34,35]. In our work we demonstrate that these diagnostics are quite conservative, i.e. they are typically unable to detect the heating to infinite temperature on system sizes which are too small in relation to the underlying heating time. If systems are too small, such that no heating is observed in real-time simulations, these measures also cannot be used to learn about heating for larger sizes [42,60]. This is illustrated in Fig. 3, where the dynamics of the energy density of |x, + is shown for various system sizes. As seen in the figure, the smallest shown sizes do not seem to heat at all and one needs to go to L ≈ 24/28 to really see a consistent heating rate. Our work based on avoided level crossings however directly focuses on the seeds of the heating processes, and is able to predict even very large heating times from rather small systems (L ≤ 12), where the circular ensemble has not yet permeated most eigenstates of the propagator. Conversely, we will also see that our method is not well suited to extract heating times in regimes where the ensemble's properties are fully expressed, but since these are the "simple" cases, where the heating happens typically very fast, this is not an important limitation. . Dynamics of the energy density starting from |x, + for various system sizes. At small sizes no heating is detectable, thus real-time simulations to extract the heating rate require large system sizes.

C. Fermi-Golden Rule
As we will see shortly, there are certain parameter regimes wherein the heating rate changes over several orders of magnitude in a small τ /ω window, rendering real-time simulations particularly expensive, since to resolve this region one needs to perform the evolution for multiple parameter values and potentially very long run times not known a priori.
In weakly driven systems, the FGR has been shown [59] to give accurate predictions at a much lower cost. The FGR is rooted in time-dependent perturbation theory/linear response theory [74,82,83] and proposes the following formulȧ for the extensive energy absorption rate (EAR)Ė(ω) as a function of the driving frequency ω. In the formula g m is the m-th Fourier component of the driving amplitude f (t) and the set of |i denotes the eigenstates of H 0 with energy E i . Subsequently, ω ji ≡ E j − E i is the energy difference between eigenstates and P i denotes the probability to find the system in eigenstate |i . The latter is needed since the heating rate is a priori state dependent, however during the evolution these probabilities of course change, which is not captured by the formula.
In a numerical FGR computation one has to calculate all the matrix elements j|V |i of the drive Hamiltonian V in the eigenstates of H 0 , amounting to one full diagonalization of H 0 . Then using specific values of ω, g m and a model for P i one can evaluate Eq. (8). Apart from providing a computation tool, the FGR also provides a way to understand the exponential increase of heating times at high frequency and sheds light on some statistical aspects of heating in many-body systems.
As a starting point one can rewrite the double sum as an integral over the density of states D(E) (details are discussed in the Supp. Mat. of Ref. [59]) Since the approach operates under the assumption that the system (i.e. average Hamiltonian) is generic, it is expected that the matrix elements of the drive are given by the ETH ansatz whereĒ is the average of E, E + mω, R is a random variable with zero mean and unit variance and f V is a smooth function independent of system size [84,85]. For local operators O, f O has been shown numerically to decay exponentially with ω for high frequencies in a variety of systems [85][86][87]. This behavior can serve as an explanation for the exponential suppression of heating and threshold behavior from a statistical perspective. Furthermore, previous results from FGR (and our results from crossing computations) suggest that small systems can already provide good estimates for this function. In evaluating the formula for such systems, one effectively cancels the density of states factors and gets an estimate for the thermodynamic limit, where the behavior is dictated by f V .

IV. HEATING RATES VIA AVOIDED CROSSING SPECTROSCOPY
As we have discussed in the previous section, the success of the FGR is rooted in a sensible separation of microscopic processes embodied in f V , from statistical factors like the density of states. We propose to construct the heating rate in a similar fashion, but using the true microscopic processes in the systems, encoded in avoided crossings, rather than the expression based on linear response theory. In later sections we will show that this method has clear advantages in several scenarios occurring in Floquet systems.

A. Avoided Crossings
Let us start by analyzing isolated avoided crossings in more detail. Examining the Floquet diagram from Fig. 1, we postulate that close to a crossing the Floquet Hamiltonian within the subspace of the two crossing states is where we note that the operators s α do not act on the physical spins but are to be understood as acting on the states within the two-dimensional subspace of crossing energy levels. This model features an avoided crossing at δ = 0, where the energy gap is ∆ c 2 . The eigenstates are (anti-)symmetric superpositions of the up and down states (in the subspace). For large δ the eigenstates are essentially the up and down states, however which one of these is higher in energy depends on the sign of δ. Going through the crossing the states switch, meaning that the up / down states have the energy of the opposite state before the crossing.
In a QMBS the first crossings are the ones between the edge states, thus after the first crossing the effective Hamiltonian is the average Hamiltonian with the two outermost edge states switched. Hence, already at this point one needs an additional many-body operator in the Floquet Hamiltonian, leading to the breakdown of perturbative expansions. This also means that if the crossings are well separated in τ , there is no energy absorption with respect to the average Hamiltonian (from their respective subspace) outside of the close vicinity of the crossing. At the crossings on the other hand, there are Rabi oscillations between the corresponding states resulting in "fast" dynamics, which might however still be very slow compared to the natural time scales of the total dynamics (see App. A).
The dynamics for the two-level system can be obtained exactly and is a classic result [88]. In case of an initial diagonal density matrix ρ(0) = diag (P 0 (0), P 1 (0)) the probability in the ground state is As we will argue later, these oscillations are the basis of FGR and our method, for which we need first to obtain some quantities for the individual crossings. Suppose we wanted to evaluate the formula in Eq. (12) for a single crossing. For this we need T, δ, ∆ c and P 0 (0) − P 1 (0), but since we are interested in the resonant oscillation we set δ = 0. This amounts to knowing the crossing time τ c , the angular gap width ∆ c and the pair of states i, j which cross. To compute these quantities in practice, we first completely diagonalize the propagator U sw (τ ) for many values of τ . The τ -resolution required depends on how narrow the gaps are in θ and how far apart they are in τ . Currently, we work with a fine uniform grid in τ and resolve the crossing locations and the minimal gaps for different grid resolutions. In future improvements this could also be done using an adaptive grid or automated root finding techniques. The determination of the energy transfer with respect to H 0 can be done in different ways. Here we chose for simplicity a scheme where we track pairs of crossing states back to their energy as τ → 0, this is done by working backwards through the crossing history of the involved states. Other possible ways to determine the energy transfer would be to measure the expectation values of H 0 in the pair of states before and after the crossing, or -for our particular protocol -to determine the slopesθ(τ ) before and after the crossing. We leave these refinements for future work however. Details of the algorithm(s) are discussed in Appendix C.
In Fig. 4 we show the gaps for the driven chain with L = 8 at two different resolutions in τ . Here, one can see that the gaps change over several orders of magnitude in a small τ window. For small values of the avoided gaps ( 10 −6 ) the resolution has a visible effect, however for the larger gaps there is not any noticeable difference.

B. Heating from Crossings
For deriving a formula for the EAR from a Rabi oscillation, we follow the linearization of the expression (12) as in the derivation of the FGR [74,82,83]. This amounts to using the identity for the delta function, which leads to a linear rate rather than an oscillation. Of course, this treatment can be valid only under certain assumptions, for example that the probability transfer is small, which are discussed in more details in the cited literature. A possible interpretation is that the dynamics can be viewed as an off-resonant (far-detuned) Rabi oscillation.
Using the linearization procedure, we obtain the following formula for the EAR (details of the derivation along with a derivation of the FGR are laid out in Appendix B) where the sum runs over all modes and the avoided crossings c attributable to the corresponding mode, with T c the period, ∆ c the gap. (∆E c ) is (the absolute value of) the energy difference of the states that cross with respect to the average Hamiltonian. In the FGR this would simply be given by the frequency ω (or multiples thereof). Generally one can use the difference in expectation values for any observable to obtain the absorption rate for that particular observable. Finally, (∆P c ) is the difference in occupation of the crossing states. This of course depends on the instantaneous state of the system during the dynamics, however we will be using a high-temperature ansatz to obtain an estimate later. Given the discussion above, the formula has an intuitive interpretation: at each avoided crossing transitions with the rate (π/2)(∆ c /T c ) 2 δ (mω − ω c ) occur and transfer an energy corresponding to the energy difference per unit time. The probabilities are a sort of "balancing" factor, such that transfer is enhanced for a large difference and vanishes in the fully mixed state (compare this to (12)).
Presumably, the gap widths cannot be related directly to matrix elements in general. Therefore, an analysis of the convergence similar to the one in Sec. III C will not be possible in general. However, in some cases, for instances at fast or strong drives, a relationship to matrix elements similar to the FGR can probably be recovered by transforming into an appropriate frame and applying the same formalism therein. Also as discussed in Sec. V B the Floquet formalism close to a discrete time crystal is very similar to the one at τ = 0, hence we expect the same convergence as in the FGR at this point.
Note as well, that even if a relationship with matrix elements is given, this does not imply that the sum over gaps in a finite system yields an accurate estimate for the thermodynamic limit. For example driven many-body localized systems have been reported to not heat in certain frequency windows [36][37][38]77]. Clearly though, the naive evaluation of the formula for a finite, small to intermediate size, system would yield an observable rate. The lack of heating thus has to stem from a different scaling with system size of the number/width of the gaps compared to the FGR case and extracting this behavior would require a more in depth analysis than just the evaluation of the formula. Nevertheless, for ergodic systems or systems, possibly transformed to a suitable frame, the simple treatment can be justified and convergence with system size is expected.
This formula, along with the automated resolution of crossings, is the central result of our work. In the next paragraph, we will show that for weak drives it is equivalent to the FGR, but in subsequent sections it will also become clear that it has a much larger region of validity, since here the actual crossing in the concrete system are used instead of perturbative approximations.

C. Comparison with FGR
Comparing formulas (8) and (13) (note that the double sum in the FGR is actually a single sum due to the delta functions as well), we recognize that the FGR is a special case of the crossing based formula, wherein the crossing time as well as the energy transfer are given by the energy difference ω ij be-tween the states. This would be the case if the lines in the Floquet diagram were perfectly straight lines, which is reasonable for weak drives. Furthermore, the matrix elements and gap widths have to be related by which corresponds to the gap width one would obtain in the Rabi model in the two-state subspace as discussed in Appendix B.
We investigate this relation numerically by introducing a factor g for the drive strength, changing V → gV in the protocol. We then compare the exact matrix elements with the appropriate expressions from our computed gaps for the first mode (g 1 = 4/π for the square drive). The results are shown in Fig. 5 for L = 8 and drive strengths g = 0.01 and g = 1, where the latter corresponds to the original model.
Here we observe that for g = 0.01 the correspondence is very good, apart from a region at very small τ where the gaps are limited by the resolution and some very small crossings which are likely due to "multi-photon resonances" i.e. levels with energy difference ω ij meeting at frequency ω = mω ij in the Floquet diagram (one can see this visually in Fig. 1   levels with the largest slope meet a second time (ω ij = 2ω) within the τ −window). At g = 1 the expressions are still qualitatively similar, but especially the locations of the gaps are noticeably different. This means that strictly speaking the assumptions of the FGR are not valid anymore for the model i.e. states do not cross at a τ given by the energy difference. However, for larger systems the density of crossings will increase and these small corrections will be washed out allowing the FGR to still make a good prediction. In this sense, in the rough region of validity of the FGR we expect our method to not improve predictions significantly. However, in the following sections we will show two scenarios in which our method, has clear advantages.

V. DRIVEN SPIN CHAIN
Having described the method in general, let us now use the driven spin chain as a concrete example to demonstrate the power of our approach. The gap widths and energy transfers obtained as explained in Sec. IV A are shown in Fig. 6 for different even system sizes, ranging from L = 6 to L = 12. On the y-axis we plot the half-period τ of the protocol. The points denote identified crossings, while the color scale of the points encodes the absolute value of the energy difference with respect to H 0 between the two Floquet eigenstates involved in the avoided crossing. Furthermore, we have indicated a dashed line at the mean angular spacing∆ = 2π/ dim H dictated by the dimension dim H of the relevant Hilbert space sector H.

A. High-Frequency Region
Let us first focus on smaller τ values τ 2: we observe that the magnitude of the gap widths increases over 4-5 orders of magnitude in this τ window. As the system size increases gaps at increasingly small values of τ appear as the spectrum now contains states with the corresponding frequency difference. For the larger system sizes some of these gaps are limited by the resolution resulting in blob like structures, which we color light gray. Furthermore, the magnitude at fixed τ remains roughly constant unless it is would be larger than the mean level spacing, which then acts as a cutoff for the magnitude of the gaps. We color the region where this is the case in darker gray.
Finally, on a technical note, there seem to be some crossings with a magnitude and energy transfer that do not fit the overall picture. These occur for two reasons: first the algorithm as outlined in the previous section is very sensitive even to small wiggles between the distance of adjacent levels and therefore detects some "ghost" crossings even between levels that seem to evolve mostly straight. Since these "ghost" crossings are not accompanied by an actual swap of two states, these cause the ordering of switched levels during the algorithm to become inaccurate after a while. For a further discussion and illustration of the wrong order introduced by "ghost" crossings see Appendix F. As a second reason, it turns out that many of the seemingly wrong crossings at τ ≈ π are actually genuine crossings with a switching. The non-fitting magnitude and energy transfer here result from the fact, that levels originating at τ = π cross within the same subspace (take the two subspaces in Fig. 1 as an example) and therefore have very similar energies. In principle one could correct for those by including some sort of curvature check in the algorithm to determine if a switching really took place or by simply discarding crossings where magnitude and energy transfer do not fit together. However, we found no noticeable effect of the non-fitting crossings, since they have a small energy transfer and there are relatively few of them (compared to "standard" crossings). Therefore, we move forward using the most straightforward scheme of the algorithm.
Let us finish the discussion by trying to understand the significance of the region, wherein the average angular level spacing due to the Hilbert space dimension is smaller than the gap width for small sizes (gray region in the figure) for our method, which is tied to some fundamental questions about heating in Floquet systems, particularly to the thermodynamic limit. To our best knowledge some of these questions have no definite answer yet, hence we give our best attempt at an interpretation of the results in the literature related to these questions in the context of our method. Clearly, in the thermodynamic limit the average level spacing vanishes and thus the information about microscopic processes as encoded in the gap widths is somehow hidden. The Floquet propagator then "has properties of matrices of the COE of random matrix theory" in the words of [34]. However, we still expect some structure depending on the frequency based on the results in [47][48][49], wherein the exponential timescale in the heating time at high frequency was established for many-body systems. This, along with the real-time simulations in current and other works (see references in Sec. I) suggest that the EAR converges in the thermodynamic limit (as also discussed in Sec. IV B). Hence, some of the structure visible at small sizes survives. In what form the information about the timescales enters the Floquet propagator for large systems is unclear to us. It might be that there are traces hidden in the spectrum, for instance there is a mechanism in the ETH leading to a "shrinking" of matrix elements with system size through the factor 1/ D(E). However, it is doubtful, whether the ETH formalism can be applied at finite frequency for large sizes, because the eigenstates of the Floquet propagator are likely to be fully mixed in the basis of the average Hamiltonian at those sizes. Thus, the spectrum might also be (statistically) equal at all frequencies for large enough sizes. This latter scenario seems to be consistent with the results in [34,35]. Therefore, we can only operate under the assumption that the information we extract is indeed relevant for large sizes without proof. Staying within this assumption though, we see that smaller sizes have a larger frequency window, wherein the gap width is separated from the mean level spacing. However, having chosen a frequency, one should strive for the largest possible sizes, for which the gap width is still unaffected, because larger sizes lead to a much more accurate estimate for the density of states and a finer frequency resolution due to more available gaps, which are needed to compute a smooth curve for the energy  Figure 6. Gap widths and energy transfer (color code) with respect to H0 for different system sizes. The kink at small τ for larger sizes (light gray) is due to finite resolution in τ (5 · 10 −8 for τ ≤ 1 and π · 10 −6 otherwise). The results from expanding the states from τ = 0 and τ = π are shown in color, while in the central region (gray) the true energy transfers cannot be obtained by the expansion. absorption.

B. Commensurate Points
Let us now focus on the upper half of the τ window from τ ≈ 2 up to τ = π. Due to the discrete nature of the protocol and the commensurability of the coupling strengths in the Hamiltonian, the Floquet Hamiltonian is not simply chaotic for all small frequencies. Instead, at some frequencies the propagator resp. Floquet Hamiltonian take on simpler forms, which shows in the spectrum as the appearance of degenerate subspaces (in our case one or two -depending on system size). This effect can clearly not be captured by a perturbation theory based on the average Hamiltonian and thus is not captured by the FGR. A detailed analysis, carried out in Appendix D, reveals that at all integer multiples of π the propagator takes on "simple forms" (but not always in the same way). At τ = π the system features a discrete time crystal [24,25] (albeit a fine-tuned one), since U 2 sw (π) = I and therefore there is no heating, but instead completely periodic dynamics with a doubled period.
In general such non-perturbative points at τ , where T H Fl (τ ) = H , can be integrated into the general formalism of Floquet expansions by expanding around τ . Writing dτ = τ − τ and expressing the propagator as the approximate Floquet Hamiltonian can again be obtained through the BCH series The radius of convergence is certainly more questionable for this expansion, and it might be more appropriate to transform into the rotating frame of H here. However, we do not make explicit use of the expansion and its only virtue is to show, how the average Hamiltonian appears away from the highfrequency limit. In fact the commutator term does not have matrix elements in the degenerate subspaces, therefore H 0 is responsible for the splitting to first order, irrespective of what H actually is. This can be observed in Fig. 1 since H Fl at θ equal to π and 2π are different from one another. For the dynamics this means that close to the time crystal the dynamics is a combination of the fast (period 2 cycles) dynamics and the much slower heating. The point here is that our method detects this, as exemplified by the vanishing of the gaps in Fig. 6 close to τ = π, and therefore the gap widths and the crossing locations can be obtained without any changes to the algorithm. For our chosen heuristic to determine the energy transfer by tracing the crossing states back to their initial energy at τ = 0, we need to alter the reference point to τ = π in the regime close to τ = π. This is however only a limitation of our simplistic heuristic, and a more robust determination of the energy transfer using previously mentioned ideas would not require a reference point to start with.

C. Heating Rates
We are now in a position to benchmark our avoided crossing spectroscopy method with large-scale real-time simulations as well as the FGR predictions for the driven spin chain. For the evaluation of Eq. (13) and Eq. (8) we follow [59] in using a high-temperature thermal state (usual Boltzmann distribution expanded to first order in β) as a model for P i , as we expect the evolved state to be sufficiently mixed in the Hilbert space for large parts of the dynamics, and using a broadened delta function, for example a normalized Gaussian with width dE, mimicking the density of states in the thermodynamic limit (see also Appendix C).
From the EAR the heating rate Γ = 1/t h is obtained where the subscripts indicate the energy evaluated at high-and infinite-temperature. Using the high-temperature expansion, the resulting heating rate is independent of temperature and should therefore give rise to a mono-exponential decay of the energy density towards zero. A more careful treatment would take into account the concrete occupations, which might be incorporated into a sort of Boltzmann equation using the ideas developed in this work, however we restrain from this here since our goal is to get a feeling for the time scales involved and especially to identify the region in τ , wherein the heating time changes drastically as discussed in previous sections. The heating rates obtained with the different methods (for dE = 0.1) are shown in Fig. 7, which features heating times (measured in cycles) extracted from real-time simulations for three product states, the prediction based on FGR and the predictions based on avoided crossing spectroscopy for different system sizes. For the latter method we estimate a range of validity following the discussion in Sec. V A.
The figure summarizes the earlier arguments, so let us also go through the main features again: for high frequencies the heating time increase rapidly and changes by several orders of magnitude, which is captured by both our method and the FGR. For lower frequencies, the FGR predicts a continuous decrease of heating times, while the observed times increase again due to the commensurate point at τ = π. This is cap-  The dotted parts indicate the estimated range of validity for each system size. As discussed, for smaller τ the FGR and gaps agree well, while the FGR is unable to detect the commensurate point at τ = π. Overall, both methods resolve a variation in the heating time of about five orders of magnitude, over a small change in half-period / frequency. tured by the computed gap widths and using our method corresponding heating times can be extracted. The range of validity decreases with the system size, hence smaller (and therefore computationally cheaper systems) can provide a more accurate estimate of the timescales involved. This comes however at the cost of featuring a lower density of crossings making it more difficult to obtain smooth curves for the heating times, if relying on broadening the delta function in the computation.
Overall, the proposed avoided crossing spectroscopy coincides with the FGR at high frequency, but also is accurate in resolving the temporal stability of the discrete time crystal as τ → π. It is impressive that the computation based on minimal assumptions such as the high-temperature ansatz and rather small systems sizes, ranging from L = 6 to L = 12, provide heating rates ranging over several orders of magnitude and capturing the regimes of rapid changes in the heating timescale very well. To increase accuracy one would need to improve on the model for occupations (the differences here are likely responsible for the different rates depending on the initial state) and to use different ways to introduce a density of states than to broaden the delta function for smaller sizes.

VI. SYSTEMS WITH FREQUENCY DEPENDENT COUPLINGS
In the cases discussed above the average Hamiltonian played an important role, which could be understood within the expansion. If the coupling strength depends on frequency itself though, specifically if it diverges with frequency, the Floquet Hamiltonian is not necessarily given by the average in lowest order. This allows to simulate dynamics (within a given time scale) with a Hamiltonian that may otherwise be inaccessible and thus is an important tool in modern experiments (see references in the Introduction, Sec. I). Note that oftentimes in the analysis of Floquet systems, no specific functional dependence of the frequency is specified a priori. Rather it turns out, that naturally a coupling strength ∝ ω results in a sensible high-frequency limit. A well known example is the modification of tunneling in Bose-Einstein condensates [4,6]. More recently, setups with strong couplings have been studied outside of the high-frequency limit and shown similar features [43,89]. As in the frequency independent case, different methods can be used to formulate high-frequency expansions for the effective Hamiltonian (see references in Sec. II A), which however often result in infinite series that cannot be summed analytically [65]. Thus, we refer to the literature for the full details and content ourselves with a sketch of the argument using the BCH series here.
We consider the switched setup from earlier, but make the drive strength proportional to the frequency V → (1/τ )V . Therefore, the expressions H ± τ appearing in the propagator are given by H 0 τ ± V . The BCH series consists of nested commutators, including commutators of the form [. . . [H 0 τ, V ], V ], . . . V ] (or other orderings), which, different from the independent case, are all O(τ ) and thus contribute to the Floquet Hamiltonian. Hence, H 0 is only one of the (typically infinitely many) terms at the lowest order. Also, the terms can introduce interactions of all ranges and lead to very complex Hamiltonians, even from basic ingredients. In the remainder of this section, we will show that our method can provide useful results even in this setup, when neither the effective Hamiltonian nor the effective drive is available 3 .

A. Driven Chain with Frequency Dependent Couplings
In order to illustrate the effect of frequency dependence, we stay with the driven spin chain from previous sections and modify it slightly V → (4/τ )V , where the factor of four is chosen such that the additional terms have visible effects, but are not strong enough to change the overall scales significantly, such that we can operate in the same τ windows as before. The main conclusions concerning the applicability of our method are however not dependent on this choice, as will be apparent from the discussion. In Fig. 8 we show the eigenangles colored by 0 , where we can see that the slopes are not governed by 0 and also display a stronger curvature overall. Also, the commensurate Floquet points at τ = π, 2π vanish as expected.   Figure 9. Energy transfer of the frequency dependent driven chain with L = 10. Compared to the frequency independent case the transfers do not correspond to the frequency and the shape of the decay region appears somewhat changed.
In Fig. 9 we show the energy transfers with respect to H 0 obtained again by tracking the switchings between states. The overall behavior of the gap widths at high frequency looks similar to the frequency independent case, therefore we verify additionally that the new Floquet Hamiltonian is in fact significantly different from the average in Appendix E. This is also visible from the mismatch between the energy transfers and the frequency of the drive, in contrast to the frequency independent case.

B. Heating Rates
In Fig. 10 we finally show the estimated heating rates, evaluated with dE = 0.3, as well as extracted rates from real-time simulations. Again, the dotted parts indicate the estimated range of validity for a given size.
As in the frequency independent case the agreement is reasonable overall, while being better at lower frequencies for the smaller sizes and better at higher frequencies for the larger sizes. Furthermore, in the Floquet diagram shown in Fig. 8 we observe that the energy (color code) seems to change during the evolution, as at the bottom there are no saturated levels while at intermediate τ there is some saturation. Hence, the extraction of energy transfers based on the original values also potentially leads to a lower accuracy. (symbols) and the gap data for different system sizes (lines). The dotted parts indicate the estimated range of validity for each system size. The agreement overall is reasonable, with different sizes showing good agreement in different parts as discussed in the text.

VII. CONCLUSION AND OUTLOOK
In this article we have shown how to analyze isolated avoided crossings in Floquet systems and how one can construct a versatile and accurate estimate for heating times based on those crossings. We have discussed that this method is closely related to the FGR, but with the demonstrated potential to go beyond it, since the crossings include nonperturbative effects. For this we have given two concrete examples using a driven spin chain: a discrete time crystal for commensurate points in a digital Floquet setting and a Floquet Hamiltonian beyond the average Hamiltonian due to frequency dependent couplings. In Appendix F we have also shown that the setup can be used to detect non-generic behavior (here weak symmetry breaking) in a seemingly generic system. Furthermore, throughout the paper we have discussed how the method combines microscopic and statistical aspects and how this understanding can be used to understand why the method performs well in small systems and to obtain estimates for the region of validity at at given system size.
The approach introduced in this paper has the potential to address and potentially solve long-standing issues, such as the detailed heating dynamics in driven Bose and Fermi-Hubbard systems, where multiply occupied sites seem to have slow dynamics, and we also believe avoided level-crossing spectroscopy in an adapted form is able to shed light on the intricate relaxation and thermalization dynamics of non-integrable quantum many body systems, such as the quenched Bose-Hubbard model [90,91].
For most of the numerical computations and the creation of the figures we use Python [92] with the packages [93][94][95][96][97]. The data for all figures, as well as corresponding plot scripts, are freely accessible online [98].

ACKNOWLEDGMENTS
The authors acknowledge support by the Austrian Science Fund FWF within the DK-ALM (W1259-N27). The computational results presented here have been achieved (in part) using the LEO HPC infrastructure of the University of Innsbruck.

Appendix A: Rabi-Oscillation Example
In the main paper we discussed that at individual crossings the eigenstates of the Floquet Hamiltonian in the subspace are the (anti-)symmetric superpositions of the original states, which therefore perform a Rabi oscillation, with the probabilities oscillating as where ∆ c is the gap width at the crossing. We verify this for two distinct crossings, by identifying the gap width, gap position and the crossing states using the methods discussed in the previous section. As we will see, the resonance region is very narrow, hence for these specific crossings we manually improve on the exact values.
The results for the first (highest frequency) crossing can be seen in Fig. 11. In the figure we show the dynamics of the lowest energy eigenstate of H 0 for three different values of τ (before, at and after the resonance). One can see clearly that the dynamics is restricted to the subspace of the lowest and highest energy state to a large degree. It is not fully in the subspace, since outside of the vicinity of the crossing, the eigenstates of the Floquet Hamiltonian are eigenstates of H 0 with perturbative corrections. Given that to obtain the heating rate one needs to sum over all states in a given window, this supports the argument that energy absorption is governed by off-resonant oscillations.

Appendix B: Fermi Golden Rule and Energy Absorption
The Fermi Golden Rule can be derived from timedependent perturbation theory. In this Appendix we present the main steps in the derivation, mostly following [82], and then show, how these ideas can be used to derive Eq. (13). Finally, the relationship between the gap widths and the matrix elements discussed in Sec. IV C will be justified.
In the derivation of the FGR, we are first concerned in the transitions between eigenstates |n of H 0 under an evolution generated by H(t) = H 0 + m>0 g m sin(mωt)V , with g m ∈ R.
In the interaction picture, the propagator can be approximated by the first term in the Dyson series and our goal is to compute the transition probability P nk (t) = | k|U (t)|n | 2 . Some steps can be performed exactly where ω kn = E k − E n and V kn = k|V |n . If we were to expand the absolute value, we would get a double sum over modes with "mixed" and "diagonal" terms. Within the diagonal terms, there are also "mixed" terms stemming from different denominators. It can be argued [82] that the contributions from the "mixed" terms can be neglected, and the remaining expression is where the identity |e ixt − 1| 2 = 4 sin 2 (xt/2) was used. The expression can now be "linearized" using the representation for the Delta function. Inserting this yields (B5) and we can define the transition rate Γ nk =Ṗ nk Γ nk = π 2 |V kn | 2 m>0 g 2 m (δ(ω kn + mω) + δ(ω kn − mω)) .

(B6)
We now consider a state with occupations P n and determine the change in energy due to transitions with rates Γ nk . Each transition has an energy transfer rate (E k − E n )Γ nk = ω kn Γ nk and the total energy absorption rate is given bẏ This is precisely the FGR as stated in Eq. (B). We now apply some further manipulations to get the formula to a form closer to Eq. (13). For this we first consider the exchange of indices n, k in the sum: P n → P k , Γ nk → Γ kn = Γ nk , ω kn → ω nk = −ω kn . We therefore can rewrite the sum aṡ E = π 2 m g 2 m n>k ω kn (P n − P k )|V kn | 2 × (δ(ω kn + mω) + δ(ω kn − mω)) .
Finally, we recognize, that due to the delta functions only terms with exactly matching energies contribute, hence the double sum is a single sum in disguise, and we can write this asĖ Comparing this to Eq. (13), we recognize that both coincide, given that ∆E c = m|ω kn |, ∆P c = P n − P k and ∆ 2 T 2 = g 2 m |V kn | 2 . Note that here we order the states such that ω kn > 0 and thus only one Delta function is included. Furthermore, the sum over all avoided crossings implicitly includes the sum over modes, since in the weak drive (FGR) regime the levels with ω kn = mω meet at T = 2π/ω as discussed in Sec. II B. Therefore, the avoided crossings include contributions from all modes.
Having derived the FGR from perturbation theory, we now consider deriving Eq. (13) from the dynamics of isolated avoided crossings. We begin though, by briefly recalling the main results from the Rabi model [74,83], using which we can make a connection between the effective Hamiltonian in the subspace (Eq. (11)) and matrix elements in a weak drive limit. For this we look at the dynamics of a two-level system under a single mode drive described by the Hamiltonian This Hamiltonian can be solved exactly within the "rotating wave approximation". The solution for the transition probability is [83] P nk (t) = From this expression, Eq. (B2) can be obtained by taking the high detuning limit (ω − ω nk ) We now compare the full solution to the dynamics in the static model H = δ T s z + ∆ T s x , where s x could be replaced by a combination of s x and s y with the same spectrum. The solution reads [88] Comparison to Eq. (B11) shows that in the weak drive limit ∆ T = g m |V nk | and δ T = mω − ω nk , provided that the resonance of the m−th mode is targeted. Of course the entire derivation assuming one mode and a fully decoupled subspace is not strictly valid, however in the regime with very small gaps the levels are well isolated and the energy differences are reasonably large for the resonances from different modes to be well separated.
Finally, let us consider the off-resonant (high detuning) limit off an avoided crossing characterized by ω c and ∆, where as above we identify δ/T = mω − ω c . Note that here we dot necessarily associate ω c with the energy difference between the states, allowing for more general scenarios, such as the discrete time crystal point discussed in Sec. V B. The transition probability then reduces to Following exactly the same procedure as described before, we arrive at Eq. (13). the result is well in the regime independent of β and thus use this value as default. A discussion of the ansatz can also be found in the Supp. Mat. of Ref. [59]. For larger sizes this procedure is more stable and leads to similar result over a range of widths, while smaller sizes are more sensitive due to the low density of crossings.

Appendix D: Propagator at Commensurate Points
We want to understand why the Floquet propagator takes on simple forms at integer multiples of π. Remember that the propagator is given by The single particle terms can be evaluated using the rotation formula x/z i and two particle terms by writing out the matrix in the computational basis , where i, i + 1 denote the spins upon which the matrix acts (note that at the boundary i = L this has to be understood rather formally). Using this one can evaluate the propagator at multiples of 2π easily • U sw (8π) = I (thus the entire angle diagram repeats between 8nπ intervals) • Usw(4π + 8nπ) = (−1) L I.
The situation at π is more difficult, since the single particle terms do not vanish but instead combine to i L i σ y i . This means that the full propagator is given by We are unable to find a full expression for the Floquet Hamiltonian, but can prove that this propagator has a period 2 dynamics for even system sizes, meaning that U 2 sw (π) = I, thus it is a toy example of a (stable but fine-tuned) discrete time crystal. For this let us work in the computational basis (product states along z-axis): i σ y i is completely anti-diagonal and the other term completely diagonal. One can readily verify that U 2 sw is then a diagonal matrix with entries The matrix elements of i σ z i σ z i+1 are −2 + 2k, with L = 2 and k being the number of kinks on top of the fully polarized state. Substituting this and some straightforward algebra leads to the claimed result that U 2 sw (π) = I. On a final note we would like to point out that the Floquet Hamiltonian at π is not a simple single particle operator, since for some product states we observe an increase in the bipartite entanglement entropy upon action with the propagator.  Table I. Energy density of product states with respect to the average Hamiltonian ( 0) and the full lowest order term . The results are based on L = 18 and τ = 10 −3 .
As one can see the energy densities of the products states along the x-and z-axis change significantly and in the case of x-states do not have opposite signs anymore, hence there must be additional operators with even number of X terms. The most striking change however is seen from the y-states, which go from a vanishing energy density to the largest/smallest one. Since they are almost the negative of another, the largest contribution comes from operators with odd number of Y . Finally, we also observe that the y-states are almost at the very edge of the spectrum, therefore we conclude that the Floquet Hamiltonian is significantly different from the average and that the heating rates, although looking similar are the result of a truly different dynamics. Additionally, to the driven spin chain in the main part, we consider a model of a spin chain, which has been studied in the context of heating in Floquet systems [52,53] and also can be regarded as an example system with weak symmetry breaking. The average and drive are given by with coefficients h x = 0.42, h y = 0.34, h z = 0.26, J x = 3, J z = 4, where the letters denote spin one-half operators analogous to the main text. Different from the cited works, we use periodic boundary conditions, which however does not seem to change the observed heating rates as well as further conclusions in this section. The average Hamiltonian is then invariant under translations, spatial reflection and spin flips about the z-axis generated by the flip operator The drive has the same spatial symmetries but is not invariant under spin flips, which leads to a (weak) breaking of this symmetry, whose effects we will observe in the spectroscopic approach. Due to the spatial symmetries we can again focus on the zero momentum and positive parity sector. Anticipating a near-conservation of spin flip parity ( F ) we color the levels in Fig. 13 by F instead of 0 and can observe clearly in Fig. 13a that the eigenstates are (almost) fixed parity states and in Fig. 13b that the matrix elements of states with the same parity are significantly smaller than the ones between different parity (this result can be explained by the oddness of the drive under spin flips [53]). In fact, the small gaps seem to be limited by resolution and might be much smaller in actuality, however they will not affect the heating rate in any case. Looking closely at the Floquet diagrams, it seems that the smaller gaps should be attributed to same parity states for all τ values and not be mixed as in the figure. The observed mixing should rather be understood as illustrating the effects of ghost gaps, as discussed in the main part.
Furthermore, upon inspection of the diagrams, we observe phenomena that we call "nested crossings": here there is a switching between two states which are not neighbors in the Floquet diagram, but rather separated by a middle state which seems unaffected. We have not observed such crossings in the driven chain and suspect that they are related to the (nearly) spin flip symmetry. Anyway, this seems to not affect the heating times displayed in Fig. 13c and further supports the robustness of the method, while also suggesting the study of crossings with near symmetries as a potential future direction.