# Adaptive Trotterization for time-dependent Hamiltonian quantum dynamics
using piecewise conservation laws

Hongzheng Zhao<EMAIL_ADDRESS>School of Physics, Peking University, 100871
Beijing, China  Max Planck Institute for the Physics of Complex Systems,
Nöthnitzer Str. 38, 01187 Dresden, Germany  Institute for Quantum Optics and
Quantum Information, Technikerstraße 21a 6020 Innsbruck, Austria  Marin Bukov
Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38,
01187 Dresden, Germany  Markus Heyl  Max Planck Institute for the Physics of
Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany  Theoretical
Physics III, Center for Electronic Correlations and Magnetism, Institute of
Physics, University of Augsburg, 86135 Augsburg, Germany  Roderich Moessner
Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38,
01187 Dresden, Germany

###### Abstract

Digital quantum simulation relies on Trotterization to discretize time
evolution into elementary quantum gates. On current quantum processors with
notable gate imperfections, there is a critical tradeoff between improved
accuracy for finer timesteps, and increased error rate on account of the
larger circuit depth. We present an adaptive Trotterization algorithm to cope
with time-dependent Hamiltonians, where we propose a concept of piecewise
“conserved” quantities to estimate errors in the time evolution between two
(nearby) points in time; these allow us to bound the errors accumulated over
the full simulation period. They reduce to standard conservation laws in the
case of time-independent Hamiltonians, for which we first developed an
adaptive Trotterization scheme [PRX Quantum 4, 030319]. We validate the
algorithm for a time-dependent quantum spin chain, demonstrating that it can
outperform the conventional Trotter algorithm with a fixed step size at a
controlled error.

Introduction.— Simulating the time evolution of non-equilibrium quantum many-
body systems poses a significant challenge for classical computers due to the
exponentially large Hilbert space dimension [1, 2]. The rapid development of
quantum processors holds the promise of resolving this key problem through
digital quantum simulation (DQS) [3, 4, 5, 6, 7, 8, 9, 10, 11].

In DQS, the continuous time evolution operator is discretized into elementary
few-body quantum gates, a procedure known as Trotterization [12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. However, due to their
noncommutativity, Trotterization introduces errors, which can accumulate over
longer simulation times. While a finer Trotter time step size $\delta t$
improves simulation precision, it also leads to increased circuit depth. In
the current era of noisy intermediate-scale quantum (NISQ) processors, gate
imperfections are inevitable, posing a significant challenge in improving the
accuracy of DQS [2], especially in the absence of experimentally efficient
error-correction schemes [30, 31, 32]. Therefore, it is crucial to identify
strategies for minimizing circuit depth while keeping the simulation error
under control.

In previous work [33], we introduced a quantum algorithm, ADA-Trotter,
allowing for adaptive step sizes $\delta t$ to optimize the usage of quantum
gates for time-independent Hamiltonians. By measuring the expectation values
of energy and energy variance, $\delta t$ is maximized as long as errors in
these conserved quantities remain bounded. According to the central limit
theorem, ADA-Trotter ensures a correct energy distribution for generic non-
integrable many-body systems. However, extending this formalism to time-
dependent Hamiltonians $H(t)$ is a demanding challenge, since: (i) Energy
conservation is absent and hence it is a priori unclear how to define a
criterion to adapt $\delta t$; (ii) without a static reference Hamiltonian the
energy distribution is difficult to define, and hence the implication of the
central limit theorem, mentioned above, is now elusive; (iii) generic many-
body systems absorb energy from time-dependent modulations and may heat up (in
the sense of the Eigenstate Thermalization Hypothesis [34, 35]) and eventually
approach states with trivial correlations [36, 37, 38]. It is not clear how to
control the additional heating generated by the piecewise constant time-
dependence of a Trotterized Hamiltonian compared to that of $H(t)$.

In this work, we propose tADA-Trotter – an adaptive algorithm for time-
dependent quantum systems [39, 40, 41, 42, 43, 44, 45, 46]. To achieve this,
we first discretize the time evolution into small time intervals
$[t,t{+}\delta t]$. In each interval, the time evolution can be generated by
an effective Hamiltonian $H_{[\infty]}^{t,\delta t}$. Note that this
Hamiltonian has only a parametric dependence on $t,\delta t$; it is time-
independent once $t$ and $\delta t$ are fixed. Consequently, the expectation
values of $H_{[\infty]}^{t,\delta t}$ and its higher moments coincide at the
boundaries of this time window for perfect time evolution, a feature we will
refer to as piecewise conservation laws, cf. Fig. 1 (a). In practice,
$H_{[\infty]}^{t,\delta t}$ can be approximated using a perturbative Magnus
expansion in the small time step $\delta t$ [47, 48]. The Magnus and
Trotterization approximations introduce errors in the piecewise conservation
laws. Our key finding is that, by constraining these errors, the step size can
be adapted to reduce circuit depth while maintaining a given simulation
accuracy for generic non-integrable many-body systems.

Figure 1: Schematics of tADA-Trotter for a time-dependent Hamiltonian. (a) The
expectation values of the piecewise conserved Hamiltonian
$H_{[\infty]}^{t,\delta t}$ coincide at times $t$ and $t+\delta t$. The exact
expectation value of $H(t)$ is depicted as black line. We use a Magnus
expansion to approximate this conservation law $H_{[\infty]}^{t,\delta
t}{\approx}H_{[k]}^{t,\delta t}$(blue) and a Trotter decomposition (green) for
the time evolution. Both approximations introduce errors to the time
evolution. We maximize $\delta t$ as long as errors in this conservation law
are bounded, i.e., deviations in the expectation value before
($\mathcal{E}_{i}$) and after ($\mathcal{E}_{f}$) the Trotterized evolution
(green) should be small. (b) Schematic depiction of state space, showing the
exact time evolution (black line) as well as the approximate state vector
evolution calculated using tADA-Trotter under either time-local (gray/black
circles) or time-global (orange/red circles) constraints.Trotter errors
accumulate with the time-local constraint (grey),
$\left|\mathcal{E}_{f}\left(t_{n},\delta
t_{n}\right)-\mathcal{E}_{i}\left(t_{n},\delta
t_{n}\right)\right|<d_{\mathcal{E}}$. This can be suppressed by using the
time-global (accumulated) constraint (red),
$\left|\sum_{n}\left[\mathcal{E}_{f}\left(t_{n},\delta
t_{n}\right)-\mathcal{E}_{i}\left(t_{n},\delta
t_{n}\right)\right]\right|<d^{\prime}_{\mathcal{E}}$, leading to more reliable
simulation with adaptive step sizes.

The adaptive step size is determined by measuring piecewise conservation laws
using Trotterized time evolution. Subsequently, we introduce the concept of a
time-global error, which represents the accumulation of time-local errors in
the piecewise conservation laws, cf. Fig. 1. To adapt the step size, we
propose a feedback procedure that can bound both the local and global errors.

Notably, Trotter-induced heating effects present in a local control scheme,
where simulation errors significantly accumulate over time by making
suboptimal choices of $\delta t$, may be efficiently suppressed by imposing
constraints on the global errors, cf. Fig. 1 (b). For time-independent
systems, this globally constrained scheme reduces to the algorithm proposed in
Ref. [33], enabling strict error bounds throughout the entire time evolution.
To determine the advantages of a globally constrained error, we perform
numerical simulations of a quantum spin chain with a time-dependent field. We
also demonstrate that tADA-Trotter outperforms the conventional fixed-step
Trotter, as depicted in Fig. 4. These findings highlight the superior
potential of tADA-Trotter with a global constraint in minimizing circuit depth
for DQS of time-dependent systems.

Piecewise conservation laws.—The time evolution operator $U(t,t^{\prime})$
follows the equation $\partial_{t}U(t,t^{\prime}){=}-iH(t)U(t,t^{\prime})$,
where $H(t)$ represents the time-dependent Hamiltonian. Its solution is given
by the time-ordered exponential $U(t+\delta
t,t){=}\mathcal{T}\exp\left(-i\int_{t}^{t{+}\delta t}H(s)\mathrm{d}s\right)$,
and the exact state evolves as $\left|\phi(t{+}\delta
t)\right\rangle{=}U(t{+}\delta t,t)\left|\phi(t)\right\rangle$.

Formally rewriting $U(t{+}\delta t,t){=}\exp\left(-iH_{[\infty]}\delta
t\right)$ indicates that the same time evolution can be generated by the
static effective Hamiltonian $H_{[\infty]}$, where we drop its parametric
dependence on $t,\delta t$ for simplicity. Hence, when $t$ and $\delta t$ are
fixed, the expectation value of $H_{[\infty]}$, and its higher-order moments,
coincide for the states $\left|\phi(t{+}\delta t)\right\rangle$ and
$\left|\phi(t)\right\rangle$. We use these piecewise conservation laws to
adapt the Trotter step size $\delta t$.

The piecewise conserved Hamiltonian can be obtained through a Magnus expansion
given by $H_{[\infty]}=i\delta t^{-1}\sum_{n=1}^{\infty}\Omega_{n}$, where the
operator $\Omega_{n}\propto\delta t^{n}$. The explicit form of $H_{[\infty]}$
can be complicated as higher-order contributions typically involve nested
commutators. To eliminate the time-ordered integral in the Magnus expansion,
we expand the time-dependence in Legendre polynomials, obtaining the concise
expression for terms of lowest orders [48]: $\Omega_{2m}{=}0$ for all even
orders $2m$, and

$\displaystyle\begin{aligned} &\Omega_{1}{=}A_{1},\
\Omega_{3}{=}-\frac{1}{6}\left[A_{1},A_{2}\right],\Omega_{5}{=}\frac{1}{60}\left[A_{1},\left[A_{1},A_{3}\right]\right]{-}\\\
&\frac{1}{60}\left[A_{2},\left[A_{1},A_{2}\right]\right]{+}\frac{1}{360}\left[A_{1},\left[A_{1},\left[A_{1},A_{2}\right]\right]\right]{-}\frac{1}{30}\left[A_{2},A_{3}\right],\end{aligned}$
(1)

where each operator $A_{n}^{t,\delta t}$ is defined as $A_{n}^{t,\delta
t}=-i(2n{-}1)\delta t\int_{0}^{1}H(t{+}x\delta t)P_{n-1}(x)dx$. Here,
$P_{n-1}$ denotes the shifted Legendre polynomials normalized to
$(2n+1)\int_{0}^{1}dsP_{m}(s)P_{n}(s){=}\delta_{mn}$, and $A_{n}^{t,\delta
t}\propto\delta t^{n}$. For a sufficiently small time interval $\delta t$,
this Magnus expansion can be truncated at a finite order $k$, resulting in an
approximation of $H_{[\infty]}$ as $H_{[k]}{=}i\delta
t^{-1}\sum_{n=1}^{k}\Omega_{n}$.

Trotterization.— On real digital quantum devices, the exact time evolution
operator $U(t{+}\delta t,t)$ for a smoothly varying Hamiltonian $H(t)$ is
usually inaccessible. Thus, one needs to decompose the former into some
elementary quantum gates using, e.g., Trotterization. For simplicity, we focus
on the time-dependence $H(t){=}g(t)G{+}f(t)F$ with smooth functions $g(t)$ and
$f(t)$, and two generic non-commuting hermitian operators $F$ and $G$. Let us
assume that quantum devices admit the exact implementation of unitaries of the
form $\exp(-iC_{1}G)$ or $\exp(-iC_{2}F)$ where $C_{1,2}$ are arbitrary real
numbers; while the implementation of linear combinations of $G$ and $F$, such
as $\exp(-iC_{1}G{-}iC_{2}F)$, is not feasible.

We aim to approximate the target unitary operator up to a given order
$\lambda$, such that $U(t{+}\delta t,t)=U_{[\lambda]}(t{+}\delta
t,t){+}\mathcal{O}(\delta t^{\lambda}).$ A larger $\lambda$ leads to more
accurate time evolution with smaller Trotter errors, but it also increases the
circuit depth. The number of gates generally scales exponentially in
$\lambda$, but better decompositions with fewer number of exponentials may
exist [43]. In this work, we use $\lambda{=}3$ and the approximation can be
obtained by the second-order Trotter formula, also known as the mid-point
rule:

$\displaystyle U_{[3]}(t{+}\delta t,t){=}\exp[{-ig(t{+}\delta t/2)G\delta
t/2}]\times$ (2) $\displaystyle\exp[{-if(t{+}\delta t/2)F\delta
t}]\exp[{-ig(t{+}\delta t/2)G\delta t/2}].$

Adaptive algorithm.— The central concept of the adaptive algorithm is to
maximize $\delta t$ while ensuring that the measurement outcome of the
expectation value and variance of $H_{[\infty]}$ remain within pre-set
tolerances. However, $H_{[\infty]}$ is generically increasingly non-local as
contributions of increasing order $\Omega_{n}$ are involved, introducing a
significant measurement overhead. Therefore, depending on the measurement
accuracy and efficiency, one may need to truncate $H_{[\infty]}$ to a finite
order $k$ to make the measurement procedure feasible on quantum computers [49,
50, 51]. Here, we consider a sufficiently large value for $k$, ensuring that
errors in the piecewise conservation law are subdominant compared to the
Trotter error. In the limit $\delta t\to 0$, this can be satisfied as long as
$k\geq\lambda$; we prove it using perturbation theory, see Sec. SM 1 for
details.

Below, we first introduce a time-local control scheme to adapt $\delta t$,
which we find can involve severe heating effects. Then we propose a global
control scheme to bound the accumulated errors and suppress heating.

In contrast to time-independent systems where conserved energies solely depend
on the initial state $|\psi(0)\rangle$, the expectation values of $H_{[k]}$
for generic time-dependent systems rely on the time-evolved state
$|\psi(t)\rangle$. As a result, there are no universal (in time) reference
expectation values known a priori for the piecewise conserved Hamiltonians.
Nonetheless, we can leverage the capability of quantum processors to measure
the expectation values, $\mathcal{E}_{i}(t,\delta
t){=}\langle\psi(t)|H_{[k]}|\psi(t)\rangle/L$ with the system size $L$, as a
reasonable approximation to the conserved quantities for the true quantum
state $|\phi(t)\rangle$. We maximize $\delta t$ such that the time-local error
in the expectation value of $H_{[k]}$ remains below a threshold
$d_{\mathcal{E}}$ (cf. Fig. 1(a)), i.e.,

$|\mathcal{E}_{f}(t,\delta t)-\mathcal{E}_{i}(t,\delta t)|<d_{\mathcal{E}},$
(3)

where $\mathcal{E}_{f}(t,\delta t){=}\langle\psi(t{+}\delta
t)|{H_{[k]}}|\psi(t{+}\delta t)\rangle/L$ represents the expectation value of
$H_{[k]}$ after the Trotterized evolution, given by $|\psi(t{+}\delta
t)\rangle{=}U_{[\lambda]}(t,\delta t)|\psi(t)\rangle$. In the ideal case of
$\lambda{\to}\infty$ and $k{\to}\infty$, we have $\mathcal{E}_{f}(t,\delta
t){-}\mathcal{E}_{i}(t,\delta t){=}0$ by definition. However, for any finite
value of $\lambda$ and $k$, this error does not vanish.

In addition, we also require the error in the variance to be bounded by the
tolerance $d_{\delta\mathcal{E}^{2}}$, i.e.,
$|\delta\mathcal{E}^{2}_{f}(t,\delta t){-}\delta\mathcal{E}^{2}_{i}(t,\delta
t)|{<}d_{\delta\mathcal{E}^{2}}$, with

$\displaystyle\begin{aligned} &\delta\mathcal{E}^{2}_{i}(t,\delta
t){=}L^{-1}\langle\psi(t)|H_{[k]}^{2}(t,\delta
t)|\psi(t)\rangle{-}L\mathcal{E}^{2}_{i},\\\
&\delta\mathcal{E}^{2}_{f}(t,\delta t){=}L^{-1}\langle\psi(t+\delta
t)|H_{[k]}^{2}(t,\delta t)|\psi(t+\delta
t)\rangle{-}L\mathcal{E}^{2}_{f}.\end{aligned}$ (4)

According to the central limit theorem, constraining the errors in the lowest
two moments of $H_{[k]}$ is sufficient to ensure the approximate conservation
of its higher moments [52, 33]. Therefore, the piecewise conservation of the
Hamiltonian $H_{[k]}$ can be satisfied reasonably well, enabling reliable DQS
of dynamics from time $t$ to $t{+}\delta t$.

Note that Trotter errors can accumulate in the time-evolved $|\psi(t)\rangle$,
leading to deviations of $\mathcal{E}_{i}$ and $\delta\mathcal{E}_{i}^{2}$,
from the exact piecewise conservation laws. This effect is also present for
time-independent systems, where $H_{[k]}$ simplifies to a static Hamiltonian
$H$. The energy constraint reduces to
$|\langle\psi(t)|H|\psi(t)\rangle-\langle\psi(t{+}\delta t)|H|\psi(t{+}\delta
t)\rangle|/L{<}d_{\mathcal{E}},$ which accumulates and cannot be bounded for
long simulation times. Consequently, many-body systems tend to heat up and
eventually become featureless. This Trotter-induced heating can be much more
pronounced in time-dependent systems, leading to unstable DQS of time
evolution over long periods.

To address this challenge, we propose restrictions on the time-global errors,
representing the accumulation of all time-local errors from previous steps:

$\displaystyle\left|\sum\nolimits_{n=1}^{m}[\mathcal{E}_{f}(t_{n},\delta
t_{n})-\mathcal{E}_{i}(t_{n},\delta t_{n})]\right|<d^{\prime}_{\mathcal{E}},$
(5)
$\displaystyle\left|\sum\nolimits_{n=1}^{m}[\delta\mathcal{E}^{2}_{f}(t_{n},\delta
t_{n})-\delta\mathcal{E}^{2}_{i}(t_{n},\delta
t_{n})]\right|<d^{\prime}_{\delta\mathcal{E}^{2}}.$

These conditions imply the time-local constraints, e.g.,
$|\mathcal{E}_{f}(t_{m},\delta t_{m}){-}\mathcal{E}_{i}(t_{m},\delta
t_{m})|{<}2d^{\prime}_{\mathcal{E}}$ 111To show this, we define
$\Delta_{m-1}{=}\sum\nolimits_{n=1}^{m-1}\mathcal{E}_{f}(t_{n},\delta
t_{n})-\mathcal{E}_{i}(t_{n},\delta t_{n}),$ which satisfies
$-d^{\prime}_{\mathcal{E}}{<}\Delta_{m-1}{<}d^{\prime}_{\mathcal{E}}$ with a
positive $d^{\prime}_{\mathcal{E}}$. For the next step, one has
$-d^{\prime}_{\mathcal{E}}<\Delta_{m-1}+\mathcal{E}_{f}(t_{m},\delta
t_{m})-\mathcal{E}_{i}(t_{m},\delta t_{m})<d^{\prime}_{\mathcal{E}}$, leading
to $-2d^{\prime}_{\mathcal{E}}<\mathcal{E}_{f}(t_{m},\delta
t_{m})-\mathcal{E}_{i}(t_{m},\delta t_{m})<2d^{\prime}_{\mathcal{E}}$, but the
converse is not true. Therefore, information from the past time steps is used
to select the current step size, such that the algorithm is capable of
automatically counteracting any accumulating Trotter-induced heating effects.
This global control is necessary to handle the lack of energy conservation in
time-dependent systems.

We enforce these constraints via a feedback loop: initially, a large time step
$\delta t_{m}$ is chosen. We then measure $\mathcal{E}_{i}$ and
$\delta\mathcal{E}_{i}^{2}$ for the current quantum state
$|\psi(t_{m})\rangle$ and for the selected $\delta t_{m}$, as a prediction of
the piecewise conserved quantities. We then implement the time evolution
$U_{[\lambda]}(t_{m},\delta t_{m})$ on the quantum processor, yielding a
candidate state $|\tilde{\psi}(t_{m}{+}\delta
t_{m})\rangle{=}U_{[\lambda]}(t_{m},\delta t_{m})|\psi(t_{m})\rangle$. For
this, we measure $\tilde{\mathcal{E}}_{f}$ and
$\delta\tilde{\mathcal{E}}_{f}^{2}$. In case the measurement outcome violates
the conditions of Eq. (LABEL:eq.GlobalConstraints), a new smaller step size is
proposed and the procedure restarts.

We use the bisection search method to find a new suitable $\delta t_{m}$. This
can be efficiently implemented with a few trials whose number does not scale
with system size and the truncation order $k$, cf. Ref. [33] and Sec. SM 3.
The extra measurement cost only depends polynomially on the system size and
can be further improved to logarithmic dependence by using classical shadows
[54, 55, 56]. Once a suitable $\delta t_{m}$ has been found, we obtain the
state $|\psi(t_{m}{+}\delta t_{m})\rangle$ at the next time, and repeat the
procedure.

Figure 2: Comparison between the time-local and global constraint. (a) Local
constraint with tolerances $d_{\mathcal{E}}=0.01$,
$d_{\delta\mathcal{E}^{2}}=0.02,$ and
$d^{\prime}_{\mathcal{E}}=d^{\prime}_{\delta\mathcal{E}^{2}}=\infty$. (b)
Time-global errors are bounded, with the tolerance
$d_{\mathcal{E}}=d_{\delta\mathcal{E}^{2}}=\infty,d^{\prime}_{\mathcal{E}}=0.01,d^{\prime}_{\delta\mathcal{E}^{2}}=0.02$.
The following Hamiltonian parameters are used for numerical simulation,
$J_{z}=1,h_{x}=1,h_{z}=0.5,\tau=1$,$\omega=4,L=16,\theta=2$.

Numerical simulation.— We next numerically compare local and global constraint
schemes and find better performance of tADA-Trotter with time-global control.

Although this algorithm is applicable to various models and initial states,
for concreteness, we start from a product state
$\left|\psi(0)\right\rangle{=}\exp(-i\theta\sum_{j}\sigma_{j}^{x})\left|\downarrow\dots\downarrow\right\rangle$
and a non-integrable quantum Ising model, with Hamiltonian
$H(t){=}g(t)H_{x}{+}f(t)H_{z}$ with
$H_{z}{=}J_{z}\sum\nolimits_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}+h_{z}\sum\nolimits_{j}\sigma_{j}^{z},H_{x}{=}h_{x}\sum\nolimits_{j}\sigma_{j}^{x},$
where $\sigma_{j}^{x}$ and $\sigma_{j}^{z}$ are Pauli matrices acting on site
$j$ of a chain consisting of $L$ sites with periodic boundary conditions. We
consider a uniform coupling $J_{z}$, and transverse and longitudinal fields
$h_{x}$ and $h_{z}$, respectively.

We choose a static longitudinal field $f(t){=}1$ and an oscillating transverse
field $g(t){=}\cos(\omega t)\exp(-t/\tau){+}1$ with nonzero mean, frequency
$\omega$, and exponentially decaying amplitude. For $t{\gg}\tau$, the system
becomes effectively time-independent. Hence, this protocol contains different
timescales and thus provides an ideal testbed for our algorithm. We employ Eq.
(LABEL:eq.trotterU) to implement the Trotterized dynamics and truncate the
piecewise conserved Hamiltonian to $H_{[k]}$ with $k=5$.

In Fig. 2, we depict the expectation value of $H_{[k]}$ with the local and
global control schemes. The exact solution varies at early times and becomes
static at later times as expected. The predicted conserved value
$\mathcal{E}_{i}$ at early times closely follows the exact solution. The
expectation value $\mathcal{E}_{f}$ after implementing the Trotterized
dynamics deviates from $\mathcal{E}_{i}$ very weakly in both cases.

A crucial difference occurs at later times. In Fig. 2 (a), for $t{>}5$ in
units of the Ising coupling $J_{z}$, the predicted value $\mathcal{E}_{i}$
exhibits a noticeable drift towards zero, indicating the accumulation of
Trotter errors. Note, this Trotter-induced heating is not the same as the
energy-non-conservation that goes along with $H(t)$ [36, 37, 38]. It happens
because statistically, a stepsize increasing the system’s entropy is more
likely. To emphasize this point, we plot the time-local error
$\mathcal{E}_{f}(t_{n},\delta t_{n}){-}\mathcal{E}_{i}(t_{n},\delta t_{n})$ at
each time, and clearly, $\delta t$ is chosen in a way that negative values
appear more frequently. By contrast, when constraining the global errors
according to Eq. (LABEL:eq.GlobalConstraints), this deviation approximately
centers around zero, indicating better controlled heating. Consequently, the
overall drift in $\mathcal{E}_{i}$ is suppressed. For a specific time $t_{0}$,
as long as $t_{0}{\gg}\tau$ where the system becomes effectively time-
independent, one can show that $\mathcal{E}_{f}(t_{m},\delta
t_{m}){=}\mathcal{E}_{i}(t_{m+1},\delta t_{m+1})$ for any $t_{m}{\geq}t_{0}$.
Hence, Eq. (LABEL:eq.GlobalConstraints) reduces to
$\left|\mathcal{E}_{i}(t_{m},\delta t_{m}){-}\mathcal{E}_{i}(t_{0},\delta
t_{0})\right|{<}d^{\prime}_{\mathcal{E}}$, strictly prohibiting the overall
drift at long times as shown in Fig. 2 (b).

Figure 3: Global constraint leads to more stable simulation results than local
constraint. Magnetization in $x$ (a) and $z$ (b) direction. We use the same
parameters as in Fig. 2. The inset shows details of the dynamics in a small
time window.

Controlling errors in the expectation values of $H_{[k]}$ ensures accurate DQS
of local observables. Fig. 3 shows the magnetizations
$M_{\alpha}=\sum_{j}\sigma_{j}^{\alpha}/L$ for $\alpha=x,z$, demonstrating
that global constraint generally yields a more accurate simulation compared to
the local constraint. In particular, in Fig. 3 (a) with the local constraint,
significant errors arise for $t{>}5$. It corresponds to the time when notable
deviation arises in the piecewise conserved quantities as shown in Fig. 2 (a).
In contrast, the globally constrained data closely follows the exact solution
for an extended period.

We now demonstrate that by constraining the time-global errors, tADA-Trotter
achieves superior simulation precision compared to the fixed-step Trotter when
the same total simulation time is reached. In the field $x(t)$ we select a
driving frequency $\omega{=}0.8$ that is comparable to other local energy
scales in the system. The characteristic decay timescale is chosen as
$\tau{=}30$. With a total number of Trotter steps $N{=}100$, the achievable
simulation time is approximately $t{\sim}20$, during which the significant
time-dependence in the Hamiltonian is still present.

Figure 4:  Comparison between tADA-Trotter and fixed-step Trotter algorithms.
Inset depicts the stepsize that varies in time. It takes larger values when
external driving fields are weak. We use the following Hamiltonian parameters
for numerical simulation,
$J_{z}=1,h_{x}=3,h_{z}=0.5,\tau=30,\omega=0.8,L=18,\theta=2,d^{\prime}_{\mathcal{E}}=0.03,d^{\prime}_{\delta\mathcal{E}^{2}}=0.1$.

In Fig. 4, $M_{x}{=}\sum_{j}\sigma_{j}^{x}/L$ is depicted with orange circles,
which closely reproduces the exact solution (black). Simulation errors only
become visible at later times, e.g., $t_{m}{>}11$. In contrast, the fixed-step
Trotter ($\delta t_{m}{=}0.2$) already introduces substantial errors in the
magnetization within a short time. The Trotter step size (Fig. 4 inset)
fluctuates within approximately one order of magnitude $\delta
t_{m}{\in}[0.1,0.7]$, highlighting the advantage and the flexibility of tADA-
Trotter. Particularly, at early times, when the quantum state undergoes rapid
changes under a strong driving field, smaller step sizes are employed ($\delta
t_{m}{\approx}0.1$). Conversely, when the driving $g(t)$ has relatively
smaller values around $t_{m}{\approx}4$ and $12$, the step size automatically
increases to $\delta t{\approx}0.7$ and $0.4$, respectively.

The time-global errors, Fig. 5, remain bounded below the specified thresholds
for the majority of the time evolution. However, it should be noted that due
to the tight tolerances at early times, the accumulated errors in the
piecewise conserved quantities may occasionally exceed the bounds. This
phenomenon can also cause the tADA-Trotter to “freeze”, wherein it tends to
select the smallest possible step size, which in this case is set to $0.1$.

Discussion.— We have introduced the concept of piecewise conserved quantities
and combined it with a time-global error constraint to devise an adaptive
Trotterisation scheme to enable reliable DQS of generic time-dependent
Hamiltonians. In the presence of weak hardware noise, we expect that tADA-
Trotter remains robust and can even self-correct the dissipation-induced
errors 222We have explicitly demonstrated this for a time-independent
Hamiltonian in Ref. [33]. We estimate that the shot overhead for a single step
is around $10^{4}$, which can be efficiently implemented on current quantum
processors with fast gate implementation, cf. details in Sec. SM 3.

Figure 5: Time-global errors under the global constraint in the expectation
value and variance of $H_{[k]}$. Errors are bounded for most of the time
evolution. Same parameters as in Fig. 4.

In a specific example, we demonstrate the superior performance of tADA-Trotter
compared to the fixed-step Trotter method. The intricate interplay between
external driving and quantum thermalization can result in highly complex many-
body dynamics [58, 59, 60, 61]. Therefore, for future investigations,
conducting a systematic benchmark of various algorithms across different
models, initial states, and time-dependence would be of great value, as would
be a comparison against higher-order truncation schemes.

Furthermore, the extension of piecewise conservation laws to open quantum
systems to enable adaptive Trotter step sizes represents an intriguing open
question [62, 63, 64, 65]. Additionally, considering the widespread use of
Trotterization in classical numerical algorithms such as the time-evolving
block decimation method, the application of tADA-Trotter to enhance the
efficiency and accuracy of these methods holds significant potential.

Note added.— During the completion of this work, we became aware of a relevant
work exploring another adaptive algorithm for DQS of time evolution [45].

Acknowledgments We thank T. Ikeda for enlightening discussions. This work is
in part supported by “The Fundamental Research Funds for the Central
Universities, Peking University”, and by “High-performance Computing Platform
of Peking University” and the Deutsche Forschungsgemeinschaft under cluster of
excellence ct.qmat (EXC 2147, project-id 390858490). This project has received
funding from the European Research Council (ERC) under the European Unions
Horizon 2020 research and innovation programme (grant agreement No. 853443).
MB was supported in part by the International Centre for Theoretical Sciences
(ICTS) for participating in the program - Periodically and quasi-periodically
driven complex systems (code: ICTS/pdcs2023/6).

## References

  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Reviews of Modern Physics 86, 153 (2014).
  * Preskill [2018] J. Preskill, Quantum computing in the nisq era and beyond, Quantum 2, 79 (2018).
  * Blatt and Roos [2012] R. Blatt and C. F. Roos, Quantum simulations with trapped ions, Nature Physics 8, 277 (2012).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, _et al._ , Programmable quantum simulations of spin systems with trapped ions, Reviews of Modern Physics 93, 025001 (2021).
  * Dumitrescu  _et al._ [2022] P. T. Dumitrescu, J. G. Bohnet, J. P. Gaebler, A. Hankin, D. Hayes, A. Kumar, B. Neyenhuis, R. Vasseur, and A. C. Potter, Dynamical topological phase realized in a trapped-ion quantum simulator, Nature 607, 463 (2022).
  * Salathé  _et al._ [2015] Y. Salathé, M. Mondal, M. Oppliger, J. Heinsoo, P. Kurpiers, A. Potočnik, A. Mezzacapo, U. Las Heras, L. Lamata, E. Solano, _et al._ , Digital quantum simulation of spin models with circuit quantum electrodynamics, Physical Review X 5, 021027 (2015).
  * Satzinger  _et al._ [2021] K. Satzinger, Y.-J. Liu, A. Smith, C. Knapp, M. Newman, C. Jones, Z. Chen, C. Quintana, X. Mi, A. Dunsworth, _et al._ , Realizing topologically ordered states on a quantum processor, Science 374, 1237 (2021).
  * Dborin  _et al._ [2022] J. Dborin, V. Wimalaweera, F. Barratt, E. Ostby, T. E. O’Brien, and A. G. Green, Simulating groundstate and dynamical quantum phase transitions on a superconducting quantum computer, Nature Communications 13, 5977 (2022).
  * Moses  _et al._ [2023] S. Moses, C. Baldwin, M. Allman, R. Ancona, L. Ascarrunz, C. Barnes, J. Bartolotta, B. Bjork, P. Blanchard, M. Bohn, _et al._ , A race track trapped-ion quantum processor, arXiv preprint arXiv:2305.03828 (2023).
  * Chen  _et al._ [2023] C. Chen, G. Bornet, M. Bintz, G. Emperauger, L. Leclerc, V. S. Liu, P. Scholl, D. Barredo, J. Hauschild, S. Chatterjee, _et al._ , Continuous symmetry breaking in a two-dimensional rydberg array, Nature 616, 691 (2023).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, _et al._ , Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * Suzuki [1991] M. Suzuki, General theory of fractal path integrals with applications to many-body theories and statistical physics, Journal of Mathematical Physics 32, 400 (1991).
  * Berry  _et al._ [2007] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Efficient quantum algorithms for simulating sparse hamiltonians, Communications in Mathematical Physics 270, 359 (2007).
  * Poulin  _et al._ [2014] D. Poulin, M. B. Hastings, D. Wecker, N. Wiebe, A. C. Doherty, and M. Troyer, The trotter step size required for accurate quantum simulation of quantum chemistry, arXiv preprint arXiv:1406.4920 (2014).
  * Babbush  _et al._ [2016] R. Babbush, D. W. Berry, I. D. Kivlichan, A. Y. Wei, P. J. Love, and A. Aspuru-Guzik, Exponentially more precise quantum simulation of fermions in second quantization, New Journal of Physics 18, 033032 (2016).
  * 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).
  * Tranter  _et al._ [2019] A. Tranter, P. J. Love, F. Mintert, N. Wiebe, and P. V. Coveney, Ordering of trotterization: Impact on errors in quantum simulation of electronic structure, Entropy 21, 1218 (2019).
  * Cirstoiu  _et al._ [2020] C. Cirstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger, Variational fast forwarding for quantum simulation beyond the coherence time, npj Quantum Information 6, 1 (2020).
  * Bolens and Heyl [2021] A. Bolens and M. Heyl, Reinforcement learning for digital quantum simulation, Physical Review Letters 127, 110502 (2021).
  * Yao  _et al._ [2021] Y.-X. Yao, N. Gomes, F. Zhang, C.-Z. Wang, K.-M. Ho, T. Iadecola, and P. P. Orth, Adaptive variational quantum dynamics simulations, PRX Quantum 2, 030307 (2021).
  * Lin  _et al._ [2021] S.-H. Lin, R. Dilip, A. G. Green, A. Smith, and F. Pollmann, Real-and imaginary-time evolution with compressed quantum circuits, PRX Quantum 2, 010342 (2021).
  * Richter and Pal [2021] J. Richter and A. Pal, Simulating hydrodynamics on noisy intermediate-scale quantum devices with random circuits, Physical Review Letters 126, 230501 (2021).
  * Mansuroglu  _et al._ [2021] R. Mansuroglu, T. Eckstein, L. Nützel, S. A. Wilkinson, and M. J. Hartmann, Variational hamiltonian simulation for translational invariant systems via classical pre-processing, arXiv e-prints , arXiv (2021).
  * Keever and Lubasch [2022] C. M. Keever and M. Lubasch, Classically optimized hamiltonian simulation, arXiv preprint arXiv:2205.11427 (2022).
  * 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, arXiv preprint arXiv:2203.15846 (2022).
  * Tepaske  _et al._ [2022] M. S. Tepaske, D. Hahn, and D. J. Luitz, Optimal compression of quantum many-body time evolution operators into brickwall circuits, arXiv preprint arXiv:2205.03445 (2022).
  * Zhang  _et al._ [2023] Z.-J. Zhang, J. Sun, X. Yuan, and M.-H. Yung, Low-depth hamiltonian simulation by an adaptive product formula, Physical Review Letters 130, 040601 (2023).
  * Mansuroglu  _et al._ [2023] R. Mansuroglu, F. Fischer, and M. J. Hartmann, Problem specific classical optimization of hamiltonian simulation, arXiv preprint arXiv:2306.07208 (2023).
  * Granet and Dreyer [2023] E. Granet and H. Dreyer, Continuous hamiltonian dynamics on noisy digital quantum computers without trotter error, arXiv preprint arXiv:2308.03694 (2023).
  * Postler  _et al._ [2022] L. Postler, S. Heu$\beta$en, I. Pogorelov, M. Rispler, T. Feldker, M. Meth, C. D. Marciniak, R. Stricker, M. Ringbauer, R. Blatt, _et al._ , Demonstration of fault-tolerant universal quantum gate operations, Nature 605, 675 (2022).
  * Krinner  _et al._ [2022] S. Krinner, N. Lacroix, A. Remm, A. Di Paolo, E. Genois, C. Leroux, C. Hellings, S. Lazar, F. Swiadek, J. Herrmann, _et al._ , Realizing repeated quantum error correction in a distance-three surface code, Nature 605, 669 (2022).
  * Chen  _et al._ [2022] E. H. Chen, T. J. Yoder, Y. Kim, N. Sundaresan, S. Srinivasan, M. Li, A. D. Córcoles, A. W. Cross, and M. Takita, Calibrated decoders for experimental quantum error correction, Physical Review Letters 128, 110504 (2022).
  * Zhao  _et al._ [2023] H. Zhao, M. Bukov, M. Heyl, and R. Moessner, Making trotterization adaptive and energy-self-correcting for nisq devices and beyond, PRX Quantum 4, 030319 (2023).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Reviews of Modern Physics 91, 021001 (2019).
  * Lazarides  _et al._ [2014] A. Lazarides, A. Das, and R. Moessner, Equilibrium states of generic quantum systems subject to periodic driving, Physical Review E 90, 012110 (2014).
  * Abanin  _et al._ [2015] D. A. Abanin, W. De Roeck, and F. Huveneers, Exponentially slow heating in periodically driven many-body systems, Physical review letters 115, 256803 (2015).
  * 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).
  * Poulin  _et al._ [2011] D. Poulin, A. Qarry, R. Somma, and F. Verstraete, Quantum simulation of time-dependent hamiltonians and the convenient illusion of hilbert space, Physical review letters 106, 170501 (2011).
  * Chen  _et al._ [2021] Y.-H. Chen, A. Kalev, and I. Hen, Quantum algorithm for time-dependent hamiltonian simulation by permutation expansion, PRX Quantum 2, 030342 (2021).
  * Lau  _et al._ [2021] J. W. Z. Lau, K. Bharti, T. Haug, and L. C. Kwek, Noisy intermediate scale quantum simulation of time dependent hamiltonians, arXiv preprint arXiv:2101.07677 (2021).
  * Watkins  _et al._ [2022] J. Watkins, N. Wiebe, A. Roggero, and D. Lee, Time-dependent hamiltonian simulation using discrete clock constructions, arXiv preprint arXiv:2203.11353 (2022).
  * Ikeda  _et al._ [2022] T. N. Ikeda, A. Abrar, I. L. Chuang, and S. Sugiura, Minimum fourth-order trotterization formula for a time-dependent hamiltonian, arXiv preprint arXiv:2212.06788 (2022).
  * An  _et al._ [2022] D. An, D. Fang, and L. Lin, Time-dependent hamiltonian simulation of highly oscillatory dynamics and superconvergence for schrödinger equation, Quantum 6, 690 (2022).
  * Ikeda and Fujii [2023] T. N. Ikeda and K. Fujii, Trotter24: A precision-guaranteed adaptive stepsize trotterization for hamiltonian simulations (2023), arXiv:2307.05406 [quant-ph] .
  * Kovalsky  _et al._ [2023] L. K. Kovalsky, F. A. Calderon-Vargas, M. D. Grace, A. B. Magann, J. B. Larsen, A. D. Baczewski, and M. Sarovar, Self-healing of trotter error in digital adiabatic state preparation, Physical Review Letters 131, 060602 (2023).
  * Hochbruck and Lubich [2003] M. Hochbruck and C. Lubich, On magnus integrators for time-dependent schrödinger equations, SIAM journal on numerical analysis 41, 945 (2003).
  * Alvermann and Fehske [2011] A. Alvermann and H. Fehske, High-order commutator-free exponential time-propagation of driven quantum systems, Journal of Computational Physics 230, 5930 (2011).
  * Kokail  _et al._ [2019] C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, _et al._ , Self-verifying variational quantum simulation of lattice models, Nature 569, 355 (2019).
  * Mi  _et al._ [2021] X. Mi, P. Roushan, C. Quintana, S. Mandra, J. Marshall, C. Neill, F. Arute, K. Arya, J. Atalaya, R. Babbush, _et al._ , Information scrambling in quantum circuits, Science 374, 1479 (2021).
  * Naldesi  _et al._ [2022] P. Naldesi, A. Elben, A. Minguzzi, D. Clément, P. Zoller, and B. Vermersch, Fermionic correlation functions from randomized measurements in programmable atomic quantum devices, arXiv preprint arXiv:2205.00981 (2022).
  * Hartmann  _et al._ [2004] M. Hartmann, O. Hess, _et al._ , Gaussian quantum fluctuations in interacting many particle systems, Letters in Mathematical Physics 68, 103 (2004).
  * Note [1] To show this, we define $\Delta_{m-1}{=}\sum@\slimits@\nolimits_{n=1}^{m-1}\mathcal{E}_{f}(t_{n},\delta t_{n})-\mathcal{E}_{i}(t_{n},\delta t_{n}),$ which satisfies $-d^{\prime}_{\mathcal{E}}{<}\Delta_{m-1}{<}d^{\prime}_{\mathcal{E}}$ with a positive $d^{\prime}_{\mathcal{E}}$. For the next step, one has $-d^{\prime}_{\mathcal{E}}<\Delta_{m-1}+\mathcal{E}_{f}(t_{m},\delta t_{m})-\mathcal{E}_{i}(t_{m},\delta t_{m})<d^{\prime}_{\mathcal{E}}$, leading to $-2d^{\prime}_{\mathcal{E}}<\mathcal{E}_{f}(t_{m},\delta t_{m})-\mathcal{E}_{i}(t_{m},\delta t_{m})<2d^{\prime}_{\mathcal{E}}$.
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics 16, 1050 (2020).
  * Huang  _et al._ [2021] H.-Y. Huang, R. Kueng, and J. Preskill, Efficient estimation of pauli observables by derandomization, Physical review letters 127, 030503 (2021).
  * Elben  _et al._ [2023] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, The randomized measurement toolbox, Nature Reviews Physics 5, 9 (2023).
  * Note [2] We have explicitly demonstrated this for a time-independent Hamiltonian in Ref. [33].
  * Kolodrubetz  _et al._ [2017] M. Kolodrubetz, D. Sels, P. Mehta, and A. Polkovnikov, Geometry and non-adiabatic response in quantum and classical systems, Physics Reports 697, 1 (2017).
  * Sugiura  _et al._ [2021] S. Sugiura, P. W. Claeys, A. Dymarsky, and A. Polkovnikov, Adiabatic landscape and optimal paths in ergodic systems, Physical Review Research 3, 013102 (2021).
  * De La Torre  _et al._ [2021] A. De La Torre, D. M. Kennes, M. Claassen, S. Gerber, J. W. McIver, and M. A. Sentef, Colloquium: Nonthermal pathways to ultrafast control in quantum materials, Reviews of Modern Physics 93, 041002 (2021).
  * 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).
  * Han  _et al._ [2021] J. Han, W. Cai, L. Hu, X. Mu, Y. Ma, Y. Xu, W. Wang, H. Wang, Y. Song, C.-L. Zou, _et al._ , Experimental simulation of open quantum system dynamics via trotterization, Physical Review Letters 127, 020504 (2021).
  * de Jong  _et al._ [2021] W. A. de Jong, M. Metcalf, J. Mulligan, M. Płoskoń, F. Ringer, and X. Yao, Quantum simulation of open quantum systems in heavy-ion collisions, Physical Review D 104, L051501 (2021).
  * Kamakari  _et al._ [2022] H. Kamakari, S.-N. Sun, M. Motta, and A. J. Minnich, Digital quantum simulation of open quantum systems using quantum imaginary–time evolution, PRX Quantum 3, 010320 (2022).
  * Vovk and Pichler [2022] T. Vovk and H. Pichler, Entanglement-optimal trajectories of many-body quantum markov processes, Physical Review Letters 128, 243601 (2022).

Supplementary Material
Adaptive Trotterization for time-dependent Hamiltonian quantum dynamics using
piecewise conservation laws

###### Contents

  1. SM 1 Perturbative analysis of errors
  2. SM 2 Numerical simulation with different $k$
  3. SM 3 Costs
    1. SM 3.1 Measurement of conservation laws
    2. SM 3.2 Number of trials
    3. SM 3.3 Simulation time versus Trotter step

##  SM 1 Perturbative analysis of errors

In this section, by using a perturbative expansion, we show that errors in the
piecewise conservation law are subdominant compared to the Trotter error if
$k\geq\lambda.$

At a given time $t$ and a Trotter step size $\delta t$, the exact time
evolution operator can be generated by the piecewise conserved quantity
$H_{[\infty]}$, $U(t{+}\delta t,t){=}\exp\left(-iH_{[\infty]}\delta t\right)$.
The target unitary operator can be Trotterized up to a given order $\lambda$,
such that $U(t{+}\delta t,t)=U_{[\lambda]}(t{+}\delta
t,t){+}\mathcal{O}(\delta t^{\lambda}).$ One can formally rewrite the
Trotterized unitary $U_{[\lambda]}(t{+}\delta
t,t){=}\exp\left(-i\tilde{H}_{[\lambda]}\delta t\right)$ with

$H_{[\infty]}=\tilde{H}_{[\lambda]}+\mathcal{O}(\delta t^{\lambda-1}),$ (S.1)

where we use $\sim$ to distinguish $\tilde{H}_{[\lambda]}$ with the
approximated piecewise conserved quantity $H_{[k]}$, and we have

$\displaystyle H_{[\infty]}=H_{[k]}+\mathcal{O}(\delta t^{k}).$ (S.2)

$H_{[k]}$ will be measured to adapt the step size. Suppose
$\mathcal{E}_{i,k}=\langle\psi|{H_{[k]}}|\psi\rangle/L$ for a given state
$\left|\psi\right\rangle$ and
$\mathcal{E}_{f,k}=\langle\psi|U_{[\lambda]}^{\dagger}{H_{[k]}}U_{[\lambda]}|\psi\rangle/L$
represents the expectation value of ${H_{[k]}}$ after the Trotterized
evolution. For $k\to\infty$, we have

$\displaystyle\mathcal{E}_{f,\infty}=L^{-1}\langle\psi|U_{[\lambda]}^{\dagger}{H_{[\infty]}}U_{[\lambda]}|\psi\rangle$
$\displaystyle=L^{-1}\langle{\psi}|\exp[i\tilde{H}_{[\lambda]}\delta
t]H_{[\infty]}\exp[-i\tilde{H}_{[\lambda]}\delta t]|{\psi}\rangle$ (S.3)
$\displaystyle=L^{-1}\langle{\psi}|H_{[\infty]}+i\delta
t[\tilde{H}_{[\lambda]},H_{[\infty]}]+\frac{1}{2}(i\delta
t)^{2}[\tilde{H}_{[\lambda]},[\tilde{H}_{[\lambda]},H_{[\infty]}]]|{{\psi}}\rangle+\dots$
$\displaystyle=\mathcal{E}_{i,\infty}+L^{-1}\langle{\psi}|i\delta
t[\tilde{H}_{[\lambda]},H_{[\infty]}]+\frac{1}{2}(i\delta
t)^{2}[\tilde{H}_{[\lambda]},[\tilde{H}_{[\lambda]},H_{[\infty]}]]|{{\psi}}\rangle+\dots$

where we use the formula

$\displaystyle e^{A}Be^{-A}=B+[A,B]+\frac{1}{2}[A,[A,B]]+\ldots,$ (S.4)

to perform the perturbative expansion in $\delta t$ and neglect contributions
from higher orders. By using Eq. (S.1), one has
$[\tilde{H}_{[\lambda]},H_{[\infty]}]\sim\delta t^{\lambda-1}$, and
$[\tilde{H}_{[\lambda]},[\tilde{H}_{[\lambda]},H_{[\infty]}]]\sim\delta
t^{\lambda-1}$. Therefore the leading order contribution to the change in the
piecewise conservation law reads

$\displaystyle\mathcal{E}_{f,\infty}-\mathcal{E}_{i,\infty}=L^{-1}\langle{\psi}|i\delta
t[\tilde{H}_{[\lambda]},H_{[\infty]}]|{{\psi}}\rangle+\mathcal{O}(\delta
t^{\lambda+1})\sim\delta t^{\lambda}.$ (S.5)

Ideally, if one can measure the expectation value of $H_{[\infty]}$,
$\mathcal{E}_{f,\infty}-\mathcal{E}_{i,\infty}$ will be used to adapt the step
size according to tADA-Trotter algorithm. However, for generic many-body
systems, $H_{\infty}$ is a highly non-local and complicated operator. We now
truncate it to a finite order $k$ such that $H_{[k]}$ can have a local
structure and can be measured in practice. However, such truncation introduces
errors in addition to the Trotter error (its dependence on $\delta t$ is
purely determined by $\lambda$). Now we define this additional error for a
finite $k$ as
$\Delta=(\mathcal{E}_{f,\infty}-\mathcal{E}_{i,\infty})-(\mathcal{E}_{f,k}-\mathcal{E}_{i,k})$,
and we aim to find the condition such that $\Delta$ is subdominant than
$\mathcal{E}_{f,\infty}-\mathcal{E}_{i,\infty}$. Such a comparison is
generally a challenging task since these quantities highly depend on the state
and the time-dependent Hamiltonian. Instead, we analyse the scaling of
$\Delta$ versus $\delta t$, and compare its asymptotic behavior between
$\Delta$ and $\mathcal{E}_{f,\infty}-\mathcal{E}_{i,\infty}$ for $\delta t\to
0$. By definition,

$\displaystyle\begin{aligned}
\Delta&=(\mathcal{E}_{f,\infty}-\mathcal{E}_{i,\infty})-(\mathcal{E}_{f,k}-\mathcal{E}_{i,k})=\langle\psi|U_{[\lambda]}^{\dagger}(H_{[\infty]}-H_{[k]})U_{[\lambda]}|\psi\rangle-\langle\psi|(H_{[\infty]}-H_{[k]})|\psi\rangle.\end{aligned}$
(S.6)

By using Eq. (S.4), we expand it perturbatively in $\delta t$,

$\displaystyle\begin{aligned} \Delta=\langle{\psi}|i\delta
t[\tilde{H}_{[\lambda]},(H_{[\infty]}-H_{[k]})]+\frac{1}{2}(i\delta
t)^{2}[\tilde{H}_{[\lambda]},[\tilde{H}_{[\lambda]},(H_{[\infty]}-H_{[k]})]]|{{\psi}}\rangle+\dots,\end{aligned}$
(S.7)

where the leading order contribution scales as $\delta t^{k+1}$ according to
Eq. (S.2). Hence, to ensure that that errors in the expectation value of
piecewise conserved quantity are subdominant compared to the Trotter error,
one needs $\mathcal{O}(\delta t^{k+1})<\mathcal{O}(\delta t^{\lambda})$, or
equivalently, $k\geq\lambda.$

Now we argue the same condition suffices to ensure the variance of piecewise
conserved quantity is dominated by Trotter error. The variance is defined as

$\displaystyle\begin{aligned}
\delta\mathcal{E}^{2}_{i,k}{=}L^{-1}\langle\psi|H_{[k]}^{2}|\psi\rangle{-}L\mathcal{E}^{2}_{i,k},\delta\mathcal{E}^{2}_{f,k}{=}L^{-1}\langle\psi|U_{[\lambda]}^{\dagger}H_{[k]}^{2}U_{[\lambda]}|\psi\rangle{-}L\mathcal{E}^{2}_{f,k}.\end{aligned}$
(S.8)

For $k\to\infty$, one has

$\displaystyle\begin{aligned}
\delta\mathcal{E}^{2}_{f,\infty}-\delta\mathcal{E}^{2}_{i,\infty}&=L^{-1}\left(\langle\psi|U_{[\lambda]}^{\dagger}H_{[\infty]}^{2}U_{[\lambda]}|\psi\rangle-\langle\psi|H_{[\infty]}^{2}|\psi\rangle\right)-L\left(\mathcal{E}^{2}_{f,\infty}-\mathcal{E}^{2}_{i,\infty}\right)\\\
&\approx L^{-1}\langle{\psi}|i\delta
t[\tilde{H}_{[\lambda]},H^{2}_{[\infty]}]|{{\psi}}\rangle-L\left(\mathcal{E}_{f,\infty}-\mathcal{E}_{i,\infty}\right)\left(\mathcal{E}_{f,\infty}+\mathcal{E}_{i,\infty}\right).\end{aligned}$
(S.9)

By using Eq. (S.1), Eq. (S.5) and the fact
$\left(\mathcal{E}_{f,\infty}+\mathcal{E}_{i,\infty}\right)\sim\delta t^{0}$,
one obtains
$\delta\mathcal{E}^{2}_{f,\infty}-\delta\mathcal{E}^{2}_{i,\infty}\sim\delta
t^{\lambda}$, which exhibits the same scaling as Eq. (S.5). Once we truncate
$k$ to a finite number, one obtains the additional error for the variance

$\displaystyle\begin{aligned}
&(\delta\mathcal{E}^{2}_{f,\infty}-\delta\mathcal{E}^{2}_{i,\infty})-(\delta\mathcal{E}^{2}_{f,k}-\delta\mathcal{E}^{2}_{i,k})\\\
&=L^{-1}\left(\langle{\psi}|U_{[\lambda]}^{\dagger}(H^{2}_{[\infty]}-H^{2}_{[k]})U_{[\lambda]}|{{\psi}}\rangle-\langle{\psi}|(H^{2}_{[\infty]}-H^{2}_{[k]})|{{\psi}}\rangle\right)-L\left[\left(\mathcal{E}^{2}_{f,\infty}-\mathcal{E}^{2}_{i,\infty}\right)-\left(\mathcal{E}^{2}_{f,k}-\mathcal{E}^{2}_{i,k}\right)\right]\end{aligned}$
(S.10)

With Eq. (S.2) and Eq. (S.4) we know the leading order contribution in the
first term scales as $\delta t^{k+1}$. To obtain the scaling for the second
term, we first use Eq. (S.3) for a finite $k$ and approximate
$\mathcal{E}_{f,k}$ as
$\mathcal{E}_{f,k}\approx\mathcal{E}_{i,k}+\langle{\psi}|i\delta
t[\tilde{H}_{[\lambda]},H_{[k]}]\left|\psi\right\rangle/L,$ hence we have
$\mathcal{E}^{2}_{f,k}\approx\mathcal{E}^{2}_{i,k}+2i\delta
t\mathcal{E}_{i,k}\langle{\psi}|[\tilde{H}_{[\lambda]},H_{[k]}]\left|\psi\right\rangle/L,$
Therefore the second term in Eq. (LABEL:eq.var_error) becomes

$\displaystyle\begin{aligned}
\left(\mathcal{E}^{2}_{f,\infty}-\mathcal{E}^{2}_{i,\infty}\right)-\left(\mathcal{E}^{2}_{f,k}-\mathcal{E}^{2}_{i,k}\right)&={2i\delta
t}{L^{-1}}\left(\mathcal{E}_{i,\infty}\langle\psi|[\tilde{H}_{\lambda},H_{[\infty]}]|\psi\rangle-\mathcal{E}_{i,k}\langle\psi|[\tilde{H}_{\lambda},H_{[k]}]|\psi\rangle\right)\\\
&={2i\delta
t}{L^{-1}}\left(\mathcal{E}_{i,k}\langle\psi|[\tilde{H}_{\lambda},H_{[\infty]}-H_{[k]}]|\psi\rangle+(\mathcal{E}_{i,\infty}-\mathcal{E}_{i,k})\langle\psi|[\tilde{H}_{\lambda},H_{[\infty]}]|\psi\rangle\right).\end{aligned}$
(S.11)

By noticing $\mathcal{E}_{i,k}$ scales as $\delta t^{0}$,
$\mathcal{E}_{i,\infty}-\mathcal{E}_{i,k}\sim\delta t^{k}$ and also using Eq.
(S.2) and Eq. (S.1), one can show that the leading order contribution to the
equation above also scales as $\delta t^{k+1}$. Consequently, to ensure that
that errors in the expectation value and variance of the piecewise
conservation law are subdominant compared to the Trotter error asymptotically
for a small $\delta t$, one needs $k\geq\lambda.$

It is worth noting that the argument above requires a small $\delta t$ for
Trotterized dynamics, implying that the circuit depth to reach certain target
time can be large and unrealistic in the NISQ era. For any practical purpose,
the optimal $k$ and $\lambda$ may highly depend on the target time-dependent
systems, the measurement resources and gate fidelity in specific experimental
setups.

##  SM 2 Numerical simulation with different $k$

Sec. SM 1 suggests that $k\geq\lambda$ would be a suitable choice to ensure
Trotter error is the only dominant error. For the Trotter decomposition with
$\lambda=3$ as considered in the main text, the minimal $k$ to satisfy this
condition would be $k=3$. In the main text, we fix $k=5$ to ensure that
$\Delta$ defined in Eq. (S.6) can be negligible. However, according to Eq.
(LABEL:eq.approximation), a large $k$ introduces extra measurement overheads.
In this section we numerically show that a smaller $k$ can also generate
sufficiently accurate Trotterized time-evolution with the global constraints,
therefore the measurement overheads can be significantly reduced.

Figure S1:  Comparison between tADA-Trotter with different truncation order
$k$ for the piecewise conservation law. (a) and (d) depict the dynamics of the
magnetization: Both $k=3$ and $5$ lead to accurate simulation results compared
with the exact solution(black). Panels (b) and (e) show the step sizes
selected from tADA-Trotter: (b) uses a smaller maximum step size ($\delta
t_{\text{max}}=0.2$) when the bisection search is implemented. Hence, errors
in the piecewise conservation law $\Delta$ is subdominant and $k=3$ and $k=5$
lead to almost the same simulation outcomes. A larger $\delta t_{\text{max}}$
is used for (e), and differences in the selection of the step sizes are
clearly visible for $k=3$ and $5$. (c) and (f) depict the number of trials
that have been implemented for the bisection search. On average this algorithm
requires 10 trials before determining the optimal step size. We use the
following Hamiltonian parameters for numerical simulation,
$J_{z}=1,h_{x}=3,h_{z}=0.5,\tau=30,\omega=0.8,L=16,\theta=2,d^{\prime}_{\mathcal{E}}=0.01,d^{\prime}_{\delta\mathcal{E}^{2}}=0.02$.

In Fig. S1, we first compare the simulation results with $k=3$ (blue) and $5$
(orange). For panels (a) (b) (c) we use $\delta t_{\text{max}}=0.2$ as the
largest allowed step size for tADA-Trotter, whereas for (d)(e)(f) we use
$\delta t_{\text{max}}=0.5$. Since the overall step sizes shown in panel (b)
are relatively small, at each step, $H_{[k]}$ is almost the same for $k=3$ and
$k=5$. Therefore the decision of $\delta t$ in panel (b) and the resulting
dynamics in the magnetization are also almost the same, especially at early
times. For (d)(e)(f), because larger step sizes are allowed, the difference in
$H_{[k]}$ becomes visible and step sizes may be chosen differently for $k=3$
and $5$. Both simulation results follow closely with the exact solution, but
the measurement overhead can be substantially smaller for $k=3$.

Figure S2: $k=1$ can also generate accurate Trotterized dynamics although the
interplay between the error in the approximated piecewise conservation law
$\Delta$ and the Trotter error becomes difficult to control. (a): dynamics of
the magnetization. (b): step sizes selected from tADA-Trotter. (c): the number
of trials that have been implemented for the bisection search. On average this
algorithm requires 10 trials before determining the optimal step size. We use
the following Hamiltonian parameters for numerical simulation,
$J_{z}=1,h_{x}=3,h_{z}=0.5,\tau=30,\omega=0.8,L=16,\theta=2,d^{\prime}_{\mathcal{E}}=0.01,d^{\prime}_{\delta\mathcal{E}^{2}}=0.02$.

The asymptotic limit $\delta t\to 0$ implies a divergent circuit depth to
reach certain physical simulation time $t$. Hence, in practice, one normally
uses larger step sizes on NISQ devices. In this case, even if the condition
$k\geq\lambda$ is not satisfied, tADA-Trotter may still be able to generate
reliable simulation of many-body dynamics. We demonstrate it in Fig. S2 where
$k=1$ is used. However, this observation is not universal and the performance
of tADA-Trotter can highly rely on specific time-dependence of the system.

##  SM 3 Costs

###  SM 3.1 Measurement of conservation laws

tADA-Trotter introduces a measurement overhead to obtain expectation values of
the piecewise conserved quantities. However, for a time-dependent and local
Hamiltonian $H(t)$, $H_{[k]}$ with a finite $k$ will always remain local.
Hence, measuring the expectation value of $H_{[k]}$ involves terms linear in
systems size $L$ and measuring $H_{[k]}^{2}$ only involves terms quadratic in
$L$. Indeed, a small value of $k$ can already generate reliable tADA-Trotter
simulation results as shown in Sec. SM 2. Furthermore, the system size
dependence can be improved based on randomized measurements and classical
shadows as recently introduced [54, 55, 56]. For instance, the measurement
cost for the $H_{[k]}$ can be $L$-independent and the cost for the Hamiltonian
variance has been improved to depend only logarithmically on the system size
[54]. The measurement shot noise normally scales as $1/\sqrt{N}$ where $N$
denotes the number of measurement shots. Therefore, to reach certain accuracy
$d^{\prime}_{\delta\mathcal{E}}$ , $N$ roughly scales as
$1/{d^{\prime}}^{2}_{\delta\mathcal{E}}$. In principle, one can always perform
more measurement shots to increase this accuracy. However, in practice,
depending on the specific experimental setup, one would estimate the total
measurement budgets and roughly distribute them over each Trotter step to
preset $d^{\prime}_{\delta\mathcal{E}}$.

In terms of actual numbers, to reach the simulation accuracy of Fig. 4, we
estimate that the shot overhead would be $10^{3}$ in a state-of-the-art
quantum simulators for a single trial. This estimate is based on recent work
which implemented a variational quantum simulation to identify the ground
state of a target Hamiltonian on a trapped ion quantum simulator [49]. There,
to quantify the variational state, measurement of energy and energy variance
were both performed. There, $10^{2}$ shots are needed to obtain a good
estimate of the energy for a given state, which is very efficient. Measurement
of energy variance is more costly: To ensure that the variational state is
sufficiently close to the ground state, the authors performed around
$4*10^{4}$ shots to resolve the energy gap (which is around 2 in units of
their local energy scale) between the ground state and first excited state.
$4*10^{4}$ shots are shown to be sufficient to determine an energy variance
around 1.2 for a given variational state. Energy variance density is around
$10^{-1}$, which is around the same accuracy required to obtain the results in
our Fig. 4. In Ref. [54], the authors further reduce this shot number to
$2*10^{3}$ by using derandomized classical shadow methods.

The Hamiltonians studied in those works are 2-local, i.e., they contain
single-site and two-site spin-1/2 operators. In our algorithm, if we choose
the order $k=1$ or 3 to approximate conservation laws, the effective
Hamiltonian is also 2-local. Hence, the shot overhead for energy and energy
variance is roughly the same for a given trial state in our algorithm. With
the bisection search method, the number of trials for a single Trotter step is
5-10 (see discussion in the next section), hence, the shot overhead for a
single step is around $10^{4}$.

On more established quantum processors with efficient gate implementation, for
instance superconducting qubits, the extra runtime for these overheads is
truly negligible. One layer of CNOT gates roughly costs $10^{-7}$ seconds on
IBM quantum processors, and a circuit involving hundreds of layers of CNOT for
Tottered dynamics only costs $10^{-5}$ seconds. Even with $10^{4}$ shot
overheads, a single step of tADA-Trotter only takes less than 1 second.

Therefore, our algorithm is definitely within the current reach of quantum
processors and it provides the great potential to substantially improve the
accuracy of Tottered dynamics.

###  SM 3.2 Number of trials

To determine the optimal step size, we currently use a bisection search
algorithm, originally proposed in our previous work [33]. To reach certain
search accuracy $\epsilon$, it introduces extra $\log(\epsilon)$ overheads.
For tolerances considered in this work, tADA-Trotter uses around 10 trials on
average before identifying the optimal step size, see Fig. S1 and Fig. S2.
This trail number does not increase for larger system sizes and a larger
truncation order $k$ (cf. panels (c) and (f) in Fig. S1), permitting the
implementation of tADA-Trotter on real quantum devices that may have access to
a much larger number of qubits.

We are currently developing a technique dubbed “Energy-Landscape Estimation”
to further reduce the number of trials. By exploiting the continuity and
analyticity of the Magnus expansion of the conserved quantity as a function of
Trotter step size, the dependence of the energy and its variance can be
efficiently estimated by only using five or even smaller trials. This method
can be applied to both time-dependent and -independent systems and will
further improve the efficiency of the tADA-Trotter.

###  SM 3.3 Simulation time versus Trotter step

Whenever a destructive measurement is performed, one needs to restart the
experiment from the beginning and run the circuit again. Hence, the total
simulation time on actual devices scales as $n^{2}$ where $n$ is the number of
Trotter steps. Note, this is indeed the same scaling as the fixed step Trotter
algorithm, because a measurement is performed after each Trotter step to
obtain the transient dynamics during the time-evolution. Therefore, this does
not introduce extra overhead.