# Reconstructing effective Hamiltonians
from nonequilibrium (pre-)thermal steady states

Sourav Nandy  Jožef Stefan Institute, 1000 Ljubljana, Slovenia  Markus Schmitt
Forschungszentrum Jülich GmbH, Peter Grünberg Institute, Quantum Control
(PGI-8), 52425 Jülich, Germany  University of Regensburg, 93053 Regensburg,
Germany  Marin Bukov  Max Planck Institute for the Physics of Complex Systems,
Nöthnitzer Str. 38, 01187 Dresden, Germany  Zala Lenarčič  Jožef Stefan
Institute, 1000 Ljubljana, Slovenia

###### Abstract

Reconstructing Hamiltonians from local measurements is key to enabling
reliable quantum simulation: both validating the implemented model, and
identifying any left-over terms with sufficient precision is a problem of
increasing importance. Here we propose a deep-learning-assisted variational
algorithm for Hamiltonian reconstruction by pre-processing a dataset that is
diagnosed to contain thermal measurements of local operators. We demonstrate
the efficient and precise reconstruction of local Hamiltonians, while long-
range interacting Hamiltonians are reconstructed approximately. Away from
equilibrium, for periodically and random multipolar driven systems, we
reconstruct the effective Hamiltonian widely used for Floquet engineering of
metastable steady states. Moreover, our approach allows us to reconstruct an
effective quasilocal Hamiltonian even in the heating regime beyond the
validity of the prethermal plateau, where perturbative expansions fail.

##  I Introduction

A central idea at the core of modern quantum simulation is the utilization of
highly controllable quantum degrees of freedom to emulate the behavior of
complex quantum systems [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Besides
developing the corresponding state preparation and stabilization techniques,
this requires the implementation of the Hamiltonian that governs the physics
of the system to be simulated, in the first place. To do this, different
approaches have been developed, which leverage specific properties of the
underlying quantum system. For instance, analog simulators, such as ultracold
atoms [9, 14, 15, 16, 11] and trapped ions [3, 5, 4], are natural platforms of
short- and long-range interacting Hubbard and Ising Hamiltonians. By contrast,
the digital approach to quantum simulation [17, 18] implemented in
superconducting qubits, Rydberg atom [19, 20, 21, 22, 23] or trapped ion [24,
23, 25] platforms, relies on constructing a time-evolution operator using
consecutively applied discrete unitary gates.

Whenever the desired Hamiltonian cannot be implemented using the available
building blocks in a straightforward way, nonequilibrium techniques, such as
Floquet engineering [26, 27, 28] or dynamical decoupling [29, 30, 31, 32],
become handy. Recently, periodic drives have enabled the realization of non-
equilibrium phase transitions [33, 34], topologically non-trivial systems [35,
36, 37, 38, 39], artificial gauge fields [40, 41, 42] and discrete time
crystals [43, 44, 45, 46, 47, 48, 49, 50]. Note that having an effective
Hamiltonian is not guaranteed for generic nonequilibrium drives due to the
lack of energy conservation; the necessary ingredients for the existence of
effective Hamiltonians are currently the subject of intensive on-going
research.

Figure 1: Schematic presentation of the Hamiltonian reconstruction algorithm
proposed in this work. (a) The input for the Hamiltonian reconstruction are
measurements of local operators performed in (locally) (pre-)thermal (steady)
states at different temperatures, for example, obtained from a quantum
simulator. (b) These serve as input data for the autoencoder, which (c)
confirms that the data is indeed parametrized with a single parameter
(temperature) by checking that the decoding ‘test error’ drops at a single
neuron in the latent space, and (c, inset) to generate the latent
representation of the dataset, which is one-dimensional for thermal data. (d)
In the next step, the potential candidate terms $O_{j}(\vec{\alpha})$ in the
Hamiltonian $H=\sum_{j,\vec{\alpha}}a_{\vec{\alpha}}O_{j}(\vec{\alpha})$ are
pre-selected as the operators with the largest average gradient of expectation
values $\overline{\partial\langle O_{j}(\vec{\alpha})\rangle}$ along the
latent representation. (e) In the last step, the Hamiltonian coefficients
$a_{\vec{\alpha}}$ are fixed. Here, we show actual data from effective
Hamiltonian reconstruction in the prethermal plateau for our Floquet protocol,
Eq. (9), presented in Fig. 3.

Whichever the implementation method may be, engineered effective Hamiltonians
often come with additional, unwanted terms. These can be residual (longer-
range) interactions, or leftover single-particle field terms; they can also
arise due to coupling to existing higher-energy levels or additional cavity
modes [51]. Moreover, nonequilibrium protocols used for Hamiltonian
engineering rely on the existence of prethermal (meta-)stable steady states,
in which heating (i.e., energy absorption and entanglement production) is
suppressed up to controllably long times [52, 53, 54, 55, 56]. Nevertheless,
residual terms in the effective Hamiltonian (e.g., higher-order corrections)
eventually lead to accumulating detrimental heating. In parallel, verifying
that a desired Hamiltonian is implemented to a given precision, is a crucial
requirement for a reliable quantum simulation; moreover, certification of the
realized Hamiltonians becomes indispensable for benchmarking the performance
of present-day noisy intermediate scale quantum (NISQ) devices [57, 14, 58,
59, 60]. For all these reasons, developing a procedure, based on measurements
of local observables, and capable of identifying and quantifying the
Hamiltonian implemented in quantum simulators to a high precision, is a
problem of ever-increasing importance.

During the last few years, several methods for Hamiltonian reconstruction have
been put forward: from measurements on a single eigenstate [61, 62, 63, 64]
using trusted quantum simulator [65, 66, 67]; or from eigenstate dynamics [68,
69, 70, 71, 72, 73, 74], via compressed sensing [75], using restricted
Boltzmann machine tomography [76] and neural networks [77]. One common feature
of these techniques is that they explicitly assume the underlying Hamiltonian
to be both local and a translationally invariant sum of few-body terms.
Furthermore, only a few of the above-mentioned methods are applicable to
Floquet systems. Although reasonable progress has been achieved recently in
this direction [78, 79], the methods employed in these works demand an apriori
intuition from the Floquet-Magnus expansion for the selection of candidate
terms. However, a typical challenge in any Hamiltonian reconstruction
algorithm often lies in developing a systematic method to select candidate
terms. It is, therefore, highly desirable to propose an algorithm that relies
on an unbiased approach to select the candidate terms of the target
Hamiltonian.

In this work, we propose a machine learning-assisted algorithm for the
reconstruction of Hamiltonians, which uses as the input the expectation values
of local observables that are diagnosed to be thermal. We show that our
algorithm can precisely reconstruct local Hamiltonians, while long-ranged
Hamiltonians are approximately reproduced with certain limitations. It is
worth mentioning that the presence of only few-body terms in the underlying
Hamiltonian is advantageous for our method but it is not a necessary
requirement. We apply our algorithm to Floquet and random multipolar driven
systems with the purpose of reconstructing the effective Hamiltonian,
responsible for the stroboscopic expectation values in the prethermal plateau
and/or in the follow-up heating regime. For constructing the effective
Hamiltonian in the prethermal plateau, our algorithm does not require any a
priori intuition from the Floquet-Magnus expansion or any other similar
perturbative technique when the selection of candidate terms is concerned.
Furthermore, we go beyond prethermal steady states and analyze the heating
regime, where perturbative expansions are known to fail. We extract an
effective, quasistatic Hamiltonian, which reproduces the thermal measurements
of local operators as the system approaches the infinite-temperature state. We
directly observe that this effective Hamiltonian becomes less and less local
as the system heats up, verifying a hypothesis laid out in previous work [52,
80].

The rest of the paper is organized as follows. In Sec. II, we discuss in
detail the algorithm of Hamiltonian reconstruction. Sec. III deals with
illustrative examples where our method is used to reconstruct both short and
long-range interacting static Hamiltonian. This is followed by the central
Sec. IV, which addresses Hamiltonian reconstruction in dynamical Floquet
setups, starting with reconstruction in the prethermal regime. Thereafter, we
discuss the insight our algorithm provides into the physics of unconstrained
heating for both Floquet systems and non-Floquet systems that exhibit eventual
thermal death. Finally, we conclude and present future perspectives in the
Sec. V.

##  II Algorithm for Hamiltonian reconstruction using autoencoders

In our Hamiltonian reconstruction algorithm, we rely on a deep autoencoder
neural network [81, 82], Fig. 1(b), which is used for the dimensional
reduction of the data; in our algorithm this amounts to a preselection of the
Hamiltonian candidate terms. In our case, each element of dataset contains
measurements of different local operators in a given quantum state, see Fig.
1(a). Our approach is applicable if these states are (effectively) thermal.
The autoencoder is used to generate a latent representation of the dataset and
to verify that the dataset indeed contains measurements on thermal states;
whenever this is the case, the latent representation is a one-dimensional
manifold, Fig. 1(c). Such compressed representation of measurement helps us
single out the candidates terms of the Hamiltonian, Fig. 1(d), significantly
reducing the complexity of the Hamiltonian reconstruction, Fig. 1(e),
performed in the subsequent step. Previously, some of us proposed part of the
algorithm as a way to reconstruct spatially local Hamiltonians [83]. Here we
extend the procedure to deal with long-range Hamiltonians as well. In the
following, we explain the algorithm in more details.

The autoencoder is a deep neural network, where the input and output layers
have the same dimension, while at least one of the middle layers has a
considerably smaller number of neurons $\mathcal{N}_{L}$ and represents a
bottleneck, Fig. 1(b). The part prior to the bottleneck is called the encoder,
which maps the input data to a lower-dimensional latent representation. The
decoder then maps the latent representation to the output. The network is set
by minimizing the difference between the input $\vec{x}$ and the output
$f_{\vec{\theta}}({\vec{x}})$ with respect to the network parameters
$\vec{\theta}$, averaged over a minibatch of data elements from the dataset
$\mathcal{D}_{T}$,

$\mathcal{L}_{\mathcal{D}_{T}}(\vec{\theta})=\frac{1}{|\mathcal{D}_{T}|}\sum_{\vec{x}\in\mathcal{D}_{T}}\big{(}f_{\vec{\theta}}(\vec{x})-\vec{x}\big{)}^{2}.$
(1)

The latent representation’s dimension, i.e., the intrinsic dimension of the
data, is extracted from the same function (called ‘test error’), evaluated for
the optimized neural network parameters $\vec{\theta}$ on the rest of the
unseen data elements. The intrinsic dimension of data is identified as the
number of neurons in the bottleneck at which the ‘test error’ drops
significantly and is flat upon further increasing $\mathcal{N}_{L}$, Fig.
1(c). Values at neurons in the bottleneck give the latent representation of
data. In the inset of Fig. 1(c), an example of latent represetation is shown
for a network with $\mathcal{N}_{L}=2$ neurons in the bottleneck. Further
technical details on our network are given in Appendix A.

Our Hamiltonian reconstruction algorithm relies on precondition (0) and
consists of Steps (1-5):

  * (0)

Work under the assumption that we are given thermal expectation values of all
local operators with support $\ell$ that is smaller than or equal to some pre-
set number $\mathcal{S}$. These operators are measured in states of different
temperatures: the set of all measurements for a given state represents one
data element, while different states yield the dataset.

  * (1)

Use autoencoder to check that the data elements indeed contain measurements
with respect to (effectively) thermal states. In that case, the latent
representation of the dataset is one-dimensional, with data elements within it
ordered with respect to the temperature, or equivalently, the energy of the
corresponding states.

  * (2)

Since data elements in the latent representation are ordered with respect to
the Hamiltonian expectation value, other operators that correspond to the
individual terms in the Hamiltonian also show a pronounced variation along
that one-dimensional manifold. Therefore, one can isolate candidate
Hamiltonian terms $O_{j}(\vec{\alpha})$,
$H=\sum_{j,\vec{\alpha}}a_{\vec{\alpha}}O_{j}(\vec{\alpha})$, via singling out
those operators
$O_{j}(\vec{\alpha})=\sigma_{j}^{\alpha_{j}}...\sigma_{j-1+|\ell|}^{\alpha_{j-1+|\ell|}}$
that have the largest average gradient $\overline{\partial\langle
O_{j}(\vec{\alpha})\rangle}$ in the expectation value along the latent
representation. We will assume periodic boundary conditions and a
translationally invariant system with $O(\vec{\alpha})\equiv
O_{j=1}(\vec{\alpha})$, but the procedure can in principle be performed in the
absence of that as well.

  * (3)

Step (2) will, in addition to the actual terms, filter out also some spurious
ones, which are not part of $H$ but are, for example, products of actual
Hamiltonian terms on neighboring sites. To fix the prefactors
$a_{\vec{\alpha}}$, and to get rid of the spurious terms, we compare the
‘trial’ thermal expectation values (with respect to the trial Hamiltonian) to
the actual measurements for one data element at a given (but unknown)
temperature $1/\beta$, and find the solution of the following equation

$\text{Tr}\left[O(\vec{\alpha}^{\prime})\frac{e^{-\beta\sum_{j,\vec{\alpha}}a_{\vec{\alpha}}O_{j}({\vec{\alpha}})}}{Z}\right]-\langle
O(\vec{\alpha}^{\prime})\rangle_{\beta}=0.$ (2)

Since one does not know the temperature of the state, $\beta$ remains
undetermined; what is actually extracted is the product $\beta
a_{\vec{\alpha}}$; $\beta$ can be determined separately, e.g., from dynamical
measurements.

For strictly local Hamiltonians and measurements that include all Hamiltonian
terms (i.e., the maximal support of measured operators $\mathcal{S}$ is larger
or equal to the support of $H$), Steps (1-3) will reconstruct $H$ (up to a
prefactor) with high precision (see examples below), as already proposed and
tested in Ref. [83]. In the next steps, we extend the algorithm from Ref. [83]
to be applicable to data from long-range Hamiltonians as well. For long-range
interacting Hamiltonians, some error in the reconstruction is inevitable,
stemming from the fact that we have measurements of only the most local terms
in the Hamiltonian and that we can at most aim at reconstructing a local
approximation of a long-range Hamiltonian. In that case, Steps (1-3) might
still generate finite prefactors for “ghost” terms within the measured
support. To remove those, we add here additional steps:

  * (4)

Compare the solutions of Step (3) at a few different $\beta$, i.e., for
different data elements. Ghost terms’ coefficients will have a large relative
variance across solutions for different states; hence, we drop terms with

$\frac{\mathrm{Var}_{\beta}(a_{\vec{\alpha}}/a_{\vec{\alpha}_{0}})}{\mathbb{E}_{\beta}(a_{\vec{\alpha}}/a_{\vec{\alpha}_{0}})}>O(1).$
(3)

where $a_{\vec{\alpha}_{0}}$ is the coefficient at the dominant term.

  * (5)

Repeat Step (3) for the Hamiltonian ansatz without the ghost terms. This step
just slightly corrects the result of Step (3).

The most expensive step of our reconstruction is Step (3); even after
performing Step (2), potentially a large number of trial terms remain for Step
(3). This can be particularly problematic for a Hamiltonian $H$ with large
support. In order to ameliorate this issue to a manageable extent, one can use
an iterative scheme:

First, perform Steps (2-4), selecting contributing operators $O(\vec{\alpha})$
with maximal support $\mathcal{S}_{k}=\ell_{k},k=0$:

  * (k1)

When increasing the maximal support to
$\mathcal{S}_{k+1}=\ell_{k+1}=\ell_{k}+1$, previously selected
$\mathcal{S}_{k}$ support terms are retained by default, while Step (2) is
performed only for operators with support $\ell_{k+1}$.

  * (k2)

Steps (3-4) are performed, including previously selected
$\mathcal{S}_{k}=\ell_{k}$ and new $\ell_{k+1}$ support trial terms from (k1).

  * (k3)

Terms surviving in (k2) are the new input for the next iteration.

In the end, Step (5) is performed.

##  III Reconstructing short- and long-range interacting static Hamiltonians

We now benchmark our approach for two classes of $H$: local and long-range
interacting Hamiltonians. The extent of applicability to the latter type of
Hamiltonians has not been discussed in Ref. [83]. In any case, each data
element consists of measurements in a thermal state (see appendix B for
details of how we prepare the initial thermal states) with temperature
$1/\beta$ for all local operators with support smaller than or equal to
$\mathcal{S}=4,6$. We use inverse temperatures from the interval
$\beta\in[0.05,0.5]$.

First, we consider strictly local Hamiltonians with maximal support
$\mathcal{R}$, for example,

$\displaystyle
H_{L}=\sum_{i}\left(h_{z}\sigma^{z}_{i}+\sum_{r=1}^{\mathcal{R}-1}(J_{1}\,\sigma^{x}_{i}\sigma^{x}_{i+r}+J_{2}\,\sigma^{y}_{i}\sigma^{y}_{i+r})\right).$
(4)

We find that the relative reconstruction error for coefficients
$h_{z},J_{1},J_{2}$ highly depends on the maximal support of the measured
operators $\mathcal{S}$:

$\displaystyle\mathcal{S}\geq\mathcal{R}:\text{relative error}\sim 10^{-6},$
(5) $\displaystyle\mathcal{S}<\mathcal{R}:\text{relative
error}\sim\mathcal{O}(1).$

Next, we consider long-range interacting Hamiltonians with power-law decaying
interactions,

$\displaystyle
H_{LR}=\sum_{i}\left(h_{z}\sigma^{z}_{i}+\sum_{r=1}^{\frac{N}{2}-1}\frac{J}{r^{\delta}}(\sigma^{x}_{i}\sigma^{x}_{i+r}+\sigma^{y}_{i}\sigma^{y}_{i+r})\right)$
(6)

with periodic boundary conditions.

Figure 2:  Evaluation of Hamiltonian reconstruction
$H=\sum_{j,\vec{\alpha}}a_{\vec{\alpha}}O_{j}(\vec{\alpha})$ for the long-
range interacting Hamiltonian from Eq. (6). (a) The relative error ${|\delta
a_{\ell}|}/{\bar{a}_{\ell}}$ of all reconstructed coefficients at a given
support $\ell$. It is seen clearly that for a fixed $\delta$ (measure of long-
rangedness of the Hamiltonian), terms with larger support have larger relative
error. Also, for a given support $\ell$, relative error decreases with
increasing $\delta$, i.e., considering more short-ranged underlying
Hamiltonian. (b) Time evolution of two different observables with respect to
(i) the original long-range $H_{LR}$ with $\delta=3$ from Eq (6), and (ii) its
reconstruction $H_{\text{ML}}$ using our method. This figure clearly shows
that even though Hamiltonian reconstruction is not perfect for long-ranged
systems, the reconstructed Hamiltonian nonetheless captures local dynamics
quite well. Parameters: $N=16,J=1.0,h_{z}=0.2,\delta=3.0$.

In this case, the reconstruction error highly depends on the maximal support
of the measured operators $\mathcal{S}$ and the interaction strength decay
exponent $\delta$. Fig. 2(a) shows the relative error ${|\delta
a_{\ell}|}/{\bar{a}_{\ell}}$, for compact representation averaged over all
reconstructed Hamiltonian terms with support $\ell$, i.e.,

$\displaystyle\frac{|\delta a_{\ell}|}{\bar{a}_{\ell}}$
$\displaystyle=\frac{|h_{z}-a_{z}|}{h_{z}},\hskip 11.38109pt\ell=1$ (7)
$\displaystyle\frac{|\delta a_{\ell}|}{\bar{a}_{\ell}}$
$\displaystyle=\frac{|2J/r^{\delta}-a_{xI\dots Ix}+a_{yI\dots
Iy}|}{2J/r^{\delta}},\hskip 11.38109pt\ell>1$ (8)

We see that terms with larger support have larger relative error. In general,
the latter increases for less local interactions (smaller $\delta$). However,
even if the Hamiltonian is not exactly reconstructed, its effect at the level
of dynamics is captured in practice. To support that, we show in Fig. 2(b)
time evolution from a domain wall initial state for the original and the
reconstructed Hamiltonian for a few local observables. For the sake of
comparison, we use the value of $\beta$ that the data was measured at to
divide the reconstructed coefficients $\beta a_{\vec{\alpha}}$ so as to obtain
the actual Hamiltonian. In a temperature-agnostic protocol, $\beta$ would be
determined as the one at which the comparison is best, as we demonstrate in
the next section.

##  IV Reconstruction of effective Hamiltonians engineered by nonequilibrium
drives

Having benchmarked our approach, we now turn to our central example: effective
Hamiltonian reconstruction in dynamical systems, which is of interest for
verification of quantum simulators, where desired Hamiltonians are often
engineered using periodic protocols. In these Floquet setups, the stroboscopic
evolution is generated by the so-called Floquet Hamiltonian; while such closed
systems are (in the absence of symmetries [84]) bound to eventually heat up to
a featureless infinite temperature state [85, 86, 87], the heating process is
delayed exponentially with increasing the frequency of the drive [52, 54].
When observed stroboscopically, at intermediate times, in the so-called
prethermal plateau [88], the system is locally effectively described by a
thermal state with respect to an effective Hamiltonian that can feature
nontrivial transient states [52, 53, 54, 55]. Floquet protocols have thus been
one of the central approaches in engineering non-equilibrium phenomena using
quantum simulators and materials.

In the inverse frequency expansion or iterative rotating frame transformation,
the effective Hamiltonian can be reconstructed up to some order [89]. However,
this is generally a daunting task [89, 53] and the former cannot be computed
in closed form for generic systems and protocols. Here we propose an unbiased
reconstruction, which does not rely on any expansion technique; instead, we
reconstruct the effective Hamiltonian up to support $\mathcal{S}$, effectively
resumming contributions from higher orders in the inverse frequency expansion
within this support. As the inverse frequency series is formally divergent,
the system eventually heats up; it is still not fully settled how this
divergence is related to heating dynamics featured in the exact stroboscopic
dynamics. Our approach reveals that stroboscopically, throughout the whole
heating regime, the system locally appears as if in a thermal state with
respect to an effective quasi-static Hamiltonian; moreover, we demonstrate
directly how its support grows in time.

###  IV.1 Floquet driving

Prethermal regime. In the following, we reconstruct the effective Hamiltonian,
characterizing the stroboscopic pre-thermal state in the prethermal plateau,
for the square pulse Floquet protocol

$\displaystyle U$
$\displaystyle=\exp\left(-iH_{1}T/2\right)\exp\left(-iH_{2}T/2\right)$ (9)
$\displaystyle H_{1}$
$\displaystyle=\sum_{j=1}^{N}J\sigma^{z}_{j}\sigma^{z}_{j+1}+h_{x}\sigma^{x}_{j}+h_{z}\sigma^{z}_{j}$
(10) $\displaystyle H_{2}$ $\displaystyle=\gamma\sum_{j=1}^{N}\sigma^{x}_{j}$
(11)

where $U$ is the unitary matrix for one drive cycle which acts repeatedly on
the system during the dynamics; $T$ denotes the time period of the drive. Fig.
3(a) shows time-dependent expectation values $\langle\sigma^{x}\rangle$ for
$J=1.0,h_{x}=0.9045,h_{z}=0.809,\gamma=0.4,\omega=2\pi/T=10.0$, for a time
evolution from an initial thermal state with respect to $H_{1}$ and inverse
temperature $\beta_{0}=0.1$.

For drive cycle numbers $n\gtrsim 5$, we observe an extended prethermal
plateau. Visible oscillations are a consequence of finite size ($N=12$) exact
diagonalization propagation. In the large frequency regime $\omega\gg 1$
considered, the Floquet Hamiltonian can be perturbatively approximated using
Baker-Campbell-Hausdorff (BCH) series expansion, as detailed in Appendix C.

Figure 3: (a) Stroboscopic measurements of $\langle\sigma^{x}_{N/2}\rangle$ as
a function of cycle number $n$, obtained from: (i) the exact Floquet protocol,
(ii) time evolution with respect to the BCH Hamiltonian up to second order in
$T$, and (iii) time evolution with respect to our reconstruction
$H_{\text{ML}}$. Initial states are thermal w.r.t. $H_{1}$ at $\beta_{0}=0.1$.
The inset demonstrates how we extract the overall multiplicative factor
$\beta$ of the coefficients obtained from our Hamiltonian reconstruction
algorithm, as discussed in the main text. (b) Comparison of Hamiltonian
coefficients computed from the BCH expansion and our Hamiltonian
reconstruction algorithm. Clearly, our reconstruction algorithm predicts the
coefficients (particularly the leading ones) quite accurately. Parameters:
$J=1.0,h_{x}=0.9045,h_{z}=0.809,\gamma=0.4,\omega=2\pi/T=10.0,N=12$.

Our Hamiltonian reconstruction is performed on a dataset with elements
containing measurement of local operators with maximal support
$\mathcal{S}=3$, in states that are evolved up to some fixed time within the
prethermal plateau. As initial states for the time evolution we consider
thermal states with respect to $H_{1}$, with inverse temperatures within the
range $\beta_{0}\in[0.05,0.40]$. In principle, one could start from non-
thermal states as well since these are expected to relax to thermal states in
the plateau anyhow, due to the nonintegrability of the exact Floquet
Hamiltonian. While Step (3) in the reconstruction algorithm [cf. Sec. II] can
in principle be performed for a single data element (a single state), we
average the relative coefficients $a_{\vec{\alpha}}/a_{x}$ ($a_{x}$ being the
dominant one) for each $\vec{\alpha}$ over a few $\beta_{0}$ and a few
measurement times in the prethermal plateau. As discussed in Sec. II, our
algorithm gives the coefficients up to an overall multiplicative factor
$\beta$. This overall multiplicative constant $\beta$ can be extracted from
dynamics in the initial time window before the system reaches the prethermal
plateau; as illustrated in the inset of Fig. 3(a). More precisely, we
determine $\beta$ by minimizing the difference
$\sum_{n}|\text{Tr}[OU^{n}\rho(0)(U^{\dagger})^{n}]-\text{Tr}[Oe^{-inTH_{\text{ML}}(\beta)}\rho(0)e^{inTH_{\text{ML}}(\beta)}]|$
between the exact and the reconstruction $H_{\text{ML}}$ based time-dependent
measurements of a dominant operator at early cycles.

In Fig. 3, our Hamiltonian reconstruction $H_{\text{ML}}$ is compared to a
second order (in $T$) BCH expansion, which contains terms with the same
maximal support $\mathcal{S}=3$. Fig. 3(a) compares the dynamics generated
from (i) the exact unitary, Eq. (9), (ii) the Hamiltonian constructed from the
BCH approximation calculated up to second order in $T$, as discussed in the
Appendix C, and (iii) our reconstructed Hamiltonian $H_{\text{ML}}$, for an
initial state with $\beta_{0}=0.1$. We emphasize that we have checked the
cases also for other observables for a range of $\beta_{0}$-values
($\beta_{0}\in[0.05,0.3])$ and all such results show similar features as shown
in Fig. 3(a).

Figure 3(b) compares the coefficients of different terms computed from the BCH
expansion to that obtained via our machine learning-assisted Hamiltonian
reconstruction algorithm. The matching is quite good for most (particularly
for the leading) coefficients. The slightly worse performance of our
reconstruction scheme in Fig. 3(a), compared to second-order BCH, is explained
as a consequence of finite size ($N=12$) effects, causing fluctuations on the
prethermal plateau and other deviations from a strictly thermal state,
affecting the reconstruction input data.

Figure 4: (a) Test error as a function of the latent space dimension
$\mathcal{N}_{L}$ (number of bottleck neurons) at different cycle number $n$
for Floquet driven systems, Eq. (12). At smaller cycle numbers e.g.,
$n=2,\;60$, the influence of initial states, prepared by applying random
single-site rotations to thermal states, is still visible. For larger cycle
numbers, the effectively one-dimensional latent representation signals that
the system is in a thermal state throughout the entire heating regime. (b)
Latent representation of dataset consisting of measurement at $n=180$ in the
heating regime. The $x$\- and $y$-axis show values of the first and the second
neuron in the bottleneck layer with $\mathcal{N}_{L}=2$ neurons, giving a one-
dimensional representation of data elements from different states. Parameters:
$N=18,\epsilon=0.08,J=1.0,h_{x}=0.9045,h_{z}=0.809,T=12.7263,\beta_{0}\in[0.05,0.45].$
Figure 5:  The considered (a) Floquet protocol, Eq. (12), with relevant
parameters $N=18,\epsilon=0.08$ and (b) random multipolar driving with
$N=16,T=1/25$, Eq. (13), exhibit heating to an infinite-temperature steady
state at easily accessible times. Reconstruction of the instantaneous
quasistatic effective Hamiltonian, $H_{\text{eff}}(n)$, responsible for
locally thermal time-dependent measurements in the heating regime, reveals a
growth of averaged normalized weight at less local terms with support
$\ell=3-5$ as a function of cycle number $n$, both, for (c) Floquet and (d)
random multipolar driving.

Heating regime.  We now turn our attention to the post-prethermal heating
regime; we use another Floquet protocol, where the length of the prethermal
plateau and the rate at which the system heats up are more easily tunable. To
this end, we follow Ref. [90] and consider the following protocol

$\displaystyle
U=\exp\left(-i\tilde{H}_{1}T/4\right)\,\exp\left(-i\tilde{H}_{2}\right)\,\exp\left(-i\tilde{H}_{1}T/4\right),$
$\displaystyle\tilde{H}_{1}=\sum_{i=1}^{N}J\sigma^{z}_{i}\sigma^{z}_{i+1}+h_{z}\sigma^{z}_{i}+h_{x}\sigma^{x}_{i},\
\tilde{H}_{2}=\epsilon\sum_{i=1}^{N}\sigma^{x}_{i}$ (12)

where $U$ is the unitary matrix for one drive cycle, which acts repeatedly on
the system during the dynamics; $T$ denotes the time period of the drive. As
argued in Ref. [90], the heating rate is regulated by $\epsilon$; the larger
the $\epsilon$, the faster the system heats up. Here we choose
$\epsilon=0.08$; Fig. 5(a), showing $\langle\tilde{H}_{1}/2\rangle$ at
different cycle numbers $n$, confirms that heating towards the infinite
temperature steady state is observed at easily accessible drive cycle numbers.
We note that the period $T$ is set to $T=2(\pi k+\epsilon)$, with $k=2$.

It has already been argued in Refs. [90, 91] that the system remains locally
describable with a Gibbs ensemble in the whole heating regime. Here, we
provide an independent check using the autoencoder analysis; it reveals that
the latent representation is one-dimensional and thus is parametrized by a
single parameter – the temperature – throughout the whole heating regime. To
confirm this, we take as initial states the thermal states with respect to the
zeroth order (in $\epsilon$) Floquet Hamiltonian $\tilde{H}_{1}/2$ with
$\beta_{0}\in[0.05,0.45]$; in addition, we perform weak, random,
translationally invariant single-site rotations so that the initial states
manifold is no longer one-dimensional. Then we evolve such states with Floquet
time evolution and probe their intrinsic dimension at chosen times. Namely, if
we considered initial states from the thermal manifold and mapped those with a
time-evolution deterministic operator to states at a fixed time, from the
perspective of an autoencoder, dataset remains on a one-dimensional manifold
since mapping due to time propagation can be absorbed in the encoder/decoder
part of the network. Therefore, to test the real, generic intrinsic dimension
of measurement in the heating regime, we must start from a manifold of initial
states with higher dimension that does not limit what we are probing. In Fig.
4(a) we show the ‘test error’, Eq.(1), as a function of the latent space
dimension $\mathcal{N}_{L}$ (i.e., number of neurons in the bottleneck). After
a few cycle numbers $n$, the ‘test error’ drops at one latent neuron and
flattens upon adding any further latent neurons, signaling that the latent
representation is one dimensional. This is confirmed also if we look directly
at the latent space representation of the dataset at a fixed time, cf. Fig.
4(b).

Next, we attempt to reconstruct the effective quasistatic Hamiltonian
$H_{\text{eff}}(n)$, responsible for thermal measurements in the heating
regime at different cycle numbers $n$, $\langle
O(\vec{\alpha})\rangle_{n}=\langle\psi(nT)|O(\vec{\alpha})|\psi(nT)\rangle\equiv\text{Tr}[O(\vec{\alpha})e^{-H_{\text{eff}}(n)}/Z]$.
This is challenging, yet interesting, because the inverse frequency expansions
do not converge in the heating regime. Note that for the reconstruction of
$H_{\text{eff}}(n)$, we resort to the iterative scheme described in Sec. II,
first selecting the leading 60 contributing terms with maximal support
$\mathcal{S}_{0}$ = 4 and then performing one iteration to increase the
maximal support to $S_{1}$ = 5, including 160 terms in total. The average
weight of all terms in the reconstructed $H_{\text{ML}}(n)$ with support
$\ell$, $\bar{a}_{\ell}$, relative to $\bar{a}_{1}$, are shown in Fig. 5(c)
plotted as a function of cycle number $n$. The relative growth of the weights
of the less local terms with increasing $n$ indicates that $H_{\text{eff}}(n)$
becomes less and less local as the system absorbs energy from the drive. On
the other hand, the Floquet Hamiltonian, defined from the unitary for one
period evolution, is a time-independent object. Our interpretation is that
$H_{\text{eff}}(n)$ does not necessarily capture the stroboscopic dynamics in
the system, but captures time-dependent thermal measurements. One way of
motivating the growing support of $H_{\text{eff}}(n)$ is in terms of operator
growth [92]. Let us say we initialized the system in a thermal state with
respect to $\tilde{H}_{1}$, $H_{\text{eff}}(0)=\tilde{H}_{1}$. Time evolution
of the corresponding thermal ensemble would be given by
$U^{n}e^{-\tilde{H}_{1}}(U^{\dagger})^{n}=e^{-U^{n}\tilde{H}_{1}(U^{\dagger})^{n}}=e^{-H_{\text{eff}}(n)}$,
with $U$ defined in Eq. (12). For parametrically small $\epsilon$, this should
naturally lead to the parametrically regulated growth of the
$H_{\text{eff}}(n)$ operator. Our $H_{\text{ML}}(n)$ reconstruction, shown in
Fig. 5(c), is consistent with an exponential growth with $n$.

The growing support of $H_{\text{eff}}(n)$ renders the reconstruction rather
difficult, since one needs to take into consideration the contributions from
terms with ever-growing support $\ell$.

###  IV.2 Aperiodic random multipolar drives

Periodically driven systems obey Floquet’s theorem, which guarantees the
existence of a static generator of stroboscopic dynamics. Away from the
strictly periodic limit, this additional structure of the time-ordered
propagator is lost, and the identification of an effective Hamiltonian becomes
far less obvious.

It was recently found that random multipolar driven (RMD) systems, which are
subject to an aperiodic energy-nonconserving drive, display similar prethermal
properties [93]: such systems exhibit a prethermal regime (albeit with a
weaker than exponential dependence of the duration of the prethermal plateau
on drive frequency), followed by an eventual relaxation to an infinite
temperature state. This behavior was proven mathematically by constructing an
effective static Hamiltonian for the prethermal plateau using a generalization
of the BCH expansion; like in Floquet systems, the expansion was shown to
diverge in the heating regime [53, 94]. To the best of our knowledge, it is
unclear whether an effective Hamiltonian exists, that describes the heating
dynamics of RMD systems beyond the prethermal plateau, nor what its locality
properties may be. This makes RMD systems an excellent candidate to test the
capabilities of our algorithm in so far uncharted territory.

To investigate further the phenomenon of heating beyond Floquet, we consider
many-body systems subject to a RMD drive. The protocol involves applying a
random sequence of unitaries $U_{\pm}$, generated from two Hamiltonians
$H_{\pm}$, each acting for a duration $T$. The level of multipolar correlation
incorporated in the drive (via the structure of the random sequence of
$U_{\pm}$) is characterised by an integer $m$. For instance, $m=0$ corresponds
to the case where random sequence of $U_{\pm}$ is considered, while $m=1$
represents the case for a random sequence of dipoles having form $U_{+}U_{-}$
or $U_{-}U_{+}$. Following the same logic, $m=2$ describes a random sequence
of quadrupoles having form $U_{-}U_{+}U_{+}U_{-}$ or $U_{+}U_{-}U_{-}U_{+}$.
It was demonstrated in [93] that RMDs, characterized by a finite value of $m$,
exhibit a prethermal regime with lifetime scaling as $(1/T)^{2m+1}$. For
$m\rightarrow\infty$, the limit itself corresponds to the Thue-Morse sequence
constructed out of $U_{\pm}$, leads to a sub-exponential (but faster than
polynomial) in $1/T$ long-lived prethermal regime.

We focus on the following Hamiltonian generators

$\displaystyle
H_{\pm}=\sum_{i}J_{x}\sigma_{i}^{x}\sigma_{i+1}^{x}+J_{z}\sigma_{i}^{z}\sigma_{i+1}^{z}+B_{\pm}\sigma_{i}^{x}+B_{z}\sigma_{i}^{z}$
(13)

where $B_{\pm}=B_{0}\pm B_{x}$. We construct the basic unitary operators
$U_{\pm}$ from $H_{\pm}$ as $U_{\pm}=\exp(-iH_{\pm}T)$.

Figure 5(b) shows the behaviour of the the mean energy
$H_{\text{RMD}}^{(0)}=(H_{+}+H_{-})/2$ for the system initialized in a thermal
state with respect to $H_{\text{RMD}}^{(0)}$. As we see, $\langle
H_{\text{RMD}}^{(0)}\rangle$ eventually reaches zero, indicating heating of
the system to an infinite temperature state at late times. In much the same
way as in the Floquet driving case, we can verify that the system remains
locally thermal throughout the whole heating regime. To the best of our
knowledge, this was not reported before.

In analogy to the Floquet case, we reconstruct the effective quasistatic
$H_{\text{eff}}(n)$ that captures the behavior of local stroboscopic thermal
expectation values in the heating regime. We use the iterative scheme
described in Sec. II, first selecting the leading 40 contributing terms with
maximal support $\mathcal{S}_{0}=4$ and then performing one further iteration
to increase the maximal support to $\mathcal{S}_{1}=5$, including 120 terms in
total. In Fig. 5(d) we show the average weights $\overline{a}_{\ell}$
(relative to $\overline{a}_{1}$) of terms with support $\ell$ for our
reconstructed $H_{\text{ML}}(n)$.

The above analysis shows that our algorithm is capable of identifying
effective Hamiltonians in new setups where no analytical theory is yet
available. Moreover, our results suggest the existence of simplified
descriptions for the heating dynamics beyond prethermal times, where BCH-like
expansions break down. Thus, the variational character of the Hamiltonian
reconstruction method we propose can prove useful both in experimental
nonequilibrium setups (e.g, to analyze dynamical decoupling sequences [29, 30,
31], or in the digital simulation of quantum dynamics [95, 96, 97]), and to
facilitate the theoretical analysis of nonequilibrium drives beyond Floquet
systems.

##  V Conclusions and Future Outlook

We propose a new algorithm for the Hamiltonian reconstruction from the
measurements of local operators. Reconstruction is possible under the
assumption that we have access to the local measurements in different thermal
states. We use autoencoder pre-processing of the data to (i) verify that the
measurements are thermal and (ii) single out the Hamiltonian candidate terms
in order to reduce the reconstruction complexity. The subsequent
reconstruction is accurate for local Hamiltonians as long as the dataset
contains measurement of all Hamiltonian terms (of course, without knowing,
which of the measured operators are the Hamiltonian terms), while it is only
approximate for long-range interacting Hamiltonians.

Beyond benchmarking, we apply Hamiltonian reconstruction to Floquet and random
multipolar driven systems. For the former, we first reconstruct the effective
Hamiltonian, responsible for the thermal local measurements in the prethermal
plateau, and compare that to its BCH approximation. Using an autoencoder, we
also confirm that the system remains locally thermal throughout the whole
heating regime, and reconstruct the effective quasistatic Hamiltonian,
corresponding to the Gibbs ensemble associated with the measurements at
different times. We interpreted this effective quasistatic Hamiltonian’s
increasing support via operator growth that is parametrically regulated by the
heating rate. Similar behaviour was observed also for the heating regime of
the random multipolar driving. Despite the lack of an analytical theory, we
confirmed that for the random multipolar driving, the heating process occurs
via thermal states with growing temperature and increasingly non-local
effective quasistatic Hamiltonian.

We envision our Hamiltonian reconstruction technique as a strong tool for the
verification of quantum simulators that strive for the realization of exotic
Hamiltonians. In our analysis, we focused on Floquet engineering, one of the
most popular nonequilibrium techniques to achieve this goal, for which,
however, Hamiltonians typically come with additional unwanted terms. Using
autoencoder analysis to select the Hamiltonian candidate terms, we demonstrate
how effective Hamiltonians can be reconstructed directly from the measurements
in the prethermal plateau or the heating regime, without relying on any
assumption on the Hamiltonian terms or the high-frequency expansion. Our
operator-based Hamiltonian reconstruction is directly suitable for the
verification of Hamiltonians in Rydberg atoms or trapped ion setups. However,
our approach would have to be extended to work with partial information stored
in state snapshots in cold atoms experiments. Hamiltonian reconstruction is
relevant also for condensed matter settings and quantum optics, where
extracting the minimal model, responsible for the observed physics, is
notoriously hard. Similarly, due to the limited types of measurement that are
accessible in these systems, the envisioned desirable generalizations of our
approach call for further studies.

Another non-trivial example where our scheme could be applied, includes
entanglement Hamiltonian reconstruction, at least in the intermediate to late-
time regime of excited systems, where entanglement Hamiltonian approaches the
actual Hamiltonian [98]. In future studies, it would also be desirable to
develop similar strategies to reconstruct other conservation laws of
integrable models that generically relax to the states that are locally
describable by the generalized Gibbs ensembles.

##  VI acknowledgements

It is our pleasure to acknowledge Wojciech De Roeck and Lev Vidmar for
stimulating discussions. Z.L. and S.N. acknowledge the support by the QuantERA
grant QuSiED, QuantERA grant T-NiSQ, J1-2463 project, and P1-0044 program of
the Slovenian Research Agency. M.S. was supported through the Helmholtz
Initiative and Networking Fund, grant no. VH-NG-1711. All the computations
were performed at the supercomputer Vega at the Institute of Information
Science (IZUM) in Maribor, Slovenia.

## References

  * Lewenstein  _et al._ [2007] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen(De), and U. Sen, Ultracold atomic gases in optical lattices: mimicking condensed matter physics and beyond, Advances in Physics 56, 243 (2007), https://doi.org/10.1080/00018730701223200 .
  * Jaksch and Zoller [2005] D. Jaksch and P. Zoller, The cold atom hubbard toolbox, Annals of Physics 315, 52 (2005), special Issue.
  * Blatt and Roos [2012] R. Blatt and C. F. Roos, Quantum simulations with trapped ions, Nature Physics 8, 277 (2012).
  * Schneider  _et al._ [2012] C. Schneider, D. Porras, and T. Schaetz, Experimental quantum simulations of many-body physics with trapped ions, Reports on Progress in Physics 75, 10.1088/0034-4885/75/2/024401 (2012).
  * Korenblit  _et al._ [2012] S. Korenblit, D. Kafri, W. C. Campbell, R. Islam, E. E. Edwards, Z. X. Gong, G. D. Lin, L. M. Duan, J. Kim, K. Kim, and C. Monroe, Quantum simulation of spin models on an arbitrary lattice with trapped ions, New Journal of Physics 14, 10.1088/1367-2630/14/9/095024 (2012), arXiv:1201.0776 .
  * Islam  _et al._ [2013] R. Islam, C. Senko, W. C. Campbell, S. Korenblit, J. Smith, A. Lee, E. E. Edwards, C. C. Wang, J. K. Freericks, and C. Monroe, Emergence and frustration of magnetism with variable-range interactions in a quantum simulator, Science 340, 583 (2013).
  * Simon  _et al._ [2011] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, and M. Greiner, Quantum simulation of antiferromagnetic spin chains in an optical lattice, Nature 472, 307 (2011), arXiv:1103.1372 .
  * Gambetta  _et al._ [2017] J. M. Gambetta, J. M. Chow, and M. Steffen, Building logical qubits in a superconducting quantum computing system, npj Quantum Information 3, 0 (2017), arXiv:1510.04375 .
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Jördens  _et al._ [2008] R. Jördens, N. Strohmaier, K. Günter, H. Moritz, and T. Esslinger, A Mott insulator of fermionic atoms in an optical lattice, Nature 455, 204 (2008), arXiv:0804.4009 .
  * Gross and Bloch [2017] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357, 995 (2017).
  * Mazurenko  _et al._ [2017] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, A cold-atom Fermi-Hubbard antiferromagnet, Nature 545, 462 (2017).
  * Zohar  _et al._ [2015] E. Zohar, J. I. Cirac, and B. Reznik, Quantum simulations of lattice gauge theories using ultracold atoms in optical lattices, Reports on Progress in Physics 79, 014401 (2015).
  * Bloch  _et al._ [2012] I. Bloch, J. Dalibard, and S. Nascimbène, Quantum simulations with ultracold quantum gases, Nature Physics 8, 267 (2012).
  * Aidelsburger  _et al._ [2013] M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Realization of the hofstadter hamiltonian with ultracold atoms in optical lattices, Phys. Rev. Lett. 111, 185301 (2013).
  * Bloch and Greiner [2005] I. Bloch and M. Greiner, Exploring quantum matter with ultracold atoms in optical lattices (Academic Press, 2005) pp. 1–47.
  * Feynman [1982] R. P. Feynman, Simulating physics with computers, International Journal of Theoretical Physics 21, 467 (1982).
  * Lloyd [1996] S. Lloyd, Universal quantum simulators, Science 273, 1073 (1996), https://www.science.org/doi/pdf/10.1126/science.273.5278.1073 .
  * Heras  _et al._ [2014] U. L. Heras, A. Mezzacapo, L. Lamata, S. Filipp, A. Wallraff, and E. Solano, Digital quantum simulation of spin systems in superconducting circuits, Phys. Rev. Lett. 112, 200501 (2014).
  * Barends  _et al._ [2015] R. Barends, L. Lamata, J. Kelly, L. García-Álvarez, A. G. Fowler, A. Megrant, E. Jeffrey, T. C. White, D. Sank, J. Y. Mutus, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, I. C. Hoi, C. Neill, P. J. O’Malley, C. Quintana, P. Roushan, A. Vainsencher, J. Wenner, E. Solano, and J. M. Martinis, Digital quantum simulation of fermionic models with a superconducting circuit, Nature Communications 6, 1 (2015), arXiv:1501.07703 .
  * Lamata  _et al._ [2018] L. Lamata, A. Parra-Rodriguez, M. Sanz, and E. Solano, Digital-analog quantum simulations with superconducting circuits, Advances in Physics: X 3, 519 (2018), arXiv:1711.09810 .
  * Weimer  _et al._ [2011] H. Weimer, M. Müller, H. P. Büchler, and I. Lesanovsky, Digital quantum simulation with Rydberg atoms, Quantum Information Processing 10, 885 (2011), arXiv:1104.3081 .
  * Müller  _et al._ [2012] M. Müller, S. Diehl, G. Pupillo, and P. Zoller, Engineered Open Systems and Quantum Simulations with Atoms and Ions, Advances in Atomic, Molecular and Optical Physics 61, 1 (2012), arXiv:1203.6595 .
  * Lanyon  _et al._ [2011] B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Gerritsma, F. Zähringer, P. Schindler, J. T. Barreiro, M. Rambach, G. Kirchmair, M. Hennrich, P. Zoller, R. Blatt, and C. F. Roos, Universal digital quantum simulation with trapped ions, Science 334, 57 (2011), https://www.science.org/doi/pdf/10.1126/science.1208001 .
  * Barreiro  _et al._ [2011] J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, An open-system quantum simulator with trapped ions, Nature 470, 486 (2011), arXiv:1104.1146 .
  * Goldman and Dalibard [2014] N. Goldman and J. Dalibard, Periodically driven quantum systems: Effective hamiltonians and engineered gauge fields, Phys. Rev. X 4, 031027 (2014).
  * Bukov  _et al._ [2015a] M. Bukov, L. D’Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to floquet engineering, Advances in Physics 64, 139 (2015a).
  * Eckardt [2017] A. Eckardt, Colloquium: Atomic quantum gases in periodically driven optical lattices, Rev. Mod. Phys. 89, 011004 (2017).
  * Viola and Lloyd [1998] L. Viola and S. Lloyd, Dynamical suppression of decoherence in two-state quantum systems, Phys. Rev. A 58, 2733 (1998).
  * Viola  _et al._ [1999] L. Viola, E. Knill, and S. Lloyd, Dynamical decoupling of open quantum systems, Phys. Rev. Lett. 82, 2417 (1999).
  * Bylander  _et al._ [2011] J. Bylander, S. Gustavsson, F. Yan, F. Yoshihara, K. Harrabi, G. Fitch, D. G. Cory, Y. Nakamura, J.-S. Tsai, and W. D. Oliver, Noise spectroscopy through dynamical decoupling with a superconducting flux qubit, Nature Physics 7, 565 (2011).
  * Choi  _et al._ [2020] J. Choi, H. Zhou, H. S. Knowles, R. Landig, S. Choi, and M. D. Lukin, Robust dynamic hamiltonian engineering of many-body spin systems, Phys. Rev. X 10, 031002 (2020).
  * Eckardt  _et al._ [2005] A. Eckardt, C. Weiss, and M. Holthaus, Superfluid-insulator transition in a periodically driven optical lattice, Phys. Rev. Lett. 95, 260404 (2005).
  * Zenesini  _et al._ [2009] A. Zenesini, H. Lignier, D. Ciampini, O. Morsch, and E. Arimondo, Coherent control of dressed matter waves, Phys. Rev. Lett. 102, 100403 (2009).
  * Oka and Aoki [2009] T. Oka and H. Aoki, Photovoltaic hall effect in graphene, Phys. Rev. B 79, 081406 (2009).
  * Kitagawa  _et al._ [2010] T. Kitagawa, E. Berg, M. Rudner, and E. Demler, Topological characterization of periodically driven quantum systems, Phys. Rev. B 82, 235114 (2010).
  * Jotzu  _et al._ [2014] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Experimental realization of the topological Haldane model with ultracold fermions, Nature 515, 237 (2014).
  * Lindner  _et al._ [2011] N. H. Lindner, G. Refael, and V. Galitski, Floquet topological insulator in semiconductor quantum wells, Nature Physics 7, 490 (2011), arXiv:1008.1792 .
  * Sun  _et al._ [2023] B.-Y. Sun, N. Goldman, M. Aidelsburger, and M. Bukov, Engineering and probing non-abelian chiral spin liquids using periodically driven ultracold atoms, PRX Quantum 4, 020329 (2023).
  * Aidelsburger  _et al._ [2011] M. Aidelsburger, M. Atala, S. Nascimbène, S. Trotzky, Y.-A. Chen, and I. Bloch, Experimental realization of strong effective magnetic fields in an optical lattice, Phys. Rev. Lett. 107, 255301 (2011).
  * Struck  _et al._ [2012] J. Struck, C. Ölschläger, M. Weinberg, P. Hauke, J. Simonet, A. Eckardt, M. Lewenstein, K. Sengstock, and P. Windpassinger, Tunable gauge potential for neutral and spinless particles in driven optical lattices, Phys. Rev. Lett. 108, 225304 (2012).
  * Bermudez  _et al._ [2011] A. Bermudez, T. Schaetz, and D. Porras, Synthetic gauge fields for vibrational excitations of trapped ions, Phys. Rev. Lett. 107, 150501 (2011).
  * Else  _et al._ [2020] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Discrete time crystals, Annual Review of Condensed Matter Physics 11, 467 (2020), https://doi.org/10.1146/annurev-conmatphys-031119-050658 .
  * Khemani  _et al._ [2019] V. Khemani, R. Moessner, and S. L. Sondhi, A Brief History of Time Crystals,  (2019), arXiv:1910.10745 .
  * Sacha and Zakrzewski [2017] K. Sacha and J. Zakrzewski, Time crystals: a review, Reports on Progress in Physics 81, 016401 (2017).
  * Mi  _et al._ [2022] X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, R. Babbush, J. C. Bardin, J. Basso, A. Bengtsson, A. Bilmes, A. Bourassa, L. Brill, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, B. Chiaro, R. Collins, W. Courtney, D. Debroy, S. Demura, A. R. Derk, A. Dunsworth, D. Eppens, C. Erickson, E. Farhi, A. G. Fowler, B. Foxen, C. Gidney, M. Giustina, M. P. Harrigan, S. D. Harrington, J. Hilton, A. Ho, S. Hong, T. Huang, A. Huff, W. J. Huggins, L. B. Ioffe, S. V. Isakov, J. Iveland, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, T. Khattar, S. Kim, A. Kitaev, P. V. Klimov, A. N. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, J. Lee, K. Lee, A. Locharla, E. Lucero, O. Martin, J. R. McClean, T. McCourt, M. McEwen, K. C. Miao, M. Mohseni, S. Montazeri, W. Mruczkiewicz, O. Naaman, M. Neeley, C. Neill, M. Newman, M. Y. Niu, T. E. O’Brien, A. Opremcak, E. Ostby, B. Pato, A. Petukhov, N. C. Rubin, D. Sank, K. J. Satzinger, V. Shvarts, Y. Su, D. Strain, M. Szalay, M. D. Trevithick, B. Villalonga, T. White, Z. J. Yao, P. Yeh, J. Yoo, A. Zalcman, H. Neven, S. Boixo, V. Smelyanskiy, A. Megrant, J. Kelly, Y. Chen, S. L. Sondhi, R. Moessner, K. Kechedzhi, V. Khemani, and P. Roushan, Time-crystalline eigenstate order on a quantum processor, Nature 601, 531 (2022).
  * Zhang  _et al._ [2017] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I. D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Observation of a discrete time crystal, Nature 543, 217 (2017), arXiv:1609.08684 .
  * Choi  _et al._ [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. Von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Observation of discrete time-crystalline order in a disordered dipolar many-body system, Nature 543, 221 (2017), arXiv:1610.08057 .
  * Gambetta  _et al._ [2019] F. M. Gambetta, F. Carollo, M. Marcuzzi, J. P. Garrahan, and I. Lesanovsky, Discrete time crystals in the absence of manifest symmetries or disorder in open quantum systems, Phys. Rev. Lett. 122, 015701 (2019).
  * Beatrez  _et al._ [2023] W. Beatrez, C. Fleckenstein, A. Pillai, E. de Leon Sanchez, A. Akkiraju, J. Diaz Alcala, S. Conti, P. Reshetikhin, E. Druga, M. Bukov, _et al._ , Critical prethermal discrete time crystal created by two-frequency driving, Nature Physics 19, 407 (2023).
  * Johansen  _et al._ [2022] C. H. Johansen, J. Lang, A. Morales, A. Baumgärtner, T. Donner, and F. Piazza, Multimode-polariton superradiance via floquet engineering, SciPost Physics 12, 094 (2022).
  * Mori  _et al._ [2016] T. Mori, T. Kuwahara, and K. Saito, Rigorous bound on energy absorption and generic relaxation in periodically driven quantum systems, Phys. Rev. Lett. 116, 120401 (2016).
  * Kuwahara  _et al._ [2016] T. Kuwahara, T. Mori, and K. Saito, Floquet–magnus theory and generic transient dynamics in periodically driven many-body quantum systems, Annals of Physics 367, 96 (2016).
  * Abanin  _et al._ [2017] D. A. Abanin, W. De Roeck, W. W. Ho, and F. m. c. Huveneers, Effective hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems, Phys. Rev. B 95, 014112 (2017).
  * Bukov  _et al._ [2015b] M. Bukov, S. Gopalakrishnan, M. Knap, and E. Demler, Prethermal floquet steady states and instabilities in the periodically driven, weakly interacting bose-hubbard model, Phys. Rev. Lett. 115, 205301 (2015b).
  * Morningstar  _et al._ [2022] A. Morningstar, D. A. Huse, and V. Khemani, Universality classes of thermalization for mesoscopic floquet systems, arXiv preprint arXiv:2210.13444 (2022).
  * Altman  _et al._ [2021] E. Altman, K. R. Brown, G. Carleo, L. D. Carr, E. Demler, C. Chin, B. DeMarco, S. E. Economou, M. A. Eriksson, K.-M. C. Fu, M. Greiner, K. R. Hazzard, R. G. Hulet, A. J. Kollár, B. L. Lev, M. D. Lukin, R. Ma, X. Mi, S. Misra, C. Monroe, K. Murch, Z. Nazario, K.-K. Ni, A. C. Potter, P. Roushan, M. Saffman, M. Schleier-Smith, I. Siddiqi, R. Simmonds, M. Singh, I. Spielman, K. Temme, D. S. Weiss, J. Vučković, V. Vuletić, J. Ye, and M. Zwierlein, Quantum simulators: Architectures and opportunities, PRX Quantum 2, 017003 (2021).
  * Daley  _et al._ [2022] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Practical quantum advantage in quantum simulation, Nature 607, 667 (2022).
  * Dür  _et al._ [2008] W. Dür, M. J. Bremner, and H. J. Briegel, Quantum simulation of interacting high-dimensional systems: The influence of noise, Phys. Rev. A 78, 052325 (2008).
  * Guerreschi  _et al._ [2020] G. G. Guerreschi, J. Hogaboam, F. Baruffa, and N. P. D. Sawaya, Intel quantum simulator: a cloud-ready high-performance simulator of quantum circuits, Quantum Science and Technology 5, 034007 (2020).
  * Qi and Ranard [2019] X. L. Qi and D. Ranard, Determining a local Hamiltonian from a single eigenstate, Quantum 3, 1 (2019), arXiv:1712.01850 .
  * Chertkov and Clark [2018] E. Chertkov and B. K. Clark, Computational inverse method for constructing spaces of quantum models from wave functions, Phys. Rev. X 8, 031029 (2018).
  * Greiter  _et al._ [2018] M. Greiter, V. Schnells, and R. Thomale, Method to identify parent hamiltonians for trial states, Phys. Rev. B 98, 081113 (2018).
  * Bairey  _et al._ [2019] E. Bairey, I. Arad, and N. H. Lindner, Learning a local hamiltonian from local measurements, Phys. Rev. Lett. 122, 020504 (2019).
  * Granade  _et al._ [2012] C. E. Granade, C. Ferrie, N. Wiebe, and D. G. Cory, Robust online hamiltonian learning, New Journal of Physics 14, 103013 (2012).
  * Wiebe  _et al._ [2014] N. Wiebe, C. Granade, C. Ferrie, and D. G. Cory, Hamiltonian learning and certification using quantum resources, Phys. Rev. Lett. 112, 190501 (2014).
  * Wang  _et al._ [2017] J. Wang, S. Paesani, R. Santagati, S. Knauer, A. A. Gentile, N. Wiebe, M. Petruzzella, J. L. O’brien, J. G. Rarity, A. Laing, and M. G. Thompson, Experimental quantum Hamiltonian learning, Nature Physics 13, 551 (2017), arXiv:1703.05402 .
  * Burgarth  _et al._ [2009] D. Burgarth, K. Maruyama, and F. Nori, Coupling strength estimation for spin chains despite restricted access, Phys. Rev. A 79, 020305 (2009).
  * Di Franco  _et al._ [2009] C. Di Franco, M. Paternostro, and M. S. Kim, Hamiltonian tomography in an access-limited setting without state initialization, Phys. Rev. Lett. 102, 187203 (2009).
  * Zhang and Sarovar [2014] J. Zhang and M. Sarovar, Quantum hamiltonian identification from measurement time traces, Phys. Rev. Lett. 113, 080401 (2014).
  * Kim  _et al._ [2016] J. Kim, J. Kwon, M. Kim, J. Do, D. Lee, and H. Han, Low-dielectric-constant polyimide aerogel composite films with low water uptake, Polymer Journal 48, 829 (2016).
  * Sone and Cappellaro [2017] A. Sone and P. Cappellaro, Hamiltonian identifiability assisted by a single-probe measurement, Phys. Rev. A 95, 022335 (2017).
  * Wang  _et al._ [2018] Y. Wang, D. Dong, B. Qi, J. Zhang, I. R. Petersen, and H. Yonezawa, A quantum hamiltonian identification algorithm: Computational complexity and error analysis, IEEE Transactions on Automatic Control 63, 1388 (2018).
  * Li  _et al._ [2020] Z. Li, L. Zou, and T. H. Hsieh, Hamiltonian Tomography via Quantum Quench, Physical Review Letters 124, 160502 (2020), arXiv:1912.09492 .
  * Rudinger and Joynt [2015] K. Rudinger and R. Joynt, Compressed sensing for hamiltonian reconstruction, Phys. Rev. A 92, 052322 (2015).
  * Kieferová and Wiebe [2017] M. Kieferová and N. Wiebe, Tomography and generative training with quantum boltzmann machines, Phys. Rev. A 96, 062327 (2017).
  * Valenti  _et al._ [2019] A. Valenti, E. Van Nieuwenburg, S. Huber, and E. Greplova, Hamiltonian learning for quantum error correction, Physical Review Research 1, 1 (2019), arXiv:1907.02540 .
  * Pastori  _et al._ [2022] L. Pastori, T. Olsacher, C. Kokail, and P. Zoller, Characterization and verification of trotterized digital quantum simulation via hamiltonian and liouvillian learning, PRX Quantum 3, 030324 (2022).
  * Olsacher  _et al._ [2022] T. Olsacher, L. Pastori, C. Kokail, L. M. Sieberer, and P. Zoller, Digital quantum simulation, learning of the floquet hamiltonian, and quantum chaos of the kicked top, Journal of Physics A: Mathematical and Theoretical 55, 334003 (2022).
  * Bukov  _et al._ [2016] M. Bukov, M. Heyl, D. A. Huse, and A. Polkovnikov, Heating and many-body resonances in a periodically driven two-band system, Phys. Rev. B 93, 155132 (2016).
  * Bourlard and Kamp [2004] H. Bourlard and Y. Kamp, Auto-association by multilayer perceptrons and singular value decomposjournaition, Biological Cybernetics 59, 291 (2004).
  * Tschannen  _et al._ [2018] M. Tschannen, O. Bachem, and M. Lucic, Recent Advances in Autoencoder-Based Representation Learning,  , 1 (2018), arXiv:1812.05069 .
  * Schmitt and Lenarčič [2022] M. Schmitt and Z. Lenarčič, From observations to complexity of quantum states via unsupervised learning, Phys. Rev. B 106, L041110 (2022).
  * Lazarides  _et al._ [2014a] A. Lazarides, A. Das, and R. Moessner, Equilibrium states of generic quantum systems subject to periodic driving, Phys. Rev. E 90, 012110 (2014a).
  * Lazarides  _et al._ [2014b] A. Lazarides, A. Das, and R. Moessner, Equilibrium states of generic quantum systems subject to periodic driving, Phys. Rev. E 90, 012110 (2014b).
  * D’Alessio and Rigol [2014] L. D’Alessio and M. Rigol, Long-time behavior of isolated periodically driven interacting lattice systems, Phys. Rev. X 4, 041048 (2014).
  * Weidinger and Knap [2017] S. A. Weidinger and M. Knap, Floquet prethermalization and regimes of heating in a periodically driven, interacting quantum system, Scientific Reports 7, 1 (2017), arXiv:1609.09089 .
  * Ho  _et al._ [2023] W. W. Ho, T. Mori, D. A. Abanin, and E. G. Dalla Torre, Quantum and classical floquet prethermalization, Annals of Physics , 169297 (2023).
  * Bukov  _et al._ [2015c] M. Bukov, L. D’Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to floquet engineering, Advances in Physics 64, 139 (2015c), https://doi.org/10.1080/00018732.2015.1055918 .
  * Fleckenstein and Bukov [2021] C. Fleckenstein and M. Bukov, Thermalization and prethermalization in periodically kicked quantum spin chains, Phys. Rev. B 103, 144307 (2021).
  * Ikeda and Polkovnikov [2021] T. N. Ikeda and A. Polkovnikov, Fermi’s golden rule for heating in strongly driven floquet systems, Phys. Rev. B 104, 134308 (2021).
  * Parker  _et al._ [2019] D. E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, and E. Altman, A universal operator growth hypothesis, Phys. Rev. X 9, 041017 (2019).
  * Zhao  _et al._ [2021] H. Zhao, F. Mintert, R. Moessner, and J. Knolle, Random multipolar driving: Tunably slow heating through spectral engineering, Phys. Rev. Lett. 126, 040601 (2021).
  * Mori  _et al._ [2021] T. Mori, H. Zhao, F. Mintert, J. Knolle, and R. Moessner, Rigorous bounds on the heating rate in thue-morse quasiperiodically and randomly driven quantum many-body systems, Phys. Rev. Lett. 127, 050602 (2021).
  * Heyl  _et al._ [2019] M. Heyl, P. Hauke, and P. Zoller, Quantum localization bounds trotter errors in digital quantum simulation, Science advances 5, eaau8342 (2019).
  * Zhao  _et al._ [2022] H. Zhao, M. Bukov, M. Heyl, and R. Moessner, Making trotterization adaptive for nisq devices and beyond, arXiv preprint arXiv:2209.12653 (2022).
  * Zhao  _et al._ [2023] H. Zhao, M. Bukov, M. Heyl, and R. Moessner, Adaptive trotterization for time-dependent hamiltonian quantum dynamics using instantaneous conservation laws, arXiv preprint arXiv:2307.10327 (2023).
  * Dalmonte  _et al._ [2022] M. Dalmonte, V. Eisler, M. Falconi, and B. Vermersch, Entanglement Hamiltonians: From Field Theory to Lattice Models and Experiments, Annalen der Physik 534, 1 (2022), arXiv:2202.05045 .
  * Bartsch and Gemmer [2009] C. Bartsch and J. Gemmer, Dynamical typicality of quantum expectation values, Phys. Rev. Lett. 102, 110403 (2009).
  * Reimann [2018] P. Reimann, Dynamical typicality of isolated many-body quantum systems, Phys. Rev. E 97, 062129 (2018).
  * Reimann and Dabelow [2019] P. Reimann and L. Dabelow, Typicality of prethermalization, Phys. Rev. Lett. 122, 080603 (2019).
  * Richter and Steinigeweg [2019] J. Richter and R. Steinigeweg, Combining dynamical quantum typicality and numerical linked cluster expansions, Phys. Rev. B 99, 094419 (2019).

##  Appendix A Details of the autoencoder

Our autoencoder consists of feed-forward neural network with the bottleneck
structure. A feed-forward neural network is made up of layers, corresponding
to a sequence of alternating affine-linear and non-linear transformations. The
activation $a_{j}^{(l)}$ of the $l-$th layer is related to that of the
previous layer $a_{j}^{(l-1)}$ as

$a_{j}^{l}=\sigma\left(\sum_{k}W_{jk}a_{k}^{(l-1)}+b_{j}\right)$ (14)

where $W_{jk}$ and $b_{j}$ are variational parameters and $\sigma$ denotes the
fixed non-linear activate function. For all the results produced in this
paper, we fix $\sigma=\tanh$. Our network consists of two encoder layers
between the input and the latent space, as well as two decoder layers between
the latent space and the output. Each layer has 400 neurons. We used the Adam
optimizer for training the network and learning rate $r=0.0001$. 80% of the
data serves as the training set, while rest are used as test set.

##  Appendix B Preparing the initial thermal state

To prepare the initial thermal state for $N>12$, we follow the prescription
detailed in Ref. [90]. The procedure relies on the principle of quantum
typicality [99, 100, 101, 102]. Typicality states that the expectation value
of an operator $\mathcal{A}$, defined on a Hilbert space $\mathcal{H}$, can be
approximated as

$\frac{1}{\text{dim}\mathcal{H}}\text{Tr}[A]\approx\frac{1}{N}\sum_{n=1}^{N}\langle
r_{n}|A|r_{n}\rangle$ (15)

where $|r_{n}\rangle$ are Haar-random states. Note that the approximation
becomes exact as $N\rightarrow\infty$. Therefore, thermal expectation value of
observable $A$, with respect to some Hamiltonian $H$ at inverse temperature
$\beta$, is given as

$\langle
A\rangle_{\beta}\approx\frac{\frac{\text{dim}(\mathcal{H})}{N}\sum_{n=1}^{N}\langle
r_{n}|e^{-\beta H/2}Ae^{-\beta
H/2}|r_{n}\rangle}{\frac{\text{dim}(\mathcal{H})}{N}\sum_{n=1}^{N}\langle
r_{n}|e^{-\beta H/2}e^{-\beta H/2}|r_{n}\rangle}$ (16)

If we interpret the RHS of the above equation as an ensemble average, then the
thermal density matrix can be approximated as

$\displaystyle\rho_{\beta}$
$\displaystyle\approx\frac{1}{Z_{\beta}}\frac{\text{dim}(\mathcal{H})}{N}\sum_{n=1}^{N}e^{-\beta
H/2}|r_{n}\rangle\langle r_{n}|e^{-\beta H/2}$ (17) $\displaystyle Z_{\beta}$
$\displaystyle\approx\frac{\text{dim}(\mathcal{H})}{N}\sum_{n}^{N}\langle
r_{n}|e^{-\beta H/2}e^{-\beta H/2}|r_{n}\rangle$

To time evolve the thermal density matrix, we just need to evolve a set of
pure states defined as $|\psi_{n}\rangle=e^{-\beta H/2}$ and then compute

$\displaystyle\rho_{\beta}(t)\approx\frac{1}{Z_{\beta}}\frac{\text{dim}(\mathcal{H})}{N}\sum_{n=1}^{N}U(t)|\psi_{n}\rangle\langle\psi_{n}|[U(t)]^{\dagger}$
(18)

where $U(t)$ is the evolution operator.

##  Appendix C BCH Expansion

First of all, we note that for that square pulse Floquet protocol used in the
main text, Floquet-Magnus expansion is identical to the BCH expansion.
Mathematically, BCH formula is the solution $Z$ to the equation

$\displaystyle e^{X}e^{Y}=e^{Z}$ (19)

for possibly noncommutative $X$ and $Y$ in the Lie algebra of a Lie group. BCH
formula essentially yields an expression for $Z$ in terms of a formal series
(not necessarily convergent) in $X$, $Y$, and iterated commutators thereof.
The first few terms of the series are:

$\displaystyle
Z=X+Y+\frac{1}{2}[X,Y]+\frac{1}{12}([X,[X,Y]]-[Y,[X,Y]])+......$

To make connection with our problem, we notice that the Floquet Hamiltonian
$H_{F}$ is defined via the relation

$\displaystyle e^{-iH_{1}T/2}e^{-iH_{2}T/2}=e^{-iH_{F}T}$ (21)

It is quite clear that Eq. (21) is same as Eq. (C) with $X=-iH_{1}T/2$ and
$Y=-iH_{2}T/2$, where $H_{1}$ and $H_{2}$ are described by Eq. (10) and
Eq.(11) respectively, in the main text. Now, resorting to Eq. (C), one can
approximate the Floquet Hamiltonian $H_{F}$ in different orders of $T$, such
that $H_{F}=\sum_{n}H_{F}^{(n)}$, where $H_{F}^{(n)}$ is proportional to
$T^{n}$. We perform such calculation up to order 2 and the results are charted
down below:

$\displaystyle H_{F}^{(0)}=$ $\displaystyle\frac{1}{2}(H_{1}+H_{2})$ (22)
$\displaystyle H_{F}^{(1)}=$ $\displaystyle\frac{\gamma
T}{4}\sum_{i}\sigma_{i}^{y}\sigma_{i+1}^{z}+\sigma_{i}^{z}\sigma_{i+1}^{y}+h_{z}\sigma_{i}^{y}$
(23) $\displaystyle H_{F}^{(2)}=$
$\displaystyle-\frac{T^{2}}{24}\gamma^{2}h_{z}\sum_{i}\sigma_{i}^{z}-\frac{T^{2}}{24}\gamma
h_{z}^{2}\sum_{i}\sigma_{i}^{x}$ (24)
$\displaystyle-\frac{T^{2}}{12}\gamma^{2}\sum_{i}(\sigma_{i}^{z}\sigma_{i+1}^{z}-\sigma_{i}^{y}\sigma_{i+1}^{y})$
$\displaystyle-\frac{T^{2}}{24}\gamma
h_{z}\sum_{i}(\sigma_{i}^{x}\sigma_{i+1}^{z}+\sigma_{i}^{z}\sigma_{i+1}^{x})$
$\displaystyle+\frac{T^{2}}{12}\gamma
h_{x}\sum_{i}(\sigma_{i}^{z}\sigma_{i+1}^{z}-\sigma_{i}^{y}\sigma_{i+1}^{y})$
$\displaystyle-\frac{T^{2}}{24}\gamma
h_{z}\sum_{i}(\sigma_{i}^{z}\sigma_{i+1}^{x}+\sigma_{i}^{x}\sigma_{i+1}^{z})$
$\displaystyle-\frac{T^{2}}{12}\gamma\sum_{i}(\sigma_{i}^{z}\sigma_{i+1}^{z}\sigma_{i+2}^{z}+\sigma_{i}^{x})$