# Better bounds for low-energy product formulas

Kasra Hejazi  Xanadu, Toronto, ON, M5G2C8, Canada  Modjtaba Shokrian Zini
Xanadu, Toronto, ON, M5G2C8, Canada  Juan Miguel Arrazola  Xanadu, Toronto,
ON, M5G2C8, Canada

###### Abstract

Product formulas are one of the main approaches for quantum simulation of the
Hamiltonian dynamics of a quantum system. Their implementation cost is
computed based on error bounds which are often pessimistic, resulting in
overestimating the total runtime. In this work, we rigorously consider the
error induced by product formulas when the state undergoing time evolution
lies in the low-energy sector with respect to the Hamiltonian of the system.
We show that in such a setting, the usual error bounds based on the operator
norm of nested commutators can be replaced by those restricted to suitably
chosen low-energy subspaces, yielding tighter error bounds. Furthermore, under
some locality and positivity assumptions, we show that the simulation of
generic product formulas acting on low-energy states can be done
asymptotically more efficiently when compared with previous results.

Simulating the Hamiltonian dynamics of a quantum system is a prominent
subroutine in many quantum algorithms, ranging from energy estimation Kitaev
(1995); Cleve  _et al._ (1998); Lin and Tong (2022) and quantum dynamics
Flannigan  _et al._ (2022), to matrix inversion Harrow  _et al._ (2009) and
quantum adiabatic optimization Albash and Lidar (2018). Several methods for
quantum simulation have been developed that include product formulas Lloyd
(1996); Suzuki (1991), truncated Taylor series Berry  _et al._ (2015), and
qubitization with quantum signal processing Low and Chuang (2017, 2019).

The performance of these algorithms has been studied extensively in general
settings, as usually time evolution is the most costly component of the
algorithms. For such cost estimation, characteristics of the Hamiltonian of
the problem are often taken into account, while the specific properties of the
initial state that undergoes time evolution are usually ignored, as
illustrated by the very few existing works that actually consider the
specifics of the initial state Şahinoğlu and Somma (2021); Su  _et al._
(2021). However, it is very likely for these initial states to be highly
structured, for example by having support only on low-energy states or having
a finite correlation length.

This suggests that in practical settings, a less costly quantum simulation
could be feasible if the above-mentioned properties are utilized. As a result,
it is an important task to find novel algorithms or estimate the costs of
existing ones when more information about the initial state is available.

In this work, we focus on product formulas as the method of quantum simulation
for initial states with negligible support over high-energy eigenstates. Such
states constitute an important set of applications, e.g., in ground-state
energy estimation Kitaev (1995); Cleve  _et al._ (1998); Lin and Tong (2022),
computation of low-energy properties of many-body systems Jordan  _et al._
(2012); Tong  _et al._ (2021), and quantum adiabatic algorithms Albash and
Lidar (2018). The characterization of the cost for product formulas is often
based on the simulation error, and extensive error analysis in generic cases
of product formulas has been performed Childs  _et al._ (2021).

The problem of product formulas as applied to low-energy states has also been
considered in Ref. Şahinoğlu and Somma (2021), where it is shown that under
certain assumptions, the first order Trotter-Suzuki formulas can induce a
lower cost when applied to low-energy states. However, Ref. Şahinoğlu and
Somma (2021) does not observe cost reduction for orders higher than one; in
fact, the cost scaling with system size becomes worse when compared with the
generic bounds in Childs  _et al._ (2021).

We develop a general framework for studying the error of a product formula for
low-energy states. We show that we can modify the error bounds of Ref. Childs
_et al._ (2021), which are based on nested commutators of the Hamiltonian
terms, so that all the norms are calculated in suitably chosen low-energy
subspaces, rather than the whole Hilbert space. Based on this framework and
imposing some further assumptions similar to those in Şahinoğlu and Somma
(2021), we show that the asymptotic scaling of the cost with system size can
in fact be improved for all orders, contrary to the findings in Şahinoğlu and
Somma (2021). We obtain a better cost scaling due to our tighter error
analysis. A summary of asymptotic cost formulas comparing our work, the
general cost scaling of Ref. Childs  _et al._ (2021), and previous work on
low-energy product formulas Şahinoğlu and Somma (2021) is shown in Table 1.

The rest of the paper is organized as follows. First, we discuss the setting
in which our results hold. Then we present our findings on product formula
errors for low-energy states, and briefly describe the ideas behind the
proofs, with details in the Supplemental Material SM . A crucial step is the
calculation of the leakage from low to higher energies due to the action of
the product formulas, which can be of independent interest. We then use the
derived low-energy error bounds to estimate the cost of time evolution for
low-energy states. We conclude with a discussion of the results and future
directions.

Order | General, Ref. Childs  _et al._ (2021) | Low-energy, Ref. Şahinoğlu and Somma (2021) | Low-energy, present work
---|---|---|---
$1$ | $\mathcal{O}\left(\frac{T^{2}N}{\epsilon}\right)$ | $\tilde{\mathcal{O}}\left(\frac{T^{2}\Delta^{2}}{\epsilon}+\frac{T^{4/3}N^{2/3}}{\epsilon^{1/3}}\right)$ | $\tilde{\mathcal{O}}\left(\frac{T^{2}\Delta^{2}}{\epsilon}+\frac{T^{6/5}N^{2/5}}{\epsilon^{1/5}}\right)$
$2$ | $\mathcal{O}\left(\frac{T^{3/2}N^{1/2}}{\epsilon^{1/2}}\right)$ | $\tilde{\mathcal{O}}\left(\frac{T^{3/2}\Delta^{3/2}}{\epsilon^{1/2}}+\frac{T^{6/5}N^{3/5}}{\epsilon^{1/5}}\right)$ | $\tilde{\mathcal{O}}\left(\frac{T^{3/2}\Delta^{3/2}}{\epsilon^{1/2}}+\frac{T^{12/11}N^{3/11}}{\epsilon^{1/11}}\right)$
$3$ | $\mathcal{O}\left(\frac{T^{4/3}N^{1/3}}{\epsilon^{1/3}}\right)$ | $\tilde{\mathcal{O}}\left(\frac{T^{4/3}\Delta^{4/3}}{\epsilon^{1/3}}+\frac{T^{8/7}N^{4/7}}{\epsilon^{1/7}}\right)$ | $\tilde{\mathcal{O}}\left(\frac{T^{4/3}\Delta^{4/3}}{\epsilon^{1/3}}+\frac{T^{20/19}N^{4/19}}{\epsilon^{1/19}}\right)$
Table 1: Comparing the costs (number of Trotter steps) between the general
case bounds of Ref. Childs  _et al._ (2021), the low-energy results from Ref.
Şahinoğlu and Somma (2021), and the present work. Hamiltonians are taken to
consist of positive semidefinite terms for this comparison. The scalings in
the low-energy cases are correct up to polylogarithmic factors in the
parameters appearing in the arguments of $\tilde{\mathcal{O}}$. The parameters
$T,N,\epsilon$, and $\Delta$ denote total evolution time, the system size, the
simulation error, and the energy cutoff of the initial state, respectively.

##  I Setting

We consider time-independent qubit Hamiltonians written in the form
$H=\sum_{m=1}^{M}H_{m}$. Each term $H_{m}$ comprises a number $L_{m}$ of
$k$-local interactions that commute with each other; while $L_{m}$ can scale
with the total number of qubits $N$, we take $M$ to be $O(1)$ and let
$L=\max_{m}L_{m}$. Each qubit is involved in at most $d$ of the $k$-local
interactions. Furthermore, we introduce the parameter $J$ such that the norm
of each local interaction is upper bounded by $J/2$. 111More precisely,
$J:=2\max_{h}\min_{c}\|h+c\|$, with $h$ running over single local terms of the
Hamiltonian and $c$ a constant, see SM .

We would like to bound the error of a product formula written in the form

$\mathcal{W}(s)=e^{-is_{q}H_{m_{q}}}\ldots
e^{-is_{2}H_{m_{2}}}\,e^{-is_{1}H_{m_{1}}},$ (1)

that approximates time evolution $e^{-iHs}$ under $H$. Here,
$m_{1},\ldots,m_{q}$ is a permutation of $1,\ldots,M$, with repetitions in the
case of higher-order product formulas, with ‘order’ defined further below.
Each $s_{j}=a_{j}s$ is proportional to $s$ with $a_{j}\in\mathbb{R}$ and the
total number $q$ of fast-forwardable exponentials in Eq. (1) is taken to be
$O(1)$.

The operator $\mathcal{W}(s)$ in Eq. (1) implements the time evolution for a
single Trotter step of length $s$; evolving for total time $T$ can be
performed through a number $r=T/s$ of Trotter steps. If the error in a single
step is $\varepsilon(s)$, the total accumulated error will be
$r\varepsilon(s)$. We characterize the cost in terms of the Trotter number
$r$; therefore, a tighter error estimation is a proxy to a lower cost.

The order of a product formula is defined based on the error that is
accumulated in a single Trotter step: a product formula of order $p$ has error
scaling $\varepsilon(s)=\mathcal{O}(s^{p+1})$. More precisely, the operator
error satisfies Childs  _et al._ (2021):

$\displaystyle\left\|e^{-iHs}-\mathcal{W}(s)\right\|=\mathcal{O}\big{(}s^{p+1}\,\sum_{m=1}^{M}\sum_{\\{m_{i}\\}}f_{\\{m_{i}\\},m}$
(2)
$\displaystyle\times\left\|[H_{m_{p}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\right\|\big{)},$

where the $m$ and $\\{m_{i}\\}={m_{1},\ldots,m_{p}}$ summations run over all
terms of the Hamiltonian. The coefficients $f_{\\{m_{i}\\},m}$ are
$\mathcal{O}(1)$ positive numbers that are determined based on the product
formula, i.e., the values of $a_{j}$.

We present our general results based on the above settings and then specialize
to the case in which we take all individual terms $H_{m}$ to be positive
semidefinite. Note that with the latter assumption, the setting becomes
similar to that of Ref. Şahinoğlu and Somma (2021). Furthermore, this is the
scenario for which the results of Table 1 are computed using our method. As
discussed in Ref. Şahinoğlu and Somma (2021), the assumption of positive
semidefinite holds for certain classes of spin Hamiltonians, e.g. frustration-
free systems Bravyi and Terhal (2010).

As we will discuss later, the setting can also be extended to include
fermionic systems in which the Hamiltonian consists of a sum of k-local terms
made of an even number of creation and annihilation fermionic operators, with
each fermionic mode participating in at most $d$ interactions. Notably, this
extends the results of this paper to wide classes of fermionic systems such as
electronic structure Hamiltonians.

##  II Low-energy Error

We calculate the error for the product formula in Eq. (1) acting on a state
with negligible support beyond an energy cutoff $\Delta$. We characterize this
error by calculating the following operator norm for a single Trotter step of
length $s$. All norms are spectral norms in this work.

$\displaystyle\varepsilon_{\Delta}(s)$
$\displaystyle=\left\|\left(\mathcal{W}(s)-e^{-iHs}\right)\Pi_{\leq\Delta}\right\|.$
(3)

We use the notation $\Pi_{\leq\Delta}$ and $\Pi_{>\Delta^{\prime}}$ to denote
the projectors onto the subspaces spanned by eigenstates of $H$ with energy
$E\leq\Delta$ and $E>\Delta^{\prime}$. We occasionally abuse this notation to
also refer to the subspaces themselves. We also use the notation
$\|A\|_{\leq\Lambda}$ to represent the restricted norm
$\|\Pi_{\leq\Lambda}\;A\;\Pi_{\leq\Lambda}\|$ to such low-energy subspaces.

One of our main results for bounding $\varepsilon_{\Delta}$ is the following:
the right-hand side of Eq. (2) maintains its accuracy in error representation
for initial states with negligible support on $\Pi_{>\Delta}$, if all the
norms are replaced by restricted norms to a suitably chosen low-energy
subspace; this amounts to using norms of projected nested commutators
$\|[H_{m_{p}},\ldots,[H_{m_{1}},H_{m}]]\|_{\leq\Delta^{\prime}}$ for some
appropriate value of $\Delta^{\prime}$, in the right-hand side of Eq. (2).
This modification will result in a lower overall error and estimated cost as
projected nested commutators have smaller norms. We identify suitable
$\Delta^{\prime}$ values in the following.

Separating the error terms. We now turn to the study of
$\varepsilon_{\Delta}(s)$. Using the triangle inequality, this error is upper
bounded by:

$\displaystyle\left\|\Pi_{\leq\Delta^{\prime}}\left(\mathcal{W}(s)-e^{-iHs}\right)\Pi_{\leq\Delta}\right\|+\left\|\Pi_{>\Delta^{\prime}}\mathcal{W}(s)\Pi_{\leq\Delta}\right\|,$
(4)

with $\Delta^{\prime}>\Delta$ to be determined later. We dropped the exact
time evolution operator $e^{-iHs}$ in the second term, as it commutes with the
two projectors and $\Pi_{>\Delta^{\prime}}\Pi_{\leq\Delta}=0$.

The above equation implies that the error can be calculated as a sum of two
separate terms: a part confined to energies below $\Delta^{\prime}$, which we
dub the retained component, and the leakage by $\mathcal{W}(s)$ to
$\Pi_{>\Delta^{\prime}}$. Below, we first state the final bounds for these two
terms for an order $p$ product formula, and then discuss the ideas behind
their proofs.

Product formula leakage. We show in SM  that:

$\left\|\Pi_{>\Delta^{\prime}}\mathcal{W}(s)\Pi_{\leq\Delta}\right\|\leq
e^{-\lambda\left[\Delta^{\prime}-\Delta-2R^{\mathcal{W}}(s)\right]},$ (5)

with $\lambda=\frac{1}{2Jdk}$ and $R^{\mathcal{W}}(s)$ a positive number,
which for an order $p$ product formula takes the form (see SM  for details):

$R^{\mathcal{W}}(s)=\alpha\,s^{p+1}L+\mathcal{O}(Ls^{p+2}),$ (6)

with
$\alpha=\frac{\ell_{p}}{\lambda(p+1)}\sum_{m,\\{m_{i}\\}}|f_{\\{m_{i}\\},m}|$
and $\ell_{p}=L\,p!\left(Jkde\right)^{p}eJ$. We note here that this scaling is
better than the scaling found in Şahinoğlu and Somma (2021) which is
$R^{\mathcal{W}}=\mathcal{O}(Ls)$, as $s$ is considered small.

The retained component. The first term in Eq. (4) can be generally bounded as
follows, as we show in SM :

$\displaystyle\left\|\Pi_{\leq\Delta^{\prime}}\left(\mathcal{W}(s)-e^{-iHs}\right)\Pi_{\leq\Delta}\right\|=$
(7)
$\displaystyle\mathcal{O}\left(s^{p+1}\,\sum_{m,\\{m_{i}\\}}\left\|[H_{m_{p}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\right\|_{\leq\Delta^{\prime}}\right),$

with $\Delta^{\prime}$ chosen large enough so that
$\Delta^{\prime}-\Delta=\mathcal{O}(s^{p+1}L)$ up to polylogarithmic
corrections in $L$ and $s$.

From the above two bounds, the total error $\varepsilon_{\Delta}(s)$ can also
be bounded as follows: we note that the leakage component shown in Eq. (5) can
also become subdominant to Eq. (7), provided that
$\Delta^{\prime}-\Delta=\mathcal{O}(s^{p+1}L)$ up to polylogarithmic
corrections in $L$ and $s$; interestingly this is the same requirement for Eq.
(7) to hold. This means that with such $\Delta^{\prime}$, the whole low-energy
error $\varepsilon_{\Delta}(s)$ of Eq. (3) is dominated by the right-hand side
of Eq. (7).

Remarkably, this is showing that $\varepsilon_{\Delta}$ can be calculated with
the use of operator norms of restricted nested commutators into a subspace of
energies of at most $\Delta^{\prime}$, provided $\Delta^{\prime}$ is chosen
according to the above discussion. Note that this projection results in
tighter bounds when contrasted with taking the operator norm without
projection as in the general result of Eq. (2). We use this error analysis
below for the case of Hamiltonians with positive semidefinite terms and
calculate the resulting implementation cost later in the text.

Positive semidefinite Hamiltonian terms. The retained component bound in Eq.
(7) can be made more concrete when the $H_{m}$’s are positive semidefinite SM
:

$\left\|\Pi_{\leq\Delta^{\prime}}\left(\mathcal{W}(s)-e^{-iHs}\right)\Pi_{\leq\Delta}\right\|=\mathcal{O}(s^{p+1}\Delta_{f}^{p+1}),$
(8)

where $\Delta_{f}>\Delta^{\prime}$ is also introduced to suppress more leakage
norms, and it is required that both $\Delta^{\prime}$ and $\Delta_{f}$ equal
$\Delta+\mathcal{O}(s^{p+1}L)$ up to polylogarithmic corrections in $L$ and
$s$ ($\Delta^{\prime}$ needs to be larger than $\Delta_{f}$ only by
polylogarithmic factors in $L$ and $s$.).

Again, the leakage error in Eq. (5) becomes subdominant with respect to the
error in Eq. (8) with the same condition on $\Delta^{\prime}$. This bounds the
error in each step as:

$\varepsilon_{\Delta}(s)=\mathcal{O}(s^{p+1}\Delta_{f}^{p+1}),$ (9)

with $\Delta_{f}-\Delta=\mathcal{O}(s^{p+1}L)$ up to polylogarithmic factors
in $L$ and $s$.

##  III Ideas behind the proofs

Let us now discuss the ideas behind the proofs. First, similar to Şahinoğlu
and Somma (2021), we prominently use theorem 2.1 of Arad  _et al._ (2016),
which allows one to bound the leakage by some arbitrary $A$ from energies
below some cutoff $\Lambda$ to those higher than some other cutoff
$\Lambda^{\prime}>\Lambda$:

$\left\|\Pi_{>\Lambda^{\prime}}\,A\,\Pi_{\leq\Lambda}\right\|\leq\|A\|\
e^{-\lambda(\Lambda^{\prime}-\Lambda-2R)},$ (10)

where $\lambda=\frac{1}{2gk}$, $g$ an upper bound for the sums of strengths of
all terms in the Hamiltonian involving a qubit, and $R$ is the sum of the
norms of all local terms in the Hamiltonian that do not commute with $A$. More
concretely, given the types of Hamiltonians assumed in this work, $g=Jd$, and
for a single $k$-local term in the Hamiltonian, $R$ can be upper bounded by
$kJd$.

Ref. Şahinoğlu and Somma (2021) also uses the same leakage bound for
estimating $\varepsilon_{\Delta}(s)$ in Eq. (3): the leakage due to the
product formula is calculated through bounding the leakage of each single
exponential in the expression Eq. (1) and then summing all of the
contributions. We show that this overestimates the leakage and hence the total
error. Concretely, the leakage happens according to Eq. (5) with
$R^{\mathcal{W}}=\mathcal{O}(s^{p+1}L)$ instead of the scaling
$R^{\mathcal{W}}=\mathcal{O}(sL)$ of Şahinoğlu and Somma (2021); this directly
results in the better cost scalings shown in Table 1. To reach this result, we
take a different approach for bounding the leakage by noting that the product
of all the exponentials in Eq. (1) is close to the exact time evolution
operator $e^{-iHs}$, which does not have any leakage; in fact, the two only
differ by terms that scale as $\mathcal{O}(s^{p+1})$ for an order $p$ product
formula Childs  _et al._ (2021). Therefore, we focus on bounding directly the
leakage of the entire product instead of each single exponential, and
subsequently obtain tighter bounds for this error.

To implement this approach, we use the theory of Trotter error developed in
Childs  _et al._ (2021) and write $\mathcal{W}(s)$ as a time-ordered integral:

$\mathcal{W}(s)=\mathcal{T}e^{-i\int_{0}^{s}d\sigma\;\left(H+\mathcal{E}(\sigma)\right)},$
(11)

where $\mathcal{T}$ is the time ordering operator and $\mathcal{E}(\sigma)$ is
written solely in terms of nested commutators of the terms of the Hamiltonian.
As we show in SM , the product formula error $\|e^{-iHs}-\mathcal{W}(s)\|$ can
be upper bounded by $\int_{0}^{s}d\sigma\;\left\|\mathcal{E}(\sigma)\right\|$,
which results in the same nested commutator error bounds of Childs  _et al._
(2021). This is because for a product formula of order $p$, the nested
commutators appearing in $\mathcal{E}(\sigma)$ are of order at least $p$ and
$\mathcal{E}(\sigma)=\mathcal{O}(\sigma^{p})$; indeed, to lowest order in $s$,
$\int_{0}^{s}d\sigma\;\left\|\mathcal{E}(\sigma)\right\|$ results in Eq. (2)
(details in SM ).

Using Eq. (11), expanding the time ordered exponential as a Dyson series and
analyzing terms separately, we show that the leakage due to the product
formula is given by Eq. (5). In particular, since all the terms in
$\mathcal{E}$ are nested commutators of the terms of the Hamiltonian, and thus
all of them scale linearly with system size, we observe that
$R^{\mathcal{W}}(s)$ in Eq. (6) also scales linearly with system size. Indeed,
the fact that the lowest order contribution to $R^{\mathcal{W}}(s)$ is
$\mathcal{O}(s^{p+1})$ in Eq. (6) is also a direct consequence of
$\mathcal{E}(\sigma)$ being $\mathcal{O}(\sigma^{p})$.

For the retained component, a similar approach applies, where we show in SM
that
$\left\|\Pi_{\leq\Delta^{\prime}}\left(\mathcal{W}(s)-e^{-iHs}\right)\Pi_{\leq\Delta}\right\|$
can be upper bounded by:

$\int_{0}^{s}d\sigma\;\left\|\mathcal{E}(\sigma)\right\|_{\leq\Delta^{\prime}}+\xi_{1}\left\|\Pi_{>\Delta^{\prime}}\mathcal{W}(s)\Pi_{\leq\Delta}\right\|.$
(12)

with $\xi_{1}=\mathcal{O}(s^{p+1}L)$. Interestingly, with taking
$\Delta^{\prime}-\Delta=\mathcal{O}(s^{p+1}L)$ up to polylogarithmic factors,
the second term, which is similar to the above leakage term, becomes
subdominant with respect to
$\int_{0}^{s}d\sigma\;\left\|\mathcal{E}(\sigma)\right\|_{\leq\Delta^{\prime}}$,
reducing the sum to the right-hand side of Eq. (7) to lowest order is $s$.
Note that the same condition on $\Delta^{\prime}$ also suppresses the leakage
term in Eq. (4), thus making
$\int_{0}^{s}d\sigma\;\left\|\mathcal{E}(\sigma)\right\|_{\leq\Delta^{\prime}}$
dominate the error $\varepsilon_{\Delta}$, which is as stated before one of
our main results. This is opposed to
$\int_{0}^{s}d\sigma\;\left\|\mathcal{E}(\sigma)\right\|$, which was mentioned
above as the whole operator error.

Furthermore, when terms $H_{m}$ of the Hamiltonian are positive semidefinite,
Eq. (7) can result in Eq. (8) by resolving the nested commutator (i.e. using
$[X,Y]\to XY-YX$) step by step; this is done by breaking every resolved
commutator into a leakage term and a retained component by including a series
of cutoffs $\Delta^{\prime}<\Delta_{1}<\ldots<\Delta_{p}=\Delta_{f}$, and
choosing them in a way so that all the leakages are negligible. As we show in
SM , if $\Delta_{f}-\Delta=\mathcal{O}(s^{p+1}L)$ up to polylogarithmic
corrections, all the leakage terms will be suppressed. With all leakage terms
negligible, the norm of every order $p$ nested commutator essentially reduces
to a product of a number $p$ of restricted norms such as
$\|H_{m}\|_{\leq\Delta_{f}}$; each one of these can be bounded by $\Delta_{f}$
because of their positive semi-definiteness. A similar discussion for non-
positive-semidefinite $H_{m}$ terms is presented in SM .

One last point to stress here is that all the steps of the proofs also hold
for fermionic Hamiltonians whose single terms consist of at most $k$ fermionic
operators and each fermionic mode enters at most $d$ interactions. Crucially,
the main theorem of Arad  _et al._ (2016), i.e., Eq. (10), also holds for such
fermionic Hamiltonians, as the entire proof only depends on the commutativity
of Hamiltonian terms and the operator $A$. This generalizes the above results
to fermionic Hamiltonians as well.

##  IV The cost

Figure 1: Scaling of the exponent of the system size $N$ in the cost (Trotter
number) based on the general case Childs  _et al._ (2021), previous low-energy
results Şahinoğlu and Somma (2021), and the present work. The best scaling
corresponds to the present work.

We focus on Hamiltonians with positive semidefinite terms here. We discussed
above how the error of a single step can be bounded as in Eq. (9), with
$\Delta_{f}-\Delta=\mathcal{O}(s^{p+1}N)$ up to polylogarithmic factors in $s$
and $N$, for an order $p$ product formula; $L$ is assumed to scale as the
system size $N$. For a quantum simulation of time $T$, with total tolerated
error $\epsilon$, we break the time evolution into $r$ segments, each with
length $s=T/r$ and require that the sum of the errors of segments is equal to
$\epsilon$, i.e. $\varepsilon_{\Delta}\left(T/r\right)r=\epsilon$.

In order to find asymptotic cost scaling, we take the effect of the two terms
in $\Delta_{f}$ separately (assuming each one is the dominant source of cost
in turn), and find an asymptotic form for the required number of Trotter steps
$r$ so that the error stays at the tolerated level. Using Eq. (9), we consider
a Trotter number $r$ that satisfies both
$\left[(T/r)\,\Delta\right]^{p+1}=\mathcal{O}(\epsilon/r),\
\left[(T/r)\,(\Delta_{f}-\Delta)\right]^{p+1}=\mathcal{O}(\epsilon/r)$, which
results in:

$r=\mathcal{O}\left(\frac{(T\Delta)^{1+1/p}}{\epsilon^{1/p}}+T\left[\frac{TN^{p+1}}{\epsilon}\right]^{\frac{1}{(p+1)^{2}+p}}\right),$
(13)

up to polylogarithmic factors in the parameters.

Now, we assume that $\Delta$ has a slow enough growth with $N$, such that the
main scaling of the cost with $N$ originates from the second term in Eq. (13).
This scaling, which reads $N^{1/(p+1+\frac{p}{p+1})}$, is better than the
general cost scaling of $N^{1/p}$ Childs  _et al._ (2021), and the previous
low-energy result $N^{1/2+1/(4p+2)}$ Şahinoğlu and Somma (2021). Note that the
previous low-energy result of Ref. Şahinoğlu and Somma (2021) only beats the
general case in terms of scaling with $N$ for the first order product
formulas, while our work provides a better cost for all orders. The advantage
is due to the better leakage bound found in Eq. (5), where
$R^{\mathcal{W}}=\mathcal{O}(s^{p+1}N)$ rather the $\mathcal{O}(sN)$ scaling
of Şahinoğlu and Somma (2021).

A comparison between present results and those in Refs. Childs  _et al._
(2021) and Şahinoğlu and Somma (2021) is presented in Table 1. Exponents with
respect to $N$ for all the cases are depicted in Fig. 1.

##  V Discussion

In this work, we establish state-of-the-art error bounds for product formulas
acting on low-energy states. In particular, we showed how the previous general
nested commutator error bounds of Childs  _et al._ (2021) can be modified to
obtain tighter bounds for low-energy states, i.e., by calculating operator
norms restricted to lower-energy subspaces. Furthermore, we showed that for
Hamiltonians with positive semidefinite terms, the scaling of the Trotter
number with respect to the system size is asymptotically better than the
general scaling results Childs  _et al._ (2021) and previous results in the
same setting Şahinoğlu and Somma (2021). More precisely, the previous results
of Ref. Şahinoğlu and Somma (2021) only showed improvement in the cost for the
first order Trotter-Suzuki formulas, while the current work shows that the
improvement can be obtained for all orders; this is especially important as
second-order Trotter-Suzuki formulas can be utilized with the same resource
requirement as the first order.

Even though the error analysis for Hamiltonians with positive semidefinite
terms has been considered here, cost estimation based on the more general
error analysis in Eq. (7) remains an open problem, and we leave it to future
work to explore implications of this more general result. It would also be
interesting to consider other possible structures in the low-energy states of
interest, e.g., a finite correlation length Hastings (2004); Anshu (2016) and
similar properties. Application of the current bound to time-dependent
Hamiltonian simulation problems as those arising in, e.g., adiabatic quantum
computing, can be an interesting topic to consider.

## References

  * Kitaev (1995) A. Y. Kitaev, arXiv preprint quant-ph/9511026 (1995).
  * Cleve  _et al._ (1998) R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 339 (1998).
  * Lin and Tong (2022) L. Lin and Y. Tong, PRX Quantum 3, 010318 (2022).
  * Flannigan  _et al._ (2022) S. Flannigan, N. Pearson, G. H. Low, A. Buyskikh, I. Bloch, P. Zoller, M. Troyer, and A. J. Daley, Quantum Science and Technology 7, 045025 (2022).
  * Harrow  _et al._ (2009) A. W. Harrow, A. Hassidim, and S. Lloyd, Physical review letters 103, 150502 (2009).
  * Albash and Lidar (2018) T. Albash and D. A. Lidar, Reviews of Modern Physics 90, 015002 (2018).
  * Lloyd (1996) S. Lloyd, Science 273, 1073 (1996).
  * Suzuki (1991) M. Suzuki, Journal of Mathematical Physics 32, 400 (1991).
  * Berry  _et al._ (2015) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Physical review letters 114, 090502 (2015).
  * Low and Chuang (2017) G. H. Low and I. L. Chuang, Physical review letters 118, 010501 (2017).
  * Low and Chuang (2019) G. H. Low and I. L. Chuang, Quantum 3, 163 (2019).
  * Şahinoğlu and Somma (2021) B. Şahinoğlu and R. D. Somma, npj Quantum Information 7, 119 (2021).
  * Su  _et al._ (2021) Y. Su, H.-Y. Huang, and E. T. Campbell, Quantum 5, 495 (2021).
  * Jordan  _et al._ (2012) S. P. Jordan, K. S. Lee, and J. Preskill, Science 336, 1130 (2012).
  * Tong  _et al._ (2021) Y. Tong, D. An, N. Wiebe, and L. Lin, Physical Review A 104, 032422 (2021).
  * Childs  _et al._ (2021) A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Physical Review X 11, 011020 (2021).
  * (17) See the Supplemental material .
  * Note (1) More precisely, $J:=2\max_{h}\min_{c}\|h+c\|$, with $h$ running over single local terms of the Hamiltonian and $c$ a constant, see SM .
  * Bravyi and Terhal (2010) S. Bravyi and B. Terhal, Siam journal on computing 39, 1462 (2010).
  * Arad  _et al._ (2016) I. Arad, T. Kuwahara, and Z. Landau, Journal of Statistical Mechanics: Theory and Experiment 2016, 033301 (2016).
  * Hastings (2004) M. B. Hastings, Physical review B 69, 104431 (2004).
  * Anshu (2016) A. Anshu, New Journal of Physics 18, 083011 (2016).

Supplemental material

##  VI Setting

We consider time-independent qubit Hamiltonians written in the form
$H=\sum_{m=1}^{M}H_{m}$. Each term $H_{m}$ comprises a number $L_{m}$ of
$k$-local interactions that commute with each other; while $L_{m}$ can scale
with the total number of qubits $N$, we take $M$ to be $O(1)$ and let
$L=\max_{m}L_{m}$. Each qubit is involved in at most $d$ of the $k$-local
interactions. Furthermore, we introduce the parameter $J$ to denote the
interaction strength. More precisely, $J:=2\max_{h}\min_{c}\|h+c\|$, with $h$
running over single local terms of the Hamiltonian and $c$ a constant.

Considering the general product formula:

$\mathcal{W}(s)=e^{-is_{q}H_{m_{q}}}\ldots
e^{-is_{2}H_{m_{2}}}\,e^{-is_{1}H_{m_{1}}},$ (14)

that approximates time evolution $e^{-iHs}$ under $H$. Here,
$m_{1},\ldots,m_{q}$ is a permutation of $1,\ldots,M$, with repetitions in the
case of higher-order product formulas. Each $s_{j}=a_{j}s$ is proportional to
$s$ with $a_{j}\in\mathbb{R}$ and the total number $q$ of simple exponentials
in Eq. (14) is taken to be $O(1)$.

##  VII The error operator

In this section, we discuss general forms for the error operator of a product
formula, the discussion is based on Ref. Childs  _et al._ (2021), and we
include it here for completeness.

Considering the general product formula in (14), we can write its derivative
obtained by directly differentiating Childs  _et al._ (2021) as:

$\frac{d}{ds}\mathcal{W}(s)=-i\mathcal{F}(s)\,\mathcal{W}(s),$ (15)

where

$\displaystyle\mathcal{F}(s)=\sum_{\nu=1}^{q}a_{\nu}\
V_{\nu}H_{m_{\nu}}V_{\nu}^{\dagger},\qquad\qquad
V_{\nu}=e^{-is_{q}H_{m_{q}}}\ldots e^{-is_{\nu+1}H_{m_{\nu}+1}}.$ (16)

Note that we are using $\nu$ to enumerate terms in the permutation of
Hamiltonian terms used in Eq. (14), while $m_{i}$ is used to enumerate terms
of the Hamiltonian. Now, $\mathcal{W}$ can be written as a time ordered
integral over $\mathcal{F}$:

$\mathcal{W}(s)=\mathcal{T}e^{-i\int_{0}^{s}d\sigma\;\mathcal{F}(\sigma)}.$
(17)

Using the Hadamard formula:

$\displaystyle
e^{-isA}\,B\,e^{isA}=B+(-is)[A,B]+\frac{(-is)^{2}}{2!}[A,[A,B]]+\ldots,$ (18)

we can see that $\mathcal{F}$ is an infinite collection of nested commutators;
each order of $s^{n}$ in $\mathcal{F}$ contains nested commutators of order
$n$:

$\displaystyle\mathcal{F}(s)$ $\displaystyle=H$ (19)
$\displaystyle+(-is)\sum_{\nu=1}^{q}\sum_{\nu_{1}=\nu+1}^{q}\tilde{f}_{\nu_{1},\nu}[H_{m_{\nu_{1}}},H_{m_{\nu}}]$
$\displaystyle+(-is)^{2}\sum_{\nu=1}^{q}\sum_{\nu_{1}=\nu+1}^{q}\sum_{\nu_{2}=\nu_{1}}^{q}\tilde{f}_{\nu_{2}\nu_{1},\nu}[H_{m_{\nu_{2}}},[H_{m_{\nu_{1}}},H_{m_{\nu}}]]$
$\displaystyle+(-is)^{3}\sum_{\nu=1}^{q}\sum_{\nu_{1}=\nu+1}^{q}\sum_{\nu_{2}=\nu_{1}}^{q}\sum_{\nu_{3}=\nu_{2}}^{q}\tilde{f}_{\nu_{3}\nu_{2}\nu_{1},\nu}[H_{m_{\nu_{3}}},[H_{m_{\nu_{2}}},[H_{m_{\nu_{1}}},H_{m_{\nu}}]]]$
$\displaystyle+\ldots.$

The coefficients can be calculated as:

$\displaystyle\tilde{f}_{\nu_{1},\nu}=a_{\nu}a_{\nu_{1}};$
$\displaystyle\nu_{1}>\nu,$ (20)
$\displaystyle\tilde{f}_{\nu_{2}\nu_{1},\nu}=a_{\nu}a_{\nu_{1}}a_{\nu_{2}};$
$\displaystyle\nu_{2}>\nu_{1}>\nu,\qquad\tilde{f}_{\nu_{1}\nu_{1},\nu}=\frac{1}{2!}a_{\nu}a_{\nu_{1}}^{2},$
$\displaystyle\tilde{f}_{\nu_{3}\nu_{2}\nu_{1},\nu}=a_{\nu}a_{\nu_{1}}a_{\nu_{2}}a_{\nu_{3}};$
$\displaystyle\nu_{3}>\nu_{2}>\nu_{1}>\nu,\qquad\tilde{f}_{\nu_{1}\nu_{1}\nu_{1},\nu}=\frac{1}{3!}a_{\nu}a_{\nu_{1}}^{3},\qquad\ldots$
$\displaystyle\vdots$

Note that in all of the above equations, the summation is performed with the
indices $\nu$ and $\nu_{i}$ enumerating the permutations of the Hamiltonian
terms occurring in the product formula (14). However, the summation can also
be done over the Hamiltonian term indices $m$ resulting:

$\displaystyle\mathcal{F}(s)$
$\displaystyle=H+(-is)\sum_{m=1}^{M}\sum_{m_{1}\neq
m}f_{m_{1},m}[H_{m_{1}},H_{m}]$ (21)
$\displaystyle+(-is)^{2}\sum_{m=1}^{M}\sum_{m_{1},m_{2}\neq
m}f_{m_{2}m_{1},m}[H_{m_{2}},[H_{m_{1}},H_{m}]]+\ldots.$

with $f$ and $\tilde{f}$ coefficients related through the particular
permutation used in Eq. (14).

The lowest nonvanishing order of $s$ in the expansion of $\mathcal{F}-H$
controls the error of the product formula Childs  _et al._ (2021). For a
Trotter Suzuki product formula of order $p$, it is well-known that the error
in a single step should scale as $s^{p+1}$; this means that all the terms with
nested commutators of orders less than $p+1$ will vanish in the expansion of
$\mathcal{F}-H$ Childs  _et al._ (2021). Indeed, the error induced by product
formula can be written as:

$\left\|e^{-iHs}-\mathcal{W}(s)\right\|\leq\frac{(s)^{p+1}}{p+1}\,\sum_{m=1}^{M}\sum_{\\{m_{i}\\}}f_{\\{m_{i}\\},m}\
\left\|[H_{m_{p}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\
\right\|+\mathcal{O}(s^{p+2}).$ (22)

##  VIII Leakage of the product formula

In order to calculate the leakage due to $\mathcal{W}$ when starting in a low-
energy state, we need to switch to the interaction picture of Eq. (17) (see
Appendix C of Ref. Childs  _et al._ (2021)). First we define
$\mathcal{F}(s)=H+\mathcal{E}(s)$, in the interaction picture, Eq. (17)
becomes:

$\mathcal{W}(s)=e^{-iHs}\
\mathcal{T}e^{-i\int_{0}^{s}d\sigma\;\mathcal{E}^{I}(\sigma)},$ (23)

where $\mathcal{E}^{I}(s)=e^{iHs}\,\mathcal{E}(s)\,e^{-iHs}$, which can again
be expanded using the Hadamard formula Eq. (18) to get a form for
$\mathcal{E}^{I}$ in terms of nested commutators of Hamiltonian terms. In fact
an expansion similar to Eq. (21) can also be written for $\mathcal{E}^{I}$:

$\displaystyle\mathcal{E}^{I}(s)$
$\displaystyle=(-is)\sum_{m,m_{1}}f^{I}_{m_{1},m}\,[H_{m_{1}},H_{m}]$ (24)
$\displaystyle+(-is)^{2}\sum_{m,m_{1},m_{2}}f^{I}_{m_{2}m_{1},m}\,[H_{m_{2}},[H_{m_{1}},H_{m}]]$
$\displaystyle+\ldots.$

Each nested commutator in the above expansion is obtained either directly from
Eq. (19) or is obtained through a number of commutations of $H$ with lower
order commutators in Eq. (19). Thus each coefficient $f^{I}$ can be written in
terms of coefficients $f$ of the same and lower order. Crucially, the lowest
nonzero order in the expansions of $\mathcal{E}=\mathcal{F}-H$ and
$\mathcal{E}^{I}$ match. As a result, for a Trotter-Suzuki product formula of
order $p$, the lowest order nested commutator in Eq. (24) is of order $p$.

Note that a Dyson series expansion can be used to write Eq. (23) as:

$\displaystyle\mathcal{W}(s)=e^{-iHs}\Bigg{[}1+\sum_{j}\frac{(-i)^{j}}{j!}\int_{0}^{s}d\sigma_{1}\int_{0}^{s}d\sigma_{2}\ldots\int_{0}^{s}d\sigma_{j}\;\mathcal{T}\mathcal{E}^{I}(\sigma_{1})\mathcal{E}^{I}(\sigma_{2})\ldots\mathcal{E}^{I}(\sigma_{j})\Bigg{]}.$
(25)

Now, we need to turn our attention to leakage from low energies; like Ref.
Şahinoğlu and Somma (2021), the main result we will be using is Theorem 2.1 of
Ref. Arad  _et al._ (2016). We first explain it here using the present
notation: take $\Pi_{>\Lambda^{\prime}}$ and $\Pi_{\leq\Lambda}$ to be
projectors onto subspaces spanned by eigenstates of $H$ with energies
satisfying $E>\Lambda^{\prime}$ and $E\leq\Lambda^{\prime}$ respectively. For
an operator A, the following inequality holds for its leakage from a low-
energy subspace:

$\left\|\Pi_{>\Lambda^{\prime}}\,A\,\Pi_{\leq\Lambda}\right\|\leq\|A\|\
e^{-\lambda(\Lambda^{\prime}-\Lambda-2R)},$ (26)

where $\lambda=\frac{1}{2gk}$ with $g$ an upper bound for the sums of
strengths of all terms involving each spin, and $R$ is the sum of the norms of
all local terms in the Hamiltonian that do not commute with $A$. Given the
Hamiltonian assumed in this work (sum of at most $k$-local terms and with
maximum degree $d$), $g=Jd$. For a single $k$-local term in the Hamiltonian,
$R$ can be upper bounded by $kJd$.
Note: Ref. Arad  _et al._ (2016) assumes positive semidefinite single terms
for the Hamiltonian; this is not a crucial assumption for their theorem to
work, and in fact their leakage bound proof also hold for nonpositive single
terms in the Hamiltonian. This just means that $J$ should be chosen with care:
in the worst case for a nonpositive operator we take $J$ to be twice as large
as the norm of the largest single term of the Hamiltonian; for example, if
single terms of the Hamiltonian are Pauli strings, $J$ should be taken as
twice the value of the largest coefficient. In general $J$ is chosen as
$2\max_{h}\min_{c}\|h+c\|$, with $h$ running over single local terms of the
Hamiltonian and $c$ showing a constant.

Now consider a nested commutator of order $n$ such as
$[H_{m_{n}},\ldots,[H_{m_{1}},H_{m}]]\ldots]$, we would like to calculate its
leakage from low energies similar to Eq. (26). We first note that altogether
such a nested commutator consists of at most $L(kd)^{n}n!$ nested commutators
of $k-$local terms (see next section and table 2). Each term within a nested
commutator of $k-$local terms has a support of at most $(n+1)k$, strength of
at most $J^{n+1}$ (see next Section and Table 2). As a result of the previous
observation on the length of single terms, the value of $R$ for these nested
commutator is upper bounded by $J(n+1)kd$. We can now use Eq. (26) for each
single term in this nested commutator to bound the leakage of the whole
commutator as follows by summing all of the contributions together:

$\displaystyle\left\|\Pi_{>\Lambda^{\prime}}\
[H_{m_{n}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\
\Pi_{\leq\Lambda}\right\|\leq\ell_{n}\
e^{-\lambda(\Lambda^{\prime}-\Lambda)},$ (27)

where

$\displaystyle\ell_{n}$ $\displaystyle=(L(kd)^{n}n!)\ (J^{n+1})\ e^{2\lambda
Jkd(n+1)}$ (28) $\displaystyle=L\,n!\left(Jkde\right)^{n}eJ.$

On the right hand side of the first row, the first factor denotes the number
of single terms in the expansion of the nested commutator, the second factor
shows their strength and the third corresponds to the $e^{2\lambda R}$ in the
original leakage form Eq. (26) and we have also used $\lambda=\frac{1}{2Jkd}$
for the second row.

One can also upper bound the leakage due to a product of nested commutators in
a similar way. Such products appear in the Dyson series in Eq. (25). Again we
should focus on using Eq. (26) for single terms in the expansion. If we focus
on a product of $j$ nested commutators with orders $n_{1},n_{2},\ldots,n_{j}$,
it is easy to see that the support of single terms will be upper bounded by
the sum $[(n_{1}+1)+(n_{2}+1)+\ldots+(n_{j}+1)]k$, this determines the new $R$
for the leakage of such term. On the other hand, the number of terms appearing
will be upper bounded by the product of the number of terms in each nested
commutator. Using these observations, and Eqs. (27) and (28), we deduce that
the leakage from a product of nested commutators can be upper bounded by the
multiplicative form

$\ell_{n_{1}}\ \ell_{n_{2}}\ldots\ell_{n_{j}}\
e^{-\lambda(\Lambda^{\prime}-\Lambda)}.$ (29)

Using the above observation, it is now easy to see the the terms in the Dyson
expansion Eq. (25) can be integrated and then recast into the exponential
form.

$\left\|\Pi_{>\Lambda^{\prime}}\mathcal{W}(s)\Pi_{\leq\Lambda}\right\|\leq
e^{-\lambda\left[\Lambda^{\prime}-\Lambda-2R^{\mathcal{W}}(s)\right]},$ (30)

with $R^{\mathcal{W}}(s)$ defined based on the expansion in Eq. 24:

$\displaystyle\lambda R^{\mathcal{W}}(s)$
$\displaystyle=\ell_{1}\frac{s^{2}}{2}\sum_{m,m_{1}}|f^{I}_{m_{1},m}|$ (31)
$\displaystyle+\ell_{2}\frac{s^{3}}{3}\sum_{m,m_{1},m_{2}}|f^{I}_{m_{2}m_{1},m}|$
$\displaystyle+\ldots.$

Note that we dropped $e^{-iHs}$ from the left in the expansion of
$\mathcal{W}$ as it commutes with $\Pi_{>\Lambda^{\prime}}$ and has norm equal
to 1.

Since the expansions for $\mathcal{E}^{I}$ and $\mathcal{F}$ agree to the
lowest nonzero order, $R^{\mathcal{W}}$ to lowest order for a Trotter-Suzuki
of order $p$ can be written as:

$R^{\mathcal{W}}(s)=\frac{s^{p+1}}{p+1}\,\frac{\ell_{p}}{\lambda}\sum_{m,\\{m_{i}\\}}|f_{\\{m_{i}\\},m}|+\mathcal{O}(s^{p+2}).$
(32)

$i$ in $\\{m_{i}\\}$ ranges from $1$ to $p$. One can check that
$R^{\mathcal{W}}$ scales linearly with $L$ and system size as $\ell_{p}$ does
so.

| Number of terms | Support of terms | Strength of terms
---|---|---|---
$H_{m}$ | $L$ | $k$ | $J$
$[H_{m_{1}},H_{m}]$ | $Lkd$ | $2k$ | $J^{2}$
$[H_{m_{2}},[H_{m_{1}},H_{m}]]$ | $Lkd(2kd)$ | $3k$ | $J^{3}$
$[H_{m_{3}},[H_{m_{2}},[H_{m_{1}},H_{m}]]]$ | $Lkd(2kd)(3kd)$ | $4k$ | $J^{4}$
Table 2: Upper bounds for number of terms, their lengths and their strengths
in nested commutators.

##  IX Terms of nested commutators

The number of terms each nested commutator such as
$[H_{m_{n}},\ldots,[H_{m_{1}},H_{m}]]\ldots]$ and other properties of such
terms is discussed here. Table 2 denotes upper bounds for number of terms and
lengths of terms in low order nested commutators. Every row of the table is
deduced from the previous one; e.g. if the total number of terms in a nested
commutator is $\tilde{L}_{j}$ and each term has a support of at most
$\tilde{k}_{j}$, for the next level the total number of terms will be upper
bounded by $\tilde{L}_{j}\tilde{k}_{j}d$ and the length of each term will be
upper bound by $\tilde{k}_{j}+k$. It is easy to see that for an order $n$
nested commutator, the number of terms, lengths of terms and their strengths
can be upper bounded as $L(kd)^{n}n!$, $(n+1)k$ and $J^{n+1}$ respectively.

##  X Bounding the retained component

The error in time evolution of a low-energy state which can be quantified by
the following operator norm is discussed in the main text:

$\displaystyle\varepsilon_{\Delta}(s)=\left\|\left(\mathcal{W}(s)-e^{-iHs}\right)\Pi_{\leq\Delta}\right\|.$
(33)

There, we broke the above norm into two separate terms as follows:

$\displaystyle\varepsilon_{\Delta}(s)=\left\|\Pi_{\leq\Delta^{\prime}}\left(\mathcal{W}(s)-e^{-iHs}\right)\Pi_{\leq\Delta}\right\|+\left\|\Pi_{>\Delta^{\prime}}\mathcal{W}(s)\Pi_{\leq\Delta}\right\|,$
(34)

The second term is the leakage of the product formula that was discussed
earlier.

Here, we would like to bound the following norm, which we call the retained
component:

$\varepsilon_{<}(s)=\left\|\Pi_{\leq\Delta^{\prime}}\left(\mathcal{W}(s)-e^{-iHs}\right)\Pi_{\leq\Delta}\right\|.$
(35)

First, we define $\delta\mathcal{W}(s)$:

$\mathcal{W}(s)-e^{-iHs}=e^{-iHs}\delta\mathcal{W}(s),$ (36)

which has a derivative:

$\frac{d}{ds}\delta\mathcal{W}=\frac{d}{ds}\left(e^{iHs}\mathcal{W}(s)-1\right)=-i\,e^{iHs}(\mathcal{F}(s)-H)\mathcal{W}(s),$
(37)

where we have used Eq. (15). Noting the definition
$\mathcal{E}=\mathcal{F}-H$, this lets us rewrite Eq. (36) as:

$\displaystyle\mathcal{W}(s)-e^{-iHs}=\int_{0}^{s}d\sigma\;e^{-iH(s-\sigma)}\,\mathcal{E}(\sigma)\,\mathcal{W}(\sigma).$
(38)

This equation can be used to give the total product formula error as:

$\left\|\mathcal{W}(s)-e^{-iHs}\right\|\leq\int_{0}^{s}d\sigma\;\left\|\mathcal{E}(\sigma)\right\|,$
(39)

where the fact that $e^{-iH(s-\sigma)}$ and $\mathcal{W}(\sigma)$ are unitary
operators are used. This is equivalent to the bounds established in Ref.
Childs  _et al._ (2021).

Getting back to evaluation of Eq. (35), we use Eq. (38) with a resolution of
the identity in the form of $\Pi_{\leq\Delta^{\prime}}+\Pi_{>\Delta^{\prime}}$
inserted between $\mathcal{E}(\sigma)$ and $\mathcal{W}(\sigma)$:

$\varepsilon_{<}(s)=\left\|\int_{0}^{s}d\sigma\;\Pi_{\leq\Delta^{\prime}}e^{-iH(s-\sigma)}\,\mathcal{E}(\sigma)\,\left(\Pi_{\leq\Delta^{\prime}}+\Pi_{>\Delta^{\prime}}\right)\,\mathcal{W}(\sigma)\,\Pi_{\leq\Delta}\right\|.$
(40)

We bound the above, by first dropping $e^{-iH(s-\sigma)}$ from the left as it
commutes with the projector and has a norm equal to 1; also using the triangle
inequality:

$\displaystyle\varepsilon_{<}(s)\leq\int_{0}^{s}d\sigma\bigg{[}\left\|\Pi_{\leq\Delta^{\prime}}\,\mathcal{E}(\sigma)\,\Pi_{>\Delta^{\prime}}\,\mathcal{W}(\sigma)\,\Pi_{\leq\Delta}\right\|+\left\|\Pi_{\leq\Delta^{\prime}}\,\mathcal{E}(\sigma)\,\Pi_{\leq\Delta^{\prime}}\,\mathcal{W}(\sigma)\,\Pi_{\leq\Delta}\right\|\bigg{]}.$
(41)

We bound the two terms separately below. Beginning with the first term in Eq.
(41), we first note:

$\int_{0}^{s}d\sigma\;\left\|\left(\Pi_{\leq\Delta^{\prime}}\,\mathcal{E}(\sigma)\,\Pi_{>\Delta^{\prime}}\,\mathcal{W}(\sigma)\,\Pi_{\leq\Delta}\right)\right\|\leq\int_{0}^{s}d\sigma\;\left\|\left(\Pi_{\leq\Delta^{\prime}}\,\mathcal{E}(\sigma)\,\Pi_{>\Delta^{\prime}}\right\|\left\|\Pi_{>\Delta^{\prime}}\,\mathcal{W}(\sigma)\,\Pi_{\leq\Delta}\right)\right\|,$
(42)

for which the leakage bound for $\mathcal{W}$ along with the fact that
$\left\|\Pi_{\leq\Delta^{\prime}}\,\mathcal{E}(\sigma)\,\Pi_{>\Delta^{\prime}}\right\|$
is bounded by $\|\mathcal{E}(\sigma)\|$ can be used. For an order $p$ product
formula:

$\displaystyle\|\mathcal{E}(\sigma)\|\leq\sigma^{p}J^{p+1}L\,\gamma_{p}\sum_{m,\\{m_{i}\\}}|f_{\\{m_{i}\\},m}|+\mathcal{O}(\sigma^{p+1}),$
(43)

with $\gamma_{p}=(kd)^{p}p!$ (see Table 2). $i$ in $\\{m_{i}\\}$ ranges from
$1$ to $p$. As a result of this and using Eqs. (30) and (32), the whole first
term can be bounded as:

$\displaystyle\int_{0}^{s}d\sigma\;\left\|\left(\Pi_{\leq\Delta^{\prime}}\,\mathcal{E}(\sigma)\,\Pi_{>\Delta^{\prime}}\,\mathcal{W}(\sigma)\,\Pi_{\leq\Delta}\right)\right\|\leq\xi_{1}\,e^{-\lambda(\Delta^{\prime}-\Delta-
R_{1})},$ (44)

where $\xi_{1}=\mathcal{O}(s^{p+1}J^{p+1}L)$ and
$R_{1}=\mathcal{O}(s^{p+1}J^{p+2}L)$. The integral over $\sigma$ has been
upper bounded by multiplying the maximum value of the integrand (which happens
at $\sigma=s$) by $s$. $\xi_{1}$ can be absorbed into $R_{1}$ resulting in a
logarithmic correction in $s$ and $L$ only. This means that, as we will
discuss more below, this leakage term can become subdominant with respect to
the first term by taking $\Delta^{\prime}-\Delta=\mathcal{O}(s^{p+1}L)$ up to
polylogarithmic factors in $s$ and $L$. Note that the same condition also
guarantees suppression of the second term in Eq. (34), i.e. the leakage of the
product formula according to Eqs. (30) and (32).

Turning to the second term in Eq. (41), this time we bound
$\|\Pi_{\leq\Delta^{\prime}}\,\mathcal{W}(\sigma)\,\Pi_{\leq\Delta}\|$ by the
norm of $\mathcal{W}(\sigma)$ which is equal to 1. Thus the second row of Eq.
(41) can be upper bounded by:

$\int_{0}^{s}d\sigma\;\left\|\Pi_{\leq\Delta^{\prime}}\,\mathcal{E}(\sigma)\,\Pi_{\leq\Delta^{\prime}}\right\|$
(45)

First we discuss this general form. It is showing us that up to an error
$\vartheta$, the expression in Eq. (45) can represent
$\varepsilon_{\Delta}(s)$ in Eq. (33) if $\Delta^{\prime}$ is chosen such that
$\Delta^{\prime}-\Delta=\mathcal{O}(s^{p+1}L+\log(1/\vartheta))$ up to
polylogarithmic factors in $s$ and $L$. The condition on $\Delta^{\prime}$
ensures that both of the leakage terms, i.e. Eq. (44) and the second term on
the right hand side of (34), are subdominant error terms compared with Eq.
(45) and in fact bounded by $\vartheta$.

Note that this is a general result and it is not dependent on the terms of the
Hamiltonian being positive semidefinite. It is showing that for a product
formula the error for low energies can be computed by projecting the operator
$\mathcal{E}$ from both sides to a low-energy subspace, provided that
$\Delta^{\prime}$ is taken to be large enough as above. If this equation is
considered to lowest order in $s$, terms with nested commutators of order $p$
such as $[H_{m_{p}},\ldots,[H_{m_{1}},H_{m}]]\ldots]$ remain, and we will
obtain:

$\frac{s^{p+1}}{p+1}\,\sum_{m=1}^{M}\sum_{\\{m_{i}\\}}f_{\\{m_{i}\\},m}\
\left\|\Pi_{\leq\Delta^{\prime}}\ [H_{m_{p}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\
\Pi_{\leq\Delta^{\prime}}\right\|+\mathcal{O}(s^{p+2}).$ (46)

Compare this result with Eq. (22). In the main text the $f$ coefficients are
dropped as they are $\mathcal{O}(1)$.

The form in Eq. (46) can also be recast in a way that is useful for
Hamiltonian with terms $H_{m}$ being positive semidefinite. From here on, a
number of energy cutoffs denoted as
$\Delta^{\prime}<\Delta_{1}<\Delta_{2}<\ldots<\Delta_{p-1}=\Delta_{f}$ should
be introduced in order to use leakage bounds. Using resolutions of identity we
can recursively reduce the nested commutator by noting that:

$\displaystyle\left\|\Pi_{\leq\Delta^{\prime}}\,[H_{m_{p}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\,\Pi_{\leq\Delta^{\prime}}\right\|$
$\displaystyle\leq\left\|\Pi_{\leq\Delta^{\prime}}\,H_{m_{p}}(\Pi_{\leq\Delta_{1}}+\Pi_{>\Delta_{1}})[H_{m_{p-1}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\,\Pi_{\leq\Delta^{\prime}}\right\|$
(47)
$\displaystyle+\left\|\Pi_{\leq\Delta^{\prime}}\,[H_{m_{p-1}},\ldots,[H_{m_{1}},H_{m}]]\ldots](\Pi_{\leq\Delta_{1}}+\Pi_{>\Delta_{1}})H_{m_{p}}\,\Pi_{\leq\Delta^{\prime}}\right\|$
$\displaystyle\leq
2\left\|\Pi_{\leq\Delta_{1}}\,H_{m_{p}}\Pi_{\leq\Delta_{1}}\right\|\left\|\Pi_{\leq\Delta_{1}}\,[H_{m_{p-1}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\Pi_{\leq\Delta_{1}}\right\|$
$\displaystyle+2\left\|\Pi_{>\Delta_{1}}\,H_{m_{p}}\Pi_{\leq\Delta^{\prime}}\right\|\left\|\Pi_{>\Delta_{1}}\,[H_{m_{p-1}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\Pi_{\leq\Delta^{\prime}}\right\|$

This procedure can be continued until all nested commutators are resolved; the
terms on the row before the last are projected to a low-energy subspace, but
we also need to bound the terms on the last row of Eq. (47), which are
leakages; we can use Eq. (27) to this end.

$\left\|\Pi_{>\Delta_{1}}\,H_{m_{p}}\Pi_{\leq\Delta^{\prime}}\right\|\left\|\Pi_{>\Delta_{1}}\,[H_{m_{p-1}},\ldots,[H_{m_{1}},H_{m}]]\ldots]\Pi_{\leq\Delta^{\prime}}\right\|\leq\ell_{0}\ell_{p-1}e^{-2\lambda(\Delta_{1}-\Delta^{\prime})},$
(48)

with $\ell_{n}$ defined in Eq. (28). Note that the prefactor
$\ell_{0}\ell_{p-1}$ scales as $\mathcal{O}(L^{2})$.

This procedure should be continued until only norms of single terms of the
Hamiltonian (without commutators) projected to low energies remain:

$\displaystyle\frac{s^{p+1}}{p+1}\bigg{(}$ $\displaystyle 2\ell_{0}\ell_{p-1}\
e^{-2\lambda(\Delta_{1}-\Delta^{\prime})}$ (49)
$\displaystyle+2^{2}\left\|\Pi_{\leq\Delta_{1}}\,H_{m_{p}}\Pi_{\leq\Delta_{1}}\right\|\
\ell_{0}\ell_{p-2}\ e^{-2\lambda(\Delta_{2}-\Delta_{1})}$
$\displaystyle+\ldots$
$\displaystyle+2^{p-1}\left\|\Pi_{\leq\Delta_{1}}\,H_{m_{p}}\Pi_{\leq\Delta_{1}}\right\|\ldots\left\|\Pi_{\leq\Delta_{p-2}}\,H_{m_{1}}\Pi_{\leq\Delta_{p-2}}\right\|\
\ell_{0}\ell_{1}\ e^{-2\lambda(\Delta_{p-1}-\Delta_{p-2})}$
$\displaystyle+2^{p}\left\|\Pi_{\leq\Delta_{1}}\,H_{m_{p}}\Pi_{\leq\Delta_{1}}\right\|\ldots\left\|\Pi_{\leq\Delta_{p-2}}\,H_{m_{1}}\Pi_{\leq\Delta_{p-2}}\right\|\left\|\Pi_{\leq\Delta_{f}}\,H_{m}\Pi_{\leq\Delta_{f}}\right\|\bigg{)}.$

Ultimately we would be interested in choosing the energy cutoffs in a way that
the exponential factors make all the terms subleading with respect to the last
term in Eq. (49). This can be done by taking the differences between adjacent
cutoffs to be large enough, since all the prefactors of the exponentials scale
as powers of $L$, the difference between adjacent cutoffs suffices to grow
only logarithmically with $L$. This means that $\Delta_{f}$ and
$\Delta^{\prime}$ only differ by polylogarithmic factors and thus
$\Delta_{f}-\Delta=\mathcal{O}(s^{p+1}L)$ up to polylogarithmic factors in $s$
and $L$ will guarantee suppression of all the leakage terms. Note that it is
simple to extend this analysis to include orders beyond the leading order by
adding more energy cutoffs, and $\Delta_{f}$ will only change
polylogarithmically.

This leaves us with the following two forms for the error:

  * •

Positive semidefinite $H_{m}$ terms: each operator norm in the above relation
can be bounded by $\Delta_{f}$ for Hamiltonians with positive semidefinite
terms $H_{m}$. As a result, the error in a single step takes the form:

$\varepsilon_{\Delta}(s)=\mathcal{O}(s^{p+1}\Delta_{f}^{p+1}).$ (50)

Notice that while we can add constants to the $H_{m}$’s to make them positive
semidefinite, that would likely make $\Delta$ acquire an undesirable scaling
with $N$ (this is similar to the argument in Şahinoğlu and Somma (2021)),
cancelling the advantage of the low-energy error analysis.

  * •

Non-positive semidefinite $H_{m}$ terms: in this case, we need to define the
quantity
$\tilde{\Delta}_{f}=\max_{m}\min_{C}\|\Pi_{\leq\Delta_{f}}(H_{m}+C)\Pi_{\leq\Delta_{f}}\|$
(this is similar to the discussion in Şahinoğlu and Somma (2021)) with $C$ a
constant. Note that we have added a constant, because the above bound is found
through the use of commutators and adding constants leaves them unchanged. In
this case, the error in a single step will take the form:

$\varepsilon_{\Delta}(s)=\mathcal{O}(s^{p+1}\tilde{\Delta}_{f}^{p+1}).$ (51)

The error in Eq. (50) has been used in the main text to calculate the cost of
performing low-energy time evolution for relevant systems. However, the cost
scaling that corresponds to the error in Eq. (51) is not as simple to
calculate. The main problem is that in this case since $H_{m}$ terms are not
taken to be positive-semidefinite, it is not simple to study the behavior of
$\tilde{\Delta}_{f}$ with system properties such as its size. Characterization
of this parameter for Hamiltonians of interest can be an compelling subject of
future study.