# Hamiltonian Simulation in the Interaction Picture using the Magnus Expansion

Kunal Sharma<EMAIL_ADDRESS>IBM Quantum, IBM T.J. Watson Research Center,
Yorktown Heights, NY 10598, USA  Minh C. Tran<EMAIL_ADDRESS>IBM Quantum,
IBM T.J. Watson Research Center, Yorktown Heights, NY 10598, USA

###### Abstract

We propose an algorithm for simulating the dynamics of a geometrically local
Hamiltonian $A$ under a small geometrically local perturbation $\alpha B$. In
certain regimes, the algorithm achieves the optimal scaling and outperforms
the state-of-the-art algorithms. By moving into the interaction frame of $A$
and classically computing the Magnus expansion of the interaction-picture
Hamiltonian, our algorithm bypasses the need for ancillary qubits. In
analyzing its performance, we develop a framework to capture the quasi-
locality of the Magnus operators, leading to a tightened bound for the error
of the Magnus truncation. The Lieb-Robinson bound also guarantees the
efficiency of computing the Magnus operators and of their subsequent
decomposition into elementary quantum gates. These features make our algorithm
appealing for near-term and early-fault-tolerant simulations.

##  I Introduction

Hamiltonian simulation is one of the most exciting applications of quantum
computers. Their universality guarantees that an arbitrary unitary evolution
can be decomposed into a series of quantum gates. Leading quantum algorithms
for such decomposition include product formulas [1, 2, 3], linear combinations
of unitaries [4, 5], quantum signal processing [6, 7], and stochastic quantum
simulation algorithms [8, 9, 10]. These algorithms each have their own
advantages depending on the structure of the Hamiltonian. High-order product
formulas, for example, achieve nearly optimal scaling in gate complexity when
applied to geometrically local Hamiltonians [3, 11], while qDRIFT [9] excels
when there are multiple energy scales in the Hamiltonian.

In this paper, we consider a family of Hamiltonians $H=A+\alpha B$, where
$\alpha B$ represents a small perturbation to the unperturbed Hamiltonian $A$.
Treatment of such Hamiltonians in the perturbation theory has led to important
results, including Fermi’s golden rule and quantum chemistry beyond the Born-
Oppenheimer approximation. Promising a systematic pathway to explore
Hamiltonian dynamics under perturbation, various quantum simulation algorithms
have been introduced. A popular approach is to move into the interaction frame
of $A$ and use time-dependent quantum simulation algorithms to simulate the
resulting interaction-picture Hamiltonians [12, 13, 14, 15]. This approach
often requires ancillary qubits and controlled-evolutions, making it
challenging to implement ons near-term quantum computers. Stochastic
algorithms can bypass these resource requirements, but their asymptotic
complexity is only comparable to that of low-order deterministic algorithms.

We propose an algorithm to simulate the evolution under a geometrically local
Hamiltonian $H=A+\alpha B$. The algorithm approximates the evolution in the
interaction frame of $A$ by a truncated Magnus expansion. Each Magnus operator
is further truncated based on the Lieb-Robinson bound, allowing an efficient
decomposition into the Pauli basis. The evolution generated by the truncated
Magnus operators is then implemented using product formulas. Compared to
existing results, our algorithm displays several favorable features, including
using no ancillas, being independent of how fast $H$ changes with time, and
achieving the optimal gate complexity in certain regimes.

##  II Setup

Given a system of $n$ qubits on a one-dimensional regular lattice $\Lambda$, a
Hamiltonian $H$ is _geometrically local_ if there exists a decomposition of
$H=\sum_{i=1}^{n}H_{i}$ such that $H_{i}$ is supported entirely on a
neighborhood of diameter $\chi=O(1)$ including $i$. To set the energy unit, we
further assume that the spectral norm $\norm{H_{i}}\leq 1$ for every site $i$.
We call $\chi$ the _interaction range_ of $H$. For example, $\chi=1$ if $H$ is
a nearest-neighbor Hamiltonian.

In this work, we consider a Hamiltonian, $H=A+\alpha B$, where both $A$ and
$B$ are geometrically local Hamiltonians. Let $\chi$ be the interaction range
of $A$. Below, we illustrate our algorithm with time-independent $A$ and $B$
and later discuss the generalization to time-dependent Hamiltonians.

To set a baseline for later comparison, we briefly explain the complexity of
simulating $e^{-iHt}$ using the product formulas. For that, we divide $[0,t]$
into $r$ equal time intervals so that $U(t)=U(t/r)^{r}$ and use the first-
order product formula to approximate

$\displaystyle U(t/r)\approx\prod_{i}e^{-iA_{i}t/r}\prod_{j}e^{-i\alpha
B_{j}t/r}.$ (1)

The total error of this approximation can be upper bounded by
$\mathcal{O}\\!\left(nt^{2}/r\right)$ [3], which is independent of $\alpha$.
Given a fixed error tolerance, we choose $r\propto nt^{2}$. In each interval,
we use $\mathcal{O}\\!\left(n\right)$ gates to implement the evolutions under
$A_{i}$ and $B_{i}$, resulting in the total gate count

$\displaystyle\mathcal{G}_{\text{Trotter,1}}=\mathcal{O}\\!\left(n^{2}t^{2}\right).$
(2)

Similarly, using a $p$th-order product formula would result in a gate
complexity

$\displaystyle\mathcal{G}_{\text{Trotter,p}}=\mathcal{O}\\!\left((nt)^{1+1/p}\right),$
(3)

which approaches $\tilde{\mathcal{O}}\\!\left(nt\right)$ in the large-$p$
limit.

Throughout the paper, we use $\tilde{\mathcal{O}}\\!\left(\cdot\right)$ to
denote $\mathcal{O}\\!\left(\cdot\right)$ up to subpolynomial corrections.
This gate count scales optimally with $n$ and $t$ [11] and is independent of
$\alpha$. Meanwhile, one would expect that, if $A$ is fast-forwardable, the
gate count of simulating $H=A+\alpha B$ would scale proportionally with
$\alpha$.

In the limit of small $\alpha$, $\alpha B$ represents a perturbation to the
Hamiltonian $A$. Naturally, we can move into the interaction frame of $A$ and
write $U(t)=U_{A}(t)U_{BI}(t)$, where $U_{A}(t)=e^{-iAt}$ and $U_{BI}(t)$ is
the evolution generated by the interaction-picture Hamiltonian

$\displaystyle B_{I}(t)=\alpha U_{A}(t)^{\dagger}BU_{A}(t).$ (4)

Various quantum algorithms have been developed to simulate the interaction-
picture evolution [12, 16, 13, 15]. Notably, Ref. [12] uses the truncated
Dyson series to approximate $U_{BI}(t)$, resulting in a gate complexity
$\tilde{\mathcal{O}}\\!\left(\alpha n^{2}t\right)$ assuming the fast-
forwardability of $A$. Compared to high-order product formulas, this approach
gains a favourable dependence on $\alpha$, but at the same time scales
quadratically worst with $n$. Another avenue to approximating $B_{I}(t)$ is
the Magnus expansion [17, 18], which is defined as follows. Let
$\alpha\norm{B}t\leq 1$. Then there exists
$\Omega(t)=\sum_{j=1}^{\infty}\Omega_{j}(t)$ such that
$U_{BI}(t)=e^{\Omega(t)}$, where $\Omega_{j}(t)$ is of the $j$th-order in
$\alpha\norm{B}t$. Truncating the expansion by keeping the first $q$ terms
results in an approximation

$\displaystyle U_{BI}(t)\approx\exp\left(\sum_{j=1}^{q}\Omega_{j}(t)\right).$
(5)

For example, at $q=1$, we get $e^{\Omega_{1}(t)}$, where
$\Omega_{1}(t)=-i\int_{0}^{t}B_{I}(s)ds$ (See Appendix A for a brief review).
Through this approximation, Ref. [13] further approximates $\Omega_{1}(t)$ by
a discrete sum, allowing $e^{\Omega_{1}(t)}$ to be implemented via a technique
based on linear combinations of unitaries [4].

Here, we propose to simply evaluate the truncated Magnus expansion, e.g.,
$\Omega_{1}(t)$, _numerically_. This strategy has been previously explored in
the context of simulating $2$-local Hamiltonians [19]. In our case, the
interaction-picture Hamiltonian $B_{I}(t)$ at $t>0$ can be supported
everywhere on the lattice, presenting a challenge for both the classical
computation of the truncated Magnus expansion and its subsequent decomposition
into elementary quantum gates. We overcome these challenges by showing that,
for a small time step, the Lieb-Robinson bound guarantees an efficient Pauli
decomposition of $\Omega_{j}(t)$ for all constant $j$. In exchange for this
mild increase in the classical computational overhead, our approach saves on
additional quantum registers and achieves more favorable scaling than other
approaches in certain regimes.

##  III Main results

In this section, we illustrate the approximation of the time evolution under
$H=A+\alpha B$ using the truncated Magnus expansion. Consider that

$\displaystyle e^{-iHt}$
$\displaystyle=e^{-iAt}\mathcal{T}\exp\left(-i\alpha\int_{0}^{t}e^{iAs}Be^{-iAs}ds\right)$
$\displaystyle=e^{-iAt}e^{\Omega(t)}\approx e^{-iAt}e^{\overline{\Omega}(t)},$
(6)

where $\overline{\Omega}(t)\equiv\sum_{j=1}^{q}\Omega_{j}(t)$. For $t=O(1)$,
we truncate terms outside the Lieb-Robinson lightcone to approximate
$\overline{\Omega}(t)$ by a quasi-local operator
$\overline{\Omega}_{\text{loc}}(t)$. The locality of
$\overline{\Omega}_{\text{loc}}(t)$ allows us to efficiently estimate its
Pauli decomposition. Finally, we approximate the evolution
$e^{\overline{\Omega}_{\text{loc}}(t)}$ using a $p$th-order product formula.
In the following analysis, we use $V_{q,p}(t)$ to denote the resulting
unitary. For large $t$, we divide it into $r$ time intervals and repeat the
decomposition for each interval as follows:

$\displaystyle U(t)=\left[e^{-iAt/r}V_{q,p}(t/r)\right]^{r}.$ (7)

There are three contributions to the error of this simulation: i) the
truncation of the Magnus expansion, ii) the Lieb-Robinson light-cone
approximation of $\overline{\Omega}$ by $\overline{\Omega}_{\text{loc}}$, iii)
and the product formula for decomposing $e^{\overline{\Omega}_{\text{loc}}}$.
We first consider the case $t=O(1)$ in the following analysis.

The Magnus expansion truncation error can be bounded using standard
techniques, giving

$\displaystyle\|e^{\Omega(t)}-e^{\overline{\Omega}(t)}\|=\mathcal{O}\\!\left((\alpha
nt)^{q+1}\right).$ (8)

Exploiting the locality of the Hamiltonian, we arrive at a tighter bound.

######  Theorem 1.

Let $A$ be a geometrically local Hamiltonian and
$B=\sum_{X\subset\Lambda}B_{X}$ be a $k$-local Hamiltonian such that $B_{X}$
is supported on a subset $X$ of at most $k$ qubits. Let
$d\equiv\max_{i\in\Lambda}\sum_{X:i\in X}\norm{B_{X}}$ denote the effective
interaction degree of $B$. Let
$e^{\Omega(t)}=\mathcal{T}[e^{\int_{0}^{t}B_{I}(\tau)d\tau}]$, where
$B_{I}(\tau)=-i\alpha e^{iA\tau}Be^{-iA\tau}$. Let
$\overline{\Omega}(t)=\sum_{j=1}^{q}\Omega_{j}(t)$ be the truncated Magnus
expansion of order $q$. We have

$\displaystyle\|e^{\Omega(t)}-e^{\overline{\Omega}(t)}\|\leq\mathcal{O}(n(\alpha
dt)^{q+1})\leavevmode\nobreak\ .$ (9)

The theorem applies to a more general class of Hamiltonians $B$ and may be of
independent interest beyond quantum simulation. For the rest of the
manuscript, we consider $B$ to be geometrically local. Our bound in Eq. (9)
reduces the error dependence on the number of qubits from
$\mathcal{O}\\!\left(n^{q+1}\right)$ to $\mathcal{O}\\!\left(n\right)$, a
quadratic improvement even at $q=1$. We present the detailed proof of Theorem
1 in Appendix C and use the $q=1$ example to illustrate the idea.

If $\overline{\Omega}(t)=\Omega_{1}(t)$, the leading contribution to the
truncation error comes from

$\displaystyle\norm{\Omega_{2}(t)}$
$\displaystyle\leq\frac{1}{2}\int_{0}^{t}d\tau\int_{0}^{\tau}ds\norm{\left[B_{I}(\tau),B_{I}(s)\right]}.$
(10)

The integrand can be bounded by
$\alpha^{2}\norm{B}^{2}=\mathcal{O}\\!\left(\alpha^{2}n^{2}\right)$, resulting
in the bound in Eq. 8. Instead, if $B_{I}(\tau)$ and $B_{I}(s)$ _were_ both
geometrically local, only Pauli terms in $B_{I}(\tau)$ of total norm
$\mathcal{O}\\!\left(1\right)$ may noncommute with a given term in $B_{I}(s)$.
Therefore, the commutator norm would scale as $\mathcal{O}\\!\left(n\right)$,
resulting in a bound of $\mathcal{O}\\!\left(n(\alpha t)^{2}\right)$ for Eq.
10.

Although $B_{I}(\tau)$ and $B_{I}(s)$ are not strictly geometrically local due
to the evolution under $A$, the geometric locality of $A$ and the Lieb-
Robinson bound ensure that they inherit a similar structure. We develop a
framework to capture the structure in Appendix D, and show that the Magnus
operators admit this structure in Appendix E, and finally apply these results
to bound the Magnus truncation error in Appendix C.

Next, we use the Lieb-Robinson bound for the classical computation of
$\overline{\Omega}(t)$. For $q=1$, this operator can be expanded into a sum
over $B_{X}$:

$\displaystyle\overline{\Omega}(t)=\sum_{X}\int_{0}^{t}e^{iAs}B_{X}e^{-iAs}ds.$
(11)

For each $X$, the evolved operator $e^{iAs}B_{X}e^{-iAs}$ is mostly supported
within a light cone of radius $O(t)=O(1)$ around $X$. In particular, let
$A_{X,R}$ be the operator constructed from $A$ by dropping terms supported
outside a radius $R$ from the set $X$. The Lieb-Robinson bound implies [20,
21]

$\displaystyle\norm{e^{iAs}B_{X}e^{-iAs}-e^{iA_{X,R}s}B_{X}e^{-iA_{X,R}s}}=\mathcal{O}\\!\left(te^{-R/\chi}\right),$
(12)

where $\chi$ is the interaction range of $A$. Therefore, the local version of
the truncated Magnus operator

$\displaystyle\overline{\Omega}_{\text{loc}}(t)\equiv\sum_{X}\int_{0}^{t}e^{iA_{X,R}s}B_{X}e^{-iA_{X,R}s}ds$
(13)

only differs from $\overline{\Omega}(t)$ by

$\displaystyle\norm{\overline{\Omega}(t)-\overline{\Omega}_{\text{loc}}(t)}=\mathcal{O}\\!\left(nte^{-R/\chi}\right).$
(14)

For a large time, we apply the construction in $r$ time steps, each of size
$t/r$, resulting in the same bound. Choosing $R=\chi\log nt+\vartheta$ for
some constant $\vartheta$ ensures that this error remains a small constant.
Note that $\overline{\Omega}_{\text{loc}}(t)$ is a sum of terms that are
supported on at most $\mathcal{O}\\!\left(R\right)=\mathcal{O}(\log nt)$
qubits. Therefore, the complexity of computing
$\overline{\Omega}_{\text{loc}}(t)$ scales polynomially with $n$ and $t$. In
Appendix E, we generalize this construction to approximate $q$th-order Magnus
operators by an operator $\overline{\Omega}_{\text{loc}}(t)$ consisting of $n$
terms, each supported on at most $R+q$ qubits. The error of this approximation
is

$\displaystyle\norm{\overline{\Omega}(t)-\overline{\Omega}_{\text{loc}}(t)}=\mathcal{O}\\!\left(nte^{-R/8\chi}\right),$
(15)

which is identical to Eq. 14 except for a mildly slower decay.

Finally, we numerically find the Pauli decomposition
$\overline{\Omega}_{\text{loc}}(t)$ and approximate
$e^{\overline{\Omega}_{\text{loc}}(t)}$ using a $p$th-order product formula.
If $\alpha t=\mathcal{O}\\!\left(1\right)$, the effective interaction degree
of $\overline{\Omega}_{\text{loc}}(t)$ is $\mathcal{O}\\!\left(\alpha
t\right)$ due to the first-order Magnus operator being the leading term (See
Appendix E). Therefore, the error of using the product formula is [3]

$\displaystyle\norm{e^{\overline{\Omega}_{\text{loc}}(t)}-V_{q,p}}=\mathcal{O}\\!\left(n(\alpha
t)^{p+1}\right).$ (16)

Combining Eqs. 9, 15 and 16, the error of our algorithm for $t=O(1)$ is

$\displaystyle\mathcal{O}\\!\left(n(\alpha t)^{q+1}+nte^{-R/8\chi}+n(\alpha
t)^{p+1}\right).$ (17)

For large $t$, we repeat the approximation for $r$ time intervals, resulting
in the total error

$\displaystyle\mathcal{O}\\!\left(\frac{n(\alpha
t)^{q+1}}{r^{q}}+nte^{-R/8\chi}+\frac{n(\alpha t)^{p+1}}{r^{p}}\right).$ (18)

Given a constant error tolerance, we choose $R=8\chi\log(nt)+{\text{const.}}$
and

$\displaystyle r\propto\max\left\\{n^{1/q}(\alpha t)^{1+1/q},n^{1/p}(\alpha
t)^{1+1/p}\right\\}.$ (19)

At this point, it is obvious that the errors of the Magnus truncation and the
product formula in our algorithm are comparable given the same approximation
orders. To simplify the expression, we choose $p=q$ in the rest of the
analysis.

To estimate the gate count, $\overline{\Omega}_{\text{loc}}(t)$ contains $n$
terms, each supported on at most $R+q$ qubits. Since $R$ scales
logarithmically with $nt$, the number of Pauli terms scales at worst
polynomially with $nt$, i.e., there exists a constant $\gamma$ such that the
number of Pauli terms is $\mathcal{O}(n(nt)^{\gamma})$. Here, $\gamma\leq
8\chi\log 4$ depends on the details of the lattice $\Lambda$ and the locality
of the Hamiltonian $A$. Therefore, the total gate count for simulating $r$
instances of $V_{q,p}(t/r)$ is

$\displaystyle\mathcal{G}_{\text{Mag},p}^{B}=\tilde{\mathcal{O}}\\!\left((nt)^{\gamma}(\alpha
nt)^{1+1/p}\right).$ (20)

Compared to that of the product formula [Eq. 3], the gate count of our
algorithm scales nearly linearly with $\alpha$. We summarize the result and
compare to the state of the art in Table 1.

If $A$ is not fast-forwardable, we also have to simulate $r$ instances of
$e^{-iAt/r}$ in addition to $V_{p,p}(t/r)$. In the above analysis, we have not
accounted for the error of simulating these evolutions. To control such error,
we further divide each interval of length $t/r$ into $r^{\prime}$ equal
subintervals and use the same $p$th-order product formula to approximate
$e^{-iAt/(rr^{\prime})}$. The total error after $rr^{\prime}$ intervals is
$\mathcal{O}\\!\left(nt^{p+1}/{(rr^{\prime})^{p}}\right)$, informing the
choice of

$\displaystyle
r^{\prime}\propto\frac{n^{1/p}t^{1+1/p}}{r}\propto\frac{1}{\alpha^{1+1/p}}.$
(21)

The total gate count for simulating $r$ instances of $e^{-iAt/r}$ is given by

$\displaystyle\mathcal{G}_{\text{Mag},p}^{A}=\tilde{\mathcal{O}}\\!\left(\\!(nt)^{1+1/p}\\!\right).$
(22)

Adding Eq. 22 to Eq. 20, we arrive at the total gate count of our algorithm,
which is asymptotically identical to that of the product formulas. The
advantage of using our algorithm is overshadowed by the cost of simulating $A$
in the asymptotic limit, but the effective reduction in the circuit depth can
be significant for near-term simulations.

Algorithm | Gate count
---|---
Product formulas [3] | $\mathcal{O}({(nt)^{1+o(1)}})$
qDRIFT [9] | $\mathcal{O}\\!\left(\alpha^{2}n^{2}t^{2}\right)$
Dyson series [12] | $\tilde{\mathcal{O}}\\!\left(\alpha n^{2}t\right)$
qHOP [13] | $\tilde{\mathcal{O}}\\!\left(\alpha n^{3}t^{2}\right)$
Cont. qDRIFT [16] | $\tilde{\mathcal{O}}\\!\left(\alpha^{2}n^{2}t^{2}\right)$
This paper | $\tilde{\mathcal{O}}\\!\left((nt)^{\gamma}(\alpha nt)^{1+o(1)}\right)$
Table 1: The gate count of several algorithms for simulating the evolution
$e^{-iHt}$ generated by the Hamiltonian $H=A+\alpha B$, as considered in this
paper. Except for the $p$th-order product formulas and qDRIFT, these
algorithms are applied to the evolution $U_{BI}(t)$ in the interaction frame.
It is assumed that the cost of implementing $e^{-iAt}$ is negligible for all
$t$. In the scenario where $\gamma$ approaches 0, our algorithm’s gate count
scales as $\tilde{\mathcal{O}}\\!\left(\alpha nt\right)$, thereby
outperforming the other algorithms.

To benchmark our algorithm on near-term simulations, we consider the
disordered XY model $H=A+\alpha B$ with $A=\sum_{i=1}^{n}Z_{i}$ and
$B=\sum_{i=1}^{n-1}r_{i}Y_{i}Y_{i+1}+s_{i}X_{i}X_{i+1}-\sum_{i=1}^{n}u_{i}X_{i}$
for $n=12$ qubits. Here, $r_{i},s_{i}$, and $u_{i}$ are sampled uniformly at
random in the range $(-1,1)$ for all $i$. In Fig. 1, we plot the approximation
error in simulating $e^{-iHt}$ for a total time $t=3$ at different $\alpha$
using our algorithm, the first-order, and the second-order product formulas.
To ensure the three algorithms have the same gate count, we use time intervals
of different sizes (1, 0.5, and 0.75 respectively) for different algorithms.

Figure 1 shows that for small values of $\alpha$, our algorithm performs
significantly (one-order magnitude) better than the first- and the second-
order product formulas. The dependence on $\alpha$ of the algorithms also
matches the theoretical prediction. In particular, errors due to the first-
and second-order product formulas scale linearly in $\alpha$, whereas the
error of our algorithm scales quadratically.

Figure 1: Comparison of the errors from our algorithm (circles), first-
(“PF1”, squares) and second-order (“PF2”, diamonds) product formulas for
simulating the disordered XY model of $n=12$ qubits for time $t=3$ at
different values of $\alpha$. Points of the same color represent different
random coefficients of the Hamiltonian. The errors for the first- and second-
order product formulas scale linearly with $\alpha$, whereas the error of our
algorithm scales quadratically.

##  IV Discussion

We have introduced an algorithm for simulating the time evolution in the
interaction frame of a geometrically local Hamiltonian. The algorithm is
ancilla-free and takes advantage of the locality of the Hamiltonian,
particularly via a tightened bound on the error of truncating the Magnus
expansion.

In the limit $\chi\\!\rightarrow 0$, $A$ contains only terms supported on
single qubits. The evolution $e^{-iAt}$ is fast-forwardable and the resulting
lightcone does not grow with time ($\gamma\\!\\!\\!\rightarrow\\!0$). The
complexity of our algorithm reduces to $\tilde{\mathcal{O}}\\!\left((\alpha
nt)^{1+o(1)}\right)$, which is optimal [11] and outperforms both
$\tilde{\mathcal{O}}\\!\left(\alpha^{o(1)}(nt)^{1+o(1)}\right)$ of product
formulas and $\tilde{\mathcal{O}}\\!\left(\alpha n^{2}t\right)$ of algorithms
based on the Dyson series.

In general, the complexity of our algorithm depends on the constant $\gamma$,
which can be as large as $4\chi\log 4$. This value assumes that
$\overline{\Omega}_{\text{loc}}(t)$ contains all exponentially many possible
Pauli strings supported inside the lightcone, preventing applications of our
algorithm to more general $B$, e.g. $k$-local Hamiltonians. Practically,
$\gamma$ depends on the Hamiltonian of interest and can be much smaller than
this upper bound, especially in the presence of symmetries.

For example, consider $A\propto\sum_{i}X_{i}X_{i+1}+Y_{i}Y_{i+1}$ and $B$
contains Pauli strings of the form $\sigma_{i}Z_{i+1}\dots Z_{j-1}\sigma_{j}$,
where $\sigma_{i},\sigma_{j}$ can be either $X,$ or $Y$ supported on qubits
$i,j$. The evolution of $B$ under $A$ preserves this structure, limiting the
number of possible Pauli strings in $\overline{\Omega}_{\text{loc}}$ to scale
polynomially with the size of the lightcone, which is logarithmic in $n$ and
$t$. In this case, $\gamma$ can be chosen arbitrarily close to 0. Such $A$ and
$B$ are the products of applying the Jordan-Wigner transformation on a
nearest-neighbor and a two-local fermionic Hamiltonian, respectively. The
argument for $\gamma\approx 0$ can be similarly generalized to $k$-local
fermionic Hamiltonians $B$.

If $\Lambda$ is a general $D$-dimensional lattice, the number of Pauli strings
in the decomposition of the Magnus operators can be exponential in $R^{D}$,
which is no longer polynomial in $n,t$ for $D>1$. This number again
corresponds to the worst case and can be significantly smaller, e.g. under the
symmetry constraints discussed earlier. For such $D$-dimensional systems, our
algorithm still applies.

While we assume $A$ and $B$ are time-independent in our analysis, our
algorithm can be immediately generalized to simulate time-dependent
Hamiltonians $A(t)$ and $B(t)$. In the special case $A(t)=0$, our algorithm
can be applied to simulate a time-dependent geometrically local Hamiltonian
$B(t)$. Since the Magnus operators are computed classically, the gate count of
our algorithm is independent of how fast $B(t)$ changes with time—a feature
similar to continuous qDRIFT [16]. For reference, the complexity of product
formulas, qHOP [13], and algorithms based on the Dyson series [12] would scale
either polynomially or logarithmically with the derivative $\norm{dB/dt}$.

Note added. While finalizing our paper, we noticed two independent papers on
arXiv studying the Hamiltonian simulation using the Magnus expansion [22, 23].

###### Acknowledgements.

We acknowledge helpful discussions with Sergey Bravyi, William Kirby, Antonio
Mezzacapo, Oles Shtanko, and Kristan Temme.

Appendix for “Hamiltonian Simulation in the Interaction Picture using the
Magnus Expansion”

##  Appendix A Preliminaries

In this section, we recall definition and prior results relevant for our
paper. We first define the notion of Hamiltonian simulation in the interaction
picture. We then summarize properties and relevant formulas for the Magnus
expansion. For completeness, we also summarize Duhamel’s principle for
Hamiltonian simulation and commutator-dependent bounds for first-order and
second-order Trotter simulation. We point readers to [24] for more details on
Hamiltonian simulation methods.

Let

$\displaystyle H=A+\alpha B$ (23)

where $A$ is geometrically local Hamiltonians and $B$ is a $k$-local
Hamiltonian with an effective interaction degree $d$. Let
$B=\sum_{X\subset\Lambda}B_{X}$ be a decomposition into a sum of $L$ terms,
each supported on a subset $X$ of at most $k$ qubits. We call

$\displaystyle d\equiv\max_{i\in\Lambda}\sum_{X:i\in X}\norm{B_{X}}$ (24)

the effective interaction degree of $B$, which upper bounds the norms of the
terms supported on a qubit. For example, $d=\mathcal{O}{(1)}$ if B is also
geometrically local and $d=\mathcal{O}{(n^{k-1})}$ in general.

We recall the definition of geometrical locality as defined in Section II. In
particular, a Hamiltonian $H$ for $n$ qubits defined on a $1$-dimensional
regular lattice is geometrically local if there is a decomposition
$H=\sum_{i=1}^{n}H_{i}$ such that $H_{i}$ is supported entirely on a
neighborhood of constant radius around $i$. In our analysis, we assume that
$\alpha<1$.

In the interaction frame of $A$, the time evolution under $H$ can be expressed
as

$\displaystyle
e^{-iHt}=e^{-iAt}\mathcal{T}\left[e^{\int_{0}^{t}B_{I}(\tau)\mathrm{d}\tau}\right]\leavevmode\nobreak\
,$ (25)

where $\mathcal{T}\left[e^{\int_{0}^{t}B_{I}(\tau)\mathrm{d}\tau}\right]$ is
the time-ordered expeonential with

$\displaystyle B_{I}(\tau)=-i\alpha e^{iA\tau}Be^{-iA\tau},$ (26)

with $\norm{B_{I}(\tau)}=\norm{B}\leq nd$ which follows trivially from the
definition of $d$ and a sum over the qubit index.

Before providing an approximation of
$\mathcal{T}\left[e^{\int_{0}^{t}B_{I}(\tau)\mathrm{d}\tau}\right]$, we
summarize important properties of the Magnus expansion. The Magnus expansion
is defined as follows.

######  Lemma 1 (Magnus expansion [25, 17, 18]).

Let $\mathcal{A}(\tau)$ be a continuous operator-valued function defined for
$0\leq\tau\leq t$ such that $\|\mathcal{A}(\tau)\|\leq 1$. Then, we get that

$\displaystyle\mathcal{T}\exp\left\\{\int_{0}^{t}\mathrm{\leavevmode\nobreak\
d}\tau\mathcal{A}(\tau)\right\\}=\exp\left(\sum_{j=1}^{\infty}\Omega_{j}(t)\right)\equiv\exp(\Omega(t))\leavevmode\nobreak\
,$ (27)

where

$\displaystyle\Omega_{j}(t)=\frac{1}{j^{2}}\sum_{\sigma}(-1)^{d_{b}}\frac{1}{\left(\begin{array}[]{c}j-1\\\
d_{b}\end{array}\right)}\int_{0}^{t}d\tau_{1}\ldots\int_{0}^{d\tau_{j-1}}d\tau_{j}\left[\mathcal{A}\left(\tau_{1}\right),\ldots\left[\mathcal{A}\left(\tau_{j-1}\right),\mathcal{A}\left(\tau_{j}\right)\right]\ldots\right],$
(30)

the sum is being taken over all permutations $\sigma$ of $\\{1,\dots,j\\}$,
and $d_{b}$ is the number of descents, i.e., pairs of consecutive numbers
$\sigma_{k},\sigma_{k+1}$ for $k=\\{1,\dots,j-1\\}$ such that
$\sigma_{k}>\sigma_{k+1}$ in the permutation.

From Lemma 1, the first two $\Omega_{j}(t)$ is given by:

$\displaystyle\Omega_{1}(t)$
$\displaystyle=\int_{0}^{t}\mathrm{\leavevmode\nobreak\
d}\tau\mathcal{A}(\tau)\leavevmode\nobreak\ ,$ (31)
$\displaystyle\Omega_{2}(t)$
$\displaystyle=\frac{1}{2}\int_{0}^{t}\mathrm{\leavevmode\nobreak\
d}\tau_{1}\int_{0}^{\tau_{1}}\mathrm{\leavevmode\nobreak\
d}\tau_{2}\left[\mathcal{A}(\tau_{1}),\mathcal{A}(\tau_{2})\right]\leavevmode\nobreak\
.$ (32)

For the problem of interest in this paper, $\mathcal{A}(\tau)$ gets replaced
by $B_{I}(\tau)$ as defined in Eq. (26) for the Hamiltonian considered in
Eq.(23).

We now recall an alternate representation for the Magnus expansion which we
employ in our proofs. In particular, we will invoke the following lemma in
Section C.

######  Lemma 2 (Magnus expansion [26]).

Let $Y^{\prime}(t)=A(t)Y(t)$ with $Y(0)=I$. Then $Y(t)$ can be written as
follows:

$\displaystyle Y(t)=\exp(\Omega(t))\leavevmode\nobreak\ ,$ (33)

where

$\displaystyle\frac{d\Omega}{dt}=\sum_{n=0}^{\infty}\frac{B_{n}}{n!}\operatorname{ad}_{\Omega}^{n}\leavevmode\nobreak\
,$ (34)

with $B_{n}$ denoting the Bernoulli numbers. Then the integration over
$d\Omega/dt$ leads to

$\displaystyle\Omega(t)=\int_{0}^{t}A\left(t_{1}\right)dt_{1}-\frac{1}{2}\int_{0}^{t}\left[\int_{0}^{t_{1}}A\left(t_{2}\right)dt_{2},A\left(t_{1}\right)\right]dt_{1}+\cdots\leavevmode\nobreak\
.$ (35)

Moreover, derivatives of $\Omega_{j}$ have the following form:

$\displaystyle\Omega_{1}^{\prime}$ $\displaystyle=A\leavevmode\nobreak\ ,$
(36) $\displaystyle\Omega_{2}^{\prime}$
$\displaystyle=-\frac{1}{2}[\Omega_{1},A]\leavevmode\nobreak\ ,$ (37)
$\displaystyle\Omega_{3}^{\prime}$
$\displaystyle=-\frac{1}{2}[\Omega_{2},A]+\frac{1}{12}[\Omega_{1},[\Omega_{1},A]]\leavevmode\nobreak\
,$ (38) $\displaystyle\Omega_{4}^{\prime}$
$\displaystyle=-\frac{1}{2}[\Omega_{3},A]+\frac{1}{12}[\Omega_{2},[\Omega_{1},A]]+\frac{1}{12}[\Omega_{1},[\Omega_{2},A]]\leavevmode\nobreak\
,$ (39) $\displaystyle\Omega_{j}^{\prime}$
$\displaystyle=\sum_{k=1}^{j-1}\frac{B_{k}}{k!}S_{j}^{(k)},\quad j\geq
2\leavevmode\nobreak\ ,$ (40)

where

$\displaystyle
S_{j}^{(k)}=\sum\left[[\Omega_{1},[\dots[\Omega_{i_{k}},A]\dots]\right]\qquad(i_{1}+\dots+i_{k}=j-1)\leavevmode\nobreak\
.$ (41)

Duhalmel’s principle. We recall Duhamel’s principle for Hamiltonian simulation
which is useful in the error analysis. Let $U(t),V(t)\in\mathbb{C}^{d\times
d}$. Let $H\in\mathbb{C}^{d\times d}$ denote a Hermitian operator and
$E(t)\in\mathbb{C}^{d\times d}$ denote an arbitrary matrix. If

$\displaystyle i\partial_{t}U(t)$ $\displaystyle=HU(t)\leavevmode\nobreak\ ,$
(42) $\displaystyle i\partial_{t}V(t)$
$\displaystyle=HV(t)+E(t)\leavevmode\nobreak\ ,$ (43) $\displaystyle U(0)$
$\displaystyle=V(0)=I\leavevmode\nobreak\ ,$ (44)

then

$\displaystyle V(t)=U(t)-i\int_{0}^{t}U(t-s)E(s)ds,$ (45)

which implies that

$\displaystyle\|U(t)-V(t)\|\leq\int_{0}^{t}ds\|E(s)\|\leavevmode\nobreak\ .$
(46)

Here, we used the data-processing inequality and the triangle inequality.

Product formulas and Trotter error. Using the Duhamel’s principle stated
above, one can derive Trotter error bounds that depends on the nested
commutators between the terms of a Hamiltonian $H$. In particular, let
$H=\sum_{i=1}^{L}H_{i}$. Let $V_{p}$ denote the $p$-th order product formula,
which is defined recursively as follows:

$\displaystyle V_{1}(t)$
$\displaystyle=\prod_{j=1}^{L}\exp\left(-iH_{j}t\right)\leavevmode\nobreak\ ,$
(47) $\displaystyle V_{2}(t)$
$\displaystyle=\prod_{j=1}^{L}\exp\left(-iH_{j}/2\right)\prod_{j=L}^{1}\exp\left(-iH_{j}/2\right)\leavevmode\nobreak\
,$ (48) $\displaystyle V_{p}(t)$
$\displaystyle=V_{2p-2}(m_{p}t)^{2}V_{2p-2}((1-4m_{p})t)V_{2p-2}(m_{p}t)^{2}\leavevmode\nobreak\
,$ (49)

where $m_{p}=\frac{1}{4-4^{1/(2p-1)}}$.

The number of $r$ Trotter steps for a $p$-th order product formula implies
$(V_{p}(t/r))^{r}$, i.e., one implements $V_{p}(t/r)$ for time $t/r$ and
repeats that $r$ times.

Let $V_{1}$ and $V_{2}$ denote the first-order and second-order product
formula, respectively, and let $U=e^{-iHt}$. Then we get that [27]

$\displaystyle\|U-V_{1}\|$
$\displaystyle\leq\frac{t^{2}}{2}\sum_{k=1}^{L}\left\|\left[\sum_{j=k+1}^{L}H_{j},H_{k}\right]\right\|\leavevmode\nobreak\
$ (50) $\displaystyle\|U-V_{2}\|$
$\displaystyle\leq\frac{t^{3}}{12}\sum_{i=1}^{L}\left\|\left[\sum_{k=i+1}^{L}H_{k},\left[\sum_{j=i+1}^{L}H_{j},H_{i}\right]\right]\right\|+\frac{t^{3}}{24}\sum_{i=1}^{L}\left\|\left[H_{i},\left[H_{i},\sum_{j=i+1}^{L}H_{j}\right]\right]\right\|$
(51)

For the case when $H=A+\alpha B$, then we have

$\displaystyle\|U-V_{1}\|$ $\displaystyle\leq\frac{\alpha
t^{2}}{2}\|[A,B]\|\leavevmode\nobreak\ ,$ (52) $\displaystyle\|U-V_{2}\|$
$\displaystyle\leq\alpha\left(\frac{t^{3}}{12}\|[A[A,B]]\|+\frac{\alpha
t^{3}}{24}\|[B,[B,A]]\|\right)\leavevmode\nobreak\ .$ (53)

Lieb-Robinson bound. Let $O_{X}$ and $O_{Y}$ operators on $n$ qubits supported
on $X$ and $Y$ sets of qubits, and $H$ denote a geometrically local
Hamiltonian with an interaction range $\chi$. Then the following inequality
holds [20]:

$\displaystyle\|[e^{iHt}O_{Y}e^{-iHt},O_{X}]\|\leq
c\min\\{|X|,|Y|\\}\|O_{X}\|\|O_{Y}\|\left(e^{v\absolutevalue{t}}-1\right)\exp\left(-\frac{R}{\chi}\right),$
(54)

where $R$ is the minimum distance between $X$ and $Y$, i.e., $\min_{i\in
X,j\in Y}d_{ij}$, and $c$ and $v$ are constants that depend on $\chi$ and the
effective interaction degree of $H$ which is also a constant.

In the next section, we derive a naive bound on the truncated Magnus
approximation.

##  Appendix B Proof of Eq. (8)

In this section we establish a bound on the truncation error in the Magnus
expansion. In particular, we highlight why exploiting the commutator structure
is necessary to establish a tight bound. As a first step, we approximate
$\Omega$ with $\overline{\Omega}=\sum_{k=1}^{q}\Omega_{k}$. Note that in our
problem, $\Omega_{k}$ are anti-Hermitian. For simplicity, we do not write the
time dependence of $\Omega$ explicitly.

Since $\norm{\Omega_{k}}\leq 2^{k}\alpha^{k}\norm{B}^{k}t^{k}\leq(2n\alpha
dt)^{k}$, we have

$\displaystyle\|e^{\Omega}-e^{\overline{\Omega}}\|$
$\displaystyle\leq\|\Omega-\overline{\Omega}\|=\sum_{k=q+1}^{\infty}\|\Omega_{k}\|\leq\sum_{k=q+1}^{\infty}(2n\alpha
dt)^{k}\leq 2(2n\alpha dt)^{q+1},$ (55)

assuming $n\alpha dt\leq 1/4$. The first inequality follows from the Duhamel’s
principle as discussed in Eq. (46). We recall it for completeness. Let $\tau$
denote a dummy variable. Note that $\partial_{\tau}(e^{\Omega\tau})=\Omega
e^{\Omega\tau}$ and $\partial_{\tau}e^{\overline{\Omega}\tau}=\Omega
e^{\overline{\Omega}\tau}+(\overline{\Omega}-\Omega)e^{\overline{\Omega}\tau}$.
Then we get

$\displaystyle e^{\overline{\Omega}\tau}$
$\displaystyle=e^{\Omega\tau}+\int_{0}^{1}e^{\Omega(\tau-s)}(\overline{\Omega}-\Omega)e^{\overline{\Omega}\tau}d\tau$
(56)

which implies that

$\displaystyle\|e^{\Omega}-e^{\overline{\Omega}}\|$
$\displaystyle=\left\|\int_{0}^{1}d\tau
e^{\Omega(\tau-s)}(\Omega-\overline{\Omega})e^{\overline{\Omega}\tau}\right\|\leq\int_{0}^{1}d\tau\left\|\Omega-\overline{\Omega}\right\|=\|\Omega-\overline{\Omega}\|\leavevmode\nobreak\
.$ (57)

We conclude from Eq. 55 that the simulation error is of order $\alpha^{q+1}$.
However, the $n$ dependence makes the bound loose as it grows with $n^{q+1}$.
In the next sections, we show how to improve the bound to a linear dependence
on $n$ by exploiting the locality $B_{I}(\tau)$.

##  Appendix C Magnus truncation error

In this section, we prove Theorem 1 by analyzing the simulation error that
takes account of the commutator between different $\Omega_{i}$. First, we note
that

$\displaystyle\norm{e^{-\overline{\Omega}}e^{\Omega}-I}=\norm{\int_{0}^{1}e^{-\overline{\Omega}\tau}(\Omega-\overline{\Omega})e^{\Omega\tau}d\tau}\leq\int_{0}^{1}\norm{\Omega-\overline{\Omega}}d\tau\leq
2(2n\alpha dt)^{q+1},$ (58)

where the last inequality follows from Eq. 55. Equation 58 implies that
$\partial_{t}(e^{-\overline{\Omega}}e^{\Omega})=\mathcal{O}\\!\left(\alpha^{q+1}\right)$
and, specifically,

$\displaystyle\partial_{\alpha}^{k}\partial_{t}(e^{-\overline{\Omega}}e^{\Omega})=0$
(59)

for all $k\leq q$.

In the following discussion, we shall prove a bound on
$\norm{\partial_{t}(e^{-\overline{\Omega}}e^{\Omega})}$ that scales linearly
with $n$ and use that to bound the Magnus truncation error.

The derivative of $e^{\Omega}$ is $\partial_{t}e^{\Omega}=B_{I}e^{\Omega}$.
Recall that for simplifying the notation, we did not write the dependence of
$\Omega$ on $t$ explicitly. For $\overline{\Omega}(t)$ we write the $t$
dependence and compute the derivative of $e^{-\overline{\Omega}(t)}$ as
follows:

$\displaystyle\partial_{t}e^{-\overline{\Omega}(t)}$
$\displaystyle=\lim_{dt\to
0}\frac{e^{-\overline{\Omega}(t+dt)}-e^{-\overline{\Omega}(t)}}{dt}$ (60)
$\displaystyle=\lim_{dt\to
0}\frac{e^{-\overline{\Omega}(t)-\overline{\Omega}(t)^{\prime}dt+O(dt^{2})}-e^{-\overline{\Omega}(t)}}{dt}$
(61) $\displaystyle=\lim_{dt\to 0}\lim_{\tau\to
1}\frac{e^{-\overline{\Omega}(t)\tau}\left(\mathcal{T}\left[\exp\left(-i\int_{0}^{\tau}dxe^{\overline{\Omega}(t)x}(\overline{\Omega}(t)^{\prime}dt)e^{-\overline{\Omega}(t)x}\right)\right]-I\right)}{dt}$
(62)
$\displaystyle=-e^{-\overline{\Omega}(t)}\left(\int_{0}^{1}dxe^{\overline{\Omega}(t)x}\overline{\Omega}(t)^{\prime}e^{-\overline{\Omega}(t)x}\right)\leavevmode\nobreak\
,$ (63)

which follows from the first-order term in the expansion of the time-ordered
exponential, and also by taking the limit $dt\to 0$. Using this identity, we
have

$\displaystyle\partial_{t}(e^{-\overline{\Omega}}e^{\Omega})$
$\displaystyle=(\partial_{t}e^{-\overline{\Omega}})e^{\Omega}+e^{-\overline{\Omega}}\partial_{t}e^{\Omega}$
(64)
$\displaystyle=-e^{-\overline{\Omega}(t)}\left(\int_{0}^{1}dxe^{\overline{\Omega}(t)x}\overline{\Omega}(t)^{\prime}e^{-\overline{\Omega}(t)x}\right)e^{\Omega}+\alpha
e^{-\overline{\Omega}}B_{I}e^{\Omega}$ (65)
$\displaystyle=-e^{-\overline{\Omega}(t)}\left(\sum_{j=1}^{q}\int_{0}^{1}dxe^{\overline{\Omega}(t)x}{\Omega}_{j}(t)^{\prime}e^{-\overline{\Omega}(t)x}-\int_{0}^{1}dx\Omega_{1}(t)^{\prime}\right)e^{\Omega(t)},$
(66)

where we used the fact that
$e^{\Omega}=\mathcal{T}\left[e^{\int_{0}^{t}B_{I}(\tau)\mathrm{d}\tau}\right]$
and therefore, $\partial_{t}e^{\Omega}=B_{I}(t)e^{\Omega}$. The second term in
the last equality follows from the fact that $B_{I}(t)=\Omega_{1}(t)^{\prime}$
and we trivially added an integration over $x$.

Terms inside the bracket in the aforementioned equation can be simplified as
follows:

$\displaystyle\sum_{j=1}^{q}\int_{0}^{1}dxe^{\overline{\Omega}(t)x}{\Omega}_{j}(t)^{\prime}e^{-\overline{\Omega}(t)x}-\int_{0}^{1}dx\Omega_{1}(t)^{\prime}$
$\displaystyle=\sum_{j=2}^{q}\int_{0}^{1}dxe^{\overline{\Omega}(t)x}{\Omega}_{j}(t)^{\prime}e^{-\overline{\Omega}(t)x}+\left(\int_{0}^{1}dxe^{\overline{\Omega}(t)x}{\Omega}_{1}(t)^{\prime}e^{-\overline{\Omega}(t)x}-\int_{0}^{1}dx\Omega_{1}(t)^{\prime}\right)$
(67)
$\displaystyle=\sum_{j=2}^{q}\int_{0}^{1}dxe^{\overline{\Omega}x}{\Omega}_{j}^{\prime}e^{-\overline{\Omega}x}+\int_{0}^{1}dx\int_{0}^{x}dye^{\overline{\Omega}y}[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]e^{-\overline{\Omega}y}\leavevmode\nobreak\
.$ (68)

Here, we used the following argument:

$\displaystyle
e^{\overline{\Omega}x}\Omega_{1}^{{}^{\prime}}e^{-\overline{\Omega}x}-\Omega_{1}^{{}^{\prime}}$
$\displaystyle=\int_{0}^{x}dy\partial_{y}e^{\overline{\Omega}y}\Omega_{1}^{{}^{\prime}}e^{-\overline{\Omega}y}=\int_{0}^{x}dye^{\overline{\Omega}y}[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]e^{-\overline{\Omega}y}$
(69)

We analyze
$\int_{0}^{x}dye^{\overline{\Omega}y}[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]e^{-\overline{\Omega}y}$
further by adding and subtracting
$\int_{0}^{x}dy[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]$. Then

$\displaystyle\int_{0}^{x}dye^{\overline{\Omega}y}[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]e^{-\overline{\Omega}y}$
$\displaystyle=\int_{0}^{x}dye^{\overline{\Omega}y}[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]e^{-\overline{\Omega}y}-\int_{0}^{x}dy[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]+\int_{0}^{x}dy[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]$
(70)
$\displaystyle=\int_{0}^{x}dy\int_{0}^{y}dze^{\overline{\Omega}z}[\overline{\Omega},[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]]e^{-\overline{\Omega}z}+\int_{0}^{x}dy[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]\leavevmode\nobreak\
.$ (71)

By repeating the same procedure multiple times, we get that

$\displaystyle\int_{0}^{x}dye^{\overline{\Omega}y}[\overline{\Omega},\Omega_{1}^{{}^{\prime}}]e^{-\overline{\Omega}y}$
$\displaystyle=\underbrace{\int_{0}^{x}dv_{q}\cdots\int_{0}^{v_{2}}dv_{1}}_{q\leavevmode\nobreak\
\operatorname{times}}e^{\overline{\Omega}v_{1}}\underbrace{[\overline{\Omega},[\dots,[\overline{\Omega},\Omega_{1}^{\prime}]]\dots]}_{q\leavevmode\nobreak\
\operatorname{nested\leavevmode\nobreak\
commutators}}e^{-\overline{\Omega}v_{1}}$
$\displaystyle\qquad\qquad+\sum_{j=1}^{q-1}\int_{0}^{x}dv_{j}\cdots\int_{0}^{v_{2}}dv_{1}\underbrace{[\overline{\Omega},[\dots,[\overline{\Omega},\Omega_{1}^{\prime}]]\dots]}_{j\leavevmode\nobreak\
\operatorname{nested\leavevmode\nobreak\ commutators}}$ (72)
$\displaystyle=\int_{0}^{x}dv_{q}\cdots\int_{0}^{v_{2}}dv_{1}e^{\overline{\Omega}v_{1}}\operatorname{adj}^{q}_{\overline{\Omega}}\Omega_{1}^{\prime}e^{-\overline{\Omega}v_{1}}+\sum_{j=1}^{q-1}\int_{0}^{x}dv_{j}\cdots\int_{0}^{v_{2}}dv_{1}\operatorname{adj}^{j}_{\overline{\Omega}}\Omega_{1}^{\prime}\leavevmode\nobreak\
$ (73)
$\displaystyle=\int_{0}^{x}dv_{q}\cdots\int_{0}^{v_{2}}dv_{1}e^{\overline{\Omega}v_{1}}\operatorname{adj}^{q}_{\overline{\Omega}}\Omega_{1}^{\prime}e^{-\overline{\Omega}v_{1}}+\mathcal{P}_{q}(\alpha)\leavevmode\nobreak\
.$ (74)

where we use $\mathcal{P}_{q}(\alpha)$ to denote polynomials of degree at most
$q$ in $\alpha$. Using similar techniques, we have

$\displaystyle\int_{0}^{1}dxe^{\overline{\Omega}x}{\Omega}_{j}^{\prime}e^{-\overline{\Omega}x}=\int_{0}^{1}dx\int_{0}^{x}dv_{q-j+1}\dots\int_{0}^{v_{2}}dv_{1}e^{\overline{\Omega}v_{1}}\operatorname{adj}_{\overline{\Omega}}^{q-j+1}\Omega_{j}^{\prime}e^{-\overline{\Omega}v_{1}}+\mathcal{P}_{q}(\alpha).$
(75)

Substituting Eqs. 74 and 75 into Eq. 68, we get

$\displaystyle\partial_{t}(e^{-\overline{\Omega}}e^{\Omega})=\underbrace{\sum_{j=1}^{q}\int_{0}^{1}dx\int_{0}^{x}dv_{q-j+1}\dots\int_{0}^{v_{2}}dv_{1}e^{\overline{\Omega}v_{1}}\operatorname{adj}_{\overline{\Omega}}^{q-j+1}\Omega_{j}^{\prime}e^{-\overline{\Omega}v_{1}}}_{=\mathcal{O}\\!\left(\alpha^{q+1}\right)}+\mathcal{P}_{q}(\alpha).$
(76)

Recall that
$\partial_{t}(e^{-\overline{\Omega}}e^{\Omega})=\mathcal{O}\\!\left(\alpha^{q+1}\right)$.
Therefore, $\mathcal{P}_{q}(\alpha)$ must vanish in Eq. 76, leaving

$\displaystyle\partial_{t}(e^{-\overline{\Omega}}e^{\Omega})=\sum_{j=1}^{q}\int_{0}^{1}dx\int_{0}^{x}dv_{q-j+1}\dots\int_{0}^{v_{2}}dv_{1}e^{\overline{\Omega}v_{1}}\operatorname{adj}_{\overline{\Omega}}^{q-j+1}\Omega_{j}^{\prime}e^{-\overline{\Omega}v_{1}}.$
(77)

We can now use Eq. 77 to bound the Magnus truncation error:

$\displaystyle\norm{e^{\overline{\Omega}}-e^{\Omega}}$
$\displaystyle\leq\norm{e^{-\overline{\Omega}}e^{\Omega}-I}\leq\int_{0}^{t}dt\norm{\partial_{t}(e^{-\overline{\Omega}}e^{\Omega})}$
(78)
$\displaystyle\leq\sum_{j=1}^{q}\int_{0}^{t}dt\int_{0}^{1}dx\int_{0}^{x}dv_{q-j+1}\dots\int_{0}^{v_{2}}dv_{1}\norm{\operatorname{adj}_{\overline{\Omega}}^{q-j+1}\Omega_{j}^{\prime}}$
(79)
$\displaystyle=\sum_{j=1}^{q}\frac{t}{(q-j+2)!}\norm{\operatorname{adj}_{\overline{\Omega}}^{q-j+1}\Omega_{j}^{\prime}}.$
(80)

Using Corollaries 2 and 6, we have that

$\displaystyle\operatorname{adj}_{\overline{\Omega}}^{q-j+1}\Omega_{j}^{\prime}\in\mathcal{H}(c_{j},\kappa_{j},k_{j},\alpha^{q+1}d^{q+1}t^{q})$
(81)

for some constants $c_{j},\kappa_{j},k_{j}$ that depend on $q,k$ but are
independent of $n,\alpha,d,t$. Applying Lemma 5 to bound their norms, we
arrive at

$\displaystyle\norm{e^{\overline{\Omega}}-e^{\Omega}}=\mathcal{O}\\!\left(n(\alpha
td)^{q+1}\right).$ (82)

##  Appendix D Concentrated $k$-local Hamiltonians

Under the evolution generated by $A$, the Hamiltonian $B$ is no longer
$k$-local. However, because $A$ is geometrically local, the interaction-
picture Hamiltonian and the Magnus operator inherits many properties of a
$k$-local Hamiltonian. In this section, we develop tools to capture these
properties.

######  Definition 1 (Concentrated Operators).

An operator $A$ is considered concentrated on site $i$ if it admits a
decomposition $A=\sum_{r=0}^{\infty}A_{i,r}$ such that

  * •

$A_{i,r}$ is supported entirely within a radius $r$ from $i$, and

  * •

there exists a constant $c$ and $\kappa$ such that $\norm{A_{i,r}}\leq
ce^{-r/\kappa}$ for all $i,r$.

We use $\mathcal{C}_{i}(c,\kappa)$ to denote the set of such operators $A$ and
call $i$ the center of $A$.

If $A\in\mathcal{C}_{i}(c,\kappa)$, it follows that $\norm{A}\leq
c(1+\kappa)$. Note that if $A\in\mathcal{C}_{i}(c,\kappa)$, then
$A\in\mathcal{C}_{i}(c^{\prime},\kappa^{\prime})$ for all
$c^{\prime}>c,\kappa^{\prime}>\kappa$.

The following lemma provides a structure for the commutator between two
concentrated operators.

######  Lemma 3.

If $A\in\mathcal{C}_{i}(c,\kappa)$ and $B\in\mathcal{C}_{j}(c,\kappa)$, there
exists a set of concentrated operators $\\{\hat{A}_{\ell}:\hat{A}_{\ell}\in
C_{i}(c,2\kappa)\\}$ and $\\{\hat{B}_{\ell}:\hat{B}_{\ell}\in
C_{j}(c,2\kappa)\\}$ such that

$\displaystyle\commutator{A}{B}=\sum_{\ell}a_{\ell}\hat{A}_{\ell}\hat{B}_{\ell},$
(83)

and $\sum_{\ell}\absolutevalue{a_{\ell}}\leq 4e^{-d_{ij}/4\kappa}$, where
$d_{ij}$ is the distance between $i$ and $j$.

###### Proof.

Let $A=\sum_{r}A_{i,r}$ and $B=\sum_{r}B_{i,r}$ be the decompositions of $A$
and $B$ according to Definition 1. We divide $A$ into a sum of
$A_{<}=\sum_{r<d_{ij}/2}A_{i,r}$ and $A_{\geq}=\sum_{r\geq d_{ij}/2}A_{i,r}$
based on the radius $r$ and divide $B=B_{<}+B_{\geq}$ similarly. The
commutator can then be expanded as

$\displaystyle\commutator{A}{B}=\commutator{A_{<}}{B_{\geq}}+\commutator{A_{\geq}}{B}.$
(84)

Note that $\commutator{A_{<}}{B_{<}}$ vanishes because the operators are
supported on distinct regions of the lattice.

Note that $A_{<},A_{\geq}$ are both in $\mathcal{C}_{i}(c,\kappa)\subseteq
C_{i}(c,2\kappa)$. We can, however, prove a stronger statement for $A_{\geq}$.
Since $r\geq d_{ij}/2$, the identity $e^{-r/\kappa}\leq
e^{-d_{ij}/4\kappa}e^{-r/2\kappa}$ implies $\tilde{A}_{\geq}\equiv
e^{d_{ij}/4\kappa}A_{\geq}\in C_{i}(c,2\kappa)$. Repeating the same argument
for $B_{\geq}$, we have

$\displaystyle\commutator{A}{B}=e^{-d_{ij}/4\kappa}\left(\commutator{A_{<}}{\tilde{B}_{\geq}}+\commutator{\tilde{A}_{\geq}}{B}\right).$
(85)

Expanding the commutators, we arrive at Lemma 3 with $a_{\ell}=\pm
e^{-d_{ij}/4\kappa}$, and $\hat{A}_{1}\equiv A_{<}$,
$\hat{A}_{2}\equiv\tilde{A}_{\geq}$ , $\hat{B}_{1}\equiv B_{\geq}$, and
$\hat{B}_{2}\equiv B$. ∎

######  Definition 2 (Concentrated $k$-local Strings).

An operator $A$ is a concentrated $k$-local string if there exist constants
$c,\kappa$ such that

$\displaystyle A=A^{(1)}_{i_{1}}A^{(2)}_{i_{2}}\dots A^{(k)}_{i_{k}},$ (86)

where $A^{(j)}_{i_{j}}\in\mathcal{C}_{i_{j}}(c,\kappa)$ for all $j=1,\dots,k$.
We use $S(A)=\\{i_{1},\dots,i_{k}\\}$ to denote the set of centers.

######  Definition 3 (Concentrated $k$-local Hamiltonians).

A Hamiltonian $H=\sum_{\mu}a_{\mu}{H_{\mu}}$ is a concentrated $k$-local
Hamiltonian with an effective interaction degree $d$ if

  1. 1.

$H_{\mu}=A^{(\mu,1)}_{i_{\mu,1}}A^{(\mu,2)}_{i_{\mu,2}}\dots
A^{(\mu,k)}_{i_{\mu,k}}$ are concentrated $k$-local strings for all $\mu$,
where $A^{(\mu,j)}_{i_{\mu,j}}\in\mathcal{C}_{i_{\mu,j}}(c,\kappa)$, and

  2. 2.

for every site $l$ on the lattice,

$\displaystyle\sum_{\mu:l\in S(H_{\mu})}\absolutevalue{a_{\mu}}\leq d.$ (87)

We use $\mathcal{H}(c,\kappa,k,d)$ to denote the set of such Hamiltonians.
Note that the definition of $H$ is invariant under the following rescaling:
$c\to gc$ and $d\to d/g^{k}$ for some constant $g$.

Note that the definition of the effective interaction degree in Definition 3
is slightly different from the one defined in Section II.

We also note that a geometrically local Hamiltonian is also a concentrated
$1$-local Hamiltonian:

######  Lemma 4.

Let $B$ be a geometrically local with an interaction range $\chi$. It follows
that $B\in\mathcal{H}(c,\kappa,1,1)$ for $c=e^{\chi/\kappa}$ and $\kappa$ is
any constant.

Concentrated $k$-local Hamiltonians are generalization of $k$-local
Hamiltonians. The following lemma upper bounds the norm of a concentrated
$k$-local Hamiltonian.

######  Lemma 5.

If $H\in\mathcal{H}(c,\kappa,k,d)$ on $n$ qubits, then $\norm{H}\leq
c^{k}(\kappa+1)^{k}nd$.

###### Proof.

We have

$\displaystyle\norm{H}\leq\sum_{\mu}\absolutevalue{a_{\mu}}\norm{H_{\mu}}\leq
c^{k}(\kappa+1)^{k}\sum_{\mu}\absolutevalue{a_{\mu}}\leq
c^{k}(\kappa+1)^{k}\sum_{i}\sum_{\mu:i\in
S(H_{\mu})}\absolutevalue{a_{\mu}}\leq c^{k}(\kappa+1)^{k}nd,$ (88)

where in the second inequality we used that $\|A_{i_{\mu,l}}^{(\mu,l)}\|\leq
c(1+\kappa),\forall A_{i_{\mu,l}}^{(\mu,l)}\in C_{i_{\mu,j}}(c,\kappa)$, and
in the final inequality we overcounted the sum over $\mu$ by summing over
index $i$ and $\mu:i\in S(H_{\mu})$. ∎

######  Lemma 6.

If $H\in\mathcal{H}(c,\kappa,k_{1},d_{1})$ and
$K\in\mathcal{H}(c,\kappa,k_{2},d_{2})$, then
$\commutator{H}{K}\in\mathcal{H}(c^{\prime},\kappa^{\prime},k^{\prime},d^{\prime})$
for

$\displaystyle c^{\prime}$
$\displaystyle=c(8(4\kappa+1)k_{1}k_{2})^{1/(k_{1}+k_{2})}\leavevmode\nobreak\
,$ (89) $\displaystyle\kappa^{\prime}$
$\displaystyle=2\kappa\leavevmode\nobreak\ ,$ (90) $\displaystyle k^{\prime}$
$\displaystyle=k_{1}+k_{2}\leavevmode\nobreak\ ,$ (91) $\displaystyle
d^{\prime}$ $\displaystyle=d_{1}d_{2}\leavevmode\nobreak\ .$ (92)

###### Proof.

Let $H=\sum_{\mu}a_{\mu}H_{\mu}$ and $K=\sum_{\nu}b_{\nu}K_{\nu}$ be
decompositions into concentrated $k_{1}$-local and $k_{2}$-local strings. Let
$V_{\mu,\nu}=a_{\mu}b_{\nu}\commutator{H_{\mu}}{K_{\nu}}=\commutator{A^{(\mu,1)}_{i_{\mu,1}}A^{(\mu,2)}_{i_{\mu,2}}\dots
A^{(\mu,k_{1})}_{i_{\mu,k_{1}}}}{B^{(\nu,1)}_{j_{\nu,1}}B^{(\nu,2)}_{j_{\nu,2}}\dots
B^{(\nu,k_{2})}_{j_{\nu,k_{2}}}}$. Using the chain rule for commutators, we
have

$\displaystyle
V_{\mu,\nu}=a_{\mu}b_{\nu}\sum_{u,v}A^{(\mu,1)}_{i_{\mu,1}}\dots
A^{(\mu,u-1)}_{i_{\mu,u-1}}B^{(\nu,1)}_{j_{\nu,1}}\dots
B^{(\nu,v-1)}_{j_{\nu,v-1}}\commutator{A^{(\mu,u)}_{i_{\mu,u}}}{B^{(\nu,v)}_{j_{\nu,v}}}B^{(\nu,v+1)}_{j_{\nu,v+1}}\dots
B^{(\nu,k_{2})}_{j_{\nu,k_{2}}}A^{(\mu,u+1)}_{i_{\mu,u+1}}\dots
A^{(\mu,k_{1})}_{i_{\mu,k_{1}}}.$ (93)

Since $A^{(\mu,u)}_{i_{\mu,u}}$ and $B^{(\nu,v)}_{j_{\nu,v}}$ are concentrated
operators, Lemma 3 implies a decomposition

$\displaystyle\commutator{A^{(\mu,u)}_{i_{u}}}{B^{(\nu,v)}_{j_{v}}}=\sum_{\ell_{\mu,\nu,u,v}}\beta_{\ell_{\mu,\nu,u,v}}\hat{A}^{(\ell_{\mu,\nu,u,v})}_{i_{\mu}}\hat{B}^{(\ell_{\mu,\nu,u,v})}_{j_{\nu}}$
(94)

such that
$\hat{A}^{(\ell_{\mu,\nu,u,v})}_{i_{\mu}}\in\mathcal{C}_{i_{\mu}}(c,2\kappa)$,
$\hat{B}^{(\ell_{\mu,\nu,u,v})}_{j_{\nu}}\in\mathcal{C}_{j_{\nu}}(c,2\kappa)$,
and
$\sum_{\ell_{\mu,\nu,u,v}}\absolutevalue{\beta_{\ell_{\mu,\nu,u,v}}}=4e^{-d_{i_{\mu}j_{\nu}}/4\kappa}$.
Therefore, the commutator $\commutator{H}{K}$ can be decomposed as

$\displaystyle\commutator{H}{K}=\sum_{\mu\nu}a_{\mu}b_{\nu}\sum_{u,v}\sum_{\ell_{\mu,\nu,u,v}}\beta_{\ell_{\mu,\nu,u,v}}V_{\mu,\nu,u,v,\ell_{\mu,\nu,u,v}}=\sum_{\tau}c_{\tau}V_{\tau}$
(95)

for some concentrated $(k_{1}+k_{2})$-local strings
$V_{\tau}=V_{\mu,\nu,u,v,\ell_{\mu,\nu,u,v}}$. Here, we use $\tau$ to denote a
tuple of $\mu,\nu,u,v,\ell_{\mu,\nu,u,v}$.

Next, we prove the second condition of Definition 3 for some $d^{\prime}$ to
be determined. A fixed site $i$ is in $S(V_{\mu,\nu,u,v,\ell_{\mu,\nu,u,v}})$
only if it is in either $S(H_{\mu})$ or $S(K_{\nu})$. Suppose $i\in
S(H_{\mu})=\\{i_{1},\dots,i_{k_{1}}\\}$. For each $\mu,u,v$, we have

$\displaystyle\sum_{\nu}\sum_{\ell_{\mu,\nu,u,v}}\absolutevalue{a_{\mu}}\absolutevalue{b_{\nu}}\absolutevalue{\beta_{\ell_{\mu,\nu,u,v}}}=4\sum_{\nu}\absolutevalue{a_{\mu}}\absolutevalue{b_{\nu}}e^{-d_{i_{\mu}j_{\nu}}/4\kappa}\leq
4\absolutevalue{a_{\mu}}\sum_{j\in\Lambda}e^{-d_{i_{\mu}j}/4\kappa}\sum_{\nu:j_{\nu}=j}\absolutevalue{b_{\nu}}\leq
4(4\kappa+1)\absolutevalue{a_{\mu}}d_{2}.$ (96)

We then upper bound the sums over $u,v$ by $k_{1},k_{2}$ respectively, and
noting that the sum over $\mu$ is restricted to $\mu$ satisfying $i\in
S(H_{\mu})$, we have

$\displaystyle\sum_{\tau:i\in
S(V_{\tau})}\absolutevalue{c_{\tau}}\leq\left(\sum_{\tau:i\in
S(H_{\mu})}\absolutevalue{c_{\tau}}+\sum_{\tau:i\in
S(K_{\mu})}\absolutevalue{c_{\tau}}\right)\leq
8(4\kappa+1)k_{1}k_{2}d_{1}d_{2},$ (97)

where a factor of two comes from the case when $i\in S(K_{\nu})$, and a factor
of $d$ comes from the fact that $\sum_{\mu}|a_{\mu}|\leq d_{1}$.

Therefore, by combining Eq. 95 and Eq. 97 and using the invariance of
$\mathcal{H}$ under rescaling, we get that

$\displaystyle[H,K]\in\mathcal{H}(c,2\kappa,k_{1}+k_{2},8(4\kappa+1)k_{1}k_{2}d_{1}d_{2})=\mathcal{H}(c(8(4\kappa+1)k_{1}k_{2})^{1/(k_{1}+k_{2})},2\kappa,k_{1}+k_{2},d_{1}d_{2})\leavevmode\nobreak\
.$ (98)

∎

The following lemma is also useful for counting the support of a nested
commutator between concentrated 1-local Hamiltonians.

######  Lemma 7.

Let $H^{(1)},\dots,H^{(q)}\in\mathcal{H}(c,\kappa,1,1)$ be $q$ concentrated
$1$-local Hamiltonians. For any integers $R\geq 1$, there exists
$V=\sum_{i\in\Lambda}V_{i}$ such that

  * •

$V\in\mathcal{H}(c_{q},2^{q}\kappa,q,1)$, where $c_{q}$ is a constant
depending on only $q,\kappa$,

  * •

$V_{i}$ is supported on at most $q+R$ qubits, and

  * •

$\norm{V-\commutator{H^{(q)}}{\dots\commutator{H^{(2)}}{H^{(1)}}}}=\mathcal{O}\\!\left(ne^{-R/4\kappa}\right)$.

###### Proof.

Let $H^{(k)}=\sum_{i\in\Lambda}H^{(k)}_{i}$ be the decomposition of $H^{(k)}$
into concentrated operators, where $H^{(k)}_{i}\in\mathcal{C}_{i}(c,\kappa)$.
Define

$\displaystyle
W_{i}=\sum_{i}\commutator{H^{(q)}}{\dots\commutator{H^{(2)}}{H^{(1)}_{i}}}=\sum_{i}\sum_{i_{2},\dots,i_{q}\in\Lambda}\sum_{r_{1},\dots,r_{q}=0}^{\infty}\commutator{H^{(q)}_{i_{q},r_{q}}}{\dots\commutator{H^{(2)}_{i_{2},r_{2}}}{H^{(1)}_{i,r_{1}}}},$
(99)

where we have further decomposed concentrated operators using their
definitions. Let $g(i,i_{2},\dots,i_{q},r_{1},\dots,r_{q})$ be the number of
qubits in the joint support of
$H^{(q)}_{i_{q},r_{q}},\dots,H^{(2)}_{i_{2},r_{2}},H^{(1)}_{i,r_{1}}$. From
$W_{i}$, we construct $V_{i}$ by discarding summands where
$g(i,i_{2},\dots,i_{q},r_{1},\dots,r_{q})>q+R$. Lemma 6 implies
$V=\sum_{i}V_{i}\in\mathcal{H}(c_{q},2^{q}\kappa,q,1)$ for a constant $c_{q}$
that may depend on $q$. By construction, $V_{i}$ is supported on at most $q+R$
qubits.

To bound the error $\norm{V_{i}-W_{i}}$, we recall $H^{(1)}_{i,r_{1}}$ is
supported on at most $2r_{1}$ qubits around $i$. Similarly,
$H^{(2)}_{i_{2},r_{2}}$ is supported on at most $2r_{2}$ qubits around
$i_{2}$. Let $m+q=g(i,i_{2},\dots,i_{q},r_{1},\dots,r_{q})$. Besides the $q$
centers, the other $m$ qubits come from the balls of radii
$r_{1},\dots,r_{q}$. It follows that

$\displaystyle m\leq 2(r_{1}+r_{2}+\dots+r_{q}).$ (100)

Therefore,

$\displaystyle\norm{\commutator{H^{(q)}_{i_{q},r_{q}}}{\dots\commutator{H^{(2)}_{i_{2},r_{2}}}{H^{(1)}_{i,r_{1}}}}}\leq
2^{q}c^{q}e^{-r_{1}/\kappa}\dots e^{-r_{q}/\kappa}\leq
2^{q}c^{q}e^{-m/2\kappa}.$ (101)

So each summand that we dropped will be bounded by $2^{q}c^{q}e^{-m/2\kappa}$.

Next, we count the number of such summand, i.e. the number of the choices of
$r_{1},\dots,r_{q}$ and $i_{2},\dots,i_{q}$ such that
$g(i,i_{2},\dots,i_{q},r_{1},\dots,r_{q})=m+q$. Denote this number by $f(m)$.
Note that $m\geq 2r_{1}$, so $0\leq r_{1}\leq m/2$ and the number of choices
for $r_{1}$ is at most $m/2+1\leq m$, assuming $m\geq 2$. Using a similar
argument for $r_{2},\dots,r_{q}$, the number of choices of these radii is at
most $m^{q}$.

To count the number of centers $i_{2},\dots,i_{q}$, we note that
$d_{i,i_{2}}\leq r_{1}+r_{2}\leq m$ (otherwise, the commutator would vanish).
Similarly, $d_{i,i_{3}}\leq r_{1}+2r_{2}+r_{3}\leq 2m$. Therefore, the number
of choices for $i_{2},\dots,i_{q}$ is at most $(2m)(4m)(6m)\dots(2(q-1)m)\leq
2^{q}q!m^{q}$. Combining with the counting argument for $r_{1},\dots,r_{q}$,
we have

$\displaystyle f(m)\leq 2^{q}q!m^{2q}.$ (102)

Combining Eqs. 102 and 101, the error of approximating $W_{i}$ with $V_{i}$ is

$\displaystyle\norm{V_{i}-W_{i}}\leq\sum_{m=R+1}^{\infty}4^{q}q!c^{q}m^{2q}e^{-m/2\kappa}\leq\sum_{m=R+1}^{\infty}4^{q}q!(2q!)(4\kappa)^{2q}c^{q}e^{-m/4\kappa}=\mathcal{O}\\!\left(e^{-R/4\kappa}\right),$
(103)

where we have used $m^{2q}\leq(2q)!(4\kappa)^{2q}e^{m/4\kappa}$ for $m\geq 2$.
Since
$\sum_{i\in\Lambda}W_{i}=\commutator{H^{(q)}}{\dots\commutator{H^{(2)}}{H^{(1)}}}$,
we arrive at the lemma. ∎

##  Appendix E The Lieb-Robinson bound and the locality of Magnus operators

In this section, we use the Lieb-Robinson bound and the machinery developed in
Appendix D to bound the error of the Magnus operator. We begin by showing that
a $k$-local operator evolved under a geometrically local Hamiltonian can be
written as a concentrated operator.

######  Lemma 8.

Let $B$ be a unit-norm operator supported on a single site $i$ and $A$ is a
geometrically local Hamiltonian. There exists constant $c,\kappa$ such that

$\displaystyle e^{iAt}Be^{-iAt}\in\mathcal{C}_{i}(c,\kappa).$ (104)

###### Proof.

Given an operator $O$, we define a superoperator

$\displaystyle\mathbb{P}_{i,r}O\equiv\int_{U\sim\mathcal{U}_{r}}dUU^{\dagger}OU,$
(105)

where the integral is over the Haar measure of unitaries supported entirely
outside a ball of radius $r$ from $i$. Consequently, $\mathbb{P}_{i,r}O$ is
supported entirely within this ball.

Let $B(t)=e^{iAt}Be^{-iAt}$. The Lieb-Robinson bound implies that
$\mathbb{P}_{i,r}B(t)$ is a good approximation to $B(t)$ in the limit of large
$r$ Eq. 54. Indeed, we have

$\displaystyle\norm{B(t)-\mathbb{P}_{i,r}B(t)}\leq\int_{U}dU\leavevmode\nobreak\
\norm{\commutator{U}{B(t)}}\leq\int_{U}dU\leavevmode\nobreak\
Ke^{vt-r/\chi}=Ke^{vt-r/\chi},$ (106)

for some constant $K,v,\chi$.

We define

$\displaystyle B_{i,r}=\mathbb{P}_{i,r}B(t)-\mathbb{P}_{i,r-1}B(t)$ (107)

for $r=0,1,\dots$. Here, we define $\mathbb{P}_{i,-1}B(t)=0$ for convenience.
Then, we have

$\displaystyle\norm{B_{i,r}}\leq\norm{B(t)-\mathbb{P}_{i,r}B(t)}+\norm{B(t)-\mathbb{P}_{i,r-1}B(t)}\leq
2Ke^{vt-1/\chi}e^{-r/\chi}.$ (108)

Therefore, $B(t)\in\mathcal{C}_{i}(c,\kappa)$ with $c=2Ke^{vt-1/\chi}$ and
$\kappa=\chi$. ∎

######  Lemma 9.

Let $A$ be a geometrically local Hamiltonian with an interaction range
$\chi\geq 1$ and $B\in\mathcal{H}(c,\kappa,1,1)$ be a concentrated $1$-local
Hamiltonian with $\kappa\leq\chi$. We have

$\displaystyle e^{iAt}Be^{-iAt}\in\mathcal{H}(c^{\prime},2\chi,1,1),$ (109)

for $c^{\prime}=16Kc\chi^{2}e^{vt-1/\chi}$ for a constant $K$. In particular,
$c^{\prime}$ is a constant if $t=\mathcal{O}(1)$.

###### Proof.

Let $B=\sum_{i\in\Lambda}B_{i}=\sum_{i\in\Lambda}\sum_{r=0}^{\infty}B_{i,r}$
be a decomposition of $B$ into concentrated operators, where
$B_{i}\in\mathcal{C}_{i}(c,\kappa)$. For each $i$, let
$B_{i,r}(t)=e^{iAt}B_{i,r}e^{-iAt}$. We have

$\displaystyle
e^{iAt}B_{i}e^{-iAt}=\sum_{r=0}^{\infty}\left[\mathbb{P}_{i,r}B_{i,r}(t)+\sum_{s=1}^{\infty}(\mathbb{P}_{i,r+s}-\mathbb{P}_{i,r+s-1})B_{i,r}(t)\right],$
(110)

where $\mathbb{P}_{i,s}$ is defined in Eq. 105. Let
$V_{i,r,s}=(\mathbb{P}_{i,r+s}-\mathbb{P}_{i,r+s-1})B_{i,r}(t)$ for $s\geq 1$
and $V_{i,r,0}=\mathbb{P}_{i,r}B_{i,r}(t)$. It follows from the Lieb-Robinson
bound that

$\displaystyle\norm{V_{i,r,s}}\leq
2Ke^{vt-1/\chi}(2r+1)e^{-s/\chi}\norm{B_{i,r}}\leq
2Kce^{vt-1/\chi}(2r+1)e^{-(r+s)/\chi},$ (111)

where the factor $2r+1$ comes from the size of the support of $B_{i,r}(0)$.
Under this notation,

$\displaystyle B_{i}^{{}^{\prime}}\equiv
e^{iAt}B_{i}e^{-iAt}=\sum_{r=0}^{\infty}\sum_{s=0}^{\infty}V_{i,r,s},$ (112)

where $V_{i,r,s}$ is supported entirely within a ball of radius $r+s$ around
$i$. Next, we group $V_{i,r,s}$ based on $t=r+s$. Define

$\displaystyle B^{\prime}_{i,t}=\sum_{r=0}^{t}V_{i,r,t-r}.$ (113)

It follows that $B^{\prime}_{i,t}$ is supported within a ball of radius $t$
from $i$ and

$\displaystyle\norm{B^{\prime}_{i,t}}\leq
2Kce^{vt-1/\chi}\sum_{r=0}^{t}(2r+1)e^{-t/\chi}=16Kc\chi^{2}e^{vt-1/\chi}e^{-t/2\chi}$
(114)

where we have used $\sum_{r=0}^{t}(2r+1)\leq 8\chi^{2}e^{-t/2\chi}$ for all
$\chi\geq 1$. Therefore $B^{\prime}_{i}\in\mathcal{C}_{i}(c^{\prime},2\chi)$
with $c^{\prime}=16Kc\chi^{2}e^{vt-1/\chi}$ and the lemma follows. ∎

Several corollaries follow immediately:

######  Corollary 1.

If $B$ is a $k$-local Hamiltonian with an effective interaction degree $d$ and
$A$ is a geometrically local Hamiltonian with an interaction range $\chi$,
then

$\displaystyle e^{iAt}Be^{-iAt}\in\mathcal{H}(c,\kappa,k,d),$ (115)

for some $c$ that may depend on $t$ and a constant $\kappa$. Moreover, if
$t=O(1)$, then $c$ is a constant.

Corollaries 1 and 6 imply

######  Corollary 2.

Let $B$ be a $k$-local Hamiltonian with an effective interaction degree $d$
and $A$ be a geometrically local Hamiltonian. Let $\Omega_{q}(t)$ be the
$q$th-order Magnus operator generating the interaction-picture evolutions of
$H=A+\alpha B$. For $t=O(1)$, we have

$\displaystyle\Omega_{q}(t)\in\mathcal{H}(c_{q},\kappa_{q},qk,\alpha^{q}d^{q}t^{q}),\text{
and}$ (116)
$\displaystyle\Omega_{q}^{\prime}(t)\in\mathcal{H}(c_{q},\kappa_{q},qk,\alpha^{q}d^{q}t^{q-1})$
(117)

for some constants $c_{q}$ and $\kappa_{q}$. In addition, Eq. 116 implies

$\displaystyle\sum_{j=1}^{q}\Omega_{j}(t)\in\mathcal{H}(c_{q}^{\prime},\kappa_{q}^{\prime},qk,\alpha
td)$ (118)

for $\alpha dt\leq 1$.

Lemmas 7, 1 and 9 imply a more robust structure of the Magnus operators:

######  Corollary 3.

Let $B$ and $A$ be geometrically local Hamiltonians and $\chi$ be the
interaction range of $A$. Let $\Omega_{q}(t)$ be the $q$th-order Magnus
operator generating the interaction-picture evolutions of $H=A+\alpha B$. For
$t=O(1)$ and for any $R\geq 1$, there exists an operator
$\Omega_{\text{loc},q}=\sum_{i\in\Lambda}\Omega_{\text{loc},q,i}$ such that

  * •

$\Omega_{\text{loc},q}\in\mathcal{H}(c_{q},2^{q+1}\chi,q,1)$ for a constant
$c_{q}$ depending on only $q,\chi$,

  * •

$\Omega_{\text{loc},q,i}$ is supported on at most $q+R$, and

  * •

$\norm{\Omega_{q}-\Omega_{\text{loc},q}}=\mathcal{O}\\!\left(nt^{q}\alpha^{q}e^{-R/8\chi}\right)$.

## References

  * [1] Masuo Suzuki.  General theory of fractal path integrals with applications to many-body theories and statistical physics.  Journal of Mathematical Physics, 32(2):400–407, February 1991.
  * [2] Andrew M. Childs and Yuan Su.  Nearly optimal lattice simulation by product formulas.  Physical Review Letters, 123:050503, Aug 2019.  arXiv:1901.00564.
  * [3] Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe, and Shuchen Zhu.  Theory of Trotter Error with Commutator Scaling.  Physical Review X, 11(1):011020, February 2021.  arXiv:1912.08854.
  * [4] Andrew M Childs and Nathan Wiebe.  Hamiltonian simulation using linear combinations of unitary operations.  Quantum Information and Computation, 12:901–924, Nov 2012.  arXiv:1202.5822.
  * [5] Dominic Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma.  Simulating Hamiltonian dynamics with a truncated taylor series.  Physical Review Letters, 114(9):090502, March 2015.  arXiv:1412.4687.
  * [6] Guang Hao Low and Isaac L Chuang.  Optimal Hamiltonian simulation by quantum signal processing.  Physical Review Letters, 118(1):010501, Jan 2017.  arXiv:1606.02685.
  * [7] Guang Hao Low and Isaac L Chuang.  Hamiltonian simulation by qubitization.  Quantum, 3:163, July 2019.  arXiv:1610.06546.
  * [8] Andrew M. Childs, Aaron Ostrander, and Yuan Su.  Faster quantum simulation by randomization.  Quantum, 3:182, September 2019.  arXiv:1805.08385.
  * [9] Earl Campbell.  Random compiler for fast hamiltonian simulation.  Physical Review Letters, 123:070503, Aug 2019.  arXiv:1811.08017.
  * [10] Yingkai Ouyang, David R. White, and Earl T. Campbell.  Compilation by stochastic Hamiltonian sparsification.  Quantum, 4:235, February 2020.  arXiv:1910.06255.
  * [11] Jeongwan Haah, Matthew B Hastings, Robin Kothari, and Guang Hao Low.  Quantum algorithm for simulating real time evolution of lattice hamiltonians.  SIAM Journal on Computing, 52(6):FOCS18–250, 2021.  arXiv:1801.03922.
  * [12] Guang Hao Low and Nathan Wiebe.  Hamiltonian simulation in the interaction picture.  arXiv preprint arXiv:1805.00675, 2018.  arXiv:1805.00675.
  * [13] Dong An, Di Fang, and Lin Lin.  Time-dependent Hamiltonian Simulation of Highly Oscillatory Dynamics and Superconvergence for Schrödinger Equation.  Quantum, 6:690, April 2022.  arXiv:2111.03103.
  * [14] Abhishek Rajput, Alessandro Roggero, and Nathan Wiebe.  Hybridized Methods for Quantum Simulation in the Interaction Picture.  Quantum, 6:780, August 2022.  arXiv:2109.03308.
  * [15] Yi-Hsiang Chen, Amir Kalev, and Itay Hen.  Quantum algorithm for time-dependent Hamiltonian simulation by permutation expansion.  PRX Quantum, 2:030342, Sep 2021.  arXiv:2103.15334.
  * [16] Dominic W. Berry, Andrew M. Childs, Yuan Su, Xin Wang, and Nathan Wiebe.  Time-dependent Hamiltonian simulation with $L^{1}$-norm scaling.  Quantum, 4:254, April 2020.  arXiv:1906.07115.
  * [17] Sergio Blanes, Fernando Casas, Jose-Angel Oteo, and José Ros.  The Magnus expansion and some of its applications.  Physics Reports, 470(5-6):151–238, 2009.  arXiv:0810.5488.
  * [18] Ana Arnal, Fernando Casas, and Cristina Chiralt.  A general formula for the Magnus expansion in terms of iterated integrals of right-nested commutators.  Journal of Physics Communications, 2(3):035024, 2018.  arXiv:1710.10851.
  * [19] Guannan Chen, Mohammadali Foroozandeh, Chris Budd, and Pranav Singh.  Quantum simulation of highly-oscillatory many-body hamiltonians for near-term devices.  arXiv preprint arXiv:2312.08310, 2023.  arXiv:2312.08310.
  * [20] Elliott H Lieb and Derek W Robinson.  The finite group velocity of quantum spin systems.  Communications in Mathematical Physics, 28(3):251–257, 1972.
  * [21] S. Bravyi, M. B. Hastings, and F. Verstraete.  Lieb-robinson bounds and the generation of correlations and topological quantum order.  Physical Review Letters, 97:050401, Jul 2006.  arXiv:quant-ph/0603121.
  * [22] Jan Lukas Bosse, Andrew M Childs, Charles Derby, Filippo Maria Gambetta, Ashley Montanaro, and Raul A Santos.  Efficient and practical Hamiltonian simulation from time-dependent product formulas.  arXiv preprint arXiv:2403.08729, 2024.  arXiv:2403.08729.
  * [23] Pablo Antonio Moreno Casares, Modjtaba Shokrian Zini, and Juan Miguel Arrazola.  Quantum simulation of time-dependent hamiltonians via commutator-free quasi-Magnus operators.  arXiv preprint arXiv:2403.13889, 2024.  arXiv:2403.13889.
  * [24] Lin Lin.  Lecture notes on quantum algorithms for scientific computation.  arXiv preprint arXiv:2201.08309, 2022.  arXiv:2201.08309.
  * [25] Per Christian Moan and Jitse Niesen.  Convergence of the Magnus series.  Foundations of Computational Mathematics, 8:291–301, 2008.  arXiv:math/0609198.
  * [26] Wilhelm Magnus.  On the exponential solution of differential equations for a linear operator.  Communications on Pure and Applied Mathematics, 7(4):649–673, 1954\.
  * [27] Andrew M. Childs, Jiaqi Leng, Tongyang Li, Jin-Peng Liu, and Chenyi Zhang.  Quantum simulation of real-space dynamics.  Quantum, 6:860, November 2022.  arXiv:2203.17006.