Provably accurate simulation of gauge theories and bosonic systems

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 the truncation error, 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, 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 an error at most $\epsilon$ can be achieved by choosing $\Lambda$ to scale polylogarithmically with $\epsilon^{-1}$, an exponential improvement over previous bounds based on energy conservation. For the Hubbard-Holstein model, we numerically compute a bound on $\Lambda$ that achieves accuracy $\epsilon$, obtaining significantly improved 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 result on the truncation error in time evolution, we also prove that spectrally isolated energy eigenstates can be approximated with accuracy $\epsilon$ by truncating local quantum numbers at $\Lambda=\textrm{polylog}(\epsilon^{-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 condensedmatter systems with bosonic degrees of freedom [21,27,43,51,52,62,64,67,76], lattice gauge theories (LGTs) [4,5,6,7,13,15,20,38,41,42,45,53,56,60,61,69,73,75,77,78], 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 de-note it by Λ.
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 using a threshold Λ scaling polynomially with −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 Λ is far too pessimistic -a truncation threshold scaling as polylog( −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 φ 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 Λ on the precision also implies exponentially improved scaling of Λ 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 [40,51]. 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 Ω(T 2 ) 1 , in stark contrast to the 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 can be achieved with truncation threshold Λ = polylog( −1 ), in contrast with the more naive estimate Λ = poly( −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 dif- 1 For functions of real variables f, g, we write f = O(g) if there exist c, t0 > 0 such that |f (T )| ≤ c|g(T )| for all |T | ≥ t0. When there is no ambiguity, we will use f = O(g) to also represent that |f (τ )| ≤ c|g(τ )| holds for all τ ∈ R. We then extend the definition of O to functions of positive integers and multivariate functions. For example, we use f (N, T, 1/ ) = O((N T ) 2 / ) to mean that |f (N, T, 1/ )| ≤ c(N |T |) 2 / for some c, n0, t0, 0 > 0 and all |T | ≥ t0, 0 < < 0, and integers N ≥ n0. We write f = Ω(g) if g = O(f ) and we use O to suppress logarithmic factors in the asymptotic expression fer 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 [12, 35,54,65], 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 Λ; 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.

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   (1). Model parameters are from [43], with ω 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 Λ obtained in this work. See Section K.2 in the Appendix for details. 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 where H f is the Hamiltonian of the Fermi-Hubbard model [34] acting on only the fermionic modes, and (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,σ is the fermionic number operator for site x and spin σ.
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 Λ on the occupation number b † x b x ("number of particles") in each bosonic mode. Then each bosonic local Hilbert space has dimension Λ+1, and is spanned by the particle-number eigenstates {|λ : λ = 0, 1, 2, . . . , Λ}. This imposed upper limit results in a truncation. For a fixed site x, we define the projection operator Π (x) [0,Λ] = Λ λ=0 |λ λ| x ⊗ I, where I is the identity operator acting on the rest of the system. Imposing the upper limit truncates a quantum state |φ to be Π One may ask how large Λ should be for the resulting truncation error to be smaller than . 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 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 [−Λ, Λ]; hence the truncated local Hilbert space at each link is 2Λ + 1 dimensional. More precisely, for a fixed gauge link ν, we define the projection operator Π  Λ] . 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 Λ/2, then the link Hilbert space has dimension (Λ + 1)(Λ + 2)(2Λ + 3)/6, and is spanned by the angular momentum eigenstates |jmm where j is a half integer less than or equal to Λ/2 and −j ≤ m, m ≤ j. Again one may ask how large Λ 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 Λ 0 . After the state evolves forward for time T , how should the truncation threshold Λ be chosen for that site or link so that the resulting error is at most ? We show that suffices, where r = 1/2 for bosons and r = 0 for LGTs, and χ 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, −1 in (3) is replaced by N −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 † where b is the annihilation operator of a bosonic mode. Then e −iT H |0 , where |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 Ω(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 −1 . To truncate a system of O(N ) sites, we scale down by a factor of N , resulting in a threshold Λ scaling polynomially with N −1 . In contrast, our bound has only polylogarithmic dependence on N −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 Λ 0 at all sites or at all links, and suppose the state evolves forward for time T using a truncated Hamiltonian H rather than the ideal untruncated Hamiltonian H (we will define H later). How should the truncation threshold Λ be chosen so that the truncated evolved state matches the ideal evolved state up to time T with an error at most ? We show that suffices, where r = 1/2 for bosons and r = 0 for LGTs, N is the system size, and χ 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 [−Λ, Λ] 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 ν, we denote by Π (ν) S the projection operator imposing the condition that the local quantum number takes values from the set S. We also denote Π all S ; this is the projection operator imposing the condition on all bosonic modes or gauge links. For any projection operator Π, we write its complement as Π = I − Π. The truncated Hamiltonian mentioned in Question 2 is H = Π all S HΠ 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 ν, we may write the full Hamiltonian H of the model as where H (ν) W is the part of the Hamiltonian that can change the value of the local quantum num-ber at ν, and H (ν) R contains all the terms in the Hamiltonian that preserve the value of the local quantum number at ν. These two parts satisfy the conditions Here Π (ν) λ projects onto the eigenspace with local quantum number λ, χ and 0 ≤ r < 1 are parameters that depend on the model, and · is the spectral norm. (The notation Π (ν) [−Λ,Λ] 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 to change the local quantum number by at most ±1. The second condition requires that the rate at which the maximal local quantum number Λ changes is sublinear in Λ. The third condition requires H (ν) R 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 onsite terms. Choosing H (x) changes the local bosonic particle number by at most ±1, and that H (x) R preserves the local bosonic particle number. Moreover, using Π (x) λ to denote the projector onto the subspace with λ bosonic particles on site x, we see that is satisfied with χ = 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 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 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 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 [38,61,69]. 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.

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 Π [−Λ,Λ] |ψ(T ) is upper bounded by the leakage. Therefore, to ensure that the truncation er-ror is at most , we only need to keep the leakage below .
As mentioned before, we assume H has the structure (5) with H (ν) W and H (ν) 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 |ψ R preserves the local quantum number, |ψ I (t) and |ψ(t) induce the same local quantum number distribution. In the interaction picture, |ψ I (t) evolves with a time-dependent Hamiltonian H We then apply the Dyson series expansion to the unitary operator generated by H (ν) W (t). In the proof of Lemma 1 in the Appendix, we show that if 0 ≤ T ≤ 1/(2χ(Λ 0 + 1) r ), the truncated Dyson series with ∆ terms approximates the exact evolution up to an error e −Ω(∆) . Moreover, such a truncated Dyson series can change the local quantum number by at most ±(∆ − 1) due to (6a). Therefore we have the short-time leakage bound Using this short-time leakage bound, we can derive the long-time bound in (3). Specifically, for any choice of Λ 0 < Λ 1 < · · · < Λ J = Λ, 0 = T 0 < T 1 < · · · < T J = T , the total leakage is at most the sum of J short-time leakages (see Lemma 2 in the Appendix).

(8)
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 Λ j are chosen to satisfy 0 ≤ T j − T j−1 ≤ 1/(2χ(Λ j−1 + 1) r ). This establishes the scaling in (3) and provides an answer to Ques-tion 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 R satisfies (6) with parameters χ and r for a fixed mode or link ν. For any t ≥ 0 and integers Λ ≥ Λ 0 ≥ 0, 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 applies truncation on all sites or links. the truncation threshold Λ is chosen large enough so that evolution governed by H is a good approximation to the exact evolution. The approximation error is upper bounded by Therefore our goal is to choose Λ to ensure that this error is at most . This is accomplished by the following theorem which we establish in Section D in the Appendix and preview here.
We now briefly explain how we upper bound the Hamiltonian truncation error for sufficiently large Λ. We prove this by expanding the target quantity using the formula for Trotter error [72, Eq.

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 evolu-  (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 nearlinear 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 . 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 O(N T polylog(Λ 0 −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 [−Λ 0 , Λ 0 ]. Thus we achieve an almost linear dependence of the gate complexity on the spacetime volume N T . 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 H into time evolution of blocks. Each block, denoted by B, has size D = O(polylog(N T −1 )) and we only need to implement its evolution for time τ = O(1). There are O(N ) such blocks and the entire time evolution is divided into 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 poly(Λ). Fortunately, these terms with Λ-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 Λ-independent constant as [28] requires.
When simulating each block B we use the interaction picture Hamiltonian simulation technique suggested in [70] and developed in [50], and the gate complexity for simulation up to ). For Λ we use the scaling (4). There are in total O(N T ) such simulations that need to be performed, leading to a total gate complexity of O(N T polylog(Λ 0 −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 H which have Λ-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 Λ-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 B is polylogarithmic in Λ, and the cost of simulating evolution of N sites for time T is nearly linear in the spacetime volume N T .

Previous work on the quantum simulation of
LGTs such as [38,61,69] does not explain how to choose the truncation threshold Λ 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 O(N 3/2 T 5/2 ) from [38,69] and O(N 2 T 2 ) from [61].
Simulating bosonic systems and an Ω(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 B of the Hubbard-Holstein model, we cannot get a polylogarithmic dependence on Λ. Rather, the gate complexity to simulate a block is O( √ Λpolylog(ΛN T −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 O(N T ) blocks to be simulated, and the scaling of Λ is given by (4), the total gate complexity is 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 [17,19]. For the Hubbard-Holstein model we use the canonical commutation relation between the bosonic position and momentum operators [X α , P α ] = 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 com- 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 [17,28]. 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 Ω(N T 2 ) gates in general, for

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 ν, the Hamiltonian H can be expressed as in (5) and satisfies (6). Consider a nondegenerate eigenvalue ε of H, with corresponding eigenstate |Ψ , where ε is separated from the rest of the spectrum of H by a gap δ, and suppose that the expectation value of the absolute value of the local quantum number in the state |Ψ is ≤ . We show that this truncation threshold can be chosen to scale with , δ, andλ according to where χ 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 Λ on the truncation error arises because the distribution of local quantum numbers in the eigenstate |Ψ 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 approx-imate eigenstate projection operator [30] When σ δ and T σ −1 , this operator is close to the eigenstate projection operator P ε = |Ψ Ψ|. We derive (12) by applying the approximate projector P ε 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 Π [−2λ,2λ] |Ψ 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 δ may be exponentially small in the system size.

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 Λ. In particular, we showed that for a fixed evolution time T , a precision can be achieved by choosing Λ scaling polylogarithmically with −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 O(N T polylog(Λ 0 −1 )), where N is the system size, assuming that the initial state can be well approximated with truncation threshold Λ 0 . For a bosonic system like the Hubbard-Holstein model, our algorithm has gate complexity O(N T ( √ Λ 0 + T )polylog( −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 using a truncation threshold polylogartihmic in −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 [1,30], 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 φ 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 φ 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 φ 4 theory and other mod-els involving bosons with anharmonic couplings. However, it has the disadvantage that the truncation threshold scales polynomially rather than polylogarithmically with −1 . Under suitable conditions, can the truncation threshold scale as polylog( −1 ) in a broader class of models than those satisfying (5) and (6), and are there models in which polylog( −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 energybased 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. for helpful discussions. YT

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 α. The c i and b α denote the fermionic and bosonic annihilation operators respectively. The Hamiltonian takes the form where 2 is the position operator corresponding to the bosonic mode α, and ij ) 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 φ x . Each link consists of a planar rotor, whose configuration space of states |θ , with θ ∈ [0, 2π] 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 for k ∈ Z.
In Hilbert space of link (x, n i ) we define operators E x,n i and U x,n i through Then the Hamiltonian of the system is , where P denotes a summation over all plaquettes P . For P whose lower-left site is x, U P is defined as The trace 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 φ l x , l = 1, 2. We write φ x = (φ 1 x , φ 2 x ) . 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 |jmm , where j, m, m are simultaneously either integers or half-integers with −j ≤ m, m ≤ j. Here j is the rotator's total angular momentum, and m, m 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 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 Because φ x has two components, where each component is a fermionic mode, U x,n i is a 2 × 2 matrix, where each of the 4 matrix entries is an operator acting on the link Hilbert space An important property that we will use later is which follows from rules for the addition of angular momentum, given that U x,n i transforms as the j = 1/2 representation of SU (2). Here O denotes the spectral norm of an operator O. We also note that relative to the basis {|jmm }, U ll x,n i 's are sparse matrices because 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 − 3/2 are the change of angular momentum as a result of applying U ll x,n i . Eqs. (21) and (22) imply that the matrix representing U ll x,n i has at most three non-zero elements in each row and column.

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 H into a direct sum of subspaces V λ with quantum numbers λ ∈ Z. The projection operator onto each subspace V λ is denoted by Π λ . Then λ∈Z Π λ = I. We consider a class of Hamiltonians of the form where Model λ r Boson-fermion Bosonic particle number 1/2 U (1) LGT Electric field value 0 SU (2) LGT Total angular momentum 0 for some χ > 0, 0 ≤ r < 1. Here Π [−Λ,Λ] = |λ|≤Λ Π λ . In such a Hamiltonian, H W changes the quantum number λ in the time evolution while H R preserves it. (24) ensures that λ is not changed too quickly. The first part of (24) ensures that the local quantum number is changed by at most ±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.
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 α 0 , we can decompose the Hilbert space according to the number of particles in the bosonic mode α 0 , which we denote by m. This means we let λ = m and where |m α 0 means the |m α 0 -particle state of the mode, and I is the identity operator acting on the rest of the system. We set Π λ = 0 for all λ < 0. We define H W to be the sum of terms in (24) that change the particle number in mode α 0 : whereas the rest of the terms in H are collected into H R . Because of the fact that where Π [0,M ] = M m=0 Π m , one can see that (24) is satisfied if we choose r = 1/2 and where |A| = √ A † 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 2 α 0 and P 2 α 0 each on their own will result in r = 1, but since in the Hamiltonian they appear together, X 2 α 0 + P 2 α 0 preserves the local quantum number and therefore does not contribute to how large r is.
Tr|A| is the trace norm of A and it is used here because ij A ij c † i c j ≤ 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 ∈ {n 1 , n 2 }. Then we decompose the Hilbert space by the electric field value on this link, i.e. we let λ = k and define Π k = |k k| (x 0 ,n 0 ) ⊗ I.
Then H W should be chosen as Because of the fact that H W ≤ 4|g B | + 2|g GM |, (24) is satisfied if we choose χ = 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 λ = 2j (j takes half-integer value), and Here we require m, m 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 χ = 16|λ B | + 8|g GM | and r = 0. There is an additional factor of 4 in χ compared to the U(1) case because there are now four operators U ll x 0 ,n 0 contributing to the growth of the quantum number instead of one. More generally, we define Π S , where S is a set of integers, as 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 Π S for any integer set S. When we need to constrain the particle number on all modes, we define As a general rule, if Π is any projection operator, we define Π = I − Π. For lattice gauge theories we adopt similar notations. For example we use Π (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 ν to index both the bosonic mode and gauge links when we discuss the two scenarios together. Therefore Π

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 λ between ±Λ 0 , what is the probability that |λ| grows beyond some given Λ as the state evolves for time t? To be more concrete, we want to bound 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 −Λ 0 and Λ 0 .

C.1 The short-time leakage bound
We first establish a bound on the leakage defined in (34) for a short time t: 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 Then writing the time evolution e −itH in the interaction picture, we have where T is the time-ordering operator. Since e −itH R commutes with Π (−Λ 0 −∆,Λ 0 +∆) , we only need to bound To this end, we consider the partial sum of the Dyson series of T e −i t 0 H W (s)ds : 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 Noting that H W (t) can only change λ by ±1, we have which implies by (24) that By repeatedly applying (40), we have Then applying (41), From the above inequality, we have where in the second inequality we have used the fact that and in the third inequality we have used t ≤ 1/(2χ(Λ 0 + 1) r ). Note that because of (40). Therefore, (47) 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 , P j be projection operators such that P j + P j = I, j = 0, 1, . . . , J. Then for any 0 = t 0 < t 1 < · · · < t J = t, which implies 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 and expand the product on the right-hand side into a sum of terms, each of which is a string of P j and P j interspersed with e −i(t j +1 −t j )H . We then recombine these terms according to where the first P j appears, or if it does not show up at all. The sum of all terms for which the first P j appears in the j-th place is We then sum over j, and multiply P J and P 0 to the left and right respectively, to get the right-hand side of (48). We now state our long-time leakage bound:  for any t ∈ R and integers Λ 0 ≥ 0 and ∆ > 1, we have 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 Λ(t) needed to ensure the leakage for the boson-fermion coupling setting and the lattice gauge theory setting. We can see that Λ(t) grows quadratically with time for the former and linearly for the latter. Moreover, very small leakage can be achieved by only slightly increasing Λ(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 ≥ 0, the intermediate times where we make the partition. First we define the instantaneous quantum number and then choose T j to be From this definition we have which can be proved using the inequality To establish the long-time leakage bound for arbitrary time t, we first prove the following lemma: Proof. We choose t j = T j for j = 0, 1, . . . , J − 1, and t J = t. By Eq. (49) in Lemma 2, we have 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 ≥ t. By (55), we have The claimed bound then follows from Lemma 4.

Theorem 5. Let H = H W + H R be a Hamiltonian satisfying (24) with parameters χ and r. For any t ≥ 0 and integers
Proof. In Theorem 3 we choose ∆ so that where in the third line we have used the inequality that a p + b ≤ (a + b) p when a ≥ 0, p ≥ 1 and b (to be chosen as ∆ − 1) is a non-negative integer. Using the fact that 2 1−∆ (∆!) −1/2 = e −Ω(∆) , the claim follows immediately from Theorem 3.
If we want to ensure that truncating at a threshold Λ has an error Π [−Λ,Λ] e −itH Π [−Λ 0 ,Λ 0 ] ≤ , then by Theorem 5 we can choose This is the scaling given in (3) of the main text.

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 H such that 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 N of them in the system. To simplify the discussion, we use ν to index either bosonic modes or gauge links, replacing the indices α and (x, n). Therefore we have Note also that all projection operators Π We will establish the following bound. .

Theorem 6 (Hamiltonian truncation). Let H be a Hamiltonian such that H = H
(64) Then for any t ∈ R, where 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, A(Λ) can be bounded by a polynomial of the Hamiltonian coefficients and Λ. This is because and in all examples we discussed in Section A, the norm of HΠ all [−Λ,Λ] is bounded by a function that is linear in all coefficients, and linear or quadratic in Λ for boson-fermion coupling and lattice gauge theories respectively. 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 where coef includes all the coefficients in the models (17) which follow immediately from the fact that for each ν, H can only change the quantum number λ by ±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 ν. This is done through a union bound, as given in the following lemma: λ be projections all commuting with each other. For any operator A and set S ⊂ Z, we have where Π all λ commute with each other, they can be simultaneously diagonalized, and by the union bound we have Π which in turn leads to where the first step can be proven using the singular value decomposition, and we have used the fact S A are all positive semidefinite for the later steps.
Note that the above union bound actually holds even when the commutativity assumption about Π We then use this, along with invariance of the spectral norm under multiplication by a unitary operator, to bound the truncation error as: Now if we choose Λ ≥ Λ 0 + 1, then by (70), 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 for some Λ to be chosen. We choose Λ = Λ − 2. With this choice and (71) we have 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 where we have used the fact that s 2 ≤ t. Substituting this bound into (78) and then (76), we have In the above derivation we used the fact that 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 On the right-hand side, [ H, H] grows polynomially with the truncation threshold Λ, while decays subexponentially with Λ . 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 Λ , 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 energybased 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.

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 H defined in (64) as opposed to the unbounded H, with the truncation threshold Λ 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 . 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 O(N T polylog(Λ 0 N T −1 )), where is the allowed error, assuming that the initial state is in the span of states whose quantum number is in the range [−Λ 0 , Λ 0 ] for each gauge link.
As mentioned above, we will be simulating H instead of H, and the resulting error has been analyzed in Section D. We use H E to denote the truncated electric field part of the Hamiltonian, i.e.
, and we adopt similar notation H M , H GM , and H B for the other three parts.
Note that the Hamiltonians for lattice gauge theories, both the original H and 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 B, has size D = O(polylog(N T −1 )) and we only need to simulate its evolution for time τ = O(1). The entire time evolution is divided into O(T ) segments and there are 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 H have norm depending on the truncation threshold Λ, which scales with the system size N , time T , and allowed error . We address this issue as follows. The only local terms that are not bounded by a constant are the electric field terms in H E , i.e. g E 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 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 B of size D for short time τ . The Hamiltonian for B, which we denote by H B , includes all the local terms in H that only act on sites and links within the block 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 [24,49]) of the Hamiltonian H B , i.e. a unitary with H B as a matrix block, with a subnormalization factor O( D Λ 2 ). Using the Hamiltonian simulation algorithm for encoded Hamiltonians [49], we can then simulate the time evolution of the block B with gate-complexity O( 2D Λ 2 τ ) = O( Λ 2 polylog(N T −1 )), which results in a total gate complexity of O(N T (Λ 0 +T ) 2 polylog( −1 )). This is however not the best method in terms of asymptotic complexity.
The polynomial dependence on Λ can be improved to be poly-logarithmic using the interactionpicture simulation technique developed in [50]. We group the local terms in B into H B M , H B GM , H B E , and H B 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 Λ comes only from H B E . We note that the time evolution under H B E 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 E 2 x,n for time t can be implemented with gate complexity O(polylog( Λ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 H B E we only need to evolve these terms separately, and thus pay a cost of O( D polylog( Λt)) = O(polylog(N T Λt −1 )).
Now instead of directly simulating the Hamiltonian H B , we simulate The original time evolution e −it H B and the interaction picture evolution are related through It then suffices to implement T e −i t 0 ds H B I (s) , namely the time evolution due to the time-dependent Hamiltonian H 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 H B E discussed above. This yields a gate complexity O( 2D τ polylog(N T Λ −1 )) = O(polylog(N T Λ −1 )) for implementing the interaction picture time evolution and consequently e −it H 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 O(N −1 T −1 ) instead of . 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 O(N T ) such simulations to perform for all the O(N ) blocks and O(T ) times steps. Therefore the total gate complexity for implementing the time evolution of the entire system is O(N T polylog(N T Λ −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 is O(N T polylog(Λ 0 N T −1 )), provided that the initial state is in the support of Π all [−Λ 0 ,Λ 0 ] , i.e. the quantum numbers on each gauge link are in the interval 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,σ (with σ 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 where H f is the Hamiltonian of the Fermi-Hubbard model: and are the boson-fermion coupling part and bosonic part respectively. The lattice sites are indexed by x and x , and spins are indexed by σ. 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 with gate complexity O(N T ( √ Λ 0 + T )polylog( −1 )), assuming the initial state has no more than Λ 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 H in (64) and simulate the evolution of H. The resulting error is analyzed in Section D. We also denote different parts of the Hamiltonian after truncation by H f , H f b , and H We apply the HHKL method to decompose the entire time evolution into evolution of blocks, each of which is denoted by B, for short time τ . Here again we encounter local terms whose norms are not bounded by a constant, and in this case these terms are contained in H f b and H b . With the help of Lemma 13 however, we can still apply the HHKL decomposition because H f b and 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 B. We denote by H B f , H B f b , and H B b the fermionic, coupling, and bosonic terms respectively. Here, the terms in H B 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 H B f b can be fast-forwarded. Therefore when we simulate the interaction picture Hamiltonian the dependence on the truncation threshold is not poly-logarithmic. Rather a factor of Λ shows up in the subnormalization factor of the encoding of the Hamiltonian because 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 B for time τ is then O( 2D Λτ polylog( ΛN T −1 )) = O( Λpolylog( ΛN T −1 )) with a target accuracy of O(N −1 T −1 ) for each block. We need to perform O(N T ) such simulations and the number of required gates is therefore O(N T Λpolylog( −1 )). The truncation threshold Λ 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 is O(N T ( √ Λ 0 + T )polylog( −1 )), assuming the initial state has at most Λ 0 particles in each bosonic mode.

F.1 Sources of error
There are two sources of error that we need to deal with. The first source of error comes from the fact that we are evolving the system with H instead of H, and this is already analyzed in Theorem 6. The second source of error is the Trotter error, which will be our focus here. A simple bound for the Trotter error is readily available if we ignore the commutation relation between pairs of the Hamiltonian terms. But here we aim for the commutator scaling described in [19], which can be much tighter when many terms in the Hamiltonian commute.
There is a technical issue that prevents us from directly applying the result of [19]. After truncation, the commutation relation between the projected position and momentum operators is different from the canonical commutation relation between the exact position and momentum operators. To address this, we use the fact that the exact commutation relation is recovered when the particle number is some distance below the truncation threshold Λ, and this in turn requires carefully tracking the particle number under the exact and truncated time evolution respectively. Our proof uses the following telescoping lemma: Lemma 9. Let Π be a projection operator and U j , U j (j = 1, 2, . . . , J) be unitary operators. We have Proof. This inequality follows immediately from the identity which can be proved by induction on J.
Our goal is to simulate the dynamics up to time T . We achieve this by dividing the entire time evolution into R = T /τ steps, each of which has duration τ and is simulated by a p-th order product formula S(τ ). Then the Trotter error can be bounded as where Π By Theorem 5 and Lemma 8 we have where χ is given in (28). We now only need to bound (S(τ ) − e −iτ H )Π all [0,Λ 0 ] .

F.2 Trotter error with bounded particle number
The main result of [19, Theorem 3] is a bound on the Trotter error in terms of the spectral norm of nested commutators of Hamiltonian terms. That bound does not take into account the fact that the initial state has a finite number of particles and is thus not suitable for our purpose. Instead, we use an exact representation of Trotter error, which is provided in Theorems 3 and 5 of [19].
In [19, Theorem 3] they derive the following expression for the Trotter error: and by [19, Theorem 5], T (τ 1 ) can be written as where is a product of operator exponentials. In the above equations H γ ν ∈ { H 1 , H 2 , . . . , H 5 }, γ = (γ 1 , γ 2 , . . . , γ p+1 ) is a string of indices for Hamiltonian terms, and L γ , C γ q , c γ ν , and R γ are constants that only depend on the Trotter formula but not on the Hamiltonian or time variable τ . Also C γ q is non-zero only when γ p+1 = γ p = · · · = γ p−q+2 , but this property does not affect the asymptotic gate complexity and will thus be ignored in the subsequent analysis.
With this exact representation of Trotter error, we have where we have used the decomposition for some Λ 1 to be chosen later. Since τ = O(1) (τ should be chosen to be much smaller than constant to suppress Trotter error) and R γ ν=0 |c γ ν | is a constant, we have from Theorem 5 that Here U γ (τ 1 , τ 2 ) involves a constant number of operator exponentials each of which is generated by a term H γ from the Hamiltonian. We note that Theorem 5 applies to each operator exponential because H γ also has the structure described in Section B. Thus we can apply Theorem 5 a constant number of times to arrive at inequality (102). We now follow [19] to define There is a Λ dependence because the truncated Hamiltonian depends on the threshold Λ. Furthermore, we define Then by (100) and (102) Combining this with (94) and (96), we bound the total error from the Trotter decomposition as where we have used the relation T = Rτ . We now choose Λ 0 , Λ 1 , Λ, and R so that the right-hand side of the above inequality is at most , while simultaneously keeping the truncation error from Theorem 6 below .
There is one other constraint in our choice of parameters: we need to ensure the canonical commutation relation of X α and P α , when replaced by X α and P α , holds exactly when evaluating β comm ( Λ, Λ 1 ). Note that in β comm ( Λ, Λ 1 ) there are at most 2(p + 1) truncated position and momentum operators multiplied together because each Hamiltonian term is at most quadratic in these operators. Thus, if then we can simply treat X α and P α as if they satisfy the exact canonical commutation relation when evaluating β comm ( Λ, Λ 1 ). We recall that Λ is the particle number truncation threshold for the Hamiltonian. To be more specific, this means when (107) is satisfied. With this extra constraint and (106), we choose where coef denotes all the coefficients t, V, g, h, ω in the Hamiltonian H, and C is defined in (28). This choice of Λ will also ensure that the Hamiltonian truncation error is upper bounded by by Theorem 6.
In choosing these parameters, we have omitted the scaling with α comm ( Λ). This is because α comm ( Λ) is upper bounded by a polynomial of Λ and the Hamiltonian coefficients, and it gets absorbed into the poly-logarithmc factors.

F.3 Bounding the nested commutators
It now remains to bound β comm ( Λ, Λ 1 ). Suppose we are given a series of indices of Hamiltonian terms γ 1 , γ 2 , . . .. We will show that where At a high level, A γ 's quantify the growth of the nested commutator when the nesting layer increases by one, while B γ 's are chosen to handle the base case when there is only one operator.
Once we have established (110), we define This implies that we need Trotter steps by (109). In the above analysis we treat the order p as a constant. The gate complexity depends on how we implement each e −it Hγ . For concreteness, we analyze the gate complexity of simulating the Hubbard-Holstein model in the next section, although the approach may be extended to simulate other quantum systems within our framework. We now derive the bound (110) for an arbitrary nested commutator. We first note that a nested commutator multiplied to a projection operator [ H γq , · · · [ H γ 2 , H γ 1 ] · · · ]Π all [0,Λ 1 ] can be written as a linear combination of products of at most q fermionic operators c † i c j , and at most q projected bosonic position or momentum operators, multiplied to the projection operator at the end. This can be proved inductively. We introduce some notations to formalize this observation. For convenience we denote We first define a set of index strings: where i and j are strings of fermionic mode indices, α is a string of bosonic mode indices, and ς is a string of 0's and 1's. Then the claimed expansion is formally given by where we have used the canonical commutation relation between X α and P α on X α and P α . This is justified because the nested commutator is multiplied to the projection operator Π all [0,Λ 1 ] and we have imposed the constraint (107). The sum of the absolute value of the coefficients is at most q max α ij |g (α) ij | and the the number of bosonic operators in the product is reduced by 1. Therefore the contribution to S q+1 is at most q max α ij |g (α) ij |(2(Λ 1 + 1)) −1/2 S q . Combining our analysis for the second and third lines of (121) we have if γ q+1 = 3. The commutators with the other H γ 's can be analyzed in a similar way. The proof of (110) is now completed.

F.4 Simulating the Hubbard-Holstein model with Trotterization
We recall the definition of the Hubbard Holstein model given in Section E.2: The Hamiltonian is where H f is the Hamiltonian of the Fermi-Hubbard model: and are the boson-fermion coupling part and bosonic part respectively. The lattice sites are indexed by x and x , and spins are indexed by σ. As in Section E.2, we assume for simplicity that all model parameters except for the system size N , i.e. g, ω 0 , U , µ, are all constants. We consider the case where the time evolution starts with an initial state that has at most Λ 0 bosonic particles at each site. We note that this Hamiltonian satisfies the general form of boson-fermion coupling Hamiltonians given in (90). Therefore we can directly apply our above analysis to analyze the number of required Trotter steps. First we note that all the quantities involved in A given in (112), i.e.
are upper bounded by some constants. This follows from the sparsity of the coefficient matrices t, V , Similarly, all the quantities involved in B given in (112), are all O(N ). Therefore we have Then by (114) the number of Trotter steps required to simulate the Hubbard-Holstein model is Note that Λ 1 has the asymptotic scaling described in (109). Taking into account the fact that for the Hubbard-Holstein model Tr|g (α) |, Tr|h (α) | = O(1), which implies further that χ = O(1), we have which gives Each Trotter step can be implemented with O(N polylog(Λ 0 T −1 )) gates, and therefore the total gate complexity is For large p, this almost matches the gate complexity derived in Section E.2 based on the HHKL decomposition.

G A gate complexity lower bound for simulating bosons
In Sections E.2 and F, we have discussed the gate complexity of simulating the Hubbard-Holstein model. One distinctive feature is that the scaling with respect to time is almost quadratic, instead of being almost linear for simulating bounded Hamiltonians. In this section we construct a class of Hamiltonians acting on a single bosonic mode and a register of qubits, for which performing simulation up to time T will require at least Ω(T 2 ) gates. This shows that simulation involving bosons cannot in general be expected to have a linear dependence on time. Note that here by simulation we mean simulating only the qubit part of the boson-qubit coupled system, and as a result we only need to deal with a finite-dimensional Hilbert space. Specifically we consider the time evolution of a bosonic mode coupled to a register containing N qubits. We will label the bosonic mode by a subscript β and the qubit register by a subscript q. A product state is written as |ψ β |φ q , where |ψ β is the state of the bosonic mode, and |φ q is the state of the qubit register. We call this qubit register the q-register because later we need an additional qubit register. When simulating the time evolution of this system, we consider a unitary circuit W acting jointly on an ancilla register, which we label as anc, and the q-register. We will also denote by |λ β the ε-particle state of the bosonic mode.
Theorem 10. For any integers N and T such that 1 ≤ √ N ≤ T ≤ 2 N/2 , there exists a boson-qubit coupled Hamiltonian H = U b+b † U † , where b and b † are the bosonic annihilation and creation operators respectively, and U is a unitary acting on the bosonic mode and N qubits (the q-register) that preserves the bosonic number. If a quantum circuit W satisfies for all |φ q , then W must use at least Ω(N T 2 ) 2-qubit gates. Here, I β and I anc are the identity operator on the bosonic mode and the ancilla register respectively, and O = |0 0| ⊗ I is the projection onto the |0 state of the first qubit of the q-register. The quantum circuit W may use an arbitrarily large number of ancilla qubits, and gates in W may be non-local and come from a possibly infinite gate set.
In essence, this theorem asserts the existence of boson-qubit coupled systems whose single-qubit measurement statistics after evolving for time T require Ω(N T 2 ) gates to approximate to constant precision.
To prove Theorem 10 we need to use the following lemma: Then Therefore b and b † can be treated as a new pair of annihilation and creation operators. By the Kermac-McCrae identity we have This can then be used to prove the lemma by using the Taylor expansion and the fact that Proof of Theorem 10. First we consider a quantum circuit U circ that acts on N qubits and has depth T 2 . It can then be written as where each U λ acts on N qubits and has depth one. We also define U λ = I for all λ ≥ T 2 . Then we let the unitary U in the theorem be Note that by construction we have [U, |λ β λ| β ] = 0 and therefore U preserves the particle number in the bosonic mode. We will show that by running time evolution e −iT H for T = Θ(T ) starting from |0 β |φ q , and performing measurement on the first qubit in the q-register, we will be able to approximately sample from the distribution generated by running U circ and then measuring the first qubit (note that U circ acts only on register q). In this procedure we trace out the bosonic mode and focus only on the qubits.
By Lemma 11, we have Now note that for any summand on the right-hand side with λ ≥ T 2 , we have As a result we can write where |Ψ ⊥ βq is the sum of the first T 2 terms on the right-hand side of (141), and Note that the normalization factor A is chosen so that |ψ β = 1. In the above quantum state e −iT H |0 β |φ q , the bosonic particle number satisfies the Poisson distribution with mean T 2 . Because the Poisson distribution decays rapidly away from the mean [14], we can choose T = Θ(T ) so that and consequently where in going from the second line to the third line we have used the fact that |Ψ ⊥ βq + A ≤ 1. Therefore where O = |0 0| ⊗ I. If a circuit W as described in the theorem satisfies the inequality (135), then by the triangle inequality This means the measurement outcome generated by running the circuit U circ can be simulated by running the circuit W .
With the above setup, we then use U circ to compute Boolean functions in the sense defined in [28]: for a Boolean function f : {0, 1} N → {0, 1}, we say U computes the Boolean function with high probability if measuring the first qubit of U |x 1 x 2 · · · x N 0 · · · 0 yields f (x) with probability at least 2/3. We also say U computes the Boolean function exactly if measuring the first qubit of U |x 1 x 2 · · · x N 0 · · · 0 yields f (x) with probability 1.
By (148), we know that if U circ computes a Boolean function f exactly, then W computes the same Boolean function with high probability. By [28, Lemma 8], we can choose U circ acting on N qubits and with depth T 2 to compute 2 Ω(T 2 N ) different Boolean functions exactly. If W uses G 2-qubit gates, then by [28, Lemma 8] different W can compute at most 2 O(G log(N )) different Boolean functions with high probability. Therefore G = Ω(T 2 N ), which completes the proof.

H Quantum number distribution tail bound in eigenstates
If we would like to prepare an eigenstate of a Hamiltonian of the form (23) on a quantum computer, then we need to be able to store this eigenstate using a finite number of qubits. This reaffirms the need to truncate infinite dimensional Hilbert spaces. A natural approach is to truncate the local quantum number λ, which, as discussed in Section B, is the local bosonic particle number in the setting of boson-fermion coupling, the electric field value in the setting of U(1) lattice gauge theory, and the total angular momentum in the setting of SU(2) lattice gauge theory.
In this section, we will show that the probability of a spectrally isolated eigenstate having a local quantum number beyond a certain threshold can be bounded, and we call this the tail bound. This tail bound justifies cutting off the part of the Hilbert space with local quantum number beyond the threshold, thus enabling us to store eigenstates using a finite number of qubits. We describe the result in the following theorem: Theorem 12 (Quantum number distribution tail bound). Let H = H W + H R be a Hamiltonian satisfying (24) with parameters χ and r. Assume that |Ψ is an eigenstate of H corresponding to an eigenvalue ε with multiplicity 1, and that ε is separated from the rest of the spectrum of H by a spectral gap δ. Moreover, we assume the absolute value of the quantum number distribution has a finite expectation λ |λ| Ψ|Π λ |Ψ =λ < ∞.
Proof. We define the projection operator into the ε-eigensubspace by This projection operator, and its approximator to be introduced later, will be the main technical tool in this proof. We first apply a projection operator to truncate the eigenstate |Ψ : We then apply an approximation of P ε to |ζ . Note that P ε |ζ is exactly the eigenstate |Ψ up to a constant factor. Therefore applying an approximation of P ε will yield a quantum state that is close to the eigenstate. The approximation of P ε is constructed as We will show that P ε is close to P ε when σ is small and T is large. First we have where we have used the identity For the first term on the second line of (155), we have and for the second term we have where we have used [16, Theorem 1] for the second inequality. Denoting the sum of these two bounds by 1 , we have We choose σ and T so that 1 ≤ 1/2 √ 2. By applying the approximate projection operator we obtain a quantum state | Ψ : where | Ψ is a normalized quantum state and β = P ε |ζ > 0. We have β − P ε |ζ = P ε |ζ − P ε |ζ ≤ P ε − P ε ≤ 1 , and as a result β ≥ P ε |ζ − 1 = | Ψ|ζ | − 1 ≥ 1/2 √ 2.

I Lieb-Robinson velocity with on-site interaction
In this section we show that the Lieb-Robinson velocity is unaffected by any on-site interaction. This fact has been proved in [58,Section 2], although their result is not completely in line with what is required in this work. Therefore we provide our own theorem and proof in this section. We use the notation in [28,Lemma 5]. We consider a lattice Λ, with dist denoting the lattice distance. A where q denotes the bosonic momentum, b q is the bosonic annihilation operator,r el and P el are electron position and momentum operators respectively. In (14) we wrote the fermionic part in second quantization and the bosonic part in first quantization. We therefore rewrite (182) accordingly. We use c kσ to denote the annihilation operator for an electron with momentum k and spin σ. Then, For the bosonic part we have b q = 1 √ 2 (X q + iP q ), and therefore we can rewrite the Hamiltonian as We thus see that the Hamiltonian is of the form (14), and therefore has the structure described in Section B. The number of particles in each bosonic mode under time evolution and energy eigenstates can be analyzed using the results in Sections C, D, and H. We also observe that the ab initio Hamiltonians describing electron-phonon coupling [25], if no anharmonic terms are included, can also be analyzd within our framework due to its similarity to the Fröhlich Hamiltonian.
Besides boson-fermion coupling, spin-boson coupling can also be analyzed within the framework of this work. As an example, we consider the Dicke model which describes light-matter interaction [22,32].
The Dicke model. The model Hamiltonian can be written as where σ x i and σ z i are the Pauli-X and Z matrices respectively acting on site i, and b is the annihilation operator for a bosonic mode corresponding to photons. We note that this Hamiltonian has the structure described in Section B. We choose Then H R preserves the bosonic particle number, H W changes the bosonic particle number by ±1, and where Π [0,Λ] is the projection operator into the subspace with at most Λ bosonic particles. Therefore (24) is satisfied if we choose χ = 2g √ N and r = 1/2. We thus see that this model can also be analyzed within our framework.

K Comparison with the energy-based truncation threshold
In Ref. [36], to simulate the φ 4 theory, a truncation threshold is chosen for the field value at each lattice site based on energy conservation and Chebyshev's inequality. This is a very general method and can be applied to the systems studied in this work. Here we compare the truncation threshold obtained using that method with the one in this work in two settings. In the first setting we consider a single bosonic mode, and in the second we consider the Hubbard-Holstein model consisting of N sites. We find that the truncation threshold in this work tends to be lower than the energy-based one if the truncation is made for short-time evolution of large systems with high precision.
where E f,0 is the ground state energy of H f . Therefore Now we assume all parameters in the model, except for ω 0 , are constants, and we only consider the scaling with respect to the system size. Consequently, |E f,0 | = O(N ), which implies Therefore we can bound the particle number expectation value on site x : Again we denote the quantum state at time t by |ψ(t) . The projection operator into the subspace with at most Λ particles in bosonic mode x is denoted by Π for all x . Using (201) and Markov's inequality, we thus need to choose the truncation threshold Λ to scale as We compare this energy-based truncation threshold with the one derived in this work in (68), which for the Hubbard-Holstein model is We see that besides the advantage mentioned in the single mode setting there is also an exponentially better scaling with respect to the system size.
In Figure 3 we compare the truncation threshold Λ computed using the method of this work and the energy-based method of [36] for the Holstein model, which is a special case of the Hubbard-Holstein model with U = 0, with parameters chosen according to [43]. We assume the initial state is a tensor product between the fermionic ground state and a quantum state of the bosonic modes that has at most Λ 0 = 4 particles in each mode. We clearly see that when the system size becomes larger or when the precision requirement is higher, our method yields a lower truncation threshold than the energy-based method.  The curves below show the Λ obtained in this work and the horizontal lines above show the Λ obtained using the energy-based method. The model parameters are chosen according to [43]: ω 0 = 1, g = 0.5, U = 0, µ = 0. We assume the initial state has at most Λ 0 = 4 particles in each bosonic mode.
In the upper left panel we set = 10 −2 and in the upper right panel we set N = 100. The panel below shows the cross-over of the two truncation thresholds using the two methods for N = 5 and = 0.1. The truncation thresholds in this work are computed by using (52) and choosing the smallest integer ∆ to satisfy the precision requirement.