# Provably accurate simulation of gauge theories and bosonic systems

Yu Tong  Google Quantum AI, Venice, CA, USA  Department of Mathematics,
University of California, Berkeley, CA, USA  Victor V. Albert  Joint Center
for Quantum Information and Computer Science, NIST and University of Maryland,
College Park, MD, USA  Jarrod R. McClean  Google Quantum AI, Venice, CA, USA
John Preskill  Institute for Quantum Information and Matter, Caltech,
Pasadena, CA, USA  AWS Center for Quantum Computing, Pasadena, CA, USA  Yuan
Su  Google Quantum AI, Venice, CA, USA  Institute for Quantum Information and
Matter, Caltech, Pasadena, CA, USA

(April 4th, 2022)

###### Abstract

Quantum many-body systems involving bosonic modes or gauge fields have
infinite-dimensional local Hilbert spaces which must be truncated to perform
simulations of real-time dynamics on classical or quantum computers. To
analyze errors resulting from truncation, we develop methods for bounding the
rate of growth of local quantum numbers such as the occupation number of a
mode at a lattice site, or the electric field at a lattice link. Our approach
applies to various models of bosons interacting with spins or fermions such as
the Hubbard-Holstein, Fröhlich, and Dicke models, and also to both abelian and
non-abelian gauge theories. We show that if states in these models are
truncated by imposing an upper limit $\Lambda$ on each local quantum number,
and if the initial state has low local quantum numbers, then a truncation
error no worse than $\epsilon$ can be achieved by choosing $\Lambda$ to
increase polylogarithmically with $\epsilon^{-1}$, an exponential improvement
over previous bounds based on energy conservation. For the Hubbard-Holstein
model, we numerically compute an upper bound on the value of $\Lambda$ that
achieves accuracy $\epsilon$, finding significant improvement over previous
estimates in various parameter regimes. We also establish a criterion for
truncating the Hamiltonian with a provable guarantee on the accuracy of time
evolution. Building on that result, we formulate quantum algorithms for
dynamical simulation of lattice gauge theories and of models with bosonic
modes; the gate complexity depends almost linearly on spacetime volume in the
former case, and almost quadratically on time in the latter case. We establish
a lower bound showing that there are systems involving bosons for which this
quadratic scaling with time cannot be improved. By applying our results on the
truncation error in time evolution, we also prove that spectrally isolated
energy eigenstates can be approximated with error at most $\epsilon$ by
truncating local quantum numbers at $\Lambda=\textrm{polylog}(\epsilon^{-1})$.

##  1 Introduction

Model physical systems are often formulated on spatial lattices, where the
local Hilbert space residing on each site or link of the lattice is infinite
dimensional. Examples include condensed-matter systems with bosonic degrees of
freedom [43, 21, 27, 67, 62, 64, 76, 51, 52], lattice gauge theories (LGTs)
[60, 45, 53, 75, 4, 78, 77, 5, 7, 69, 38, 61, 73, 6, 42, 41, 56, 13, 15, 20],
and other lattice field theories [36, 37]. In such models, it is convenient to
characterize the local state of the system in terms of a local quantum number,
such as the occupation number of a bosonic mode at a particular site, or the
electric field of a gauge variable at a particular link. When simulating a
lattice model using a classical or quantum computer, it is typically necessary
to truncate the local Hilbert space, replacing it by a finite-dimensional
space in which the local quantum number has a maximum value. We call this
maximum value the _truncation threshold_ , and denote it by $\Lambda$.

Quantum states of the ideal untruncated model, if concentrated on relatively
low values of the local quantum numbers, can be accurately approximated within
the truncated model. However, in a dynamical simulation governed by a
specified Hamiltonian, local quantum numbers may increase as the system
evolves. Therefore, even if the initial state is well approximated within the
truncated model, the approximation might no longer be accurate after evolution
for a sufficiently long time. To ensure that the truncated model can
accommodate the evolved state we need to bound the rate of growth of the local
quantum numbers in the ideal model.

One way to obtain such a bound is to invoke conservation of the total energy.
However, even though the total energy is conserved, the local quantum numbers
are not, and we need to worry about whether energy which is initially
distributed among many lattice sites might become focused on a much smaller
number of sites, pushing the local quantum numbers at some sites beyond the
capacity of the truncated local Hilbert space. Using conservation of energy,
combined with the Chebyshev inequality to bound the probability of large
deviations from mean values, one may infer that (for a fixed evolution time),
quantum states can be truncated with an error at most $\epsilon$ using a
threshold $\Lambda$ scaling polynomially with $\epsilon^{-1}$ [36, 37].
However, it is unclear whether this _energy-based bound_ can be used to
truncate Hamiltonians with a provable accuracy guarantee when the local
quantum numbers are not conserved under time evolution. We will further
clarify this issue in Section 2.

In this work, we develop a unified framework that shows, for a large class of
models, this energy-based estimate of $\Lambda$ is far too pessimistic — a
truncation threshold scaling as $\polylog(\epsilon^{-1})$ actually suffices,
as previously suggested in [51, 52]. This model class includes systems
involving bosons such as the Hubbard-Holstein model [33], the Fröhlich model
[23], and the Dicke model [22, 32], as well as both U(1) and SU(2) LGTs
(although our results do not apply to interacting scalar field theories such
as $\phi^{4}$ theory). For a system with many bosonic modes or gauge links,
the truncation error scales with the total number of truncated local
variables; therefore the exponentially improved dependence of $\Lambda$ on the
precision also implies exponentially improved scaling of $\Lambda$ with the
total system size. To illustrate the improvement, Figure 1 compares our
truncation threshold with the energy-based estimate for the case of the
Hubbard-Holstein model. See Section K in the Appendix for a more detailed
comparison. We further establish a threshold for truncating the Hamiltonian
such that the time evolution is provably accurate when the initial state is
assumed to have low local quantum numbers.

Previous analytical studies of the truncation problem have been mostly
restricted to simple models, while only limited small-scale numerical results
are available for more complicated systems [51, 40]. For instance, Ref. [70]
proposed one such method for simulating a single quantum harmonic oscillator.
In Ref. [51] the authors argued via the Nyquist-Shannon sampling theorem that
for a single bosonic mode with an occupation number cutoff, a grid
discretization leads to exponentially small error. This argument was further
extended in Ref. [40] to the setting of scalar field theories. The occupation
number cutoff is justified by considering a forced Harmonic oscillator, for
which analytic solution can be obtained. However, the model of a forced
Harmonic oscillator does not cover all features of boson-fermion interaction,
because by modeling the interaction between the bosonic mode and the rest of
the system by a time-dependent force, this model ignores the entanglement
between the two parts of the system. To the best of our knowledge, the
framework we develop provides the first exponential accuracy guarantee for
truncating a wide range of unbounded quantum systems of physical interest.

The new truncation threshold enables us to more accurately analyze the
computational cost of simulating dynamical evolution in the systems mentioned
above.Although we will mainly consider applications of our result to quantum
simulation, our techniques can be used to determine truncation threshold for
classical simulation as well. Using standard estimates, the simulation cost
typically depends on norms of local terms in the Hamiltonian, which are
formally infinite in bosonic systems and LGTs. We can obtain a tighter
estimate by considering evolution governed by a truncated Hamiltonian acting
on the truncated Hilbert space. We focus specifically on digital quantum
simulation of time evolution in the Hubbard-Holstein model and the U(1) and
SU(2) LGTs. For the latter, by adapting the simulation algorithm of [28] to
our truncated Hamiltonian, we find a gate complexity that scales almost
linearly with the spacetime volume. In doing so, we establish a constant Lieb-
Robinson velocity for LGT models which is essential for the method of [28] and
may be of independent interest. We also observe that there are Hamiltonians in
the class we consider such that the gate complexity of simulation for time $T$
is $\widetilde{\Omega}(T^{2})$ 111For functions of real variables $f,g$, we
write $f=\mathcal{O}(g)$ if there exist $c,t_{0}>0$ such that $|f(T)|\leq
c|g(T)|$ for all $|T|\geq t_{0}$. When there is no ambiguity, we will use
$f=\mathcal{O}(g)$ to also represent that $|f(\tau)|\leq c|g(\tau)|$ holds for
all $\tau\in\R$. We then extend the definition of $\mathcal{O}$ to functions
of positive integers and multivariate functions. For example, we use
$f(N,T,1/\epsilon)=\mathcal{O}((NT)^{2}/\epsilon)$ to mean that
$|f(N,T,1/\epsilon)|\leq c(N|T|)^{2}/\epsilon$ for some
$c,n_{0},t_{0},\epsilon_{0}>0$ and all $|T|\geq t_{0}$,
$0<\epsilon<\epsilon_{0}$, and integers $N\geq n_{0}$. We write $f=\Omega(g)$
if $g=\mathcal{O}(f)$ and we use $\widetilde{\mathcal{O}}$ to suppress
logarithmic factors in the asymptotic expression, in stark contrast to the
$\widetilde{\mathcal{O}}(T)$ cost that applies when local Hilbert spaces are
finite dimensional [10, 11, 48]. The cost can increase quadratically with $T$
in cases where local quantum numbers rise without bound as $T$ increases.

Although our main focus here is on the cost of dynamical simulation, our
bounds on truncation error also have consequences for approximating
eigenstates of the ideal untruncated Hamiltonian within the truncated Hilbert
space. For energy eigenvalues separated from the rest of the spectrum by a
specified gap, we derive a “tail bound” showing that the corresponding
eigenstates have very little support on large values of the local quantum
numbers. It follows that, for the class of models we study, a truncation error
less than $\epsilon$ can be achieved with truncation threshold
$\Lambda=\polylog(\epsilon^{-1})$, in contrast with the more naive estimate
$\Lambda=\poly(\epsilon^{-1})$ obtained using energy-based methods.

In our analysis of the cost of simulating time evolution, we assume that in
the initial state all local quantum numbers lie within a bounded range, and
then derive bounds on how much the local quantum numbers can increase during
time evolution. Our focus is somewhat related to previous work using
conservation of energy or particle number to tighten the analysis of Trotter
product formulas [63, 71], but our techniques differ from previous works in
that we need to deal with non-conserved quantities and unbounded local terms,
the latter of which makes the main tools in [63], namely Lemmas 1 and 2, no
longer apply. Our bounds also have potential applications to error mitigation
in quantum simulations, as an unexpectedly large value of a local quantum
number might flag an error that occurred during execution of the simulation
algorithm. Similar proposals have been based on conserved quantities [54, 12,
65, 35], and here we note that the same idea can be applied to non-conserved
quantities if we can rigorously bound the growth of those quantities during a
specified time interval.

Quantum simulations of non-abelian LGTs should eventually enable us to probe
particle physics in regimes where classical simulations are intractable.
Therefore the computational cost of such simulations is of fundamental
interest. Though for the sake of concreteness we focus on SU(2) in this work,
we anticipate that similar conclusions apply for other non-abelian gauge
groups, including SU(3), the relevant case for quantum chromodynamics. We
emphasize, though, that our results apply to LGTs where the lattice spacing is
a fixed physical length; we have not studied the approach to the continuum
limit or other formulations of quantum field theories without using lattices
[46]. We also emphasize that our analysis of the cost of simulating dynamics
assumes that the initial state is well approximated by a state in which all
local quantum numbers are less than the truncation threshold $\Lambda$; for
appropriate initial states, for example when the initial state is a
superposition of low-energy eigenstates, this assumption might be justified by
our tail bounds. However, we do not consider the computational cost of the
initial state preparation [55]. Despite these important caveats, our findings
strengthen the expectation that quantum computers will become powerful
instruments for scientific discovery.

(a)  (b)

Figure 1: Truncation threshold $\Lambda$ that achieves precision $\epsilon$
for the Hubbard-Holstein model with $N$ lattice sites and evolution time $T$.
(a) $\Lambda$ as a function of $T$ for $N=100$ and three different values of
$\epsilon$. (b) $\Lambda$ as a function of $T$ for $\epsilon=.01$ and three
different values of $N$. The Hamiltonian is given in (1). Model parameters are
from [43], with $\omega_{0}=1$ and $g=0.5$. The horizonal lines are the time-
independent truncation thresholds obtained using an energy-based method as in
[36]; the other curves are values of $\Lambda$ obtained in this work. See
Section K.2 in the Appendix for details.

##  2 Framework

We begin by setting up our framework and concisely stating our results, to be
proven in subsequent sections. For a more formal introduction of the framework
see Sections A and B. To illustrate our framework in a concrete setting, we
first consider the Hubbard-Holstein model [33], a model of electron-phonon
interactions. The model is defined on a $D$-dimensional lattice with linear
size $L$ and $L^{D}=N$ sites. Each site in the lattice, indexed by $x$,
contains two fermionic modes (spin up and down) and a bosonic mode. The
Hamiltonian is

$H=H_{f}+H_{fb}+H_{b},$ (1)

where $H_{f}$ is the Hamiltonian of the Fermi-Hubbard model [34] acting on
only the fermionic modes, and

$H_{fb}=g\sum_{x=1}^{N}(b_{x}^{\dagger}+b_{x})(n_{x,\uparrow}+n_{x,\downarrow}-1),\quad
H_{b}=\omega_{0}\sum_{x=1}^{N}b_{x}^{\dagger}b_{x},$ (2)

are the boson-fermion coupling and purely bosonic parts of the Hamiltonian
respectively. Here, $b_{x}$ is the bosonic annihilation operator on site $x$,
and $n_{x,\sigma}$ is the fermionic number operator for site $x$ and spin
$\sigma$.

In this setting, the local Hilbert space of each bosonic mode is infinite-
dimensional. In order to have a finite-dimensional local Hilbert space, a
natural idea is to impose an upper limit $\Lambda$ on the occupation number
$b_{x}^{\dagger}b_{x}$ (“number of particles”) in each bosonic mode. Then each
bosonic local Hilbert space has dimension $\Lambda+1$, and is spanned by the
particle-number eigenstates
$\\{|\lambda\rangle:\lambda=0,1,2,\dots,\Lambda\\}$. This imposed upper limit
results in a truncation. For a fixed site $x$, we define the projection
operator
$\Pi^{(x)}_{[0,\Lambda]}=\sum_{\lambda=0}^{\Lambda}\ket{\lambda}\bra{\lambda}_{x}\otimes
I$, where $I$ is the identity operator acting on the rest of the system.
Imposing the upper limit truncates a quantum state $\ket{\phi}$ to be
$\Pi^{(x)}_{[0,\Lambda]}\ket{\phi}$, and truncates the Hamiltonian $H$ to be
$\Pi^{(x)}_{[0,\Lambda]}H\Pi^{(x)}_{[0,\Lambda]}$. One may ask how large
$\Lambda$ should be for the resulting truncation error to be smaller than
$\epsilon$. There is some ambiguity regarding what “truncation error” means,
and we will refine this question later.

A similar situation is encountered in LGTs, where we consider the Hamiltonian
formulation proposed in [44]. For a more detailed introduction to the LGTs we
consider, see Section A. We have a $D$-dimensional lattice consisting of $N$
total sites and $\mathcal{O}(N)$ gauge links. Each gauge link has an infinite-
dimensional local Hilbert space, and the Hamiltonian contains unbounded
operators associated with each link. We no longer have a natural notion of
particle number, but the truncation of the link Hilbert space can still be
performed according to what we call the local quantum number. We focus on two
cases: the U(1) and SU(2) LGTs. For the U(1) case, we choose the local quantum
number to be the integer-valued electric field. We retain only the part of
local Hilbert space with electric field value in the interval
$[-\Lambda,\Lambda]$; hence the truncated local Hilbert space at each link is
$2\Lambda+1$ dimensional. More precisely, for a fixed gauge link $\nu$, we
define the projection operator
$\Pi^{(\nu)}_{[0,\Lambda]}=\sum_{\lambda=-\Lambda}^{\Lambda}\ket{\lambda}\bra{\lambda}_{\nu}\otimes
I$, where $\ket{\lambda}$’s are elements of the electric eigenbasis. Similar
to the previous example of the Hubbard-Holstein model, after truncation a
quantum state $\ket{\phi}$ becomes
$\Pi^{(\nu)}_{[-\Lambda,\Lambda]}\ket{\phi}$, and the Hamiltonian $H$ becomes
$\Pi^{(\nu)}_{[-\Lambda,\Lambda]}H\Pi^{(\nu)}_{[-\Lambda,\Lambda]}$. For the
SU(2) case, we choose the local quantum number to be 2 times the total angular
momentum (the multiplication by 2 makes the local quantum number an integer).
If we retain only the part of the link Hilbert space with total angular
momentum no larger than $\Lambda/2$, then the link Hilbert space has dimension
$(\Lambda+1)(\Lambda+2)(2\Lambda+3)/6$, and is spanned by the angular momentum
eigenstates $\ket{jmm^{\prime}}$ where $j$ is a half integer less than or
equal to $\Lambda/2$ and $-j\leq m,m^{\prime}\leq j$. Again one may ask how
large $\Lambda$ should be for the resulting truncation error to be small.

When analyzing time evolution, this question can be refined into two different
but related questions.

Question 1 (Truncating an evolved quantum state): Consider an initial state
such that at some particular site or link the local quantum number is no
larger than $\Lambda_{0}$. After the state evolves forward for time $T$, how
should the truncation threshold $\Lambda$ be chosen for that site or link so
that the resulting error is at most $\epsilon$? We show that

$\Lambda^{1-r}=\Lambda_{0}^{1-r}+\widetilde{\mathcal{O}}((\chi
T+1)\polylog(\epsilon^{-1}))$ (3)

suffices, where $r=1/2$ for bosons and $r=0$ for LGTs, and $\chi$ is a
constant that only depends on the model parameters but not on the system size
or on $T$. If we want to truncate every bosonic mode or gauge link in the
model, then, to account for the accumulation of error, $\epsilon^{-1}$ in (3)
is replaced by $N\epsilon^{-1}$, where $N$ is the system size.

A simple example shows that this scaling of the truncation threshold is
optimal in certain cases, such as quadratic scaling in time for bosons
$(r=1/2)$. Suppose $H=b+b^{\dagger}$ where $b$ is the annihilation operator of
a bosonic mode. Then $e^{-iTH}\ket{0}$, where $\ket{0}$ is the vacuum state,
is a coherent state such that the particle number distribution is Poissonian
with mean $T^{2}$. Because the Poisson distribution concentrates around its
mean, a truncation threshold that achieves constant precision must scale like
$\Omega(T^{2})$, which matches (3) for $r=1/2$.

It is instructive to compare our approach with the method based on energy
conservation described in [36, 37]. That method yields a truncation threshold
for a single site with a polynomial dependence on the inverse accuracy
$\epsilon^{-1}$. To truncate a system of $\widetilde{\mathcal{O}}(N)$ sites,
we scale down $\epsilon$ by a factor of $N$, resulting in a threshold
$\Lambda$ scaling polynomially with $N\epsilon^{-1}$. In contrast, our bound
has only polylogarithmic dependence on $N\epsilon^{-1}$, an exponential
improvement compared to the truncation threshold obtained using the energy-
based method. Importantly, this advantage holds not only in the asymptotic
regime, but also when the constant prefactors are incorporated. We numerically
compare our bound with the energy-based bound for the Hubbard-Holstein model,
observing a significantly better estimate in various parameter regimes. We
illustrate this comparison in Figure 1 and discuss it in more detail in
Section K in the Appendix.

Question 2 (Truncating the Hamiltonian): Consider an initial state such that
the local quantum number is no larger than $\Lambda_{0}$ at all sites or at
all links, and suppose the state evolves forward for time $T$ using a
truncated Hamiltonian $\widetilde{H}$ rather than the ideal untruncated
Hamiltonian $H$ (we will define $\widetilde{H}$ later). How should the
truncation threshold $\Lambda$ be chosen so that the truncated evolved state
matches the ideal evolved state up to time $T$ with an error at most
$\epsilon$? We show that

$\Lambda^{1-r}=\Lambda_{0}^{1-r}+\widetilde{\mathcal{O}}((\chi
T+1)\polylog(N\epsilon^{-1}))$ (4)

suffices, where $r=1/2$ for bosons and $r=0$ for LGTs, $N$ is the system size,
and $\chi$ is again a constant that does not depend on $N$ or on $T$.

Our above two questions both concern the truncation of the local quantum
number, albeit from different perspectives: the first focuses on evolved
quantum states while the second focuses on Hamiltonians. In fact, a threshold
for truncating the Hamiltonian can be directly used to truncate evolved
quantum states, although some extra efforts are required to handle the
converse. The truncation of an evolved quantum state in Question 1 is only for
some fixed time $T$, but when we perform truncation on the Hamiltonian in
Question 2, the evolved quantum state will never have a local quantum number
beyond $[-\Lambda,\Lambda]$ throughout the evolution up to time $T$. In this
sense, our second result is stronger than the first one.

It is worth noting that while the energy-based method in Refs.[36] is enough
to establish a bound to address Question 1, it cannot be used to address
Question 2. This is because the state truncation error does not decay fast
enough as we increase the truncation threshold, and as a result the energy-
based bound is not enough for the derivation in Section D.

Before stating our results for a more general class of Hamiltonians, we first
introduce some notation. For a bosonic mode or gauge link which we denote by
$\nu$, we denote by $\Pi^{(\nu)}_{S}$ the projection operator imposing the
condition that the local quantum number takes values from the set $S$. We also
denote $\Pi^{\mathrm{all}}_{S}=\prod_{\nu}\Pi^{(\nu)}_{S}$; this is the
projection operator imposing the condition on all bosonic modes or gauge
links. For any projection operator $\Pi$, we write its complement as
$\overline{\Pi}=I-\Pi$. The truncated Hamiltonian mentioned in Question 2 is
$\widetilde{H}=\Pi^{\mathrm{all}}_{S}H\Pi^{\mathrm{all}}_{S}$, where $H$ is
the untruncated Hamiltonian, and $S$ is the set of local quantum numbers less
than or equal to the truncation threshold.

Using this notation we can readily pinpoint the common structure of the
Hamiltonians in the Hubbard-Holstein model and the U(1) and SU(2) LGTs. In all
three examples, although the Hamiltonian contains local terms with unbounded
norm, each of these terms changes the local quantum number at only a single
site or a single link; there are no unbounded terms that allow the local
quantum number to propagate from site to site or from link to link. For each
site or link, denoted by $\nu$, we may write the full Hamiltonian $H$ of the
model as

$H=H_{W}^{(\nu)}+H_{R}^{(\nu)},$ (5)

where $H_{W}^{(\nu)}$ is the part of the Hamiltonian that can change the value
of the local quantum number at $\nu$, and $H_{R}^{(\nu)}$ contains all the
terms in the Hamiltonian that preserve the value of the local quantum number
at $\nu$. These two parts satisfy the conditions

$\displaystyle\Pi^{(\nu)}_{\lambda}H_{W}^{(\nu)}\Pi^{(\nu)}_{\lambda^{\prime}}=0,\
\text{if }|\lambda-\lambda^{\prime}|>1,$ (6a)
$\displaystyle\|H_{W}^{(\nu)}\Pi^{(\nu)}_{[-\Lambda,\Lambda]}\|\leq\chi(\Lambda+1)^{r},$
(6b) $\displaystyle[H_{R}^{(\nu)},\Pi^{(\nu)}_{\lambda}]=0.$ (6c)

Here $\Pi^{(\nu)}_{\lambda}$ projects onto the eigenspace with local quantum
number $\lambda$, $\chi$ and $0\leq r<1$ are parameters that depend on the
model, and $\|\cdot\|$ is the spectral norm. (The notation
$\Pi^{(\nu)}_{[-\Lambda,\Lambda]}$ is appropriate for the U(1) gauge theory,
where the electric field can take either positive or negative integer values,
but we will use this same notation for the other models as well, even though
in those models the local quantum number takes only nonnegative values.) These
three conditions can be interpreted as follows: the first condition requires
$H_{W}^{(\nu)}$ to change the local quantum number by at most $\pm 1$. The
second condition requires that the rate at which the maximal local quantum
number $\Lambda$ changes is sublinear in $\Lambda$. The third condition
requires $H_{R}^{(\nu)}$ to preserve the local quantum number. See Section B
in the Appendix for a more detailed explanation of this framework.

Let us verify that the Hubbard-Holstein Hamiltonian in (1) fits the general
framework of (5) and (6). The bosonic mode appears only in on-site terms.
Choosing
$H_{W}^{(x)}=g(b_{x}^{\dagger}+b_{x})(n_{x,\uparrow}+n_{x,\downarrow}-1)$ and
$H_{R}^{(x)}=\sum_{x^{\prime}}\omega_{0}b_{x^{\prime}}^{\dagger}b_{x^{\prime}}+\sum_{x^{\prime}\neq
x}g(b_{x^{\prime}}^{\dagger}+b_{x^{\prime}})(n_{x^{\prime},\uparrow}+n_{x^{\prime},\downarrow}-1)+H_{f}$,
we see that $H_{W}^{(x)}$ changes the local bosonic particle number by at most
$\pm 1$, and that $H_{R}^{(x)}$ preserves the local bosonic particle number.
Moreover, using $\Pi^{(x)}_{\lambda}$ to denote the projector onto the
subspace with $\lambda$ bosonic particles on site $x$, we see that
$\|(b_{x}^{\dagger}+b_{x})\Pi^{(x)}_{[0,\Lambda]}\|\leq 2\sqrt{\Lambda+1}$,
which implies (6b) is satisfied with $\chi=2g$ and $r=1/2$. In Sections A and
J in the Appendix, we explain how other examples fit this framework, including
U(1) and SU(2) LGTs, the spin-fermion coupling in the Fröhlich model [23], and
spin-boson coupling in the Dicke model. In Section 3 we show that for
Hamiltonians with the structure indicated in (5), (6), local quantum numbers
may be truncated as specified by the answer (3) to Question 1 and the answer
(4) to Question 2.

The linear dependence on the evolution time $T$ in (3), (4) has a simple
interpretation. Specifically, for the case of a bosonic mode ($r=1/2$) where
$H_{W}^{(\nu)}$ is linear in creation and annihilation operators, the
conditions (6a), (6b) impose that in time $T$ the position of the mode in
phase space is translated by $\mathcal{O}(T)$. Since the particle number
scales like the square of the displacement from the origin of phase space, a
truncation threshold growing quadratically with $T$, as specified in (3), (4),
suffices to approximate the translated state accurately.

Given the scaling of the truncation threshold expressed in (4), we can
accurately approximate time evolution using the truncated Hamiltonian
$\widetilde{H}$, in which all local terms in the Hamiltonian have bounded
norm. In Section 4, we leverage this observation to analyze the cost of
simulating time evolution on a digital quantum computer for the Hamiltonians
characterized above. In particular, we develop algorithms for simulating the
U(1) and SU(2) LGTs that achieve an almost linear dependence on the spacetime
volume, a substantial improvement over previous estimates of the gate
complexity [69, 38, 61]. We also analyze the cost of simulating the Hubbard-
Holstein model in Section 4. In Section 5, by applying these results on time
evolution, we establish that spectrally isolated energy eigenstates can be
approximated using a local quantum number truncation threshold scaling
polylogrithmically with the allowed error.

##  3 Hilbert space truncation in time evolution

We now show how the truncation threshold scaling relations (3) and (4) are
obtained. Recall that in our two questions about the truncation threshold,
Question 1 concerns truncating the quantum state obtained from exact time
evolution for time $T$. Using the notations introduced earlier, we can clarify
this question and our result. We define a quantity
$\|\overline{\Pi}^{(\nu)}_{[-\Lambda,\Lambda]}e^{-iTH}\Pi^{(\nu)}_{[-\Lambda_{0},\Lambda_{0}]}\|$
which we call leakage. If we start with an initial state $\ket{\psi_{0}}$ with
the local quantum number on $\nu$ taking a value in
$[-\Lambda_{0},\Lambda_{0}]$, and denote the state at time $t$ by
$\ket{\psi(t)}$, then the truncation error
$\|\overline{\Pi}^{(\nu)}_{[-\Lambda,\Lambda]}\ket{\psi(T)}\|$ is upper
bounded by the leakage. Therefore, to ensure that the truncation error is at
most $\epsilon$, we only need to keep the leakage below $\epsilon$.

As mentioned before, we assume $H$ has the structure (5) with $H^{(\nu)}_{W}$
and $H^{(\nu)}_{R}$ satisfying (6). First we prove a leakage bound that holds
for relatively short evolution time governed by such $H$, and then establish
(3) by extending the short-time leakage bound to longer times. We view the
time evolution in the interaction picture, and consider the evolution of
$\ket{\psi_{I}(t)}=e^{itH^{(\nu)}_{R}}\ket{\psi(t)}$. Because we assume
$H^{(\nu)}_{R}$ preserves the local quantum number, $\ket{\psi_{I}(t)}$ and
$\ket{\psi(t)}$ induce the same local quantum number distribution. In the
interaction picture, $\ket{\psi_{I}(t)}$ evolves with a time-dependent
Hamiltonian
$H^{(\nu)}_{W}(t)=e^{itH^{(\nu)}_{R}}H^{(\nu)}_{W}e^{-itH^{(\nu)}_{R}}$. We
then apply the Dyson series expansion to the unitary operator generated by
$H^{(\nu)}_{W}(t)$. In the proof of Lemma 1 in the Appendix, we show that if
$0\leq T\leq 1/(2\chi(\Lambda_{0}+1)^{r})$, the truncated Dyson series with
$\Delta$ terms approximates the exact evolution up to an error
$e^{-\Omega(\Delta)}$. Moreover, such a truncated Dyson series can change the
local quantum number by at most $\pm(\Delta-1)$ due to (6a). Therefore we have
the short-time leakage bound

$\|\overline{\Pi}^{(\nu)}_{(-\Lambda_{0}-\Delta,\Lambda_{0}+\Delta)}e^{-iTH}\Pi^{(\nu)}_{[-\Lambda_{0},\Lambda_{0}]}\|\leq
e^{-\Omega(\Delta)}.$ (7)

Using this short-time leakage bound, we can derive the long-time bound in (3).
Specifically, for any choice of
$\Lambda_{0}<\Lambda_{1}<\cdots<\Lambda_{J}=\Lambda$,
$0=T_{0}<T_{1}<\cdots<T_{J}=T$, the total leakage is at most the sum of $J$
short-time leakages (see Lemma 2 in the Appendix).

$\displaystyle\|\overline{\Pi}^{(\nu)}_{[-\Lambda,\Lambda]}e^{-iTH}\Pi^{(\nu)}_{[-\Lambda_{0},\Lambda_{0}]}\|$
(8)
$\displaystyle\leq\sum_{j=1}^{J}\|\overline{\Pi}^{(\nu)}_{[-\Lambda_{j},\Lambda_{j}]}e^{-i(T_{j}-T_{j-1})H}\Pi^{(\nu)}_{[-\Lambda_{j-1},\Lambda_{j-1}]}\|.$

We then carefully choose $T_{j}$’s and apply the short-time leakage bound to
each segment $[T_{j-1},T_{j}]$, which gives an upper bound on the right-hand
side of (8). Since the local quantum number can potentially change as the
system evolves, we define the length of time steps adaptively based on the
instantaneous quantum number to reach the same target accuracy. Specifically,
$T_{j}$ and $\Lambda_{j}$ are chosen to satisfy $0\leq T_{j}-T_{j-1}\leq
1/(2\chi(\Lambda_{j-1}+1)^{r})$. This establishes the scaling in (3) and
provides an answer to Question 1. We summarize our result below and leave
details of the proof to Section C in the Appendix.

Theorem (State truncation (Theorem 5 in the Appendix)). Let $H$ be a
Hamiltonian such that $H=H_{W}^{(\nu)}+H_{R}^{(\nu)}$ satisfies (6) with
parameters $\chi$ and $r$ for a fixed mode or link $\nu$. For any $t\geq 0$
and integers $\Lambda\geq\Lambda_{0}\geq 0$,

$\displaystyle\left\|\overline{\Pi}_{[-\Lambda,\Lambda]}^{(\nu)}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{(\nu)}\right\|$
(9) $\displaystyle\leq\poly(\chi
t,\Lambda_{0},\Lambda)\exp\left(-\Omega\left(\frac{\Lambda^{1-r}-\Lambda_{0}^{1-r}}{\chi
t+1}\right)\right).$

We now set out to answer Question 2. First we clarify the question using our
notation for projection operators. Here we consider replacing $H$ by a
truncated Hamiltonian
$\widetilde{H}=\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}$,
where

$\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}=\prod_{\nu}\Pi^{(\nu)}_{[-\Lambda,\Lambda]}$

applies truncation on all sites or links. the truncation threshold $\Lambda$
is chosen large enough so that evolution governed by $\widetilde{H}$ is a good
approximation to the exact evolution. The approximation error is upper bounded
by

$\max_{0\leq t\leq
T}\|(e^{-itH}-e^{-it\widetilde{H}})\Pi^{\mathrm{all}}_{[-\Lambda_{0},\Lambda_{0}]}\|.$
(10)

Therefore our goal is to choose $\Lambda$ to ensure that this error is at most
$\epsilon$. This is accomplished by the following theorem which we establish
in Section D in the Appendix and preview here.

Theorem (Hamiltonian truncation (Theorem 6 in the Appendix)). Let $H$ be a
Hamiltonian with $\mathcal{O}(N)$ bosonic modes or gauge links, such that
$H=H_{W}^{(\nu)}+H_{R}^{(\nu)}$ satisfies (6) with parameters $\chi$ and $r$
for every mode or link $\nu$. Assume that all projection operators
$\Pi^{(\nu)}_{\lambda}$ commute with each other. For any integers
$\Lambda\geq\Lambda_{0}\geq 0$, define
$\widetilde{H}=\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}$
and assume $\|[\widetilde{H},H]\|=\poly(N,\Lambda)$. Then for any $t>0$,

$\displaystyle\left\|(e^{-it\widetilde{H}}-e^{-itH})\Pi^{\mathrm{all}}_{[-\Lambda_{0},\Lambda_{0}]}\right\|$
(11) $\displaystyle\leq\poly(\chi
t,\Lambda_{0},\Lambda,N)\exp\left(-\Omega\left(\frac{\Lambda^{1-r}-\Lambda_{0}^{1-r}}{\chi
t+1}\right)\right).$

We now briefly explain how we upper bound the Hamiltonian truncation error
$\|(e^{-itH}-e^{-it\widetilde{H}})\Pi^{\mathrm{all}}_{[-\Lambda_{0},\Lambda_{0}]}\|$.
Since this is equal to
$\|(I-e^{itH}e^{-it\widetilde{H}})\Pi^{\mathrm{all}}_{[-\Lambda_{0},\Lambda_{0}]}\|$,
we instead show that
$e^{itH}e^{-it\widetilde{H}}\Pi^{\mathrm{all}}_{[-\Lambda_{0},\Lambda_{0}]}$
is close to $\Pi^{\mathrm{all}}_{[-\Lambda_{0},\Lambda_{0}]}$ for sufficiently
large $\Lambda$. We prove this by expanding the target quantity using the
formula for Trotter error [72, Eq. (3.4)]. We then invoke the leakage bound
Theorem 5 in the Appendix, as well as the fact that the truncated Hamiltonian
and the original Hamiltonian act identically on a state with local quantum
numbers all in $[-\Lambda+1,\Lambda-1]$, i.e.
$(H-\widetilde{H})\Pi^{\mathrm{all}}_{[-\Lambda+1,\Lambda-1]}=0$. The
resulting bound depends on the commutator norm
$\mathcal{A}(\Lambda)=\|[H,\widetilde{H}]\|$. If this quantity scales
polynomially with the system size $N$ and the truncation threshold $\Lambda$,
we can then establish the desired scaling in (4) and answer our Question 2.
This is indeed true for the Hubbard-Holstein model and the LGTs, which can be
seen as follows. We have
$[H,\widetilde{H}]=(H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]})^{2}-(\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}H)^{2},$
where each local term in $H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}$ can be
bounded by a polynomial of $\Lambda$ (linear for Hubbard-Holstein model and
quadratic for LGTs) and there are $\mathcal{O}(N)$ of such local terms. Thus
we have $\|[H,\widetilde{H}]\|=\mathcal{O}(N^{2}\poly(\Lambda))$ in all three
examples. We discuss the Hamiltonian truncation in more detail in Section D in
the Appendix.

##  4 Application to Hamiltonian simulation

Our main results on the truncation of unbounded Hamiltonians allow us to
simulate such systems more efficiently with a provable accuracy guarantee. For
concretenesss, we consider the problem of digital Hamiltonian simulation,
wherein the dynamics of a quantum system are approximated on a quantum
computer by elementary gates, and the cost of simulation is determined by the
gate complexity. While the majority of the past work on Hamiltonian simulation
has focused on quantum systems with finite-dimensional local Hilbert spaces,
there are also systems of physical interest whose local Hilbert spaces are
infinite dimensional. In such cases, it is typically necessary to perform
truncation, so that quantum states can be represented and processed on a
digital quantum computer.

In Section 3 we established that time evolution governed by the Hamiltonian
$H$ can be accurately approximated by time evolution governed by the truncated
Hamiltonian $\widetilde{H}$, if we choose the truncation threshold according
to (4). Here we propose Hamiltonian simulation approaches that take advantage
of our results and discuss the implications for expected costs on a quantum
computer. Unlike $H$, $\widetilde{H}$ acts nontrivally on only a finite-
dimensional subspace of the infinite-dimensional Hilbert space at each site or
link. Furthermore, because each unbounded local term in $H$ is truncated
separately, the truncation does not affect the geometric locality of the
Hamiltonian. Therefore, to simulate $e^{-itH}$ we can instead simulate
$e^{-it\widetilde{H}}$, which can be done on a quantum computer using existing
simulation techniques for local Hamiltonians [47, 9, 11, 49, 48, 50, 28]. In
what follows, we consider simulations of the U(1), SU(2) LGTs, as well as the
Hubbard-Holstein model, although the quantum algorithms we present can in
principle be extended to simulate other gauge theories and bosonic systems
within our framework.

Simulating lattice gauge theories with near-linear spacetime volume scaling.
We propose an algorithm to simulate the time evolution of the U(1) and SU(2)
LGTs in $D$ spatial dimensions; Hamiltonians of these models are described in
Eq. (17) in the Appendix. The goal is to simulate a lattice with $N$ sites for
time $T$ with total error at most $\epsilon$. Our algorithm combines the Haah-
Hastings-Kothari-Low (HHKL) decomposition [28], which provides a nearly
optimal approach for geometrically local Hamiltonians, with the interaction-
picture simulation method [50], which gives further improved scaling with the
truncation threshold. We show that the simulation can be done with gate
complexity $\widetilde{\mathcal{O}}(NT\polylog(\Lambda_{0}\epsilon^{-1}))$,
assuming that in the initial state the local quantum number (electric field
value for U(1) or total angular momentum for SU(2)) on each gauge link is in
the interval $[-\Lambda_{0},\Lambda_{0}]$. Thus we achieve an almost linear
dependence of the gate complexity on the spacetime volume $NT$. We briefly
outline the algorithm here; further details are presented in Section E.1 in
the Appendix.

We first use [28, Lemma 6] to decompose the time evolution of the entire
system due to $\widetilde{H}$ into time evolution of blocks. Each block,
denoted by $\mathcal{B}$, has size
$\ell^{D}=\mathcal{O}(\polylog(NT\epsilon^{-1}))$ and we only need to
implement its evolution for time $\tau=\mathcal{O}(1)$. There are
$\mathcal{O}(N)$ such blocks and the entire time evolution is divided into
$\mathcal{O}(T)$ segments. We note that [28, Lemma 6] requires a constant
Lieb-Robinson velocity, which was guaranteed by [28, Lemma 5] since all terms
in their Hamiltonian were geometrically local with norm upper bounded by a
constant. In our case, however, there are terms in the truncated Hamiltonian
with norm $\textrm{poly}(\Lambda)$. Fortunately, these terms with
$\Lambda$-dependent norm act on either a single lattice site (in the models
with bosonic modes) or a single gauge link (in LGTs). We show in Section I in
the Appendix that for Hamiltonians of this form, the Lieb-Robinson velocity is
indeed bounded above by a $\Lambda$-independent constant as [28] requires.

When simulating each block $\mathcal{B}$ we use the interaction picture
Hamiltonian simulation technique suggested in [70] and developed in [50], and
the gate complexity for simulation up to time $\tau=\mathcal{O}(1)$ is
$\mathcal{O}(\polylog(\Lambda NT\epsilon^{-1}))$. For $\Lambda$ we use the
scaling (4). There are in total $\mathcal{O}(NT)$ such simulations that need
to be performed, leading to a total gate complexity of
$\widetilde{\mathcal{O}}(NT\polylog(\Lambda_{0}\epsilon^{-1}))$. The
interaction picture is useful because it allows us to express the time
evolution operator as a product of two operators. One factor in this product
is the evolution arising from the terms in the truncated Hamiltonian
$\widetilde{H}$ which have $\Lambda$-dependent norms, the terms involving the
electric field at each link. This evolution can be “fast-forwarded” [3, 26]
because the Hamiltonian is diagonal in a natural basis, and the evolution
operator is just the tensor product of simple unitary operators, each acting
on a single link. The other factor in the product is the interaction-picture
evolution operator generated by the time-dependent interaction-picture
Hamiltonian, in which each term has $\Lambda$-independent norm because the
evolution induced by the electric field has been “rotated away.” As a result,
the cost of simulating the evolution of a block $\mathcal{B}$ is
polylogarithmic in $\Lambda$, and the cost of simulating evolution of $N$
sites for time $T$ is nearly linear in the spacetime volume $NT$.

Previous work on the quantum simulation of LGTs such as [69, 38, 61] does not
explain how to choose the truncation threshold $\Lambda$ to perform simulation
with a provable accuracy. While this issue can be remedied by using our
Hamiltonian truncation threshold (4), our result still substantially improves
over the previous results $\widetilde{\mathcal{O}}(N^{3/2}T^{5/2})$ from [69,
38] and $\widetilde{\mathcal{O}}(N^{2}T^{2})$ from [61].

Simulating bosonic systems and an $\widetilde{\Omega}(T^{2})$ gate complexity
lower bound. Here we outline two methods for simulating bosonic systems, using
the Hubbard-Holstein model as an example model. In the first method we again
use the HHKL decomposition combined with the interaction-picture Hamiltonian
simulation; see Section E in the Appendix for a detailed discussion. The
important difference from the setting of LGTs is that, when simulating a block
$\mathcal{B}$ of the Hubbard-Holstein model, we cannot get a polylogarithmic
dependence on $\Lambda$. Rather, the gate complexity to simulate a block is
$\mathcal{O}(\sqrt{\Lambda}\polylog(\Lambda NT\epsilon^{-1}))$, because, as
explained in Section E.2 in the Appendix, the Hubbard-Holstein Hamiltonian has
multiple unbounded terms and it is not known how to fast-forward them
simultaneously. Since there are $\mathcal{O}(NT)$ blocks to be simulated, and
the scaling of $\Lambda$ is given by (4), the total gate complexity is
$\widetilde{\mathcal{O}}(NT(\sqrt{\Lambda_{0}}+T)\polylog(\epsilon^{-1}))$.

In the second method we use the $p$-th order Trotter product formula, which
can be easier to implement in practice. To obtain a tight error bound in this
case one may use the commutation relations among the Hamiltonian terms [19,
17]. For the Hubbard-Holstein model we use the canonical commutation relation
between the bosonic position and momentum operators
$[X_{\alpha},P_{\alpha}]=i$, and also invoke geometric locality to tightly
bound the error. A subtle issue with this naive analysis is that the canonical
commutation relation no longer holds when acting on arbitrary states due to
the truncation of the Hamiltonian terms. However, we recover the commutation
relation by restricting to states with low particle numbers. A detailed
discussion of all the issues involved can be found in Section F in the
Appendix. In the end we obtain a gate complexity
$\widetilde{\mathcal{O}}\left(N^{1+1/p}(\sqrt{\Lambda_{0}}+T)^{1+2/p}T^{1+1/p}\epsilon^{-1/p}\right)$,
nearly matching the complexity of the method based on HHKL decomposition for
large values of $p$.

Notice that the gate complexity of simulating the Hubbard-Holstein model has
an almost quadratic dependence on the time $T$, in stark contrast with the
almost linear dependence that applies when all local terms in the Hamiltonian
have bounded norm [28, 17]. In fact, there exist unbounded Hamiltonians which
are impossible to simulate with an almost linear scaling in $T$. In Section G
in the Appendix we construct a class of Hamiltonians acting on one bosonic
mode and $N$ qubits for which simulating the evolution of qubits for time $T$
requires $\widetilde{\Omega}(NT^{2})$ gates in general, for $\sqrt{N}\leq
T\leq 2^{N/2}$.

##  5 The eigenstate tail bound

Aside from studies of dynamics, classical or quantum computers may be used to
study the static properties of ground states or low-energy states in quantum
systems involving bosons or gauge fields. As in simulations of dynamics, we
must truncate the local quantum numbers to ensure that local Hilbert spaces at
sites or links are finite dimensional. How well can we approximate energy
eigenstates of the ideal untruncated Hamiltonian within the truncated Hilbert
space?

Suppose that for each site or link, denoted by $\nu$, the Hamiltonian $H$ can
be expressed as in (5) and satisfies (6). Consider a nondegenerate eigenvalue
$\varepsilon$ of $H$, with corresponding eigenstate $\ket{\Psi}$, where
$\varepsilon$ is separated from the rest of the spectrum of $H$ by a gap
$\delta$, and suppose that the expectation value of the absolute value of the
local quantum number in the state $\ket{\Psi}$ is finite,
$\bar{\lambda}=\sum_{\lambda}\left|\lambda\right|\braket{\Psi}{\Pi^{(\nu)}_{\lambda}}{\Psi}<\infty$.
Our goal is to find a truncation threshold $\Lambda$ such that
$\|\overline{\Pi}^{(\nu)}_{[-\Lambda,\Lambda]}\ket{\Psi}\|\leq\epsilon$. We
show that this truncation threshold can be chosen to scale with $\epsilon$,
$\delta$, and $\bar{\lambda}$ according to

$\Lambda^{1-r}=(2\bar{\lambda})^{1-r}+\mathcal{O}(\chi\delta^{-1}\log^{2}(\epsilon^{-1})+\log(\epsilon^{-1})),$
(12)

where $\chi$ is a constant independent of system size.

A detailed proof of (12) can be found in Section H in the Appendix. The
polylogarithmic dependence of the truncation threshold $\Lambda$ on the
truncation error $\epsilon$ arises because the distribution of local quantum
numbers in the eigenstate $|\Psi\rangle$ decays exponentially. This contrasts
with the polynomial decay one can derive using Markov’s or Chebyshev’s
inequality.

The main tool used in our proof is an approximate eigenstate projection
operator [30]

$\widetilde{P}_{\varepsilon}=\frac{\sigma}{\sqrt{2\pi}}\int_{-T}^{T}\mathrm{d}te^{-\frac{1}{2}\sigma^{2}t^{2}}e^{-i\varepsilon
t}e^{itH}.$ (13)

When $\sigma\ll\delta$ and $T\gg\sigma^{-1}$, this operator is close to the
eigenstate projection operator $P_{\varepsilon}=\ket{\Psi}\bra{\Psi}$. We
derive (12) by applying the approximate projector
$\widetilde{P}_{\varepsilon}$ to a suitable initial state and using properties
of the time evolution operator $e^{-iHt}$, in particular the truncation
threshold result (3). We may choose the initial state to be
$\Pi^{(\nu)}_{[-2\bar{\lambda},2\bar{\lambda}]}\ket{\Psi}$, which by Markov’s
inequality has an $\mathcal{O}(1)$ overlap with $\ket{\Psi}$. By observing
that $e^{-iHt}\Pi^{(\nu)}_{[-2\bar{\lambda},2\bar{\lambda}]}\ket{\Psi}$ can be
well approximated by a state with an appropriately chosen truncation
threshold, we obtain (12).

Note that (12) does not apply to eigenstates that are degenerate due to
symmetries of the Hamiltonian $H$. Nor is it particularly useful when applied
to generic highly excited eigenstates, for which the gap $\delta$ may be
exponentially small in the system size.

##  6 Discussion

We have studied the task of simulating Hamiltonian dynamics for quantum
systems on a lattice, where local Hilbert spaces at lattice sites or links are
infinite dimensional. In these systems, local quantum numbers on sites or
links can be arbitrarily large in principle. For a large class of such models,
we derived upper bounds on how rapidly these local quantum numbers can
increase with time, hence showing that time evolved states can be well
approximated in a truncated Hilbert space in which each local quantum number
is no larger than a truncation threshold $\Lambda$. In particular, we showed
that for a fixed evolution time $T$, a precision $\epsilon$ can be achieved by
choosing $\Lambda$ scaling polylogarithmically with $\epsilon^{-1}$, as
indicated in (3) and (4). Leveraging this finding, we established a threshold
for truncating the Hamiltonian with a provable accuracy guarantee and
developed algorithms for quantum simulation of LGTs with gate complexity
$\widetilde{\mathcal{O}}(NT\polylog(\Lambda_{0}\epsilon^{-1}))$, where $N$ is
the system size, assuming that the initial state can be well approximated with
truncation threshold $\Lambda_{0}$. For a bosonic system like the Hubbard-
Holstein model, our algorithm has gate complexity
$\widetilde{\mathcal{O}}(NT(\sqrt{\Lambda_{0}}+T)\polylog(\epsilon^{-1}))$. By
applying our bounds on the growth of local quantum numbers, we also showed
that spectrally isolated energy eigenstates can be approximated with precision
$\epsilon$ using a truncation threshold polylogartihmic in $\epsilon^{-1}$, as
indicated in (12).

Although formally the local Hilbert spaces are infinite dimensional in the
models we considered, our results show that at least for some purposes these
models can be accurately approximated by models with finite-dimensional local
Hilbert spaces of relatively modest size. Many fundamental results have been
derived for quantum spin systems with finite-dimensional spins on each lattice
site, such as the exponential clustering theorem [29, 31, 57],

the area law in one dimension [30, 1], and the connection between local and
global eigenstates [2]. Perhaps the tools we have developed can be exploited
to extend some of these results to systems with infinite-dimensional local
degrees of freedom.

There are certain models of physical interest that do not immediately fit in
our framework. These include models that involve a quadratic coupling between
bosonic modes, such as the Bose-Hubbard model ($r=1$ in (6b)) and the
discretized $\phi^{4}$ theory ($r=2$); our analysis handles the case where
$r<1$ in (6). Our framework also does not apply to boson-fermion coupling
models where anharmonicity is involved that leads to $r=2$. Nevertheless, the
method we have developed already provides a unified treatment for a wide range
of bosonic systems and lattice gauge theories, and we hope future work could
study other physical systems that have not been considered in our work.

For $\phi^{4}$ theory on a lattice, truncation thresholds were previously
analyzed using energy conservation and Chebyshev’s inequality [36], a method
that can be extended to other models as well. Our results apply only to models
that satisfy (5) and (6). For models in this class, we compare our methods
with energy-based methods in Section K in the Appendix, finding that our
methods yield a more favorable truncation threshold in the limit of short
time, high precision, or large system size.

The energy-based truncation threshold in [36] has the advantage of being time
independent, and it can also be applied to models that do not satisfy (5) and
(6), such as $\phi^{4}$ theory and other models involving bosons with
anharmonic couplings. However, it has the disadvantage that the truncation
threshold scales polynomially rather than polylogarithmically with
$\epsilon^{-1}$. Under suitable conditions, can the truncation threshold scale
as $\polylog(\epsilon^{-1})$ in a broader class of models than those
satisfying (5) and (6), and are there models in which
$\polylog(\epsilon^{-1})$ scaling can be achieved by a time-independent
truncation threshold? Moreover, the energy-based truncation threshold provides
an answer to Question 1 in the context of truncating a quantum state, but it
has not been shown, at the same level of rigor, that the energy-based method
also provides an answer to Question 2 in the context of truncating the
Hamiltonian. The latter is however necessary if we want to rigorously apply
the energy-based truncation threshold to Hamiltonian simulation. These are
open questions to be addressed in future work.

Another question that has yet to be answered is how to control the error for
observables in boson and gauge theory simulations. For bounded observables,
once we can control the error in the quantum state, we can automatically
control the error of observables. However for unbounded observables, such as
the boson occupation number and the electric field value, this simple approach
is not suitable. For local observables in lattice models with a finite speed
limit for information propagation, one intuitively expects the observable
error to also have a bound that respects this locality. We hope our approach
can be extended to address questions of this kind, ultimately leading to a
theoretical foundation for studying quantum systems with infinite degrees of
freedom.

## Acknowledgement

We thank Kunal Sharma, Mark Wilde, Minh Cong Tran, Junyu Liu, Chi-Fang
(Anthony) Chen, Ryan Babbush, Joonho Lee, Di Luo, Nathan Wiebe, Dominic Berry,
and Lin Lin for helpful discussions. YT was partly supported by the NSF
Quantum Leap Challenge Institute (QLCI) program through Grant No. OMA-2016245,
and by the Department of Energy under Grant No. FWP-NQISCCAWL. JP was partly
supported by the U.S. Department of Energy Office of Advanced Scientific
Computing Research (DE-NA0003525, DE-SC0020290) and Office of High Energy
Physics (DE-ACO2-07CH11359, DE-SC0018407), the Simons Foundation It from Qubit
Collaboration, the Air Force Office of Scientific Research (FA9550-19-1-0360),
and the National Science Foundation (PHY-1733907). YS was partly supported by
the National Science Foundation RAISE-TAQS 1839204. VVA acknowledges support
from the NSF Quantum Leap Challenge Institute (QLCI) program through Grant No.
OMA-2120757. The Institute for Quantum Information and Matter is an NSF
Physics Frontiers Center. Contributions to this work by NIST, an agency of the
US government, are not subject to US copyright. Any mention of commercial
products does not indicate endorsement by NIST.

## References

  * Arad et al. [2013] I. Arad, A. Kitaev, Z. Landau, and U. Vazirani.  An area law and sub-exponential algorithm for 1D systems.  _arXiv preprint arXiv:1301.1162_ , 2013.  doi: 10.48550/arXiv.1301.1162.
  * Arad et al. [2016] I. Arad, T. Kuwahara, and Z. Landau.  Connecting global and local energy distributions in quantum spin models on a lattice.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2016(3):033301, 2016.  doi: 10.1088/1742-5468/2016/03/033301.
  * Atia and Aharonov [2017] Y. Atia and D. Aharonov.  Fast-forwarding of Hamiltonians and exponentially precise measurements.  _Nature Communications_ , 8(1):1572, Nov 2017\.  doi: 10.1038/s41467-017-01637-7.
  * Banerjee et al. [2012] D. Banerjee, M. Dalmonte, M. Müller, E. Rico, P. Stebler, U.-J. Wiese, and P. Zoller.  Atomic quantum simulation of dynamical gauge fields coupled to fermionic matter: From string breaking to evolution after a quench.  _Physical Review Letters_ , 109(17):175302, 2012\.  doi: 10.1103/PhysRevLett.109.175302.
  * Bañuls et al. [2017] M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, and S. Kühn.  Efficient basis formulation for $(1+1)$-dimensional SU(2) lattice gauge theory: Spectral calculations with matrix product states.  _Physical Review X_ , 7(4):041046, 2017.  doi: 10.1103/PhysRevX.7.041046.
  * Banuls et al. [2020] M. C. Banuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, et al.  Simulating lattice gauge theories within quantum technologies.  _The European physical journal D_ , 74(8):1–42, 2020.  doi: 10.1140/epjd/e2020-100571-8.
  * Bender et al. [2018] J. Bender, E. Zohar, A. Farace, and J. I. Cirac.  Digital quantum simulation of lattice gauge theories in three spatial dimensions.  _New Journal of Physics_ , 20(9):093001, 2018\.  doi: 10.1088/1367-2630/aadb71.
  * Berry and Childs [2012] D. W. Berry and A. M. Childs.  Black-box hamiltonian simulation and unitary implementation.  _Quantum Information & Computation_, 12(1-2):29–62, 2012.  doi: 10.26421/QIC12.1-2.
  * Berry et al. [2006] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders.  Efficient quantum algorithms for simulating sparse Hamiltonians.  _Communications in Mathematical Physics_ , 270(2):359–371, 2006.  doi: 10.1007/s00220-006-0150-x.
  * Berry et al. [2014] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma.  Exponential improvement in precision for simulating sparse Hamiltonians.  In _Proceedings of the forty-sixth annual ACM symposium on Theory of computing_ , pages 283–292, 2014.  doi: 10.1145/2591796.2591854.
  * Berry et al. [2015] D. W. Berry, A. M. Childs, and R. Kothari.  Hamiltonian simulation with nearly optimal dependence on all parameters.  In _2015 IEEE 56th Annual Symposium on Foundations of Computer Science_ , pages 792–809, 2015.  doi: 10.1145/3313276.3316386.
  * Bonet-Monroig et al. [2018] X. Bonet-Monroig, R. Sagastizabal, M. Singh, and T. O’Brien.  Low-cost error mitigation by symmetry verification.  _Physical Review A_ , 98(6):062339, 2018.  doi: 10.1103/PhysRevA.98.062339.
  * Byrnes and Yamamoto [2006] T. Byrnes and Y. Yamamoto.  Simulating lattice gauge theories on a quantum computer.  _Physical Review A_ , 73(2):022328, 2006.  doi: 10.1103/PhysRevA.73.022328.
  * Canonne [2017] C. Canonne.  A short note on Poisson tail bounds.  2017\.  URL http://www.cs.columbia.edu/ ccanonne/files/misc/2017-poissonconcentration.pdf.
  * Chakraborty et al. [2022] B. Chakraborty, M. Honda, T. Izubuchi, Y. Kikuchi, and A. Tomiya.  Classically emulated digital quantum simulation of the schwinger model with a topological term via adiabatic state preparation.  _Phys. Rev. D_ , 105:094503, May 2022.  doi: 10.1103/PhysRevD.105.094503.  URL https://link.aps.org/doi/10.1103/PhysRevD.105.094503.
  * Chang et al. [2011] S.-H. Chang, P. C. Cosman, and L. B. Milstein.  Chernoff-type bounds for the Gaussian error function.  _IEEE Transactions on Communications_ , 59(11):2939–2944, 2011.  doi: 10.1109/TCOMM.2011.072011.100049.
  * Childs and Su [2019] A. M. Childs and Y. Su.  Nearly optimal lattice simulation by product formulas.  _Physical Review Letters_ , 123(5):050503, 2019\.  doi: 10.1103/PhysRevLett.123.050503.
  * Childs et al. [2017] A. M. Childs, R. Kothari, and R. D. Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  _SIAM J. Comput._ , 46(6):1920–1950, 2017.  doi: 10.1137/16m1087072.
  * Childs et al. [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu.  Theory of Trotter error with commutator scaling.  _Physical Review X_ , 11(1):011020, 2021.  doi: 10.1103/PhysRevX.11.011020.
  * Davoudi et al. [2021] Z. Davoudi, N. M. Linke, and G. Pagano.  Toward simulating quantum field theories with controlled phonon-ion dynamics: A hybrid analog-digital approach.  _Phys. Rev. Research_ , 3:043072, Oct 2021.  doi: 10.1103/PhysRevResearch.3.043072.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.3.043072.
  * Del Pino et al. [2018] J. Del Pino, F. A. Schröder, A. W. Chin, J. Feist, and F. J. Garcia-Vidal.  Tensor network simulation of non-Markovian dynamics in organic polaritons.  _Physical Review Letters_ , 121(22):227401, 2018\.  doi: 10.1103/PhysRevLett.121.227401.
  * Dicke [1954] R. H. Dicke.  Coherence in spontaneous radiation processes.  _Physical Review_ , 93(1):99, 1954.  doi: 10.1103/PhysRev.93.99.
  * Fröhlich [1954] H. Fröhlich.  Electrons in lattice fields.  _Advances in Physics_ , 3(11):325–361, 1954.  doi: 10.1080/00018735400101213.
  * Gilyén et al. [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ , pages 193–204, 2019.  doi: 10.1145/3313276.3316366.
  * Giustino [2017] F. Giustino.  Electron-phonon interactions from first principles.  _Reviews of Modern Physics_ , 89(1):015003, 2017\.  doi: 10.1103/RevModPhys.89.015003.
  * Gu et al. [2021] S. Gu, R. D. Somma, and B. Şahinoğlu.  Fast-forwarding quantum evolution.  _Quantum_ , 5:577, 2021.  doi: 10.22331/q-2021-11-15-577.
  * Guo et al. [2012] C. Guo, A. Weichselbaum, J. von Delft, and M. Vojta.  Critical and strong-coupling phases in one-and two-bath spin-boson models.  _Physical Review Letters_ , 108(16):160401, 2012\.  doi: 10.1103/PhysRevLett.108.160401.
  * Haah et al. [2021] J. Haah, M. B. Hastings, R. Kothari, and G. H. Low.  Quantum algorithm for simulating real time evolution of lattice Hamiltonians.  _SIAM Journal on Computing_ , (0):FOCS18–250, 2021.  doi: 10.1137/18M1231511.
  * Hastings [2004] M. B. Hastings.  Locality in quantum and Markov dynamics on lattices and networks.  _Physical Review Letters_ , 93(14):140402, 2004\.  doi: 10.1103/PhysRevLett.93.140402.
  * Hastings [2007] M. B. Hastings.  An area law for one-dimensional quantum systems.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2007(08):P08024, 2007.  doi: 10.1088/1742-5468/2007/08/p08024.
  * Hastings and Koma [2006] M. B. Hastings and T. Koma.  Spectral gap and exponential decay of correlations.  _Communications in Mathematical Physics_ , 265(3):781–804, 2006.  doi: 10.1007/s00220-006-0030-4.
  * Hepp and Lieb [1973] K. Hepp and E. H. Lieb.  On the superradiant phase transition for molecules in a quantized radiation field: the Dicke maser model.  _Annals of Physics_ , 76(2):360–404, 1973.  doi: https://doi.org/10.1016/0003-4916(73)90039-0.
  * Holstein [1959] T. Holstein.  Studies of polaron motion: Part I. the molecular-crystal model.  _Annals of Physics_ , 8(3):325–342, 1959.  doi: https://doi.org/10.1016/0003-4916(59)90002-8.
  * Hubbard [1963] J. Hubbard.  Electron correlations in narrow energy bands.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ , 276(1365):238–257, 1963\.  doi: 10.1098/rspa.1963.0204.
  * Huggins et al. [2021] W. J. Huggins, S. McArdle, T. E. O’Brien, J. Lee, N. C. Rubin, S. Boixo, K. B. Whaley, R. Babbush, and J. R. McClean.  Virtual distillation for quantum error mitigation.  _Phys. Rev. X_ , 11:041036, Nov 2021.  doi: 10.1103/PhysRevX.11.041036.  URL https://link.aps.org/doi/10.1103/PhysRevX.11.041036.
  * Jordan et al. [2012] S. P. Jordan, K. S. Lee, and J. Preskill.  Quantum algorithms for quantum field theories.  _Science_ , 336(6085):1130–1133, 2012.  doi: 10.1126/science.1217069.
  * Jordan et al. [2014] S. P. Jordan, K. S. Lee, and J. Preskill.  Quantum computation of scattering in scalar quantum field theories.  _Quantum Information & Computation_, 14(11-12):1014–1080, 2014.  doi: 10.5555/2685155.2685163.
  * Kan and Nam [2021] A. Kan and Y. Nam.  Lattice quantum chromodynamics and electrodynamics on a universal quantum computer.  _arXiv preprint arXiv:2107.12769_ , 2021.  doi: 10.48550/arXiv.2107.12769.
  * Kivlichan et al. [2018] I. D. Kivlichan, J. McClean, N. Wiebe, C. Gidney, A. Aspuru-Guzik, G. K.-L. Chan, and R. Babbush.  Quantum simulation of electronic structure with linear depth and connectivity.  _Physical Review Letters_ , 120(11):110501, 2018\.  doi: 10.1103/PhysRevLett.120.110501.
  * Klco and Savage [2019] N. Klco and M. J. Savage.  Digitization of scalar fields for quantum computing.  _Physical Review A_ , 99(5):052335, 2019.  doi: 10.1103/PhysRevA.99.052335.
  * Klco et al. [2018] N. Klco, E. F. Dumitrescu, A. J. McCaskey, T. D. Morris, R. C. Pooser, M. Sanz, E. Solano, P. Lougovski, and M. J. Savage.  Quantum-classical computation of Schwinger model dynamics using quantum computers.  _Physical Review A_ , 98(3):032331, 2018.  doi: 10.1103/PhysRevA.98.032331.
  * Klco et al. [2020] N. Klco, M. J. Savage, and J. R. Stryker.  Su(2) non-abelian gauge field theory in one dimension on digital quantum computers.  _Physical Review D_ , 101(7):074512, 2020.  doi: 10.1103/PhysRevD.101.074512.
  * Kloss et al. [2019] B. Kloss, D. R. Reichman, and R. Tempelaar.  Multiset matrix product state calculations reveal mobile Franck-Condon excitations under strong Holstein-type coupling.  _Physical Review Letters_ , 123(12):126601, 2019\.  doi: 10.1103/PhysRevLett.123.126601.
  * Kogut and Susskind [1975] J. Kogut and L. Susskind.  Hamiltonian formulation of Wilson’s lattice gauge theories.  _Physical Review D_ , 11(2):395, 1975.  doi: 10.1103/PhysRevD.11.395.
  * Kühn et al. [2015] S. Kühn, E. Zohar, J. I. Cirac, and M. C. Bañuls.  Non-Abelian string breaking phenomena with matrix product states.  _Journal of High Energy Physics_ , 2015(7):1–26, 2015.  doi: 10.1007/JHEP07(2015)130.
  * Liu and Xin [2020] J. Liu and Y. Xin.  Quantum simulation of quantum field theories as quantum chemistry.  _Journal of High Energy Physics_ , 2020(12):11, Dec 2020.  ISSN 1029-8479.  doi: 10.1007/JHEP12(2020)011.
  * Lloyd [1996] S. Lloyd.  Universal quantum simulators.  _Science_ , 273(5278):1073–1078, 1996.  doi: 10.1126/science.273.5278.1073.
  * Low and Chuang [2017] G. H. Low and I. L. Chuang.  Optimal Hamiltonian simulation by quantum signal processing.  _Physical Review Letters_ , 118(1):010501, 2017\.  doi: 10.1103/physrevlett.118.010501.
  * Low and Chuang [2019] G. H. Low and I. L. Chuang.  Hamiltonian simulation by qubitization.  _Quantum_ , 3:163, 2019.  doi: 10.22331/q-2019-07-12-163.
  * Low and Wiebe [2018] G. H. Low and N. Wiebe.  Hamiltonian simulation in the interaction picture.  _arXiv preprint arXiv:1805.00675_ , 2018.  doi: 10.48550/arXiv.1805.00675.
  * Macridin et al. [2018a] A. Macridin, P. Spentzouris, J. Amundson, and R. Harnik.  Digital quantum computation of fermion-boson interacting systems.  _Physical Review A_ , 98(4), 2018a.  doi: 10.1103/PhysRevA.98.042312.
  * Macridin et al. [2018b] A. Macridin, P. Spentzouris, J. Amundson, and R. Harnik.  Electron-phonon systems on a universal quantum computer.  _Physical Review Letters_ , 121(11), 2018b.  doi: 10.1103/PhysRevLett.121.110504.
  * Magnifico et al. [2021] G. Magnifico, T. Felser, P. Silvi, and S. Montangero.  Lattice quantum electrodynamics in $(3+1)$-dimensions at finite density with tensor networks.  _Nature Communications_ , 12(1):1–13, 2021.  doi: 10.1038/s41467-021-23646-3.
  * McArdle et al. [2019] S. McArdle, X. Yuan, and S. Benjamin.  Error-mitigated digital quantum simulation.  _Physical Review Letters_ , 122:180501, May 2019.  doi: 10.1103/PhysRevLett.122.180501.
  * Moosavian et al. [2019] A. H. Moosavian, J. R. Garrison, and S. P. Jordan.  Site-by-site quantum state preparation algorithm for preparing vacua of fermionic lattice field theories.  _arXiv preprint arXiv:1911.03505_ , 2019.  doi: 10.48550/arXiv.1911.03505.
  * Muschik et al. [2017] C. Muschik, M. Heyl, E. Martinez, T. Monz, P. Schindler, B. Vogell, M. Dalmonte, P. Hauke, R. Blatt, and P. Zoller.  U(1) Wilson lattice gauge theories in digital quantum simulators.  _New Journal of Physics_ , 19(10):103020, 2017\.  doi: 10.1088/1367-2630/aa89ab.
  * Nachtergaele and Sims [2006] B. Nachtergaele and R. Sims.  Lieb-Robinson bounds and the exponential clustering theorem.  _Communications in Mathematical Physics_ , 265(1):119–130, 2006.  doi: 10.1007/s00220-006-1556-1.
  * Nachtergaele et al. [2009] B. Nachtergaele, H. Raz, B. Schlein, and R. Sims.  Lieb-Robinson bounds for harmonic and anharmonic lattice systems.  _Communications in Mathematical Physics_ , 286(3):1073–1098, 2009.  doi: 10.1007/s00220-008-0630-2.
  * Otte [2010] P. Otte.  Boundedness properties of fermionic operators.  _Journal of Mathematical Physics_ , 51(8):083503, 2010.  doi: 10.1063/1.3464264.
  * Pichler et al. [2016] T. Pichler, M. Dalmonte, E. Rico, P. Zoller, and S. Montangero.  Real-time dynamics in U(1) lattice gauge theories with tensor networks.  _Physical Review X_ , 6(1):011023, 2016.  doi: 10.1103/PhysRevX.6.011023.
  * Rajput et al. [2022] A. Rajput, A. Roggero, and N. Wiebe.  Hybridized methods for quantum simulation in the interaction picture.  _Quantum_ , 6:780, 2022.  doi: 10.22331/q-2022-08-17-780.
  * Reinhard et al. [2019] T. E. Reinhard, U. Mordovina, C. Hubig, J. S. Kretchmer, U. Schollwöck, H. Appel, M. A. Sentef, and A. Rubio.  Density-matrix embedding theory study of the one-dimensional Hubbard-Holstein model.  _Journal of chemical theory and computation_ , 15(4):2221–2232, 2019.  doi: 10.1021/acs.jctc.8b01116.
  * Şahinoğlu and Somma [2021] B. Şahinoğlu and R. D. Somma.  Hamiltonian simulation in the low-energy subspace.  _npj Quantum Information_ , 7(1):119, Jul 2021\.  ISSN 2056-6387.  doi: 10.1038/s41534-021-00451-w.
  * Sandhoefer and Chan [2016] B. Sandhoefer and G. K.-L. Chan.  Density matrix embedding theory for interacting electron-phonon systems.  _Physical Review B_ , 94(8):085115, 2016.  doi: 10.1103/PhysRevB.94.085115.
  * Sawaya et al. [2016] N. P. D. Sawaya, M. Smelyanskiy, J. R. McClean, and A. Aspuru-Guzik.  Error sensitivity to environmental noise in quantum circuits for chemical state preparation.  _Journal of Chemical Theory and Computation_ , 12(7):3097–3108, 2016.  doi: 10.1021/acs.jctc.6b00220.
  * Sawaya et al. [2020] N. P. D. Sawaya, T. Menke, T. H. Kyaw, S. Johri, A. Aspuru-Guzik, and G. G. Guerreschi.  Resource-efficient digital quantum simulation of $d$-level systems for photonic, vibrational, and spin-$s$ Hamiltonians.  _npj Quantum Information_ , 6(1):49, Jun 2020\.  ISSN 2056-6387.  doi: 10.1038/s41534-020-0278-0.
  * Schröder and Chin [2016] F. A. Schröder and A. W. Chin.  Simulating open quantum dynamics with time-dependent variational matrix product states: Towards microscopic correlation of environment dynamics and reduced system evolution.  _Physical Review B_ , 93(7):075105, 2016.  doi: 10.1103/PhysRevB.93.075105.
  * Sen [2011] P. Sen.  Achieving the Han-Kobayashi inner bound for the quantum interference channel by sequential decoding.  _arXiv preprint arXiv:1109.0802_ , 2011.  doi: 10.48550/arXiv.1109.0802.
  * Shaw et al. [2020] A. F. Shaw, P. Lougovski, J. R. Stryker, and N. Wiebe.  Quantum algorithms for simulating the lattice Schwinger model.  _Quantum_ , 4:306, 2020.  doi: 10.22331/q-2020-08-10-306.
  * Somma [2015] R. D. Somma.  Quantum simulations of one dimensional quantum systems.  _arXiv preprint arXiv:1503.06319_ , 2015.  doi: 10.48550/arXiv.1503.06319.
  * Su et al. [2021] Y. Su, H.-Y. Huang, and E. T. Campbell.  Nearly tight Trotterization of interacting electrons.  _Quantum_ , 5:495, 2021.  doi: 10.22331/q-2021-07-05-495.
  * Suzuki [1985] M. Suzuki.  Decomposition formulas of exponential operators and Lie exponentials with some applications to quantum mechanics and statistical physics.  _Journal of Mathematical Physics_ , 26(4):601–612, 1985.  doi: 10.1063/1.526596.
  * Tran et al. [2021] M. C. Tran, Y. Su, D. Carney, and J. M. Taylor.  Faster digital quantum simulation by symmetry protection.  _PRX Quantum_ , 2:010323, Feb 2021.  doi: 10.1103/PRXQuantum.2.010323.
  * Verstraete and Cirac [2005] F. Verstraete and J. I. Cirac.  Mapping local Hamiltonians of fermions to local Hamiltonians of spins.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2005(09):P09012, 2005.  doi: 10.1088/1742-5468/2005/09/p09012.
  * Wiese [2013] U.-J. Wiese.  Ultracold quantum gases and lattice systems: quantum simulation of lattice gauge theories.  _Annalen der Physik_ , 525(10-11):777–796, 2013\.  doi: https://doi.org/10.1002/andp.201300104.
  * Woods et al. [2015] M. P. Woods, M. Cramer, and M. B. Plenio.  Simulating bosonic baths with error bars.  _Physical Review Letters_ , 115(13):130401, 2015\.  doi: 10.1103/PhysRevLett.115.130401.
  * Zohar et al. [2012] E. Zohar, J. I. Cirac, and B. Reznik.  Simulating compact quantum electrodynamics with ultracold atoms: Probing confinement and nonperturbative effects.  _Physical Review Letters_ , 109(12):125302, 2012\.  doi: 10.1103/PhysRevLett.109.125302.
  * Zohar et al. [2013] E. Zohar, J. I. Cirac, and B. Reznik.  Cold-atom quantum simulator for SU(2) Yang-Mills lattice gauge theory.  _Physical Review Letters_ , 110(12):125304, 2013\.  doi: 10.1103/PhysRevLett.110.125304.

##  Appendix A Motivating examples

We begin by introducing example quantum systems that we will analyze and
simulate. These include a general model for boson-fermion coupling, U(1)
lattice gauge theory, and SU(2) lattice gauge theory. We refer the reader to
Section J for other common models that can be analyzed within our framework.

Boson-fermion coupling. We assume that there are $N_{f}$ fermionic modes and
$N_{b}$ bosonic modes in the system. We label the fermionic modes by $i,j$ and
bosonic modes by $\alpha$. The $c_{i}$ and $b_{\alpha}$ denote the fermionic
and bosonic annihilation operators respectively. The Hamiltonian takes the
form

$\displaystyle H$ $\displaystyle=H_{f}+H_{fb}+H_{b},$ (14) $\displaystyle
H_{f}$
$\displaystyle=\sum_{ij}t_{ij}c_{i}^{\dagger}c_{j}+\sum_{ij}V_{ijkl}c_{i}^{\dagger}c_{j}^{\dagger}c_{k}c_{l},$
$\displaystyle H_{fb}$ $\displaystyle=\sum_{\alpha
ij}g_{ij}^{(\alpha)}c_{i}^{\dagger}c_{j}X_{\alpha}+\sum_{\alpha
ij}h_{ij}^{(\alpha)}c_{i}^{\dagger}c_{j}P_{\alpha},$ $\displaystyle H_{b}$
$\displaystyle=\frac{1}{2}\sum_{\alpha}\omega_{\alpha}X_{\alpha}^{2}+\frac{1}{2}\sum_{\alpha}\omega_{\alpha}P_{\alpha}^{2},$

where $X_{\alpha}=(b_{\alpha}+b_{\alpha}^{\dagger})/\sqrt{2}$ is the position
operator corresponding to the bosonic mode $\alpha$, and
$P_{\alpha}=i(b^{\dagger}_{\alpha}-b_{\alpha})/\sqrt{2}$ is the momentum
operator. $t=(t_{ij})$, $g^{(\alpha)}=(g^{(\alpha)}_{ij})$, and
$h^{(\alpha)}=(h_{ij}^{(\alpha)})$ are all Hermitian matrices, and
$V=(V_{ijkl})$ is the electron repulsion integral tensor satisfying the usual
symmetry. We remark that the commonly seen Hubbard-Holstein model [33] and the
Fröhlich model [23] both take the above form.

U(1) lattice gauge theory. For notation simplicity we consider only the
$(2+1)$-dimensional theory. Extension to the $(3+1)$-dimensional case is
straightforward. The system consists of a square lattice of $N$ sites. We
denote each site by $x$, and the lattice vector in the horizontal and vertical
directions are noted $n_{1}$ and $n_{2}$ respectively. We use $(x,n_{i})$ to
represent the link between sites $x$ and $x+n_{i}$, $i=1,2$. The links are
sometimes called gauge links.

On each site $x$ we have a fermionic mode whose annihilation operator is
denoted by $\phi_{x}$. Each link consists of a planar rotor, whose
configuration space of states $\ket{\theta}$, with $\theta\in[0,2\pi]$ being
an angle, is equivalent to that of a particle on a ring. An orthonormal basis
of the Hilbert space can be chosen to be

$\ket{k}=\frac{1}{\sqrt{2\pi}}\int_{0}^{2\pi}e^{ik\theta}\ket{\theta}\mathrm{d}\theta,$
(15)

for $k\in\mathbb{Z}$.

In Hilbert space of link $(x,n_{i})$ we define operators $E_{x,n_{i}}$ and
$U_{x,n_{i}}$ through

$E_{x,n_{i}}\ket{k}=k\ket{k},\quad U_{x,n_{i}}\ket{k}=\ket{k-1}.$ (16)

Then the Hamiltonian of the system is

$\displaystyle H$ $\displaystyle=H_{M}+H_{GM}+H_{E}+H_{B},$ (17)
$\displaystyle H_{M}$
$\displaystyle=g_{M}\sum_{x}(-1)^{x}\phi_{x}^{\dagger}\phi_{x},$
$\displaystyle H_{GM}$
$\displaystyle=g_{GM}\sum_{x,i}(\phi_{x}^{\dagger}U_{x,n_{i}}\phi_{x+n_{i}}+\mathrm{h.c.}),$
$\displaystyle H_{E}$ $\displaystyle=g_{E}\sum_{x,i}E_{x,n_{i}}^{2},$
$\displaystyle H_{B}$
$\displaystyle=g_{B}\sum_{P}(\mathrm{Tr}U_{P}+\mathrm{h.c.}),$

where $\sum_{P}$ denotes a summation over all plaquettes $P$. For $P$ whose
lower-left site is $x$, $U_{P}$ is defined as

$U_{P}=U_{x,n_{1}}U_{x+n_{1},n_{2}}U_{x+n_{2},n_{1}}^{\dagger}U_{x,n_{2}}^{\dagger}.$
(18)

The trace $\mathrm{Tr}$ in (17) is not needed here but will be required in the
setting of SU(2) lattice gauge theory. The four terms
$H_{M},H_{GM},H_{E},H_{B}$ describe the fermionic mass (using staggered
fermions [44]), the gauge-matter interaction, the electric energy, and the
magnetic energy respectively.

SU(2) lattice gauge theory. The setup of the SU(2) lattice gauge theory is
very similar to the U(1) case. Here for simplicity we only consider the theory
using the fundamental representation of SU(2). Compared to the U(1) theory,
each site $x$ now contains two fermionic modes, whose annihilation operators
are denoted by $\phi_{x}^{l}$, $l=1,2$. We write
$\phi_{x}=(\phi_{x}^{1},\phi_{x}^{2})^{\top}$. Each link consists of a rigid
rotator whose configuration is described by an element of the group SU(2)
[44]. An orthonormal basis of the link Hilbert space consists of the quantum
states $\ket{jmm^{\prime}}$, where $j,m,m^{\prime}$ are simultaneously either
integers or half-integers with $-j\leq m,m^{\prime}\leq j$. Here $j$ is the
rotator’s total angular momentum, and $m$, $m^{\prime}$ denote the components
of angular momentum along the $z$-axis in the body-fixed and space-fixed
coordinate systems.

The Hamiltonian takes the form (17), and is invariant under $\mathrm{SU}(2)$
transformations acting either from the left or from the right, which may be
interpreted as rotations of the rigid rotator with respect to space-fixed or
body-fixed axes respectively. The operators $E^{2}_{x,n_{i}}$ and
$U_{x,n_{i}}$ are different from the U(1) case. The operator $E^{2}_{x,n_{i}}$
is defined through

$E^{2}_{x,n_{i}}\ket{jmm^{\prime}}=j(j+1)\ket{jmm^{\prime}}.$ (19)

Because $\phi_{x}$ has two components, where each component is a fermionic
mode, $U_{x,n_{i}}$ is a $2\times 2$ matrix, where each of the 4 matrix
entries is an operator acting on the link Hilbert space

$U_{x,n_{i}}=\begin{pmatrix}U^{11}_{x,n_{i}}&U^{12}_{x,n_{i}}\\\
U^{21}_{x,n_{i}}&U^{22}_{x,n_{i}}\end{pmatrix}.$ (20)

An important property that we will use later is

$\braket{j_{1}m_{1}m_{1}^{\prime}}{U^{ll^{\prime}}_{x,n_{i}}}{j_{2}m_{2}m_{2}^{\prime}}=0,\
\text{if }|j_{1}-j_{2}|>1/2,\quad\|U^{ll^{\prime}}_{x,n_{i}}\|\leq 1,$ (21)

which follows from rules for the addition of angular momentum, given that
$U_{x,n_{i}}$ transforms as the $j=1/2$ representation of $\mathrm{SU}(2)$.
Here $\|O\|$ denotes the spectral norm of an operator $O$. We also note that
relative to the basis $\\{\ket{jmm^{\prime}}\\}$,
$U^{ll^{\prime}}_{x,n_{i}}$’s are sparse matrices because

$\braket{j_{1}m_{1}m_{1}^{\prime}}{U^{ll^{\prime}}_{x,n_{i}}}{j_{2}m_{2}m_{2}^{\prime}}=0,\
\text{if }m_{1}-m_{2}\neq l-3/2\ \text{or }m_{1}^{\prime}-m_{2}^{\prime}\neq
l^{\prime}-3/2,$ (22)

due to the conservation of angular momentum along the $z$-axis in the body-
fixed and space-fixed coordinate systems. Here $l-3/2$ and $l^{\prime}-3/2$
are the change of angular momentum as a result of applying
$U^{ll^{\prime}}_{x,n_{i}}$. Eqs. (21) and (22) imply that the matrix
representing $U^{ll^{\prime}}_{x,n_{i}}$ has at most three non-zero elements
in each row and column.

##  Appendix B The common structure

Here we identify a common structure in all the examples introduced in Section
A. We first decompose the entire Hilbert space $\mathcal{H}$ into a direct sum
of subspaces $\mathcal{V}_{\lambda}$ with quantum numbers
$\lambda\in\mathbb{Z}$. The projection operator onto each subspace
$\mathcal{V}_{\lambda}$ is denoted by $\Pi_{\lambda}$. Then
$\sum_{\lambda\in\mathbb{Z}}\Pi_{\lambda}=I$. We consider a class of
Hamiltonians of the form

$H=H_{W}+H_{R},$ (23)

where

$\Pi_{\lambda}H_{W}\Pi_{\lambda^{\prime}}=0,\ \text{if
}|\lambda-\lambda^{\prime}|>1,\quad\|H_{W}\Pi_{[-\Lambda,\Lambda]}\|\leq\chi(\Lambda+1)^{r},\quad[H_{R},\Pi_{\lambda}]=0,$
(24)

for some $\chi>0$, $0\leq r<1$. Here
$\Pi_{[-\Lambda,\Lambda]}=\sum_{|\lambda|\leq\Lambda}\Pi_{\lambda}$. In such a
Hamiltonian, $H_{W}$ changes the quantum number $\lambda$ in the time
evolution while $H_{R}$ preserves it. (24) ensures that $\lambda$ is not
changed too quickly. The first part of (24) ensures that the local quantum
number is changed by at most $\pm 1$ each time the Hamiltonian is applied, and
the second part ensures that the rate of the change is sublinear in the
current local quantum number. The meaning of these conditions will be made
clearer when we discuss the leakage bound in Section C.

Model | $\lambda$ | $r$
---|---|---
Boson-fermion | Bosonic particle number | $1/2$
U(1) LGT | Electric field value | $0$
SU(2) LGT | Total angular momentum | $0$
Table 1: Local quantum number $\lambda$ and the exponent $r$ for models
discussed in Section B.

We check that all the models introduced in Section A satisfy (23) and (24).
For the boson-fermion coupling Hamiltonian defined in (14), fixing a bosonic
mode $\alpha_{0}$, we can decompose the Hilbert space according to the number
of particles in the bosonic mode $\alpha_{0}$, which we denote by $m$. This
means we let $\lambda=m$ and

$\Pi_{m}=\ket{m}\\!\bra{m}_{\alpha_{0}}\otimes I,$ (25)

where $\ket{m}_{\alpha_{0}}$ means the $\ket{m}_{\alpha_{0}}$-particle state
of the mode, and $I$ is the identity operator acting on the rest of the
system. We set $\Pi_{\lambda}=0$ for all $\lambda<0$. We define $H_{W}$ to be
the sum of terms in (24) that change the particle number in mode $\alpha_{0}$:

$H_{W}=\sum_{ij}\left(g^{\alpha_{0}}_{ij}c_{i}^{\dagger}c_{j}X_{\alpha_{0}}+h^{\alpha_{0}}_{ij}c_{i}^{\dagger}c_{j}P_{\alpha_{0}}\right),$
(26)

whereas the rest of the terms in $H$ are collected into $H_{R}$. Because of
the fact that

$\|X_{\alpha_{0}}\Pi_{[0,M]}\|,\|P_{\alpha_{0}}\Pi_{[0,M]}\|<\sqrt{2(M+1)},$
(27)

where $\Pi_{[0,M]}=\sum_{m=0}^{M}\Pi_{m}$, one can see that (24) is satisfied
if we choose $r=1/2$ and

$\chi=\sqrt{2}(\max_{\alpha}\mathrm{Tr}|g^{(\alpha)}|+\max_{\alpha}\mathrm{Tr}|h^{(\alpha)}|),$
(28)

where $|A|=\sqrt{A^{\dagger}A}$ for any matrix $A$. Note that $r$ is not
determined by the highest order terms in terms of the position and momentum
operators, but rather the highest order terms that change the local quantum
number. For example, $X_{\alpha_{0}}^{2}$ and $P_{\alpha_{0}}^{2}$ each on
their own will result in $r=1$, but since in the Hamiltonian they appear
together, $X_{\alpha_{0}}^{2}+P_{\alpha_{0}}^{2}$ preserves the local quantum
number and therefore does not contribute to how large $r$ is.

$\mathrm{Tr}|A|$ is the trace norm of $A$ and it is used here because
$\|\sum_{ij}A_{ij}c_{i}^{\dagger}c_{j}\|\leq\mathrm{Tr}|A|$ for any Hermitian
matrix $A=(A_{ij})$. This can be proved for any matrix $A$ using the singular
value decomposition [59].

In the setting of U(1) lattice gauge theory, again we fix a given link indexed
by $(x_{0},n_{0})$ where $n_{0}\in\\{n_{1},n_{2}\\}$. Then we decompose the
Hilbert space by the electric field value on this link, i.e. we let
$\lambda=k$ and define

$\Pi_{k}=\ket{k}\\!\bra{k}_{(x_{0},n_{0})}\otimes I.$ (29)

Then $H_{W}$ should be chosen as

$H_{W}=g_{GM}(\phi_{x_{0}}^{\dagger}U_{x_{0},n_{0}}\phi_{x_{0}+n_{0}}+\mathrm{h.c.})+g_{B}\sum_{P\ni(x_{0},n_{0})}(\mathrm{Tr}U_{P}+\mathrm{h.c.}).$
(30)

Because of the fact that $\|H_{W}\|\leq 4|g_{B}|+2|g_{GM}|$, (24) is satisfied
if we choose $\chi=4|g_{B}|+2|g_{GM}|$ and $r=0$. Here we have $r=0$ because,
unlike the bosonic position and momentum operators, $U_{x_{0},n_{0}}$ is a
bounded operator.

In the setting of SU(2) lattice gauge theory, again we fix a given link
indexed by $(x_{0},n_{0})$. We decompose the Hilbert space according to the
total angular momentum on this link. This is to say, we let $\lambda=2j$ ($j$
takes half-integer value), and

$\Pi_{2j}=\sum_{-j\leq m,m^{\prime}\leq
j}\ket{jmm^{\prime}}\\!\bra{jmm^{\prime}}_{(x_{0},n_{0})}\otimes I.$ (31)

Here we require $m,m^{\prime}$ to be integers when $j$ is an integer and half-
integers when $j$ is a half integer. Then $H_{W}$ takes the same form as in
(30). Eq. (21) ensures that (24) is satisfied if we choose
$\chi=16|\lambda_{B}|+8|g_{GM}|$ and $r=0$. There is an additional factor of
$4$ in $\chi$ compared to the U(1) case because there are now four operators
$U_{x_{0},n_{0}}^{ll^{\prime}}$ contributing to the growth of the quantum
number instead of one.

More generally, we define $\Pi_{S}$, where $S$ is a set of integers, as

$\Pi_{S}=\sum_{\lambda\in S}\Pi_{\lambda}.$ (32)

In the examples introduced above, we have focused on the quantum numbers on a
single fixed bosonic mode or gauge link, and decomposed the Hilbert space
accordingly. In fact this procedure can be done for every mode and link.
Therefore we sometimes need to designate projection operators for each mode or
link. In the boson-fermion coupling situation, we denote by
$\Pi^{(\alpha)}_{m}$ the projection operator into the subspace with $m$
particles in the bosonic mode $\alpha$. Similarly we define
$\Pi^{(\alpha)}_{S}$ for any integer set $S$. When we need to constrain the
particle number on all modes, we define

$\Pi^{\mathrm{all}}_{S}=\prod_{\alpha}\Pi^{(\alpha)}_{S}.$ (33)

As a general rule, if $\Pi$ is any projection operator, we define
$\overline{\Pi}=I-\Pi$. For lattice gauge theories we adopt similar notations.
For example we use $\Pi^{(x,n)}_{S}$ to denote the projection operator into
the subspace with the quantum number taking value in set $S$ on gauge link
$(x,n)$. Moreover, we sometimes use $\nu$ to index both the bosonic mode and
gauge links when we discuss the two scenarios together. Therefore
$\Pi^{(\nu)}_{S}$ can mean either $\Pi^{(\alpha)}_{S}$ or $\Pi^{(x,n)}_{S}$
depending on the context.

##  Appendix C Truncating an evolved quantum state

Our first goal is to answer the following question: suppose we start from an
initial state with quantum number $\lambda$ between $\pm\Lambda_{0}$, what is
the probability that $|\lambda|$ grows beyond some given $\Lambda$ as the
state evolves for time $t$? To be more concrete, we want to bound

$\left\|\overline{\Pi}_{[-\Lambda,\Lambda]}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|,$
(34)

when $H$ has the structure (23). We call this quantity the _leakage_. As
discussed in the main article, this upper bounds the error of truncating the
quantum state at time $t$ when the initial state has a quantum number between
$-\Lambda_{0}$ and $\Lambda_{0}$.

###  C.1 The short-time leakage bound

We first establish a bound on the leakage defined in (34) for a short time
$t$:

######  Lemma 1 (Short-time leakage bound).

Given Hamiltonian $H=H_{W}+H_{R}$ satisfying (24) with parameters $\chi$ and
$r$, we have

$\left\|\overline{\Pi}_{(-\Lambda_{0}-\Delta,\Lambda_{0}+\Delta)}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|\leq\frac{1}{2^{\Delta-1}(\Delta!)^{1-r}}$

for any $0\leq|t|\leq 1/(2\chi(\Lambda_{0}+1)^{r})$ and integers
$\Lambda_{0},\Delta\geq 0$, where
$\overline{\Pi}_{(-\Lambda_{0}-\Delta,\Lambda_{0}+\Delta)}=I-\Pi_{(-\Lambda_{0}-\Delta,\Lambda_{0}+\Delta)}$.

###### Proof.

This proof is based on rewriting the time evolution using the interaction
picture, and truncating the Dyson series of the new time evolution. Below we
only consider $t>0$. The proof can readily be extended to $t<0$ because we
only need to replace $H$ by $-H$, and the structure in (24) is preserved by
this transformation. First we define

$H_{W}(t)=e^{itH_{R}}H_{W}e^{-itH_{R}}.$ (35)

Then writing the time evolution $e^{-itH}$ in the interaction picture, we have

$e^{-itH}=e^{-itH_{R}}\mathcal{T}e^{-i\int_{0}^{t}H_{W}(s)\mathrm{d}s},$ (36)

where $\mathcal{T}$ is the time-ordering operator. Since $e^{-itH_{R}}$
commutes with $\overline{\Pi}_{(-\Lambda_{0}-\Delta,\Lambda_{0}+\Delta)}$, we
only need to bound

$\left\|\overline{\Pi}_{(-\Lambda_{0}-\Delta,\Lambda_{0}+\Delta)}\mathcal{T}e^{-i\int_{0}^{t}H_{W}(s)\mathrm{d}s}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|.$
(37)

To this end, we consider the partial sum of the Dyson series of
$\mathcal{T}e^{-i\int_{0}^{t}H_{W}(s)\mathrm{d}s}$:

$S(\Delta)=\sum_{k=0}^{\Delta-1}(-i)^{k}\int_{0}^{t}\mathrm{d}t_{k}\cdots\int_{0}^{t_{3}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{1}H_{W}(t_{k})\cdots
H_{W}(t_{2})H_{W}(t_{1}).$ (38)

In order to bound the error from replacing the exact time evolution with this
truncated Dyson series, we need to estimate the norms of terms of the form

$H_{W}(t_{k})\cdots H_{W}(t_{2})H_{W}(t_{1})\Pi_{[-\Lambda,\Lambda]}.$ (39)

Noting that $H_{W}(t)$ can only change $\lambda$ by $\pm 1$, we have

$\Pi_{[-\Lambda-1,\Lambda+1]}H_{W}(t)\Pi_{[-\Lambda,\Lambda]}=H_{W}(t)\Pi_{[-\Lambda,\Lambda]},$
(40)

which implies by (24) that

$\|H_{W}(t)\Pi_{[-\Lambda,\Lambda]}\|=\|H_{W}\Pi_{[-\Lambda,\Lambda]}\|\leq\chi(\Lambda+1)^{r}.$
(41)

By repeatedly applying (40), we have

$H_{W}(t_{k})\cdots
H_{W}(t_{2})H_{W}(t_{1})\Pi_{[-\Lambda_{0},\Lambda_{0}]}=\prod_{j=1}^{k}\left(\Pi_{[-\Lambda_{0}-j,\Lambda_{0}+j]}H_{W}(t_{j})\Pi_{[-\Lambda_{0}-j+1,\Lambda_{0}+j-1]}\right).$
(42)

Then applying (41),

$\displaystyle\left\|H_{W}(t_{k})\cdots
H_{W}(t_{2})H_{W}(t_{1})\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|$
$\displaystyle\leq\prod_{j=1}^{k}\left\|\Pi_{[-\Lambda_{0}-j,\Lambda_{0}+j]}H_{W}(t_{j})\Pi_{[-\Lambda_{0}-j+1,\Lambda_{0}+j-1]}\right\|$
(43)
$\displaystyle\leq\chi^{k}\left(\frac{(\Lambda_{0}+k)!}{\Lambda_{0}!}\right)^{r}.$

From the above inequality, we have

$\displaystyle\left\|\left(\mathcal{T}e^{-i\int_{0}^{t}H_{W}(s)\mathrm{d}s}-S(\Delta)\right)\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|$
$\displaystyle\leq\sum_{k=\Delta}^{\infty}\frac{(\chi
t)^{k}}{k!}\left(\frac{(\Lambda_{0}+k)!}{\Lambda_{0}!}\right)^{r}$ (44)
$\displaystyle\leq\sum_{k=\Delta}^{\infty}\frac{(\chi
t(\Lambda_{0}+1)^{r})^{k}}{(k!)^{1-r}}$
$\displaystyle\leq\frac{1}{(\Delta!)^{1-r}}\sum_{k=\Delta}^{\infty}\frac{1}{2^{k}}$
$\displaystyle=\frac{1}{2^{\Delta-1}(\Delta!)^{1-r}},$

where in the second inequality we have used the fact that

$\frac{(\Lambda_{0}+k)!}{\Lambda_{0}!k!}=\frac{\Lambda_{0}+1}{1}\frac{\Lambda_{0}+2}{2}\cdots\frac{\Lambda_{0}+k}{k}\leq(\Lambda_{0}+1)^{k},$
(45)

and in the third inequality we have used $t\leq 1/(2\chi(\Lambda_{0}+1)^{r})$.

Note that

$\overline{\Pi}_{(-\Lambda_{0}-\Delta,\Lambda_{0}+\Delta)}S(\Delta)\Pi_{[-\Lambda_{0},\Lambda_{0}]}=0$
(46)

because of (40). Therefore,

$\displaystyle\left\|\overline{\Pi}_{(-\Lambda_{0}-\Delta,\Lambda_{0}+\Delta)}\mathcal{T}e^{-i\int_{0}^{t}H_{W}(s)\mathrm{d}s}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|$
$\displaystyle=\left\|\overline{\Pi}_{(-\Lambda_{0}-\Delta,\Lambda_{0}+\Delta)}\left(\mathcal{T}e^{-i\int_{0}^{t}H_{W}(s)\mathrm{d}s}-S(\Delta)\right)\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|$
(47) $\displaystyle\leq\frac{1}{2^{\Delta-1}(\Delta!)^{1-r}}.$

This finishes the proof. ∎

###  C.2 The long-time leakage bound

The long time bound is based on the following decomposition.

######  Lemma 2.

Let $P_{j}$, $\overline{P}_{j}$ be projection operators such that
$P_{j}+\overline{P}_{j}=I$, $j=0,1,\ldots,J$. Then for any
$0=t_{0}<t_{1}<\cdots<t_{J}=t$,

$\overline{P}_{J}e^{-itH}P_{0}=\overline{P}_{J}\sum_{j=1}^{J}e^{-i(t-t_{j})H}\overline{P}_{j}\prod_{j^{\prime}=0}^{j-1}\left(e^{-i(t_{j^{\prime}+1}-t_{j^{\prime}})H}P_{j^{\prime}}\right),$
(48)

which implies

$\|\overline{P}_{J}e^{-itH}P_{0}\|\leq\sum_{j=1}^{J}\|\overline{P}_{j}e^{-i(t_{j}-t_{j-1})H}P_{j-1}\|.$
(49)

At a high level, this lemma suggests that the total leakage (quantified by the
spectral norm) is upper bounded by the sum of the leakage in each time step.
This lemma can be easily proved by induction on $J$. A more intuitive way of
proving it is to write

$e^{-itH}=\prod_{j^{\prime}=0}^{J-1}\left((P_{j^{\prime}+1}+\overline{P}_{j^{\prime}+1})e^{-i(t_{j^{\prime}+1}-t_{j^{\prime}})H}\right),$
(50)

and expand the product on the right-hand side into a sum of terms, each of
which is a string of $P_{j^{\prime}}$ and $\overline{P}_{j^{\prime}}$
interspersed with $e^{-i(t_{j^{\prime}+1}-t_{j^{\prime}})H}$. We then
recombine these terms according to where the first $\overline{P}_{j^{\prime}}$
appears, or if it does not show up at all. The sum of all terms for which the
first $\overline{P}_{j^{\prime}}$ appears in the $j$-th place is

$e^{-i(t-t_{j})H}\overline{P}_{j}\prod_{j^{\prime}=1}^{j-1}\left(e^{-i(t_{j^{\prime}+1}-t_{j^{\prime}})H}P_{j^{\prime}}\right)e^{-i(t_{1}-t_{0})H}.$
(51)

We then sum over $j$, and multiply $\overline{P}_{J}$ and $P_{0}$ to the left
and right respectively, to get the right-hand side of (48).

Figure 2: The truncation threshold $\Lambda(t)$ needed to keep the error below
$\epsilon$ for boson-fermion coupling (left) and lattice gauge theories
(right). Note that $\Lambda(t)$ grows quadratically for the former and
linearly for the latter. The truncation threshold is obtained, according to
(52), in the boson-fermion coupling setting for the Hubbard-Holstein model
(86) with $g=1$, and in the lattice gauge theory setting for the U(1) lattice
gauge theory (17) with $g_{B}=g_{GM}=1/6$.

We now state our long-time leakage bound:

######  Theorem 3 (Long-time leakage bound).

Let $H=H_{W}+H_{R}$ be a Hamiltonian satisfying (24) with parameters $\chi$
and $r$. Defining

$\Lambda(t)=\Lambda_{0}+\Big{\lceil}\frac{1}{\Delta-1}\left((\Lambda_{0}^{1-r}+2\chi|t|(1-r)(\Delta-1))^{\frac{1}{1-r}}-\Lambda_{0}\right)\Big{\rceil}(\Delta-1)$
(52)

for any $t\in\mathbb{R}$ and integers $\Lambda_{0}\geq 0$ and $\Delta>1$, we
have

$\left\|\overline{\Pi}_{[-\Lambda(t),\Lambda(t)]}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|\leq\frac{\lceil((\Lambda_{0}^{1-r}+2\chi|t|(1-r)(\Delta-1))^{\frac{1}{1-r}}-\Lambda_{0})/(\Delta-1)\rceil}{2^{\Delta-1}(\Delta!)^{1-r}}.$

Below we will focus on the case of $t>0$. The $t<0$ case can be dealt with in
the same way because of the reason explained in the proof of Lemma 1.

In Figure 2 we plot the truncation threshold $\Lambda(t)$ needed to ensure the
leakage
$\left\|\overline{\Pi}_{[-\Lambda(t),\Lambda(t)]}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|\leq\epsilon$,
for the boson-fermion coupling setting and the lattice gauge theory setting.
We can see that $\Lambda(t)$ grows quadratically with time for the former and
linearly for the latter. Moreover, very small leakage can be achieved by only
slightly increasing $\Lambda(t)$. This follows from the exponential
suppression of leakage that we will describe in Theorem 5.

The basic idea of the proof is to partition the time evolution into small
segments, and apply the short-time bound in Lemma 1 to each segment. We denote
by $T_{j}$, for $j\geq 0$, the intermediate times where we make the partition.
First we define the instantaneous quantum number

$\Lambda_{j}=\Lambda_{0}+j(\Delta-1),$ (53)

and then choose $T_{j}$ to be

$T_{j}=T_{j-1}+\frac{1}{2\chi(\Lambda_{j-1}+1)^{r}}.$ (54)

From this definition we have

$T_{j}\geq\frac{1}{2\chi(\Delta-1)(1-r)}\left((\Lambda_{0}+j(\Delta-1))^{1-r}-\Lambda_{0}^{1-r}\right),$
(55)

which can be proved using the inequality

$\frac{1}{(\Lambda_{j-1}+1)^{r}}\geq\frac{1}{\Delta-1}\int_{\Lambda_{0}+j(\Delta-1)}^{\Lambda_{0}+(j+1)(\Delta-1)}\frac{\mathrm{d}w}{w^{r}}.$
(56)

To establish the long-time leakage bound for arbitrary time $t$, we first
prove the following lemma:

######  Lemma 4.

Let $H=H_{W}+H_{R}$ be a Hamiltonian satisfying (24) with parameters $\chi$
and $r$. Given integers $\Lambda_{0}\geq 0$ and $\Delta>0$, define
$\Lambda_{j}$ and $T_{j}$ as in (53) and (54). For any $T_{J-1}<t\leq T_{J}$
with integer $J\geq 0$,

$\displaystyle\left\|\overline{\Pi}_{[-\Lambda_{J},\Lambda_{J}]}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|\leq\frac{J}{2^{\Delta-1}(\Delta!)^{1-r}}.$
(57)

###### Proof.

We choose $t_{j}=T_{j}$ for $j=0,1,\ldots,J-1$, and $t_{J}=t$. By Eq. (49) in
Lemma 2, we have

$\displaystyle\left\|\overline{\Pi}_{[-\Lambda_{J},\Lambda_{J}]}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|\leq\sum_{j=1}^{J}\left\|\overline{\Pi}_{[-\Lambda_{j},\Lambda_{j}]}e^{-i(t_{j}-t_{j-1})H}\Pi_{[-\Lambda_{j-1},\Lambda_{j-1}]}\right\|.$
(58)

Because of (54) each $t_{j}-t_{j-1}$ is short enough for us to apply our
short-time bound Lemma 1. This completes the proof. ∎

With this lemma we can prove the theorem by appropriately choosing $J$.

###### Proof of Theorem 3.

We choose $J$ to be the first integer that makes $T_{J}\geq t$. By (55), we
have

$J=\Big{\lceil}\frac{1}{\Delta-1}\left((\Lambda_{0}^{1-r}+2\chi
t(1-r)(\Delta-1))^{\frac{1}{1-r}}-\Lambda_{0}\right)\Big{\rceil}.$ (59)

The claimed bound then follows from Lemma 4. ∎

Theorem 3 might not be very straightforward to apply in many situations
because it involves an extra parameter $\Delta$. A more versatile version
would be

######  Theorem 5.

Let $H=H_{W}+H_{R}$ be a Hamiltonian satisfying (24) with parameters $\chi$
and $r$. For any $t\geq 0$ and integers $\Lambda\geq\Lambda_{0}\geq 0$,

$\left\|\overline{\Pi}_{[-\Lambda,\Lambda]}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|\leq\poly(\chi
t,\Lambda_{0},\Lambda)\exp\left(-\Omega\left(\frac{\Lambda^{1-r}-\Lambda_{0}^{1-r}}{2(1-r)\chi
t+1}\right)\right).$

###### Proof.

In Theorem 3 we choose $\Delta$ so that

$\Delta=\Big{\lfloor}\frac{\Lambda^{1-r}-\Lambda_{0}^{1-r}}{2(1-r)\chi
t+1}\Big{\rfloor}+1.$ (60)

If $\Delta=1$ then the claimed bound holds trivially since
$\left\|\overline{\Pi}_{[-\Lambda,\Lambda]}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|\leq
1$. If $\Delta>1$ we use Theorem 3. For the $\Lambda(t)$ defined in Theorem 3
we have

$\displaystyle\Lambda(t)$
$\displaystyle=\Lambda_{0}+\Big{\lceil}\frac{1}{\Delta-1}\left((\Lambda_{0}^{1-r}+2\chi
t(1-r)(\Delta-1))^{\frac{1}{1-r}}-\Lambda_{0}\right)\Big{\rceil}(\Delta-1)$
(61) $\displaystyle\leq(\Lambda_{0}^{1-r}+2\chi
t(1-r)(\Delta-1))^{\frac{1}{1-r}}+\Delta-1$
$\displaystyle\leq(\Lambda_{0}^{1-r}+2\chi
t(1-r)(\Delta-1)+(\Delta-1))^{\frac{1}{1-r}}$ $\displaystyle\leq\Lambda,$

where in the third line we have used the inequality that
$a^{p}+b\leq(a+b)^{p}$ when $a\geq 0$, $p\geq 1$ and $b$ (to be chosen as
$\Delta-1$) is a non-negative integer. Using the fact that
$2^{1-\Delta}\left(\Delta!\right)^{-1/2}=e^{-\Omega(\Delta)}$, the claim
follows immediately from Theorem 3. ∎

If we want to ensure that truncating at a threshold $\Lambda$ has an error
$\left\|\overline{\Pi}_{[-\Lambda,\Lambda]}e^{-itH}\Pi_{[-\Lambda_{0},\Lambda_{0}]}\right\|\leq\epsilon$,
then by Theorem 5 we can choose

$\Lambda^{1-r}=\Lambda_{0}^{1-r}+\widetilde{\mathcal{O}}((\chi
t+1)\polylog(\epsilon^{-1})).$

This is the scaling given in (3) of the main text.

##  Appendix D Truncating the Hamiltonian

In this section we consider the problem of replacing an unbounded Hamiltonian
$H$, such as one describing boson-fermion interactions or lattice gauge
theories, with a bounded Hamiltonian, while keeping the error in time
evolution small. More precisely, we want to construct some bounded
$\widetilde{H}$ such that

$\left\|(e^{-it\widetilde{H}}-e^{-itH})\Pi^{\mathrm{all}}_{[-\Lambda_{0},\Lambda_{0}]}\right\|$
(62)

is sufficiently small.

In the previous section we have focused on a single bosonic mode or gauge
link, but here the truncation needs to be performed for every bosonic mode or
gauge link, and we assume there are $\mathcal{N}$ of them in the system. To
simplify the discussion, we use $\nu$ to index either bosonic modes or gauge
links, replacing the indices $\alpha$ and $(x,n)$. Therefore we have

$\Pi_{S}^{\mathrm{all}}=\prod_{\nu=1}^{\mathcal{N}}\Pi_{S}^{(\nu)}.$ (63)

Note also that all projection operators $\Pi_{S}^{(\nu)}$ commute with each
other. Then for each $\nu$, there is a decomposition of the Hamiltonian
$H=H_{W}^{(\nu)}+H_{R}^{(\nu)}$.

We will establish the following bound.

######  Theorem 6 (Hamiltonian truncation).

Let $H$ be a Hamiltonian such that $H=H_{W}^{(\nu)}+H_{R}^{(\nu)}$ satisfies
(24) with parameters $\chi$ and $r$ for every $1\leq\nu\leq\mathcal{N}$.
Assume that all projection operators $\Pi^{(\nu)}_{\lambda}$ commute with each
other. For any integers $\widetilde{\Lambda}\geq\Lambda_{0}\geq 0$, define

$\widetilde{H}=\Pi^{\mathrm{all}}_{[-\widetilde{\Lambda},\widetilde{\Lambda}]}H\Pi^{\mathrm{all}}_{[-\widetilde{\Lambda},\widetilde{\Lambda}]}.$
(64)

Then for any $t\in\mathbb{R}$,

$\left\|(e^{-it\widetilde{H}}-e^{-itH})\Pi^{\mathrm{all}}_{[-\Lambda_{0},\Lambda_{0}]}\right\|\leq\sqrt{\mathcal{N}}\mathcal{A}(\widetilde{\Lambda})\poly(\chi|t|,\Lambda_{0},\widetilde{\Lambda})\exp\left(-\Omega\left(\frac{\widetilde{\Lambda}^{1-r}-\Lambda_{0}^{1-r}}{4(1-r)\chi|t|+1}\right)\right),$
(65)

where

$\mathcal{A}(\Lambda)=\left\|\left[H,\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}\right]\right\|.$
(66)

We recall that $r=1/2$ for boson-fermion coupling and $r=0$ for lattice gauge
theories. We also note that for both boson-fermion coupling and the lattice
gauge theories, $\mathcal{A}(\Lambda)$ can be bounded by a polynomial of the
Hamiltonian coefficients and $\Lambda$. This is because

$\left[H,\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}\right]=(H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]})^{2}-(\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}H)^{2},$
(67)

and in all examples we discussed in Section A, the norm of
$H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}$ is bounded by a function that is
linear in all coefficients, and linear or quadratic in $\Lambda$ for boson-
fermion coupling and lattice gauge theories respectively.

######  Remark 7.

For the boson-fermion coupling in (14), to ensure that the error in truncating
the time evolution (62) is below $\epsilon$, Theorem 6 implies that a
threshold of

$\widetilde{\Lambda}=\left(\sqrt{\Lambda_{0}}+\widetilde{\mathcal{O}}((\chi|t|+1)\polylog(\mathrm{coef},\Lambda_{0},|t|,\epsilon^{-1}))\right)^{2}$
(68)

suffices, where $\mathrm{coef}$ includes all the coefficients in the model
(14) as well as the number of bosonic and fermionic modes.

Similarly, for U(1) and SU(2) lattice gauge theories, it suffices to truncate
the electric field and total angular momentum, for the two situations
respectively, at

$\widetilde{\Lambda}=\Lambda_{0}+\widetilde{\mathcal{O}}((\chi|t|+1)\polylog(\mathrm{coef},\Lambda_{0},|t|,\epsilon^{-1})),$
(69)

where $\mathrm{coef}$ includes all the coefficients in the models (17) as well
as the number of lattice sites. The values of $\chi$ in all cases are
discussed in Section B.

To establish Theorem 6, we use the following properties of the truncated
Hamiltonian
$\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}$:

$\displaystyle(H-\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]})\Pi^{\mathrm{all}}_{[-\Lambda+1,\Lambda-1]}=0,$
(70)
$\displaystyle\left[H,\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}H\Pi^{\mathrm{all}}_{[-\Lambda,\Lambda]}\right]\Pi^{\mathrm{all}}_{[-\Lambda+2,\Lambda-2]}=0,$
(71)

which follow immediately from the fact that for each $\nu$, $H$ can only
change the quantum number $\lambda$ by $\pm 1$.

Because we are now studying the whole system rather than a single mode or
link, we need to bound the total leakage from the leakage at each individual
$\nu$. This is done through a union bound, as given in the following lemma:

######  Lemma 8.

Let $\Pi^{(\nu)}_{\lambda}$ be projections all commuting with each other. For
any operator $A$ and set $S\subset\mathbb{Z}$, we have

$\left\|\overline{\Pi}^{\mathrm{all}}_{S}A\right\|^{2}\leq\sum_{\nu=1}^{\mathcal{N}}\left\|\overline{\Pi}^{(\nu)}_{S}A\right\|^{2},$
(72)

where $\overline{\Pi}^{\mathrm{all}}_{S}=I-{\Pi}^{\mathrm{all}}_{S}$.

###### Proof.

Because $\Pi^{(\nu)}_{\lambda}$ commute with each other, they can be
simultaneously diagonalized, and by the union bound we have

$\overline{\Pi}^{\mathrm{all}}_{S}\preceq\sum_{\nu}\overline{\Pi}^{(\nu)}_{S},$
(73)

which in turn leads to

$\left\|\overline{\Pi}^{\mathrm{all}}_{S}A\right\|^{2}=\left\|A^{\dagger}\overline{\Pi}^{\mathrm{all}}_{S}A\right\|\leq\left\|A^{\dagger}\sum_{\nu}\overline{\Pi}^{(\nu)}_{S}A\right\|\leq\sum_{\nu}\left\|\overline{\Pi}^{(\nu)}_{S}A\right\|^{2},$
(74)

where the first step can be proven using the singular value decomposition, and
we have used the fact that $A^{\dagger}\overline{\Pi}^{\mathrm{all}}_{S}A$ and
$A^{\dagger}\overline{\Pi}^{(\nu)}_{S}A$ are all positive semidefinite for the
later steps. ∎

Note that the above union bound actually holds even when the commutativity
assumption about $\Pi^{(\nu)}_{\lambda}$ is dropped; see [68, Lemma 2] for
details. With these we are ready to prove the main result of this section.

###### Proof of Theorem 6.

Using the formula for the error of the first-order Trotter decomposition [72,
Eq. (3.4)], we have

$e^{it\widetilde{H}}e^{-itH}=e^{-it(H-\widetilde{H})}+\int_{0}^{t}\mathrm{d}s_{1}\int_{0}^{s_{1}}\mathrm{d}s_{2}\
e^{is_{1}\widetilde{H}}e^{-i(s_{1}-s_{2}){H}}[\widetilde{H},H]e^{-is_{2}{H}}e^{-i(t-s_{1})(H-\widetilde{H})}.$
(75)

We then use this, along with invariance of the spectral norm under
multiplication by a unitary operator, to bound the truncation error as:

$\displaystyle\left\|\left(e^{-it\widetilde{H}}-e^{-itH}\right)\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|$
$\displaystyle=\left\|\left(I-e^{it\widetilde{H}}e^{-itH}\right)\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|$
(76)
$\displaystyle\leq\left\|\left(I-e^{-it(H-\widetilde{H})}\right)\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|$
$\displaystyle\quad+\int_{0}^{t}\mathrm{d}s_{1}\int_{0}^{s_{1}}\mathrm{d}s_{2}\left\|[\widetilde{H},H]e^{-is_{2}{H}}e^{-i(t-s_{1})(H-\widetilde{H})}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|.$

Now if we choose $\widetilde{\Lambda}\geq\Lambda_{0}+1$, then by (70),

$\left(I-e^{-it(H-\widetilde{H})}\right)\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}=0.$
(77)

As a result the second line of (76) is $0$. We now only need to bound the
integrand in the third line.

For this integrand we have

$\displaystyle\left\|[\widetilde{H},H]e^{-is_{2}{H}}e^{-i(t-s_{1})(H-\widetilde{H})}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|$
$\displaystyle=\left\|[\widetilde{H},H]e^{-is_{2}{H}}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|$
(78)
$\displaystyle\leq\left\|[\widetilde{H},H]\Pi_{[-\Lambda^{\prime},\Lambda^{\prime}]}^{\mathrm{all}}e^{-is_{2}{H}}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|$
$\displaystyle\quad+\left\|[\widetilde{H},H]\overline{\Pi}_{[-\Lambda^{\prime},\Lambda^{\prime}]}^{\mathrm{all}}e^{-is_{2}{H}}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|,$

for some $\Lambda^{\prime}$ to be chosen. We choose
$\Lambda^{\prime}=\widetilde{\Lambda}-2$. With this choice and (71) we have

$[\widetilde{H},H]\Pi_{[-\Lambda^{\prime},\Lambda^{\prime}]}^{\mathrm{all}}=0.$
(79)

This eliminates the right-hand side on the second line of (78). Therefore we
are only left with the third line of (78) to deal with.

We apply Theorem 5, as well as Lemma 8, to get

$\left\|\overline{\Pi}_{[-\Lambda^{\prime},\Lambda^{\prime}]}^{\mathrm{all}}e^{-is_{2}{H}}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|\leq\sqrt{\mathcal{N}}\poly(\chi
t,\Lambda_{0},\Lambda^{\prime})\exp\left(-\Omega\left(\frac{\Lambda^{\prime
1-r}-\Lambda_{0}^{1-r}}{2(1-r)\chi t+1}\right)\right),$ (80)

where we have used the fact that $s_{2}\leq t$. Substituting this bound into
(78) and then (76), we have

$\left\|\left(e^{-it\widetilde{H}}-e^{-itH}\right)\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|\leq\sqrt{\mathcal{N}}\mathcal{A}(\widetilde{\Lambda})\poly(\chi
t,\Lambda_{0},\widetilde{\Lambda})\exp\left(-\Omega\left(\frac{\widetilde{\Lambda}^{1-r}-\Lambda_{0}^{1-r}}{2(1-r)\chi
t+1}\right)\right).$ (81)

In the above derivation we used the fact that $\Lambda^{\prime
1-r}=(\widetilde{\Lambda}-2)^{1-r}=\widetilde{\Lambda}^{1-r}+o(1)$ as
$\widetilde{\Lambda}\to\infty$. This completes the proof of the theorem.

∎

One can ask the following question about the proof above: can the energy-based
truncation threshold proposed in [36], and discussed in detail in Section K,
be justified as a truncation threshold for Hamiltonian truncation, through a
proof that is similar to the proof above? We remark that this may require
different assumptions and the proof will need to be substantially modified. If
one were to use the above proof strategy, together with the energy-based
truncation threshold, to derive a truncation threshold for the Hamiltonian
truncation, then an important obstacle is bounding the third line in (78).
This line is bounded, in the proof above, through

$\left\|[\widetilde{H},H]\overline{\Pi}_{[-\Lambda^{\prime},\Lambda^{\prime}]}^{\mathrm{all}}e^{-is_{2}{H}}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|\leq\left\|[\widetilde{H},H]\right\|\left\|\overline{\Pi}_{[-\Lambda^{\prime},\Lambda^{\prime}]}^{\mathrm{all}}e^{-is_{2}{H}}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|.$

On the right-hand side, $\left\|[\widetilde{H},H]\right\|$ grows polynomially
with the truncation threshold $\widetilde{\Lambda}$, while
$\left\|\overline{\Pi}_{[-\Lambda^{\prime},\Lambda^{\prime}]}^{\mathrm{all}}e^{-is_{2}{H}}\Pi_{[-\Lambda_{0},\Lambda_{0}]}^{\mathrm{all}}\right\|$
decays subexponentially with $\Lambda^{\prime}$. Therefore asymptotically the
latter decays faster than the former and consequently we can reach an
arbitrarily high precision. If we could only use the energy-based truncation
threshold, then the latter term only decays polynomially with
$\Lambda^{\prime}$, and as a result a careful comparison between the rates of
growth and decay of the two terms would be needed, and we could only reach an
arbitrarily high precision when the latter decays faster than the former. This
would require further assumptions not included in our framework.

Moreover, the most appealing feature of the energy-based truncation threshold
is that it does not depend on time. However, suppose one could overcome the
above mentioned difficulty; then the energy-based quantum state truncation
threshold would lead to a time-dependent Hamiltonian truncation threshold,
because of the integration over time in (75), and thus the above appealing
feature no longer holds.

##  Appendix E Hamiltonian simulations using the HHKL decomposition

In this section we consider performing Hamiltonian simulation for U(1) and
SU(2) lattice gauge theories and boson-fermion coupling. The basic idea is to
simulate the truncated Hamiltonian $\widetilde{H}$ defined in (64) as opposed
to the unbounded $H$, with the truncation threshold $\widetilde{\Lambda}$
chosen according to (68) and (69) for boson-fermion coupling and lattice gauge
theories respectively.

###  E.1 Simulating lattice gauge theories

In this section we propose an algorithm to simulate the time evolution of the
$(D+1)$-dimensional U(1) and SU(2) lattice gauge theories whose Hamiltonians
are of the form (17). The goal is to perform simulation of a square lattice
consisting of $N$ sites up to time $T$ with an error at most $\epsilon$. This
algorithm is based on a combination of the HHKL decomposition [28] and
interaction picture Hamiltonian simulation [50]. We will show that the
simulation can be done with gate complexity
$\mathcal{O}(NT\polylog(\Lambda_{0}NT\epsilon^{-1}))$, where $\epsilon$ is the
allowed error, assuming that the initial state is in the span of states whose
quantum number is in the range $[-\Lambda_{0},\Lambda_{0}]$ for each gauge
link.

As mentioned above, we will be simulating $\widetilde{H}$ instead of $H$, and
the resulting error has been analyzed in Section D. We use $\widetilde{H}_{E}$
to denote the truncated electric field part of the Hamiltonian, i.e.
$\widetilde{H}_{E}=\Pi^{\mathrm{all}}_{[-\widetilde{\Lambda},\widetilde{\Lambda}]}H_{E}\Pi^{\mathrm{all}}_{[-\widetilde{\Lambda},\widetilde{\Lambda}]}$,
and we adopt similar notation $\widetilde{H}_{M}$, $\widetilde{H}_{GM}$, and
$\widetilde{H}_{B}$ for the other three parts. Moreover we denote

$\widetilde{E}^{2}_{x,n}=\Pi^{(x,n)}_{[-\widetilde{\Lambda},\widetilde{\Lambda}]}E^{2}_{x,n}\Pi^{(x,n)}_{[-\widetilde{\Lambda},\widetilde{\Lambda}]}.$
(82)

Note that the Hamiltonians for lattice gauge theories, both the original $H$
and $\widetilde{H}$, consist of geometrically local terms, and to achieve a
linear scaling in both system size and time we consider using the HHKL
decomposition developed in [28].

####  E.1.1 The HHKL decomposition

We first use [28, Lemma 6] to decompose the time evolution of the entire
system into evolution of blocks, each of which, denoted by $\mathcal{B}$, has
size $\ell^{D}=\mathcal{O}(\polylog(NT\epsilon^{-1}))$ and we only need to
simulate its evolution for time $\tau=\mathcal{O}(1)$. The entire time
evolution is divided into $\mathcal{O}(T)$ segments and there are
$\mathcal{O}(N)$ such blocks within each segment.

The original [28, Lemma 6] requires that all the local terms in the target
Hamiltonian have norm bounded by a constant. However, local terms in
$\widetilde{H}$ have norm depending on the truncation threshold
$\widetilde{\Lambda}$, which scales with the system size $N$, time $T$, and
allowed error $\epsilon$. We address this issue as follows. The only local
terms that are not bounded by a constant are the electric field terms in
$\widetilde{H}_{E}$, i.e. $g_{E}\widetilde{E}^{2}_{x,n}$ for each $(x,n)$, and
each of these terms only acts on a single gauge link. We call such terms, i.e.
terms that act only on a single lattice site or gauge link (which can be seen
as a lattice site as well for this purpose), _on-site_ interactions. In Lemma
13 of Section I we show that on-site interactions do not change the Lieb-
Robinson velocity. Therefore, even with the terms in $\widetilde{H}_{E}$, the
system still has a constant Lieb-Robinson velocity, and as a result we can
invoke [28, Lemma 6] to decompose the time evolution.

####  E.1.2 Simulating the blocks

We now consider simulating the dynamics of each individual block $\mathcal{B}$
of size $\ell^{D}$ for short time $\tau$. The Hamiltonian for $\mathcal{B}$,
which we denote by $\widetilde{H}^{\mathcal{B}}$, includes all the local terms
in $\widetilde{H}$ that only act on sites and links within the block
$\mathcal{B}$. As discussed in Section A, each local term can be represented
by a sparse matrix with respect to the basis discussed in Section B, and can
therefore be encoded by a quantum walk operator [8, 10, 18]. A sum of these
terms can be encoded in an unitary using the linear combination of unitaries
(LCU) method [10, 18]. In this way we have an encoding (known as “block-
encoding” or “standard-form” in [49, 24]) of the Hamiltonian
$\widetilde{H}^{\mathcal{B}}$, i.e. a unitary with
$\widetilde{H}^{\mathcal{B}}$ as a matrix block, with a subnormalization
factor $\mathcal{O}(\ell^{D}\widetilde{\Lambda}^{2})$. Using the Hamiltonian
simulation algorithm for encoded Hamiltonians [49], we can then simulate the
time evolution of the block $\mathcal{B}$ with gate-complexity
$\mathcal{O}(\ell^{2D}\widetilde{\Lambda}^{2}\tau)=\mathcal{O}(\widetilde{\Lambda}^{2}\polylog(NT\epsilon^{-1}))$,
which results in a total gate complexity of
$\widetilde{\mathcal{O}}(NT(\Lambda_{0}+T)^{2}\polylog(\epsilon^{-1}))$. This
is however not the best method in terms of asymptotic complexity.

The polynomial dependence on $\widetilde{\Lambda}$ can be improved to be poly-
logarithmic using the interaction-picture simulation technique developed in
[50]. We group the local terms in $\mathcal{B}$ into
$\widetilde{H}_{M}^{\mathcal{B}}$, $\widetilde{H}_{GM}^{\mathcal{B}}$,
$\widetilde{H}_{E}^{\mathcal{B}}$, and $\widetilde{H}_{B}^{\mathcal{B}}$
depending on whether the term describes fermionic mass energy, gauge-matter
interaction, electric field energy, or magnetic field energy. Then the
polynomial dependence on $\widetilde{\Lambda}$ comes only from
$\widetilde{H}_{E}^{\mathcal{B}}$. We note that the time evolution under
$\widetilde{H}_{E}^{\mathcal{B}}$ can be fast-forwarded, i.e. the number of
gates required to implement it has a poly-logarithmic dependence on the
evolution time multiplied by the Hamiltonian norm. To be more specific, the
time evolution due to each electric field term $g_{E}\widetilde{E}^{2}_{x,n}$
for time $t$ can be implemented with gate complexity
$\mathcal{O}(\polylog(\widetilde{\Lambda}t))$ because this term is represented
by a diagonal matrix in both U(1) and SU(2) settings (see (16) and (19) for
the two settings respectively). And all of these terms act on different gauge
links and therefore commute with each other. To implement
$e^{-it\widetilde{H}_{E}^{\mathcal{B}}}$ we only need to evolve these terms
separately, and thus pay a cost of
$\mathcal{O}(\ell^{D}\polylog(\widetilde{\Lambda}t))=\mathcal{O}(\polylog(NT\widetilde{\Lambda}t\epsilon^{-1}))$.

Now instead of directly simulating the Hamiltonian
$\widetilde{H}^{\mathcal{B}}$, we simulate

$\widetilde{H}^{\mathcal{B}}_{I}(t)=e^{it\widetilde{H}_{E}^{\mathcal{B}}}(\widetilde{H}_{M}^{\mathcal{B}}+\widetilde{H}_{GM}^{\mathcal{B}}+\widetilde{H}_{B}^{\mathcal{B}})e^{-it\widetilde{H}_{E}^{\mathcal{B}}}.$
(83)

The original time evolution $e^{-it\widetilde{H}^{\mathcal{B}}}$ and the
interaction picture evolution are related through

$e^{-it\widetilde{H}^{\mathcal{B}}}=e^{-it\widetilde{H}_{E}^{\mathcal{B}}}\mathcal{T}e^{-i\int_{0}^{t}\mathrm{d}s\widetilde{H}^{\mathcal{B}}_{I}(s)}.$
(84)

It then suffices to implement
$\mathcal{T}e^{-i\int_{0}^{t}\mathrm{d}s\widetilde{H}^{\mathcal{B}}_{I}(s)}$,
namely the time evolution due to the time-dependent Hamiltonian
$\widetilde{H}^{\mathcal{B}}_{I}(s)$. We accomplish this using the truncated
Dyson series method in [50, Corollary 4]. The time-dependent matrix encoding
in [50, Definition 2] can be constructed using the the encoding of the local
Hamiltonian terms as well as the fast-forwarding of
$\widetilde{H}_{E}^{\mathcal{B}}$ discussed above. This yields a gate
complexity
$\mathcal{O}(\ell^{2D}\tau\polylog(NT\widetilde{\Lambda}\epsilon^{-1}))=\mathcal{O}(\polylog(NT\widetilde{\Lambda}\epsilon^{-1}))$
for implementing the interaction picture time evolution and consequently
$e^{-it\widetilde{H}^{\mathcal{B}}}$ can be implemented with the same gate
complexity scaling. Note that here we want to keep the error for simulating
this block to be at most $\mathcal{O}(N^{-1}T^{-1}\epsilon)$ instead of
$\epsilon$. This however does not significantly increase the asymptotic
scaling because the scaling with respect to the allowed error is poly-
logarithmic.

There are in total $\mathcal{O}(NT)$ such simulations to perform for all the
$\mathcal{O}(N)$ blocks and $\mathcal{O}(T)$ times steps. Therefore the total
gate complexity for implementing the time evolution of the entire system is
$\mathcal{O}(NT\polylog(NT\widetilde{\Lambda}\epsilon^{-1}))$. Using the
truncation threshold given in (69), we have that the total gate complexity for
simulating the U(1) or SU(2) lattice gauge theory with $N$ sites up to time
$T$ and allowed error $\epsilon$ is
$\mathcal{O}(NT\polylog(\Lambda_{0}NT\epsilon^{-1}))$, provided that the
initial state is in the support of
$\Pi^{\mathrm{all}}_{[-\Lambda_{0},\Lambda_{0}]}$, i.e. the quantum numbers on
each gauge link are in the interval $[-\Lambda_{0},\Lambda_{0}]$.

We remark that in using the HHKL decomposition, we need to preserve the
locality of fermionic operators after the Jordan-Wigner transformation. This
can be done by introducing auxiliary fermionic modes as discussed in [28]
using the method developed in [74].

###  E.2 Simulating boson-fermion coupling

In this section we consider simulating the Hubbard-Holstein model [33], which
is the simplest model describing the electron-phonon interaction. This model
is defined on a $D$-dimensional lattice, and each side of the lattice contains
$L$ sites where $L^{D}=N$. Each site $x$ in the lattice contains two fermionic
modes $c_{x,\sigma}$ (with $\sigma$ denoting either spin up and down) and a
bosonic mode $b_{x}$. We are interested in the case where $D$ is a constant.
The Hamiltonian is

$H=H_{f}+H_{fb}+H_{b},$ (85)

where $H_{f}$ is the Hamiltonian of the Fermi-Hubbard model:

$H_{f}=-\sum_{\braket{x,x^{\prime}},\sigma}(c_{x,\sigma}^{\dagger}c_{x^{\prime},\sigma}+\mathrm{h.c.})+U\sum_{x=1}^{N}(n_{x,\uparrow}-\frac{1}{2})(n_{x,\downarrow}-\frac{1}{2})-\mu\sum_{r=1}^{N}n_{x},$
(86)

and

$H_{fb}=g\sum_{x=1}^{N}(b_{x}^{\dagger}+b_{x})(n_{x,\uparrow}+n_{x,\downarrow}-1)\quad
H_{b}=\omega_{0}\sum_{x=1}^{N}b_{x}^{\dagger}b_{x},$ (87)

are the boson-fermion coupling part and bosonic part respectively. The lattice
sites are indexed by $x$ and $x^{\prime}$, and spins are indexed by $\sigma$.
It is easy to verify that this model satisfies the general form of boson-
fermion coupling in (14).

Here, we propose an algorithm that simulates the above model up to time $T$
and error $\epsilon$ with gate complexity
$\widetilde{\mathcal{O}}(NT(\sqrt{\Lambda_{0}}+T)\polylog(\epsilon^{-1}))$,
assuming the initial state has no more than $\Lambda_{0}$ particles on each
bosonic mode. Just like in the previous section the algorithm is based on HHKL
decomposition [28] and interaction picture Hamiltonian simulation [50].

First we replace the exact Hamiltonian $H$ with the truncated Hamiltonian
$\widetilde{H}$ in (64) and simulate the evolution of $\widetilde{H}$. The
resulting error is analyzed in Section D. We also denote different parts of
the Hamiltonian after truncation by $\widetilde{H}_{f}$, $\widetilde{H}_{fb}$,
and $\widetilde{H}_{b}$.

We apply the HHKL method to decompose the entire time evolution into evolution
of blocks, each of which is denoted by $\mathcal{B}$, for short time $\tau$.
Here again we encounter local terms whose norms are not bounded by a constant,
and in this case these terms are contained in $\widetilde{H}_{fb}$ and
$\widetilde{H}_{b}$. With the help of Lemma 13 however, we can still apply the
HHKL decomposition because $\widetilde{H}_{fb}$ and $\widetilde{H}_{b}$ are
both on-site. So far the algorithm proceeds in a similar way as that for the
lattice gauge theories.

Then we apply interaction picture Hamiltonian simulation to simulate the
evolution in each block $\mathcal{B}$. We denote by
$\widetilde{H}_{f}^{\mathcal{B}}$, $\widetilde{H}_{fb}^{\mathcal{B}}$, and
$\widetilde{H}_{b}^{\mathcal{B}}$ the fermionic, coupling, and bosonic terms
respectively. Here, the terms in $\widetilde{H}_{b}^{\mathcal{B}}$ can still
be fast-forwarded the same way as the electric field terms in lattice gauge
theories. However, it is not known whether the boson-fermion coupling terms in
$\widetilde{H}_{fb}^{\mathcal{B}}$ can be fast-forwarded. Therefore when we
simulate the interaction picture Hamiltonian

$H_{I}^{\mathcal{B}}=e^{it\widetilde{H}_{b}^{\mathcal{B}}}(\widetilde{H}_{f}^{\mathcal{B}}+\widetilde{H}_{fb}^{\mathcal{B}})e^{-it\widetilde{H}_{b}^{\mathcal{B}}},$
(88)

the dependence on the truncation threshold is not poly-logarithmic. Rather a
factor of $\sqrt{\widetilde{\Lambda}}$ shows up in the subnormalization factor
of the encoding of the Hamiltonian because

$\|\Pi^{\mathrm{all}}_{[0,\widetilde{\Lambda}]}(b_{x}^{\dagger}+b_{x})(n_{x,\uparrow}+n_{x,\downarrow}-1)\Pi^{\mathrm{all}}_{[0,\widetilde{\Lambda}]}\|\leq
2\sqrt{\widetilde{\Lambda}+1}.$ (89)

There is some subtle difference between the subnormalization factor and the
spectral norm, but in the present case they have the same asymptotic scaling.
The number of gates required to simulate the time evolution of a block
$\mathcal{B}$ for time $\tau$ is then
$\mathcal{O}(\ell^{2D}\sqrt{\widetilde{\Lambda}}\tau\polylog(\widetilde{\Lambda}NT\epsilon^{-1}))=\mathcal{O}(\sqrt{\widetilde{\Lambda}}\polylog(\widetilde{\Lambda}NT\epsilon^{-1}))$
with a target accuracy of $\mathcal{O}(N^{-1}T^{-1}\epsilon)$ for each block.

We need to perform $\mathcal{O}(NT)$ such simulations and the number of
required gates is therefore
$\widetilde{\mathcal{O}}(NT\sqrt{\widetilde{\Lambda}}\polylog(\epsilon^{-1}))$.
The truncation threshold $\widetilde{\Lambda}$ can be chosen according to
(68). As a result the total gate complexity for simulating an $N$-site
Hubbard-Holstein model for time $T$ up to error $\epsilon$ is
$\widetilde{\mathcal{O}}(NT(\sqrt{\Lambda_{0}}+T)\polylog(\epsilon^{-1}))$,
assuming the initial state has at most $\Lambda_{0}$ particles in each bosonic
mode.

##  Appendix F Simulating boson-fermion coupling using Trotterization

In this section, we consider a generic Hamiltonian of a fermion-boson coupled
system. We assume there are $N_{f}$ fermionic modes (labeled by $i,j$) and
$N_{b}$ bosonic modes (labeled by $\alpha$). We have

$\displaystyle H$ $\displaystyle=\sum_{\gamma=1}^{6}H_{\gamma},$ (90)
$\displaystyle H_{1}$ $\displaystyle=\sum_{ij}t_{ij}c_{i}^{\dagger}c_{j},\quad
H_{2}=\sum_{ij}V_{ij}n_{i}n_{j},$ $\displaystyle H_{3}$
$\displaystyle=\sum_{\alpha
ij}g_{ij}^{(\alpha)}c_{i}^{\dagger}c_{j}X_{\alpha},\quad H_{4}=\sum_{\alpha
ij}h_{ij}^{(\alpha)}c_{i}^{\dagger}c_{j}P_{\alpha},$ $\displaystyle H_{5}$
$\displaystyle=\frac{1}{2}\sum_{\alpha}\omega_{\alpha}X_{\alpha}^{2},\quad
H_{6}=\frac{1}{2}\sum_{\alpha}\omega_{\alpha}P_{\alpha}^{2},$

where $X_{\alpha}=(b_{\alpha}+b_{\alpha}^{\dagger})/\sqrt{2}$ is the position
operator corresponding to the bosonic mode $\alpha$, and
$P_{\alpha}=i(b^{\dagger}_{\alpha}-b_{\alpha})/\sqrt{2}$ is the momentum
operator. Matrices $t=(t_{ij})$, $V=(V_{ij})$,
$g^{(\alpha)}=(g^{(\alpha)}_{ij})$, and $h^{(\alpha)}=(h_{ij}^{(\alpha)})$ are
all Hermitian. Such a Hamiltonian is in general not geometrically local, and
therefore the HHKL decomposition we have introduced in the previous section is
not directly applicable here. Instead, we consider performing quantum
simulation using Trotterization, which can be very efficient when many terms
in the Hamiltonian commute [19] and can be easier to realize in practice.

The Hamiltonians introduced above belong to a subclass of (14) as the fermion
interaction term has two indices instead of four. We impose this restriction
to ensure that this term can be fast-forwarded. We also require
$\sum_{ij}g_{ij}^{(\alpha)}c_{i}^{\dagger}c_{j}$ to commute with each other
for all $\alpha$, and require the same for
$\sum_{ij}h_{ij}^{(\alpha)}c_{i}^{\dagger}c_{j}$ for all $\alpha$. These
requirements are equivalent to

$[g^{(\alpha)},g^{(\alpha^{\prime})}]=[h^{(\alpha)},h^{(\alpha^{\prime})}]=0,$
(91)

for all bosonic modes $\alpha,\alpha^{\prime}$, where
$g^{(\alpha)}=(g_{ij}^{(\alpha)})$ and $h^{(\alpha)}=(h_{ij}^{(\alpha)})$.
Note in particular that this holds for the Hubbard-Holstein model, since
$(n_{x,\uparrow}+n_{x,\downarrow}-1)$ commutes with each other for all sites
$x$.

This decomposition of the Hamiltonian into $H_{\gamma}$’s and the above
requirements are to ensure that all the $H_{\gamma}$’s can be fast-forwarded.
More precisely, the evolution $e^{-itH_{\gamma}}$ due to the truncated
Hamiltonian terms can be implemented with gate complexity that depends poly-
logarithmically on $t$. The fermionic part $e^{-itH_{1}}$ can be implemented
using a series of Givens rotations [39]. The circuit implementation becomes
more challenging when a boson is involved. However, with a suitable particle
number truncation, we can still implement the time evolution of these terms.
We denote by $\widetilde{H}_{\gamma}$ the truncated Hamiltonian terms:

$\widetilde{H}_{\gamma}=\Pi_{[0,\widetilde{\Lambda}]}^{\mathrm{all}}H_{\gamma}\Pi_{[0,\widetilde{\Lambda}]}^{\mathrm{all}},$
(92)

for some $\widetilde{\Lambda}$ to be chosen. The coupling parts
$e^{-it\widetilde{H}_{2}}$ and $e^{-it\widetilde{H}_{3}}$ can be implemented
using the technique described in [51, Section IV], where the boson is
represented using the eigenbasis of $X_{\alpha}$, combined with the Givens
rotation technique. We can implement the coupling term associated with each