Time-dependent unbounded Hamiltonian simulation with vector norm scaling

The accuracy of quantum dynamics simulation is usually measured by the error of the unitary evolution operator in the operator norm, which in turn depends on certain norm of the Hamiltonian. For unbounded operators, after suitable discretization, the norm of the Hamiltonian can be very large, which significantly increases the simulation cost. However, the operator norm measures the worst-case error of the quantum simulation, while practical simulation concerns the error with respect to a given initial vector at hand. We demonstrate that under suitable assumptions of the Hamiltonian and the initial vector, if the error is measured in terms of the vector norm, the computational cost may not increase at all as the norm of the Hamiltonian increases using Trotter type methods. In this sense, our result outperforms all previous error bounds in the quantum simulation literature. Our result extends that of [Jahnke, Lubich, BIT Numer. Math. 2000] to the time-dependent setting. We also clarify the existence and the importance of commutator scalings of Trotter and generalized Trotter methods for time-dependent Hamiltonian simulations.


Introduction
Simulation of the quantum dynamics is widely viewed as one of the most important applications of a quantum computer. Let H(t) be a Hamiltonian defined on the interval [0, T ], and |ψ 0 be the initial vector, then the time-dependent Hamiltonian simulation problem aims to find |ψ(T ) , which solves the time-dependent Schrödinger equation If H(t) ≡ H is time-independent, then the solution can be expressed in closed form as |ψ(T ) = exp(−iT H) |ψ 0 . The task of creating such a universal quantum simulator was first conceptualized by Lloyd [34], and the past few years have witnessed significant progresses in the development of new quantum algorithms as well as the improvement of theoretical error bounds of existing quantum algorithms for time-independent Hamiltonian simulation [3,4,9,6,7,36,14,37,15,11,35,16,17,12,49]. In particular, for a d-sparse Hamiltonian with bounded H max (the largest element of H in absolute value), the complexity of the quantum signal processing (QSP) method [36] is O (T d H max + log(1/ )/ log log(1/ )), which matches complexity lower bounds in all parameters. Meanwhile the error bound of the high order Trotter-Suzuki scheme has also been significantly improved [17], which yields near-best asymptotic complexities for simulating problems such as k-local Hamiltonians, and improves previous error bounds for Hamiltonians with long range interactions. On the other hand, simulation with time-dependent Hamiltonians appears ubiquitously, such as in the context of quantum controls [57,38,42,21,43,39], non-adiabatic quantum dynamics [48,18], and adiabatic quantum computation [22,47,2], to name a few. Compared to the time-independent setting, there are significantly fewer quantum algorithms available [28,56,46,5,6,55,37,8]. The time-dependent setting has so far ruled out the usage of quantum signal processing [36] and quantum singular value transformation (QSVT) [24] types of techniques. To the extent of our knowledge, the best results available are given by the truncated Dyson series simulation [6,37], and the rescaled Dyson series method [8], which scales with respect to H max,∞ := sup t∈[0,T ] H(t) max , and H max,1 := T 0 H(t) max dt, respectively. In this paper, we are concerned with the simulation of a time-dependent and unbounded Hamiltonian H(t), which naturally includes the simulation of a time-independent Hamiltonian H(t) ≡ H as a special case. More precisely, we assume that there is a family of Hamiltonians H (n) (t) such that as n → ∞, the norm of H (e.g. the max of the operator norm or the L 1 norm) also increases towards infinity.
For concreteness, we will consider the bilinear quantum control Hamiltonian of the following form H (n) (t) = f 1 (t)H (n) Here H can be efficiently simulated. More specifically, if n also denotes the dimension of the Hamiltonian, we assume that the cost of the time-independent simulations depends at most poly-logarithmically in terms of n and the error . Such an assumption is standard for H has spectral norm asymptotically independent of n thus can be efficiently simulated, e.g. via the QSP technique [36]. However, the assumption on the effectiveness of simulating H (n) 1 is very strong especially when H (n) 1 grows polynomially in terms of n, and the no-fast-forwarding theorem [3,7] requires roughly Ω( H (n) 1 ) queries for generic quantum algorithms to simulate exp −iH (n) 1 . Nevertheless, for a subset of Hamiltonians with special structures, such a time-independent simulation can indeed be fast-forwarded and the query complexity is still poly-logarithmic of n. Typical examples include 1-sparse Hamiltonians [13,1,37] and thus unitarily diagonalized Hamiltonians where the diagonalization procedure can be efficiently implemented. We will show later the H (n) 1 of interest in this paper can also be fast-forwarded. The availability of the fast-forwarded time-independent Hamiltonian simulation allows us to measure the cost directly in terms of the number of Trotter steps.
When the context is clear, we will drop the superscript n and assume instead that H 1 is sufficiently large. In particular, we have H 1 H 2 . The form of Eq. (2) allows us to efficiently evaluate terms of the form t 2 t 1 H(t)dt = t 2 t 1 f 1 (t)dt H 1 + t 2 t 1 f 2 (t)dt H 2 , where the coefficients in the parentheses can be precomputed on classical computers when f 1 (t), f 2 (t) are available.
As an example, consider the following Schrödinger equation with a time-dependent effective mass M eff (t) (see e.g. [19,45,44,30,23,50]) in a domain D with proper boundary conditions as Here ω(t) is a frequency parameter. Then we set f 1 (t) = 1/(2M eff (t)), f 2 (t) = M eff (t)ω 2 (t)/2. When V (x) ≡ x 2 the system is a quantum harmonic oscillator with time-dependent effective mass. In general we assume the potential has suitable regularity conditions and is bounded on D. After proper spatial discretization using n degrees of freedom, H (n) 1 is the discretized negative Laplacian operator −∆ which is unbounded, and H (n) 2 is the discretized diagonal potential V (x) which is bounded. We notice that the simulation of H (n) 1 can be fast-forwarded since it can be diagonalized under the quantum Fourier transform procedure [41]. In order to demonstrate the behavior of the Trotter formulae for unbounded operators, we require n to grow polynomially with respect to −1 , where is the relative 2-norm error of the solution. This is the case, for instance, when the potential V (x) is of limited regularity. Throughout the paper we only require V (x) to be a C 4 function on the domain D. 1 Again for concreteness of discussion about the computational cost, unless otherwise specified, we will assume the system is one-dimensional, D = [0, 1] with the periodic boundary condition, and use the second order finite difference method with n equidistant nodes for spatial discretization.
To the extent of our knowledge, all previous results in the quantum simulation literature (for both time-independent and time-dependent Hamiltonians) measure the error of the evolution operator U (T ) − U (T ) , where U (T ) = exp T −i T 0 H(t)dt is the exact evolution operator expressed in terms of a time-ordered matrix exponential, and U (T ) is an approximate evolution operator obtained via the numerical scheme. We then directly obtain the vector norm error | ψ(T ) − |ψ(T ) ≤ U (T ) − U (T ) . However, since U (T ) − U (T ) typically depends polynomially on the operator norm H 1 , as H 1 increases, if the computational cost does not increase accordingly, then all error bounds of the operator norm U (T ) − U (T ) would increase to O(1), with the exception of the interaction picture method for time-independent Hamiltonian simulations [37]. 2 While the operator norm error provides an upper bound of the error given any initial vector, for a particular simulation instance, it is the vector norm error | ψ(T ) − |ψ(T ) that matters. It turns out that for certain unbounded operators and initial vectors, the vector norm bound can be significantly improved. The key reason is that the magnitude of terms such as H 1 |ψ , [H 1 , H 2 ] |ψ can be much smaller than the corresponding operator norm estimates. In fact the importance of the vector norm estimates has long been recognized in the numerical analysis literature, and the vector norm error bounds have been established for time-independent Hamiltonian simulation using second and higher order Trotter methods of the form H = −∆ + V (x) [29,53,20], and for time-dependent Hamiltonian simulation using Magnus integrators of the form H = −∆ + V (t, x) [27]. Under suitable discretization and choice of the initial vector, the vector norm error | ψ(T ) − |ψ(T ) obtained by the standard Trotter method remains small, even as H → ∞ and the operator norm U (T ) − U (T ) becomes O(1).

Contribution:
The first contribution of this work is to extend the vector norm estimate [29] to time-dependent unbounded Hamiltonian simulations. For concreteness we focus on the standard first and second-order Trotter methods, as well as a class of generalized Trotter methods proposed in [28], which will be introduced in Section 3. Our main result for a given control Hamiltonian Eq. (2) is Theorem 4. It states that under suitable assumptions, the vector norm error obtained from both standard and generalized Trotter methods depends mainly on sup t∈[0,T ] H 1 |ψ(t) , which can be significantly smaller than H 1 . 3 In order to simulate the Hamiltonian of the form Eq. (3), we take both the spatial and temporal discretization into account, and our complexity estimates are given in Theorem 5. Our result compared to existing results are given in Table 1, where the complexity for timeindependent simulations are obtained by treating M eff (t), ω(t) as constants. In particular, the vector norm is asymptotically independent of the spatial discretization parameter n, and complexity in terms of the error matches that of the time-independent Hamiltonian simulation obtained by [29]. Under the same second order spatial discretization, our complexity estimate for second order Trotter formulae outperforms state-of-the-art error bounds using high order Trotter and post-Trotter schemes [6,37,8] in terms of the desired level of accuracy, due to their dependence on the spectral norm of H 1 and thus on n.
The effectiveness of the vector norm bound depends on the initial vector |ψ 0 . We remark that recently [49] establishes improved error estimates of low-order Trotter methods for time-independent Hamiltonian simulation, when the initial vector is constrained to be within a low energy subspace. Another recent work [51] obtains an improved complexity estimate for simulating a system with η interacting electrons using time-independent Trotter formula, by considering the operator norm constrained on this η-electron sub-manifold. Our vector norm estimate provides a complementary perspective in understanding why such improved estimates are possible. When sup t∈[0,T ] H 1 |ψ(t) is indeed comparable to H 1 , the operator norm bound still serves as a good indicator of the error.
Given the improved error commutator scaling estimates for time-independent simulations [17], it is natural to ask whether the commutator scaling of the operator norm still holds for time-dependent simulations. The second contribution of this paper is to reveal that for time-dependent simulations, the error of standard Trotter method does not exhibit commutator scalings, while the commutator scaling holds for the generalized Trotter method (Theorem 2). Therefore in the context of time-dependent simulations, the use of the generalized Trotter method could reduce the simulation cost. Our proof of the operator norm error bounds mainly follow the procedure proposed in [17], and our results generalize the first and second order time-independent results in [17] in the sense that, when the scalar functions f 1 and f 2 are constant functions, both time-dependent standard Trotter formula and time-dependent generalized Trotter formula degenerate to the same time-independent Trotter formula, and the corresponding operator norm error bound is of commutator scaling.
Yet another twist comes when we ask the question  Section 6). Therefore the first-order generalized Trotter method outperforms the first-order standard Trotter method, but the asymptotic efficiency of the second-order generalized and standard Trotter methods are the same (Lemma 11). Table 2 summarizes the performance of Trotter and generalized Trotter methods. Though both second-order schemes share the same asymptotic scaling, the generalized Trotter formula may still be a better choice in practice due to smaller preconstants, which is observed numerically. Moreover, the p-th generalized Trotter scheme depends only on the (p − 1)-th derivatives of the control functions, while the p-th standard Trotter method on its p-th derivative. Therefore when the control functions have high frequency or limited regularity, the generalized Trotter scheme may significantly outperform the standard one. Such an advantage under first-order schemes has been demonstrated in [28]   (3) in one-dimension using second order Trotter method or higher order Trotter or post-Trotter method, with C 4 potential function V (x), and time-independent mass and frequency (top 2) or time-dependent mass and frequency (bottom 2). For all the methods, we use a second order finite difference discretization with n degrees of freedom. The third column summarizes the scaling of the cost with respect to n in order to reach constant target accuracy, and the fourth column summarizes the overall query complexity in order to achieve a desired level of relative 2-norm error . Since we assume the efficiency of time-independent simulation for both H 1 and H 2 , the query complexity is measured by the number of required  (3) in one-dimension with time-dependent effective mass and frequency. For all the methods, we use a second order finite difference discretization with n degrees of freedom. The third column summarizes the scaling of the cost with respect to n in order to reach constant target accuracy (Lemma 11), and the fourth column summarizes the overall number of required Trotter steps to achieve a desired level of relative 2-norm error estimated from error bounds in different norms (Theorem 5). The simulation time T is O(1).

Preliminaries
In this section we first introduce several notations and preliminary lemmas used in this paper. Then we briefly sketch the main ideas for proving the main theorems.

Notations
We refer to a (possibly unnormalized) vector as ψ, u or v depending on the context, and use |ψ to denote the corresponding quantum state (i.e. normalized vector under vector 2norm). We define two vector norms for a vector ψ = (ψ 0 , · · · , ψ n−1 ), namely the standard 2-norm The rescaled 2-norm is directly motivated by the discretization of the continuous L 2 norm [33,54]. Specifically, for a real-space function u(x) discretized in the real space using n equidistant nodes, we apply the trapezoidal rule and obtain Since the · simply rescales the standard vector 2-norm, the estimates that we will derive for 2-norm also hold for this rescaled 2-norm. Furthermore, the corresponding matrix norm induced by the rescaled 2-norm is still the standard matrix 2-norm without any rescaling factor, as We remark that it is equivalent to use either 2-norm or rescaled 2-norm if we wish to bound the relative error of the numerical solutions.
For two matrices A, B, define the adjoint mapping ad A as The following conjugation of matrix exponentials will be commonly used for a scale-valued function f and matrices A, B For a scalar-valued continuous function f (t) defined on time domain t ∈ [0, T ], we use f ∞ to denote the supremum of the function in this time interval, i.e.

Elementary lemmas
We review two elementary lemmas to be used in the proof of the paper. Proofs of the results can be found in, e.g. [26,32].
Lemma 2 (Variation of parameters formula). Assume U (t, s) solves the differential equa- Then 1. For any matrix-valued continuous function R(t), the solution of the differential equa- can be represented as 2. For any vector-valued continuous function r(t), the solution of the differential equation can be represented as

Main ideas
Here we discuss the main ideas for our operator and vector norm error bounds, and they are applicable to both standard and generalized Trotter formulae to be introduced in Section 3.1. For simplicity, we only discuss first order formulae here, and the ideas for second order formulae are similar. First, in Section 3, we derive the error representations between the exact evolution operator U (h, 0) and the Trotterized evolution operator U s (h, 0) or U g (h, 0), by establishing the differential equations that these unitary operators satisfy and by using variation of parameters. Such error representations are exact. Furthermore, although full error representations can be technically complicated, they are simply linear combinations of integrals with integrand of the form Here functions g j 's can be f 1  Next we focus on the situation when H 1 is very large, and we still would like to obtain a well approximated quantum state of the exact wavefunction. In this case, the operator norm error bounds do not offer useful performance guarantees. To obtain a vector norm error bound, the starting point of our approach is still the exact error representation. Notice that the error between the exact state |ψ and the approximate state | ψ obtained by Trotter formulae can be expressed as | ψ(h) − |ψ(h) ≤ U (h, 0) − U (h, 0) |ψ(0) . Therefore the vector norm error bounds should be a linear combination of the terms of the form (15) This can be obtained by applying the operator Eq. (14) to some vector v, and v is related to the initial condition as well as the exact Schrödinger wavefunctions.
To further bound Eq. (15), the key observation is as follows. When A is H 1 (or commutators involving H 1 ), although A can be very large, it is possible for A v to be small for certain vectors v. As an example, let us consider the continuous case in one dimension and we take H 1 = −∆ and H 2 being a bounded, smooth potential function V (x). The direct computation shows that Both of the terms on the right hand side depend on the spatial derivatives of the wavefunctions, of which the norm can be small if ψ is a smooth function. Then according to Eq. (15), we only need to somehow exchange the order between A and the exponentials without introducing much overhead (Lemma 7). Combining all previous arguments, we can obtain the desired vector norm error bounds.

Trotter type algorithms and error representations
In this section, we consider two different types of Trotter algorithms -the standard and generalized Trotter formulae -in simulating time dependent Hamiltonian Eq. (2), and derive their error representations explicitly. Here for simplicity we restrict ourselves to the first-order and second-order cases. We point out that such schemes and results regarding the (non-)existence of commutator scaling can be generalized to their higher order counterparts.

The standard and generalized Trotter formulae
Both the standard and the generalized Trotter formulae (proposed in [28]) belong to the class of splitting methods [25]. The first-order standard Trotter algorithm is The first-order generalized Trotter formula is The second-order standard Trotter formula is (18) The second-order generalized Trotter formula is It is clear that the difference between the standard and the generalized Trotter formulae lies in the temporal treatment of f 1 and f 2 , and the standard Trotter formula can be viewed as applying certain quadrature rules in representing the integrals of f 1 and f 2 . From now on we assume that b a f (s)ds can be accurately computed with negligible extra cost for any scalar-valued smooth function f (s). Furthermore, we remark that although in our definitions of the schemes we perform evolution governed by H 1 at first and then by H 2 , the order of H 1 and H 2 only affects the absolute preconstants in the error bounds and will not lead to any difference in the asymptotic scalings.

Error representations
For the time-independent Trotter formula, the work of [17,20] prove a commutator type of error of any order by writing down an explicit error representation via variation of parameters formula. Here we follow the procedure in [17] to write down the corresponding error representations for standard and generalized time-dependent Trotter formulae, which turn out to be the starting point for proving both the operator norm and the vector norm error bounds. Although we only present the error representation on the interval [0, h], this is just for notation simplicity and with minor modifications the results naturally hold on [t, t + h] for any t. The proofs are given in Appendix B.
Lemma 3 (Error representation of the first-order standard Trotter formula).
Lemma 4 (Error representation of the first-order generalized Trotter formula).
Lemma 5 (Error representation of the second-order standard Trotter formula).
Lemma 6 (Error representation of the second-order generalized Trotter formula).
The expressions of these exact error representations are somewhat complicated, but the structures for all the representations are the same. As introduced in Section 2.3, the error representations for both standard and generalized Trotter formulae of first and second-order are linear combinations of integrals with integrands in the form of Eq. (14), which can be expressed as the multiplication of matrix exponentials, Hamiltonians, and commutators of Hamiltonians.
Before we proceed, we remark that in the error representations of standard Trotter formulae of first and second-order, besides the O(h p+1 ) terms, we also include higher order terms O(h p+2 ) and/or O(h p+3 ). This is because we aim at writing down the exact error terms, and the Taylor expansion of terms like exp(−ihf 1 (h)H 1 ) will naturally involve higher order term even though we only expand it up to the desired lower order. For example, if we look at the first-order derivative of exp We can observe that the first term is O(1) and the second term is O(h), which are on different scales. Therefore, the same order term in the Taylor expansion of the unitaries does not necessarily have the same scaling in terms of h.

Remark 1 (Exact error representation). It is possible to derive a simpler error bound by considering only the lowest order term and discarding all the higher order terms in h.
However, such an error bound would not reveal the commutator scaling in the higher order remainder terms. For instance, for the time-independent Hamiltonian simulation, [53] deduces an error bound for the p-th order Trotter formula, in which the O(h p ) term has a commutator structure, but the higher order terms do not. This leads to complexity overhead when the spectral norms of the Hamiltonians become large. The work [17] fixes this issue by deriving an exact error representation, demonstrating the validity of the commutator scaling for high order terms as well. Therefore for time-dependent simulation, we also preserve all the terms in the error representation (at least for now). We can observe that all the terms in the exact representation, regardless of the order in h, are in the form of Eq. (14), thus no overhead will be introduced by higher order terms and it is safe to bound them by the lowest order term later in estimating complexity.

Operator norm error bounds
We first establish the operator norm error bounds. In this section we assume that H 1 and H 2 are two bounded operators. The operator norm error bounds can be directly obtained from the error representations by bounding the operator norms of all the unitaries by 1.

Theorem 1.
The error of each standard/generalized Trotter step measured in the operator norm is as follows: 1. First-order standard Trotter formula: where and 2. First-order generalized Trotter formula: where 3. Second-order standard Trotter formula: where α s,2 = 7 24 and 4. Second-order generalized Trotter formula: where Now we move on to the global error bounds. To obtain an approximation of the exact unitary evolution up to time T , we can divide the time interval [0, T ] into L equilength segments and implement Trotter discretization on each segment. Since the evolutions for both continuous and discretized cases are unitary, the global error is simply a linear accumulation of local errors at each time step. For sufficiently large L, the total error can be controlled to be arbitrarily small.

Theorem 2. Let T > 0 be the evolution time, and the dynamics Eq. (1) is discretized via standard and generalized Trotter formulae with L equidistant time steps (thus the time step size h = T /L). Then
where preconstants α and β are defined in Theorem 1.

Vector norm error bounds
Now we consider the case when H 1 is an approximation of an unbounded operator, while H 2 remains reasonably bounded. In this section, we assume that H 1 H 2 , and the functions f 1 , f 2 , as well as their first and second-order derivatives are bounded. To simplify the proof and emphasize our focus on overcoming the difficulty brought by H 1 , throughout this section we will not track the explicit dependence on H 2 , f 1 and f 2 . We use the notation C with a tilde above to denote preconstants (with possibly varying sizes in different inequalities) which can depend polynomially on ∞ but do not depend on H 1 . Furthermore, we make the following assumptions.

Assumption 1 (Bounds of commutators).
We assume H 1 is a positive semidefinite operator, and there exists an operator D 1 such that H 1 = D † 1 D 1 . Furthermore, we assume for any vector v, there exist constants C 1 , C 2 such that and Assumption 1 has been used in previous works [29,27] related to the vector norm error bounds of time-independent Trotter formula and exponential integrators for is a bounded operator. It will be helpful to understand Assumption 1 from a continuous analog, in which H 1 = D † 1 D 1 and D 1 is a first-order differential operator. A direct calculation shows that the operator is a first-order differential operator, and is a second-order differential operator, given that V (x) is a C 4 function. These are exactly what Eqs. (41) and (42) are addressing. Furthermore, if the wavefunction v is smooth enough with bounded derivatives, then the right hand sides of these inequalities are bounded, which provides the key motivation and possibility to establish vector norm error bounds and obtain improvement in complexity estimates. In the context of quantum simulation, since all matrices and vectors are finite dimensional, we omit the explicit statements of regularity assumptions of v below. As we have briefly discussed before, starting from the error representations, the vector norm error bounds are just linear combinations of the terms in the form of Eq. (15), and the key step to prove vector norm error bounds is to exchange the order of a Hamiltonian or an commutator with matrix exponentials. We find that such an exchange of order will not introduce any overhead in the error bounds. This is established by the following lemma.

Lemma 7.
Under Assumption 1, we have the following: 2. Let ξ be any real number such that ( C 1 + H 2 )|ξ| ≤ 1/2. Then for any vector v, 3. Let K be a positive integer, and H l k be either H 1 or H 2 , and ξ k be some real numbers , then for any vector v, all the following inequalities hold: for some constant C > 0 depending only on H 2 .
Proof. 1. By the definition of D 1 and the Cauchy-Schwarz inequality, 2. We start with Taylor's theorem of exp(itξH 2 ) up to first-order, then The norm can be estimated using Eq. (41) as Define M (t) := H 1 exp(itξH 2 ) v , and it follows from Eq. (44) that Thus Eq. (48) can be rewritten as Since H 1 v + 2 C 1 |ξ| v t is non-decreasing with respect to t, applying Gronwall's inequality yields the bound Finally taking t = 1 and applying the condition on ξ, the desired result is achieved 3. We first show that it suffices to only prove the first inequality in Eq. (46).
By Eqs. (41) and (47), and the fact w = v , we have Similarly Eq. (42) gives Furthermore, Therefore we only need to bound H 1 w further by C ( H 1 v + v ). (45), for H l being either H 1 or H 2 . Then we have the recursive relation By applying this estimation K times, we obtain the desired result Now we are ready to state our main theorems for the vector norm estimates.

Theorem 3. For any vector v and time step size
Proof. We start with the error representations (Lemmas 3 to 6). By multiplying a vector v on the right of the error representations, bounding all the f (k) j by its supremum, bounding all the higher order terms involving the (nested) commutators by second-order terms, and bounding all the unitaries multiplied on the left by 1, we obtain where θ s and θ g are linear combinations of the terms in the form of Eq. (46), with an exception that θ g,1 only consists terms like [H 1 , H 2 ] w . Then part 3 of Lemma 7 completes the proof.
Similar to the operator norm error bound case, from Theorem 3, we can establish the global error bound and estimate the total number of time steps we need to achieve a desired accuracy using standard and generalized Trotter formulae. The proof of the global error bound is essentially the same as the standard argument for error accumulation in quantum computing, that is, to replace the exact evolution operator by numerical evolution operator step by step and bound each step by local error bound. In the operator norm case, it does not matter whether we replace the local evolution operator in a forward or backward fashion. However, in the vector norm case, the order of the replacements indeed matters. In particular, we would like to obtain an error bound that depends on the exact, instead of the numerical solution of the dynamics. We state our vector norm global error bound in Theorem 4, and provide a complete proof for the second-order generalized Trotter formula.  Proof. Here we only present the proof for second-order generalized Trotter formula Eq. (58). The other three cases can be proved using the same approach.
For the second-order generalized Trotter formula, according to Theorem 3 and notice Now we compare the error bounds in terms of operator norm (Theorem 2) with those in terms of vector norm (Theorem 4). We notice that the scalings with respect to T and L are the same for schemes of the same order, and the difference is in the dependence of the preconstants on H 1 and H 2 . More precisely, the operator norm error bounds still depend on H 1 for standard Trotter formula and depend on norms of commutators like [H 1 , [H 1 , H 2 ]] for generalized Trotter formula. On the other hand, the vector norm error bounds only depend on H 2 , and the dependence on H 1 only appears in the form of H 1 ψ . Such a difference implies that the vector norm bounds can be much sharper than the operator norm bounds when H 1 is very large, but H 1 ψ and H 2 are relatively small. We will show later that this is indeed the case for the model of interest in Eq. (3). Furthermore, the difference in the error bounds can influence the scaling of total required Trotter steps with respect to the accuracy .

Application to Schrödinger equation with time-dependent effective mass and frequency
The model of the Schrödinger equation with a time-dependent effective mass and frequency in Eq. (3) has been studied in many works [19,45,44,30,23,50]. Our goal is to study the complexity to obtain an -approximation of the wavefunction at time T ∼ O(1), where D = [0, 1] with periodic boundary conditions. Throughout the section we make the following assumptions.
1. M eff (t) is positive function, and is uniformly bounded from below.

M eff (t)
, ω(t) are second-order continuously differentiable functions in t with uniformly bounded function and derivative values up to second order.

V (x) is a fourth-order continuously differentiable function in x with bounded function and derivative values up to fourth order.
Here the fourth order derivative of V (x) is required when estimating the errors of the second-order formulae in the operator norm. To be specific, it guarantees the nested commutator [H 1 , [H 1 , H 2 ]] in its spatial discretization with n spatial grids to have an operator norm bounded by n 2 (instead of n 4 ). Without going into details of the discretization, which will be presented in the proofs, here we provide an intuition of this requirement on the continuous level -the presence of the fourth derivative of V in [−∆, [−∆, V ]] as in Eq. (43). Since the control functions and potential are bounded, throughout we will not track explicitly the dependence on them and absorb them into the preconstant denoted by C or the big-O notation O in our estimates. We discretize the dynamics Eq. (3) as follows. First we perform spatial discretization using a central finite difference scheme with n equidistant nodes x k = k/n, 0 ≤ k ≤ n − 1.

Then the semi-discretized dynamics becomes i∂ t ψ(t) = H(t) ψ(t).
Here the k-th entry of ψ(t) will be an approximation of the exact wavefunction evaluated at t and x = (k − 1)/n.
The standard or generalized Trotter formulae are used to discretize the dynamics in time with equidistant time steps and obtain numerical approximation of the wavefunction. Furthermore, the H 1 and H 2 under central finite difference scheme using n equidistant nodes satisfy Assumption 1 with H 1 = D † 1 D 1 , and therefore the vector norm error bounds proved in Section 5 can be applied. This can be verified by straightforward but somewhat tedious matrix computations. We formally state the result in Lemma 8, and its proof is given in Appendix C.
j=0 . Therefore exp(−iH 1 ) = F exp(−iΛ)F † can be simulated efficiently, by first applying inverse quantum Fourier transform F † , then applying fast-forwarding techniques [13,1] for exp(−iΛ), and finally applying quantum Fourier transform F . For H 2 , it can be implemented via either QSP technique [36] due to the boundedness of H 2 , or fast-forwarding techniques since H 2 is a diagonal matrix as well. Due to the efficiency of simulating exp(−iH 1 ) and exp(−iH 2 ), it is reasonable to estimate the query complexity by counting the total number of Trotter steps.
Now we are ready to analyze the errors. We measure the discretization errors using rescaled 2-norm, i.e. ψ(t) − (φ(t, k/n)) n−1 k=0 , where ψ(t) is the numerical solution after spatial and time discretization at time t, and φ(t, x) denotes the exact solution. As discussed before, the reason why we use rescaled 2-norm, rather than regular 2-norm, to measure the error is because the exact solution (φ(t, k/n)) n−1 k=0 is a discrete representation of the function φ, which is normalized under continuous L 2 norm rather than discrete 2norm. Furthermore, if we encode the spatial discretized solution ψ into a quantum state, then the normalized condition requires |ψ ∼ 1 √ n ψ, thus |ψ ∼ 1 √ n ψ = ψ , that is, under correct normalization in each scenario, bounding regular 2-norm error for quantum states is equivalent to bounding rescaled 2-norm error for spatial discretized vectors. We remark that if we would like to control the relative error, then it does not matter whether the rescaled 2-norm or the 2-norm is used.
The errors are from two sources: spatial discretization of the Laplacian and potential operator, and the time discretization by Trotter formulae. We first bound the error from spatial discretization in Lemma 10, and its proof is given in Appendix D.

Let ψ(t) denote the solution of the dynamics
where H 1 and H 2 are the discretized Hamiltonian defined in Eq. (60) and Eq. (61), then for any 0 ≤ t ∈ T , This lemma shows that in order to make the vector error induced by the spatial discretization bounded by , it suffices to choose The second source of the error is the time discretization using standard or generalized Trotter formula by applying the aforementioned Theorems for the errors in the operator and vector norm. The following lemma makes explicit the scaling in n.

Lemma 11. Let ψ(t) be the solution of spatially discretized Schrödinger equation Eq. (3)
using finite difference with n grid points, and U (t, 0) be the evolution operator. Let ψ s,p (t) and ψ g,p (t) be the corresponding numerical solution from p-th order standard and generalized Trotter formula with L equidistant time steps, respectively, and U s,p (t, 0), U g,p (t, 0) be the corresponding evolution operators. Assume that L is sufficiently large, then there exists a constant C > 0, independent of T, L, H 1 , n, ψ, φ, such that 1.
2. 1] |φ(0, y)| , 1] |φ(0, y)| 1] |φ(0, y)| ,  , x) is the exact solution before any discretization. By Lemma 10, Finally, we combine both spatial and temporal errors. It is not possible to obtain an operator norm error bound between the evolution operator of an unbounded operator and that of a finite dimensional matrix. Hence the operator norm bounds below are obtained by plugging in the estimate of n that achieves the vector norm error with precision . In particular, the operator norm error bound involves the derivatives of the exact solution of interest φ. Combining Eq. (68) and Lemma 11, we obtain the total complexity estimates.

Theorem 5.
We use central finite difference for spatial discretization and Trotter formulae for time discretization to obtain an -approximation in rescaled 2-norm of the solution φ(t, x). Let L ope,s,1 and L ope,g,1 denote the total number of required time steps of first-order standard and generalized Trotter estimated from operator norm error bounds, respectively, L vec,s,1 and L vec,g,1 denote the estimates from vector norm error bounds, and L ope,s,2 , L ope,g,2 , L vec,s,2 , L vec,g,2 are the corresponding estimates for second-order schemes. Then for sufficiently small , 1] |φ(0, y)| , 1] |φ(0, y)| and Proof. The complexity can be estimated by requiring both error bounds in Lemma 10 and Lemma 11 to be smaller than . First, by requiring the right hand side of Eq. (67) to be bounded by , the scaling of n should be that in Eq. (68). Plug this back into Lemma 11 and also let the bounds in Lemma 11 to be bounded by , we obtain the complexity estimates.
Theorem 5 clearly illustrates the advantage of vector norm error bounds in terms of the desired level of error . More precisely, the total number of required Trotter steps estimated from vector norm bounds only scales O(1/ 1/p ) for p-th order schemes. This is because the operator norm error bounds depend on the spatial discretization n, where n = O(1/ 1/2 ) for second order spatial discretization, but the vector norm error bounds do not. We summarize our complexity estimates in terms of the spatial discretization as well as the error level in Table 2, where the simulation time T is O(1).
The best scaling is achieved by the second order standard and generalized Trotter formulae with the vector norm error bound, which is the result we are referring to as 'This work' in Table 1 for comparison with existing estimates. As discussed earlier, in order to demonstrate the behavior of the Trotter formulae for unbounded operators, we require n to grow as poly(1/ ). Therefore we choose V (x) to be a C 4 function so that the commutator scaling of the second order Trotter formulae are valid.
Numerical tests in Section 7 demonstrate that the complexity estimates from vector norm error bounds are sharp in terms of for all the schemes we consider.
Remark 5 (a priori estimates of the solution φ). Due to the potential growth of the derivatives of the exact solution with respect to T , a priori estimates are required if we would like to obtain the overall scalings in T . Such a priori estimates, where the spatial derivatives are bounded by polynomials of T , have been established in the literature for various special cases, such as when f 1 ≡ 1, f 2 is smooth in t and V is a real potential, smooth in x and periodic in x as considered in [10], and for strictly positive f 1 in [40]. The corresponding estimates are usually technical, while the common approach to derive them is a combination of various analytical tools and a careful capture on the resonance in the dynamics. Detailed discussions are orthogonal to the topic here and are omitted.
As we have already observed in Theorem 2, the generalized Trotter formula exhibits commutator type error bounds, while the standard Trotter formula does not. However, the commutator error bound only translates to improved asymptotic complexity for the first order generalized Trotter scheme. For second order schemes, there is no significant difference in the scaling with respect to between the standard and generalized Trotter formulae. As discussed before, this is due to the fact that H 1 and [H 1 , [H 1 , H 2 ]] have the same asymptotic scaling in n. The generalized Trotter is less restrictive on the control functions, namely, the p-th order generalized Trotter formula (p = 1, 2) only requires the boundedness of the derivatives of control functions up to the (p − 1)-th order while the p-th standard one requires the boundedness up to the p-th order. We expect the same situation for higher order Trotter formulae.

Numerical Results
To illustrate the difference between the operator norm and vector norm, we consider the following Hamiltonian with periodic boundary conditions. Here a > 0 controls the magnitude of the derivatives f 1 ∞ and f 1 ∞ . These sizes play a role in the preconstants of the errors as shown in Theorem 1. As discussed in Section 6, H 1 and H 2 correspond to the discretized matrices of −∆ and V (x), respectively. Besides the second order finite difference scheme, we also demonstrate that our estimates are equally applicable to the Fourier discretization. Though in this particular example V (x) is smooth, the scaling of n is still chosen according to Eq. (68), which only requires the regularity of V up to its fourth order derivatives and hence works for more general potentials. We first demonstrate the scaling of the vector norms and the operator norms, respectively. Consider the vector v as the discretization of the smooth function cos(x). Fig. 1 plots the operator norms and the vector norms for various number of spatial grids n using the finite difference and Fourier spatial discretization. We find that  (41) and (42), which is also proved for the finite difference scheme in Lemma 8.
We then verify the scaling of the errors with respect to n. The initial wavefunction is φ(x, 0) = cos(x). The time step size h is fixed to be 10 −4 . We run the Trotter formulae for 10 steps, which is sufficient for demonstrating the difference in scalings. The relative errors for both the operator and vector norms are plotted in Fig. 2 for a = 1 and a = 10. In terms of the operator norm, the generalized Trotter formula has a smaller error compared to the standard one: the relative error in the operator norms for the first-order standard Trotter scheme scales quadratically with respect to the number of grids while the first-order generalized Trotter schemes admits a linear scaling thanks to the commutator bounds. On the other hand, the relative errors in the vector norm do not grow with respect to n.
For second-order schemes, it can been seen that the errors measured by the operator norm for both methods grow quadratically with respect to n, while the corresponding errors in the vector norm are stable as n increases. These results agree with Lemma 11. Note that though the operator norm errors of the second-order schemes have the same asymptotic scaling in n, their preconstants may differ. When a = 1, the sizes of f 1 ∞ , f 1 ∞ , f 1 ∞ are comparable, and there is no significant difference in the preconstants. However, when a = 10, f 1 ∞ is one order of magnitude larger than f 1 ∞ , f 1 ∞ . In this case, we find from Fig. 2 that the generalized Trotter formula has a smaller preconstant, which agrees with the preconstant estimates as described in Theorem 1.

Conclusion
We have studied in detail the behavior of first and second order standard and generalized Trotter formulae for time-dependent Hamiltonian simulation with unbounded, control type Hamiltonians. We demonstrated that the error of the Hamiltonian simulation for a given initial state can often be overestimated using the standard analysis based on operator norms, which overestimates the computational cost. By taking into account the information of the initial state in the error analysis, sharper error estimates can be derived via the vector norm scaling. As a side product, we also obtained improved error bounds of the standard and generalized Trotter formulae in operator norm as well in the time-dependent setting.
As an example, we applied our results to the time-dependent Schrödinger equation with a time-dependent effective mass and frequency. While the complexities of existing quantum algorithms for time-dependent Hamiltonian simulation scale at least linearly in the spatial discretization parameter n, we demonstrate that, the error bounds in vector norm do not suffer from such overheads (for both the standard and generalized Trotter formulae). Thus in this setting, our results outperform all existing quantum algorithms, including higher order Trotter and post-Trotter methods.
The bilinear form in Eq. (2) facilitates the discussion of the error of the Trotter formulae. For the most general Hamiltonian H(t) = H 1 (t) + H 2 (t), it has been demonstrated that the error bound can be much more complicated even in the second order case [28]. Nevertheless, under suitable modifications, the main conclusion of this paper can still be applicable to more general time-dependent Hamiltonian under further assumption that ∂ k t H j (s) and ∂ k t H j (s ) commute for any j = 1, 2 and k, k , s, s (thus no essential difference is introduced in taking derivatives of unitaries and deriving error representation). This allows us to simulate e.g. Schrödinger equation with general time-dependent potential function V (x, t).
A natural extension of this work is to consider general high order time-dependent standard and generalized Trotter formulae defined by Suzuki recursion [52,56]. For the operator norm error bound, our results can be generalized to higher order case with a     control Hamiltonian Eq. (2). More specifically, let C k denote the set of the norms of all possible k-th order nested commutators of H 1 and H 2 , for example For p-th order schemes, we expect that the one-step operator norm error bounds for the standard and generalized Trotter formula scales as O(α s,p h p+1 ), O(α g,p h p+1 ), respectively. Here α s,p is a linear combination of terms in the set p k=0 C k , while α g,p is expressed as a linear combination of terms in the set p k=1 C k . Hence the difference lies in whether C 0 is included, and generalized Trotter formula allows a commutator scaling. Notice that such an error bound will improve the best existing estimate [56], which depends on the norms of the Hamiltonians as well as their high order derivatives, and does not demonstrate possible commutator scalings.
The extension of our vector norm error bounds to p-th order time-dependent Trotter formula is also possible. The corresponding assumption on the bounds of commutators (i.e. counterpart of Assumption 1 in this work) becomes for any 1 ≤ k ≤ p. Compared with the operator norm error bounds, for the Schrödinger equation with a time-dependent mass and frequency, such vector norm error bounds can still remove the dependence on the spatial discretization thus provide speedup in terms of the accuracy. However, the significance of such improvement might be subtle: in order to satisfy the assumption in Eq. (75), the potential function V (x) needs to be much smoother with bounded higher order derivatives. Hence, the dependence of n on the error may become much weaker by employing higher order discretization schemes. In such a scenario, the spectral norms H 1 and H 2 may even become comparable, and the Hamiltonian H(t) may not be regarded as an unbounded operator after all.
In this work, we mainly focus on the improvement brought by vector norm error bounds in terms of the accuracy. It is also interesting to study whether vector norm error bounds can improve the scalings of other parameters. For example, if the Schrödinger equation is in d dimension rather than one dimension considered in this paper, then a vector norm error bound may offer speedup in terms of d, since the degree of freedom for spatial discretization can scale linearly in d [31]. Another related topic is the scaling with respect to the number of the particles in quantum many-body systems. Recently [51] obtained an improved estimate in terms of the number of electrons for electronic structure problem with plane-wave basis functions in a second quantized formulation, by combining sparsity, commutator scalings and initial-state knowledge and bounding the operator norm on an η-electron sub-manifold. Although much smaller than that on the entire space, the operator norm on the η-electron sub-manifold may still overestimate the error, and a vector norm error bound might offer further improvement by taking the smoothness and low-energy property of the wavefunction into consideration. It is also an interesting question to investigate whether a vector norm error bound can provide any benefit for other applications such as spin systems.
Our work suggests that it may be of interest to explore the gap between the operator norm and vector norm error bounds in other schemes for Hamiltonian simulations with unbounded operators. Note that such a gap may not exist in all methods. For instance, for time-independent Hamiltonian simulation, the quantum signal processing (QSP) method [36] is based on the polynomial approximation to the function cos(xt) and sin(xt), and we do not expect that the error bound can be significantly improved by considering vector norms. However, it may be possible to prove vector norm error bounds for other post-Trotter methods.
A Derivation of results in Table 1 In this section we show explicitly how to derive the results in Table 1. Throughout this section we are considering the setup in Section 6 with T = O(1).
To obtain Table 1, we first restate all the complexity estimates for different methods proved in existing literature and show how they depend on as well as the scale of the Hamiltonians H 1 and H 2 . The dependence on H 1 naturally gives rise to the dependence on n, by noticing that as is discussed in Lemma 9. Then, under second order finite difference spatial discretization, Lemma 10 and Eq. (68) tell that n should be chosen as large as O( −1/2 ). Plugging this back into the complexity estimates leads to the overall scaling in terms of , as shown in the last column of Table 1.

A.1 Time-independent schemes
Time-independent second order Trotter formula [17,Proposition 16] gives an operator norm error bound for time-independent second order Trotter formula that the one-step local Trotter error is bounded by To bound this by , it suffices to choose [29, Theorem 3.2] provides a vector norm error bound for time-independent second order Trotter formula that the global Trotter error is bounded by We remark that [29] does not track explicitly the dependence on T . Noticing that H 1 v and v are independent of and n (shown in the proof of Lemma 11), the number of required Trotter steps scales Time-independent high order Trotter formula [17,Corollary 12] shows that for a p-th order time-independent Trotter formula, the number of required Trotter steps to obtain an -approximation of the exact evolution operator is Straightforward bounds for these p-th nested commutators are that which results in Notice that the scaling of is not improved by higher order Trotter formula. This is because such an estimate is made under the assumption that the potential V (x) is a C 4 function, therefore we only have better scaling for nested commutator up to second order.
In that case the complexity can be improved to L = O(n/ 1/p ), although there is still a linear dependence on n.
Truncated Taylor series [5,Theorem 1] shows that to obtain an -approximation of the exact evolution operator using truncated Taylor series, the query complexity is . We remark that the work [31] studies further the complexity of simulating time-independent many-body Hamiltonian and discusses carefully the errors from both time and space discretization. In this work, the authors use truncated Taylor series as well for time discretization, and use high order finite difference formula for spatial discretization. However, they only assume that the potential V (x) is first-order continuous differentiable thus the high order finite difference formula does not offer improved scaling of n To bound this by , it suffices to choose The second order complexity estimate from [56] is a special case of their general high order result. We will show the general case later. [ To bound this by , it suffices to choose [56,Theorem 1] proves that, to simulate a system with Hamiltonian H(t) = m j=1 H j (t) within operator spectral norm error using a 2k-th order standard Trotter formula, the total number of exponentials is

Time-dependent high order Trotter formula
We first notice that the total number of exponentials only differ from the total number of Trotter steps by a factor of 2m5 k−1 . After absorbing all the terms independent of n and into the big-O notation, in the case of the Schrödinger equation with a time-dependent effective mass, the total number of Trotter steps becomes Therefore the total number of Trotter steps becomes [37,Theorem 9] shows that to obtain an approximation of the exact evolution operator with success probability at least (1 − ) using truncated Dyson series method, the query complexity is

Truncated Dyson series
. Rescaled Dyson series [8,Theorem 10] shows that to obtain an approximation of the exact evolution operator using rescaled Dyson series method, the query complexity is Here d is the sparsity of the Hamiltonian, H max,1 = T 0 H(t) max dt where A max denotes the largest matrix element of A in absolute value. In the case of the model Eq. (3), noticing that H 1 max = O(n 2 ) because every non-zero entry of H 1 is either n 2 or (−2n 2 ), we have H max,1 = O(n 2 ). Therefore the query complexity becomes We mention that in [8] another method called continuous qDRIFT is also proposed to successfully achieve L 1 scaling of the Hamiltonian. However, continuous qDRIFT is a first order method, and its complexity dependence on H max,1 is quadratic, which is worse than that of rescaled Dyson series. Hence we only include the rescaled Dyson series method in our table for comparison.

B Proof of error representations
In this part, we derive the error representations of the first-order and second-order Trotter formulae, as presented in Lemma 3 -Lemma 6. All of the proofs consisting of the following two steps: One first compares the derivatives and its numerical analogs ∂ h U m,p (h, 0) (m = g, s and p = 1, 2), and apply the variation of parameter formula (Lemma 2); Then the Taylor theorem (Lemma 1) is applied to further simplify the terms. We first present the proof of Lemma 4, since its error representation contains fewest terms. The rest of the error representations, Lemma 3, Lemma 5 and Lemma 6, follow the exact same idea of proof, just involving more calculations.
Proof of Lemma 4. By taking derivative of U g,1 (h, 0) with respect to h, one has By applying Lemma 2 to Eq. (76) and Eq. (77), one obtains where E s,1 (h) is defined as Before proceeding, we first define the following quantities needed in the error repre-sentations of the second order standard and generalized Trotter formulae and Then We further split [H 1 , H 2 ] = D L + D R + S where and S = −diag(V (2) n−1 , V (2) 0 , V 1 , · · · , V (2) n−2 ). Notice that for any vector v = (v k ) n−1 k=0 , and similarly D R v 2 ≤ sup |V To

−2V
(2) n−1 and W = n 2 diag(V k−1 ) n−1 k=0 . Notice that H L , H R , H C are modifications from H 1 , with the center of the central formula to be on the diagonal or subdiagonal and multiplying each row by different bounded parameters. D DL and D DR are very similar to D L and D R , with higher order potential. Furthermore, we have and similarly For D DL and D DR , they can be bounded by the same way we bound D L and D R before, and thus and Finally since V has bounded fourth order derivative, the term n 2 (V Combining all the estimates together, we have  = −f 1 (t)n 2 (φ(t, x + 1/n) − 2φ(t, x) + φ(t, x − 1/n)) + f 2 (t)V (x)φ(t, x) + f 1 (t)r(t, x) where r(t, x) = n 2 (φ(t, x + 1/n) − 2φ(t, x) + φ(t, x − 1/n)) − ∆φ(t, x), the vector φ(t) = (φ(t, k/n)) n−1 k=0 satisfies the ordinary differential equation where R(t) = (r(t, k/n)) n−1 k=0 . Same as our previous notations, let U (t, s) denote the evolution operator from time s to t of the dynamics Eq. (1) with Hamiltonian Eq. (3). By the variation of parameters formula (Lemma 2), and thus It remains to bound R(s) .