# Time-dependent unbounded Hamiltonian simulation with vector norm scaling

Dong An  Department of Mathematics, University of California, Berkeley, CA
94720, USA  Di Fang  Department of Mathematics, University of California,
Berkeley, CA 94720, USA  Lin Lin  Department of Mathematics, University of
California, Berkeley, CA 94720, USA  Computational Research Division, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, USA  Challenge Institute for
Quantum Computation, University of California, Berkeley, CA 94720, USA

###### Abstract

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.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
    1. 2.1 Notations
    2. 2.2 Elementary lemmas
    3. 2.3 Main ideas
  3. 3 Trotter type algorithms and error representations
    1. 3.1 The standard and generalized Trotter formulae
    2. 3.2 Error representations
  4. 4 Operator norm error bounds
  5. 5 Vector norm error bounds
  6. 6 Application to Schrödinger equation with time-dependent effective mass and frequency
  7. 7 Numerical Results
  8. 8 Conclusion
  9. A Derivation of results in Table 1
    1. A.1 Time-independent schemes
    2. A.2 Time-dependent schemes
  10. B Proof of error representations
  11. C Proof of Lemma 8
  12. D Proof of Lemma 10

##  1 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 $\ket{\psi_{0}}$ be the initial vector,
then the time-dependent Hamiltonian simulation problem aims to find
$\ket{\psi(T)}$, which solves the time-dependent Schrödinger equation

$\mathrm{i}\partial_{t}\ket{\psi(t)}=H(t)\ket{\psi(t)},\quad\ket{\psi(0)}=\ket{\psi_{0}}.$
(1)

If $H(t)\equiv H$ is time-independent, then the solution can be expressed in
closed form as $\ket{\psi(T)}=\exp(-\mathrm{i}TH)\ket{\psi_{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 $\left\lVert
H\right\rVert_{\max}$ (the largest element of $H$ in absolute value), the
complexity of the quantum signal processing (QSP) method [36] is
$\mathcal{O}\left(Td\left\lVert
H\right\rVert_{\max}+\log(1/\epsilon)/\log\log(1/\epsilon)\right)$, 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,\infty}:=\sup_{t\in[0,T]}\|H(t)\|_{\max}$, and
$\|H\|_{\max,1}:=\int_{0}^{T}\left\lVert H(t)\right\rVert_{\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)\equiv H$ as a special case. More precisely,
we assume that there is a family of Hamiltonians $H^{(n)}(t)$ such that as
$n\to\infty$, 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)}_{1}+f_{2}(t)H^{(n)}_{2}.$ (2)

Here $H^{(n)}_{1}$ and $H^{(n)}_{2}$ are time-independent Hamiltonians, and
$f_{1}$ and $f_{2}$ are two bounded, smooth scalar functions on a time
interval $[0,T]$. Without loss of generality we assume that
$\lim_{n\to\infty}\left\lVert H_{1}^{(n)}\right\rVert=\infty$, while
$\lim_{n\to\infty}\left\lVert H_{2}^{(n)}\right\rVert<\infty$, i.e.
$H_{1}^{(n)}$ approaches an unbounded operator, while the limit of
$H_{2}^{(n)}$ is a bounded operator. We also assume that
$\exp\left(-\mathrm{i}H_{1}^{(n)}\right)$ and
$\exp\left(-\mathrm{i}H_{2}^{(n)}\right)$ 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 $\epsilon$. Such an assumption
is standard for $H_{2}^{(n)}$ because $H_{2}^{(n)}$ 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_{1}^{(n)}$ is very strong especially when $\|H_{1}^{(n)}\|$
grows polynomially in terms of $n$, and the no-fast-forwarding theorem [3, 7]
requires roughly $\Omega(\|H_{1}^{(n)}\|)$ queries for generic quantum
algorithms to simulate $\exp\left(-\mathrm{i}H_{1}^{(n)}\right)$.
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_{1}^{(n)}$ 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 $\left\lVert H_{1}\right\rVert$ is sufficiently large. In particular, we
have $\left\lVert H_{1}\right\rVert\gg\left\lVert H_{2}\right\rVert$. The form
of Eq. 2 allows us to efficiently evaluate terms of the form
$\int_{t_{1}}^{t_{2}}H(t)dt=\left(\int_{t_{1}}^{t_{2}}f_{1}(t)dt\right)H_{1}+\left(\int_{t_{1}}^{t_{2}}f_{2}(t)dt\right)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_{\text{eff}}(t)$ (see e.g. [19, 45, 44, 30, 23,
50]) in a domain $D$ with proper boundary conditions as

$H(t)=-\frac{1}{2M_{\text{eff}}(t)}\Delta+\frac{1}{2}M_{\text{eff}}(t)\omega^{2}(t)V(x),\quad
x\in D.$ (3)

Here $\omega(t)$ is a frequency parameter. Then we set
$f_{1}(t)=1/(2M_{\text{eff}}(t)),f_{2}(t)=M_{\text{eff}}(t)\omega^{2}(t)/2$.
When $V(x)\equiv 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_{1}^{(n)}$ is the discretized
negative Laplacian operator $-\Delta$ which is unbounded, and $H_{2}^{(n)}$ is
the discretized diagonal potential $V(x)$ which is bounded. We notice that the
simulation of $H_{1}^{(n)}$ 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 $\epsilon^{-1}$, where $\epsilon$ 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$.111Here the $C^{4}$
regularity is a technical assumption to bound the norm of the nested
commutators, which will be detailed in Section 6. 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
$\left\lVert\widetilde{U}(T)-U(T)\right\rVert$, where
$U(T)=\exp_{\mathcal{T}}\left(-\mathrm{i}\int_{0}^{T}H(t)dt\right)$ is the
exact evolution operator expressed in terms of a time-ordered matrix
exponential, and $\widetilde{U}(T)$ is an approximate evolution operator
obtained via the numerical scheme. We then directly obtain the vector norm
error
$\left\lVert\ket{\widetilde{\psi}(T)}-\ket{\psi(T)}\right\rVert\leq\left\lVert\widetilde{U}(T)-U(T)\right\rVert$.
However, since $\left\lVert\widetilde{U}(T)-U(T)\right\rVert$ typically
depends polynomially on the operator norm $\|H_{1}\|$, as $\left\lVert
H_{1}\right\rVert$ increases, if the computational cost does not increase
accordingly, then all error bounds of the operator norm
$\left\lVert\widetilde{U}(T)-U(T)\right\rVert$ would increase to
$\mathcal{O}(1)$, with the exception of the interaction picture method for
time-independent Hamiltonian simulations [37].222In the context of time-
independent simulation [17], the error of Trotter methods does not scale
directly with respect to the operator norm $\left\lVert H_{1}\right\rVert$,
but with respect to the norm of the (high-order) commutators
$[H_{1},H_{2}],[H_{1},[H_{1},H_{2}]],[H_{2},[H_{2},H_{1}]]$ and so on.
However, this does not change our conclusion here. In principle, the
interaction picture method can also be generalized to efficiently simulate
time-dependent Hamiltonians. However, its practical performance has not been
well understood.  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
$\left\lVert\ket{\widetilde{\psi}(T)}-\ket{\psi(T)}\right\rVert$ 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 $\left\lVert
H_{1}\ket{\psi}\right\rVert,\left\lVert[H_{1},H_{2}]\ket{\psi}\right\rVert$
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=-\Delta+V(x)$ [29, 53, 20], and
for time-dependent Hamiltonian simulation using Magnus integrators of the form
$H=-\Delta+V(t,x)$ [27]. Under suitable discretization and choice of the
initial vector, the vector norm error
$\left\lVert\ket{\widetilde{\psi}(T)}-\ket{\psi(T)}\right\rVert$ obtained by
the standard Trotter method remains small, even as $\left\lVert
H\right\rVert\to\infty$ and the operator norm
$\left\lVert\widetilde{U}(T)-U(T)\right\rVert$ becomes $\mathcal{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\in[0,T]}\|H_{1}\ket{\psi(t)}\|$, which can be significantly smaller
than $\left\lVert H_{1}\right\rVert$.333Theorem 4 shows that the number of
Trotter steps may not scale with respect to $\|H_{1}\|$. The mechanism of the
improvement is very different from that of the interaction picture approach,
where the number of the time steps is still linear in $\|H_{1}\|$.

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 time-independent simulations are obtained by
treating $M_{\text{eff}}(t),\omega(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
$\ket{\psi_{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 $\eta$ interacting electrons using time-independent
Trotter formula, by considering the operator norm constrained on this
$\eta$-electron sub-manifold. Our vector norm estimate provides a
complementary perspective in understanding why such improved estimates are
possible. When $\sup_{t\in[0,T]}\|H_{1}\ket{\psi(t)}\|$ is indeed comparable
to $\left\lVert H_{1}\right\rVert$, 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: when $H_{1}$ is unbounded,
is it clear that the norm of the commutators
$\left\lVert[H_{1},H_{2}]\right\rVert,\left\lVert[H_{1},[H_{1},H_{2}]]\right\rVert,\left\lVert[H_{2},[H_{2},H_{1}]]\right\rVert$
must be smaller than $\left\lVert H_{1}\right\rVert$? It turns out that for
the Hamiltonian Eq. 3, we may directly analyze that
$\left\lVert[H_{1},H_{2}]\right\rVert,\left\lVert[H_{2},[H_{2},H_{1}]]\right\rVert\in\mathcal{O}(n)$,
while $\left\lVert[H_{1},[H_{1},H_{2}]]\right\rVert,\left\lVert
H_{1}\right\rVert\in\mathcal{O}(n^{2})$ (see 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] as well.

All results above are confirmed by numerical experiments for the model Eq. 3
in Section 7, which verifies the sharpness of our estimates.

|  Work/Method |  Scaling w. spatial discretization |  Overall query complexity
---|---|---|---
Time-independent second order Trotter |  Childs et al. [17] |  $\mathcal{O}(n)$ |  $\mathcal{O}(\epsilon^{-1})$
Jahnke et al. [29] |  $\mathcal{O}(1)$ |  $\mathcal{O}(\epsilon^{-0.5})$
Time-independent higher order methods |  $p$-th order Trotter [17] |  $\mathcal{O}(n^{2-2/p})$ |  $\mathcal{O}(\epsilon^{-1})$
Truncated Taylor series [5, 31] |  $\widetilde{\mathcal{O}}(n^{2})$ |  $\widetilde{\mathcal{O}}(\epsilon^{-1})$
Quantum signal processing [36] |  ${\mathcal{O}}(n^{2})$ |  ${\mathcal{O}}(\epsilon^{-1})$
Interaction picture [37] |  $\mathcal{O}(\log(n))$ |  $\mathcal{O}(\operatorname{polylog}(1/\epsilon))$
Time-dependent second order Trotter |  Huyghebaert et al. [28] |  $\mathcal{O}(n)$ |  $\mathcal{O}(\epsilon^{-1})$
Wiebe et al. [56] |  $\mathcal{O}(n^{3})$ |  $\mathcal{O}(\epsilon^{-2})$
Wecker et al. [55] |  $\mathcal{O}(n)$ |  $\mathcal{O}(\epsilon^{-1})$
|  This work |  $\mathcal{O}(1)$ |  $\mathcal{O}(\epsilon^{-0.5})$
Time-dependent higher order methods |  $p$-th order Trotter [56] |  $\mathcal{O}(n^{2+2/p})$ |  $\mathcal{O}(\epsilon^{-1-2/p})$
Truncated Dyson series [6, 37] |  $\widetilde{\mathcal{O}}(n^{2})$ |  $\widetilde{\mathcal{O}}(\epsilon^{-1})$
|  Rescaled Dyson series [8] |  $\widetilde{\mathcal{O}}(n^{2})$ |  $\widetilde{\mathcal{O}}(\epsilon^{-1})$
Table 1: Comparison of complexity estimates for simulating the model Eq. 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 $\epsilon$. 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 Trotter steps for Trotter-type methods, or the query complexity under standard query model for post-Trotter methods. The simulation time $T$ is $\mathcal{O}(1)$. ‘This work’ refers to the vector norm error bound using the second order standard or generalized Trotter formula. Throughout the paper $f=\widetilde{\mathcal{O}}(g)$ if $f=\mathcal{O}(g\operatorname{polylog}(g))$. See Appendix A for details of the derivation of the scalings.  |  Method & Error type |  Scaling w. spatial discretization |  Overall number of Trotter steps
---|---|---|---
First-order Trotter |  standard, operator norm |  $\mathcal{O}(n^{2})$ |  $\mathcal{O}(\epsilon^{-2})$
generalized, operator norm |  $\mathcal{O}(n)$ |  $\mathcal{O}(\epsilon^{-1.5})$
standard, vector norm |  $\mathcal{O}(1)$ |  $\mathcal{O}(\epsilon^{-1})$
|  generalized, vector norm |  $\mathcal{O}(1)$ |  $\mathcal{O}(\epsilon^{-1})$
Second-order Trotter |  standard, operator norm |  $\mathcal{O}(n)$ |  $\mathcal{O}(\epsilon^{-1})$
generalized, operator norm |  $\mathcal{O}(n)$ |  $\mathcal{O}(\epsilon^{-1})$
standard, vector norm |  $\mathcal{O}(1)$ |  $\mathcal{O}(\epsilon^{-0.5})$
|  generalized, vector norm |  $\mathcal{O}(1)$ |  $\mathcal{O}(\epsilon^{-0.5})$
Table 2: Summary of results for first and second order Trotter formulae
applied to simulating the model Eq. 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 $\epsilon$ estimated from error bounds in different
norms (Theorem 5). The simulation time $T$ is $\mathcal{O}(1)$.

##  2 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.

###  2.1 Notations

We refer to a (possibly unnormalized) vector as $\vec{\psi},\vec{u}$ or
$\vec{v}$ depending on the context, and use $\ket{\psi}$ to denote the
corresponding quantum state (_i.e._ normalized vector under vector 2-norm). We
define two vector norms for a vector
$\vec{\psi}=(\psi_{0},\cdots,\psi_{n-1})$, namely the standard 2-norm

$\|\vec{\psi}\|=\sqrt{\sum_{k=0}^{n-1}|\psi_{k}|^{2}},$

and the rescaled 2-norm

$\|\vec{\psi}\|_{\star}=\frac{1}{\sqrt{n}}\|\vec{\psi}\|.$

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

$\int_{0}^{1}|u(x)|^{2}dx\approx\sum_{k=0}^{n-1}\left(|u(k/n)|^{2}\frac{1}{n}\right)=\frac{1}{n}\|(u(k/n))_{k=0}^{n-1}\|^{2}=\|(u(k/n))_{k=0}^{n-1}\|_{\star}^{2}.$
(4)

Since the $\|\cdot\|_{\star}$ 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

$\|A\|=\sup_{\|\vec{u}\|\neq
0}\|A\vec{u}\|/\|\vec{u}\|=\sup_{\|\vec{u}\|_{\star}\neq
0}\|A\vec{u}\|_{\star}/\|\vec{u}\|_{\star}.$

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 $\operatorname{ad}_{A}$ as

$\operatorname{ad}_{A}(B)=[A,B]=AB-BA,$ (5)

and then the conjugation of matrix exponentials of the form $\exp(A)B\exp(-A)$
can be simply expressed as

$\exp(\operatorname{ad}_{A})B=\exp(A)B\exp(-A).$ (6)

The following conjugation of matrix exponentials will be commonly used for a
scale-valued function $f$ and matrices $A,B$

$\exp\left(\operatorname{ad}_{\mathrm{i}\int_{t_{1}}^{t_{2}}f(s)dsA}\right)B=\exp\left(\mathrm{i}\int_{t_{1}}^{t_{2}}f(s)dsA\right)B\exp\left(-\mathrm{i}\int_{t_{1}}^{t_{2}}f(s)dsA\right).$
(7)

For a scalar-valued continuous function $f(t)$ defined on time domain
$t\in[0,T]$, we use $\|f\|_{\infty}$ to denote the supremum of the function in
this time interval, _i.e._

$\|f\|_{\infty}=\sup_{t\in[0,T]}|f(t)|.$

###  2.2 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 1 (Taylor’s theorem).

For any $k$-th order continuously differentiable function $f$ (scale-valued or
matrix-valued) defined on an interval $[a,t]$, we have

$f(t)=\sum_{j=0}^{k-1}\frac{f^{(j)}(a)}{j!}(t-a)^{j}+\int_{a}^{t}\frac{f^{(k)}(s)}{(k-1)!}(t-s)^{k-1}ds.$
(8)

######  Lemma 2 (Variation of parameters formula).

Assume $U(t,s)$ solves the differential equation

$\partial_{t}U(t,s)=H(t)U(t,s),\quad U(s,s)=I.$ (9)

Then

  1. 1.

For any matrix-valued continuous function $R(t)$, the solution of the
differential equation

$\partial_{t}\widetilde{U}(t,0)=H(t)\widetilde{U}(t,0)+R(t),\quad\widetilde{U}(0,0)=I$
(10)

can be represented as

$\widetilde{U}(t,0)=U(t,0)+\int_{0}^{t}U(t,s)R(s)ds.$ (11)

  2. 2.

For any vector-valued continuous function $\vec{r}(t)$, the solution of the
differential equation

$\partial_{t}\vec{\widetilde{u}}(t)=H(t)\vec{\widetilde{u}}(t)+\vec{r}(t),\quad\vec{\widetilde{u}}(0)=\vec{u}_{0}$
(12)

can be represented as

$\vec{\widetilde{u}}(t)=U(t,0)\vec{u}_{0}+\int_{0}^{t}U(t,s)\vec{r}(s)ds.$
(13)

###  2.3 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

$\left(\prod_{j=1}^{J}g_{j}\right)\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]A\left[\prod_{k^{\prime}=1}^{K^{\prime}}\exp\left(\mathrm{i}h\xi^{\prime}_{k^{\prime}}H_{l^{\prime}_{k^{\prime}}}\right)\right].$
(14)

Here functions $g_{j}$’s can be $f_{1}$, $f_{2}$, or their derivatives;
$H_{l_{k}}$ and $H_{l^{\prime}_{k^{\prime}}}$ are either $H_{1}$ or $H_{2}$;
and $\xi_{k}$ and $\xi^{\prime}_{k^{\prime}}$ are some bounded real numbers.
The matrix $A$ is in the set $\\{H_{1},H_{2},[H_{1},H_{2}]\\}$ for standard
Trotter formula and can only be $[H_{1},H_{2}]$ for generalized Trotter
formula. Therefore it suffices to focus on each term in the form of Eq. 14 to
obtain error bounds.

The operator norm error bounds directly follow the error representations.
Under the assumption that $H_{1}$ and $H_{2}$ are bounded operators, we can
simply bound all the unitaries by $1$ and the (local) operator norm error
bounds become $\alpha h^{2}$ where $\alpha$ can be expressed in terms of
$\|H_{1}\|,\|H_{2}\|$ and $\|[H_{1},H_{2}]\|$ for the standard Trotter
formula, and of $\|[H_{1},H_{2}]\|$ for the generalized Trotter formula,
respectively. Notice that the local errors for both formulae are
$\mathcal{O}(h^{2})$ for first order schemes, which agrees with the order
condition. Furthermore, the preconstant for the generalized Trotter formula
only consists of commutators, while the preconstant for standard Trotter
formula still includes norms of $H_{1}$ and $H_{2}$ themselves.

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 $\ket{\psi}$ and the approximate state
$\ket{\widetilde{\psi}}$ obtained by Trotter formulae can be expressed as
$\left\lVert\ket{\widetilde{\psi}(h)}-\ket{\psi(h)}\right\rVert\leq\left\lVert\left(\widetilde{U}(h,0)-U(h,0)\right)\ket{\psi(0)}\right\rVert$.
Therefore the vector norm error bounds should be a linear combination of the
terms of the form

$\left\|\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]A\left[\prod_{k^{\prime}=1}^{K^{\prime}}\exp\left(\mathrm{i}h\xi^{\prime}_{k^{\prime}}H_{l^{\prime}_{k^{\prime}}}\right)\right]\vec{v}\right\|=\left\|A\left[\prod_{k^{\prime}=1}^{K^{\prime}}\exp\left(\mathrm{i}h\xi^{\prime}_{k^{\prime}}H_{l^{\prime}_{k^{\prime}}}\right)\right]\vec{v}\right\|.$
(15)

This can be obtained by applying the operator Eq. 14 to some vector $\vec{v}$,
and $\vec{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 $\left\lVert
A\right\rVert$ can be very large, it is possible for $\left\lVert
A\vec{v}\right\rVert$ to be small for certain vectors $\vec{v}$. As an
example, let us consider the continuous case in one dimension and we take
$H_{1}=-\Delta$ and $H_{2}$ being a bounded, smooth potential function $V(x)$.
The direct computation shows that

$H_{1}\psi=-\partial_{x}^{2}\psi,$
$[H_{1},H_{2}]\psi=[-\Delta,V]\psi=-(\partial_{x}^{2}V)\psi-2(\partial_{x}V)\partial_{x}\psi.$

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 $\psi$ 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.

##  3 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.

###  3.1 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

$U_{s,1}(t+h,t)=\exp\left(-\mathrm{i}f_{2}(t+h)H_{2}h\right)\exp\left(-\mathrm{i}f_{1}(t+h)H_{1}h\right).$
(16)

The first-order generalized Trotter formula is

$U_{g,1}(t+h,t)=\exp\left(-\mathrm{i}\int_{t}^{t+h}f_{2}(s)dsH_{2}\right)\exp\left(-\mathrm{i}\int_{t}^{t+h}f_{1}(s)dsH_{1}\right).$
(17)

The second-order standard Trotter formula is

$U_{s,2}(t+h,t)=\exp\left(-\frac{\mathrm{i}h}{2}f_{1}(t+h/2)H_{1}\right)\exp\left(-\mathrm{i}hf_{2}(t+h/2)H_{2}\right)\exp\left(-\frac{\mathrm{i}h}{2}f_{1}(t+h/2)H_{1}\right).$
(18)

The second-order generalized Trotter formula is

$U_{g,2}(t+h,t)=\exp\left(-\mathrm{i}\int_{t+h/2}^{t+h}f_{1}(s)dsH_{1}\right)\exp\left(-\mathrm{i}\int_{t}^{t+h}f_{2}(s)dsH_{2}\right)\exp\left(-\mathrm{i}\int_{t}^{t+h/2}f_{1}(s)dsH_{1}\right).$
(19)

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 $\int_{a}^{b}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.

###  3.2 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).

$\begin{split}U_{s,1}(h,0)-U(h,0)=\int_{0}^{h}U(h,s)\exp\left(-\mathrm{i}sf_{2}(s)H_{2}\right)E_{s,1}(s)\exp\left(-\mathrm{i}sf_{1}(s)H_{1}\right)ds\end{split}$
(20)

where

$\begin{split}E_{s,1}(h)&=\int_{0}^{h}f_{1}(h)f_{2}(s)\left(\exp\left(\operatorname{ad}_{\mathrm{i}sf_{2}(s)H_{2}}\right)([H_{1},H_{2}])\right)ds-\mathrm{i}hf_{1}^{\prime}(h)H_{1}-\mathrm{i}hf_{2}^{\prime}(h)H_{2}\\\
&\quad+\int_{0}^{h}sf_{1}(h)f_{2}^{\prime}(s)\left(\exp\left(\operatorname{ad}_{\mathrm{i}sf_{2}(s)H_{2}}\right)([H_{1},H_{2}])\right)ds.\end{split}$
(21)

######  Lemma 4 (Error representation of the first-order generalized Trotter
formula).

$\begin{split}U_{g,1}(h,0)-U(h,0)=\int_{0}^{h}U(h,s)\exp\left(-\mathrm{i}\int_{0}^{s}f_{2}(s^{\prime})ds^{\prime}H_{2}\right)E_{g,1}(s)\exp\left(-\mathrm{i}\int_{0}^{s}f_{1}(s^{\prime})ds^{\prime}H_{1}\right)ds\end{split}$
(22)

where

$\begin{split}E_{g,1}(h)&=\int_{0}^{h}f_{1}(h)f_{2}(s)\left(\exp\left(\operatorname{ad}_{\mathrm{i}\int_{0}^{h}f_{2}(s^{\prime})ds^{\prime}H_{2}}\right)([H_{1},H_{2}])\right)ds.\end{split}$
(23)

######  Lemma 5 (Error representation of the second-order standard Trotter
formula).

$\begin{split}U_{s,2}(h,0)-U(h,0)=\int_{0}^{h}&U(h,s)\exp\left(-\frac{\mathrm{i}s}{2}f_{1}(s/2)H_{1}\right)E_{s,2}(s)\\\
&\exp\left(-\mathrm{i}sf_{2}(s/2)H_{2}\right)\exp\left(-\frac{\mathrm{i}s}{2}f_{1}(s/2)H_{1}\right)ds\end{split}$
(24)

where $E_{s,2}$ is defined in Eq. 86.

######  Lemma 6 (Error representation of the second-order generalized Trotter
formula).

$\begin{split}U_{g,2}(h,0)-U(h,0)=\int_{0}^{h}&U(h,s)\exp\left(-\mathrm{i}\int_{s/2}^{s}f_{1}(s^{\prime})ds^{\prime}H_{1}\right)E_{g,2}(s)\\\
&\exp\left(-\mathrm{i}\int_{0}^{s}f_{2}(s^{\prime})ds^{\prime}H_{2}\right)\exp\left(-\mathrm{i}\int_{0}^{s/2}f_{1}(s^{\prime})ds^{\prime}H_{1}\right)ds\end{split}$
(25)

where $E_{g,2}$ is defined in Eq. 87.

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 $\mathcal{O}(h^{p+1})$
terms, we also include higher order terms $\mathcal{O}(h^{p+2})$ and/or
$\mathcal{O}(h^{p+3})$. This is because we aim at writing down the exact error
terms, and the Taylor expansion of terms like
$\exp(-\mathrm{i}hf_{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(-\mathrm{i}hf_{1}(h)H_{1})$,
then

$\frac{d}{dh}(\exp(-\mathrm{i}hf_{1}(h)H_{1}))=-\mathrm{i}f_{1}(h)H_{1}\exp(-\mathrm{i}hf_{1}(h)H_{1})-\mathrm{i}hf_{1}^{\prime}(h)H_{1}\exp(-\mathrm{i}hf_{1}(h)H_{1}).$

We can observe that the first term is $\mathcal{O}(1)$ and the second term is
$\mathcal{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 $\mathcal{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.

##  4 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. 1.

First-order standard Trotter formula:

$\|U_{s,1}(h,0)-U(h,0)\|\leq\alpha_{s,1}h^{2}+\beta_{s,1}h^{3},$ (26)

where

$\alpha_{s,1}=\frac{1}{2}\|f_{1}^{\prime}\|_{\infty}\|H_{1}\|+\frac{1}{2}\|f_{2}^{\prime}\|_{\infty}\|H_{2}\|+\frac{1}{2}\|f_{1}\|_{\infty}\|f_{2}\|_{\infty}\|[H_{1},H_{2}]\|$
(27)

and

$\beta_{s,1}=\frac{1}{6}\|f_{1}\|_{\infty}\|f_{2}^{\prime}\|_{\infty}\|[H_{1},H_{2}]\|.$
(28)

  2. 2.

First-order generalized Trotter formula:

$\|U_{g,1}(h,0)-U(h,0)\|\leq\alpha_{g,1}h^{2}$ (29)

where

$\alpha_{g,1}=\frac{1}{2}\|f_{1}\|_{\infty}\|f_{2}\|_{\infty}\|[H_{1},H_{2}]\|.$
(30)

  3. 3.

Second-order standard Trotter formula:

$\|U_{s,2}(h,0)-U(h,0)\|\leq\alpha_{s,2}h^{3}+\beta_{s,2}h^{4}+\gamma_{s,2}h^{5},$
(31)

where

$\begin{split}\alpha_{s,2}&=\frac{7}{24}\|f_{1}^{\prime\prime}\|_{\infty}\|H_{1}\|+\frac{1}{12}\|f_{1}^{\prime}\|_{\infty}\|f_{2}\|_{\infty}\|H_{1}\|+\frac{7}{24}\|f_{2}^{\prime\prime}\|_{\infty}\|H_{2}\|\\\
&\quad+\frac{1}{6}(\|f_{1}^{\prime}\|_{\infty}\|f_{2}\|_{\infty}+\|f_{1}\|_{\infty}\|f_{2}^{\prime}\|_{\infty})\|[H_{1},H_{2}]\|\\\
&\quad+\frac{1}{24}\|f_{1}\|_{\infty}^{2}\|f_{2}\|_{\infty}\|[H_{1},[H_{1},H_{2}]]\|+\frac{1}{12}\|f_{1}\|_{\infty}\|f_{2}\|_{\infty}^{2}\|[H_{2},[H_{1},H_{2}]]\|,\end{split}$
(32)
$\begin{split}\beta_{s,2}&=\frac{1}{64}\|f_{1}^{\prime}\|_{\infty}\|f_{2}^{\prime}\|_{\infty}\|H_{1}\|\\\
&\quad+\left(\frac{1}{192}\|f_{1}\|_{\infty}\|f_{2}^{\prime\prime}\|_{\infty}+\frac{1}{192}\|f_{1}^{\prime\prime}\|_{\infty}\|f_{2}\|_{\infty}+\frac{1}{48}\|f_{1}^{\prime}\|_{\infty}\|f_{2}^{\prime}\|_{\infty}\right)\|[H_{1},H_{2}]\|\\\
&\quad+\frac{1}{96}\|f_{1}\|_{\infty}\|f_{1}^{\prime}\|_{\infty}\|f_{2}\|_{\infty}\|[H_{1},[H_{1},H_{2}]]\|+\frac{1}{48}\|f_{1}\|_{\infty}\|f_{2}\|_{\infty}\|f_{2}^{\prime}\|_{\infty}\|[H_{2},[H_{1},H_{2}]]\|,\end{split}$
(33)

and

$\begin{split}\gamma_{s,2}=\frac{1}{960}\|f_{1}^{\prime}\|_{\infty}^{2}\|f_{2}\|_{\infty}\|[H_{1},[H_{1},H_{2}]]\|+\frac{1}{480}\|f_{1}\|_{\infty}\|f_{2}^{\prime}\|_{\infty}^{2}\|[H_{2},[H_{1},H_{2}]]\|.\end{split}$
(34)

  4. 4.

Second-order generalized Trotter formula:

$\|U_{g,2}(h,0)-U(h,0)\|\leq\alpha_{g,2}h^{3},$ (35)

where

$\begin{split}\alpha_{g,2}&=\left(\frac{7}{12}\|f_{1}\|_{\infty}\|f_{2}^{\prime}\|_{\infty}+\frac{11}{24}\|f_{1}^{\prime}\|_{\infty}\|f_{2}\|_{\infty}\right)\|[H_{1},H_{2}]\|\\\
&\quad+\frac{3}{8}\|f_{1}\|_{\infty}^{2}\|f_{2}\|_{\infty}\|[H_{1},[H_{1},H_{2}]]\|+\frac{1}{12}\|f_{1}\|_{\infty}\|f_{2}\|_{\infty}^{2}\|[H_{2},[H_{2},H_{1}]]\|.\end{split}$
(36)

######  Remark 2 (The preconstants in Theorem 1).

First, the preconstants of the standard Trotter formula involve the norms of
both the Hamiltonians as well as their commutators, while the preconstants of
the generalized Trotter formula of the first order only involve the commutator
$[H_{1},H_{2}]$, and those of the second order involve further nested
commutators. Second, the $p$-th order standard Trotter scheme ($p=1,2$)
depends on the $p$-th order derivatives of the control functions while the
$p$-th generalized Trotter scheme depends only on $(p-1)$-th order
derivatives. In this sense, when the time derivatives of $f_{1}$ and $f_{2}$
are large (or when the regularity of $f_{1},f_{2}$ are limited), the
generalized Trotter formula can further outperform the standard Trotter
method.

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

$\displaystyle\left\|\prod_{l=1}^{L}U_{s,1}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)-U(T,0)\right\|\leq\alpha_{s,1}\frac{T^{2}}{L}+\beta_{s,1}\frac{T^{3}}{L^{2}},$
(37)
$\displaystyle\left\|\prod_{l=1}^{L}U_{g,1}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)-U(T,0)\right\|\leq\alpha_{g,1}\frac{T^{2}}{L},$
(38)
$\displaystyle\left\|\prod_{l=1}^{L}U_{s,2}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)-U(T,0)\right\|\leq\alpha_{s,2}\frac{T^{3}}{L^{2}}+\beta_{s,2}\frac{T^{4}}{L^{3}}+\gamma_{s,2}\frac{T^{5}}{L^{4}},$
(39)
$\displaystyle\left\|\prod_{l=1}^{L}U_{g,2}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)-U(T,0)\right\|\leq\alpha_{g,2}\frac{T^{3}}{L^{2}},$
(40)

where preconstants $\alpha$ and $\beta$ are defined in Theorem 1.

##  5 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 $\left\lVert H_{1}\right\rVert\gg\left\lVert H_{2}\right\rVert$, 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 $\widetilde{C}$ with a tilde above to denote preconstants (with
possibly varying sizes in different inequalities) which can depend
polynomially on $H_{2}$, $\|f_{1}^{(k)}\|_{\infty}$ and
$\|f_{2}^{(k)}\|_{\infty}$ 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}^{\dagger}D_{1}$. Furthermore, we
assume for any vector $\vec{v}$, there exist constants
$\widetilde{C}_{1},\widetilde{C}_{2}$ such that

$\|[H_{1},H_{2}]\vec{v}\|\leq\widetilde{C}_{1}(\|D_{1}\vec{v}\|+\|\vec{v}\|),$
(41)

and

$\|[H_{1},[H_{1},H_{2}]]\vec{v}\|\leq\widetilde{C}_{2}(\|H_{1}\vec{v}\|+\|\vec{v}\|).$
(42)

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
$H=-\Delta+V(x)$, where $H_{1}=-\Delta$ is positive semidefinite and
$H_{2}=V(x)$ is a bounded operator. It will be helpful to understand 1 from a
continuous analog, in which $H_{1}=D_{1}^{\dagger}D_{1}$ and $D_{1}$ is a
first-order differential operator. A direct calculation shows that the
operator

$[-\Delta,V]=-\partial_{x}^{2}V-2(\partial_{x}V)\partial_{x}$

is a first-order differential operator, and

$[-\Delta,[-\Delta,V]]=\partial_{x}^{4}V+4\partial_{x}^{3}V\partial_{x}+4\partial_{x}^{2}V\partial_{x}^{2}$
(43)

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 1, we have the following:

  1. 1.

For any vector $\vec{v}$,

$\|D_{1}\vec{v}\|\leq\|H_{1}\vec{v}\|+\|\vec{v}\|.$ (44)

  2. 2.

Let $\xi$ be any real number such that $(\widetilde{C}_{1}+\|H_{2}\|)|\xi|\leq
1/2$. Then for any vector $\vec{v}$,

$\|H_{1}\exp\left(\mathrm{i}\xi H_{2}\right)\vec{v}\|\leq
2(\|H_{1}\vec{v}\|+\|\vec{v}\|).$ (45)

  3. 3.

Let $K$ be a positive integer, and $H_{l_{k}}$ be either $H_{1}$ or $H_{2}$,
and $\xi_{k}$ be some real numbers for $1\leq k\leq K$. Assume that
$(\widetilde{C}_{1}+\|H_{2}\|)|\xi_{k}|h\leq 1/2$, then for any vector
$\vec{v}$, all the following inequalities hold:

$\begin{split}\left\|H_{1}\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|&\leq\widetilde{C}(\|H_{1}\vec{v}\|+\|\vec{v}\|),\\\
\left\|H_{2}\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|&\leq\widetilde{C}\|\vec{v}\|,\\\
\left\|[H_{1},H_{2}]\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|&\leq\widetilde{C}\left(\sqrt{\|\vec{v}\|\|H_{1}\vec{v}\|}+\|\vec{v}\|\right),\\\
\left\|[H_{1},[H_{1},H_{2}]]\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|&\leq\widetilde{C}(\|H_{1}\vec{v}\|+\|\vec{v}\|),\\\
\left\|[H_{2},[H_{2},H_{1}]]\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|&\leq\widetilde{C}(\|H_{1}\vec{v}\|+\|\vec{v}\|).\end{split}$
(46)

for some constant $\widetilde{C}>0$ depending only on $\|H_{2}\|$.

###### Proof.

1\. By the definition of $D_{1}$ and the Cauchy-Schwarz inequality,

$\|D_{1}\vec{v}\|^{2}=(D_{1}\vec{v})^{\dagger}D_{1}\vec{v}=\vec{v}^{\dagger}H_{1}\vec{v}\leq\|\vec{v}\|\|H_{1}\vec{v}\|\leq(\|H_{1}\vec{v}\|+\|\vec{v}\|)^{2}.$
(47)

2\. We start with Taylor’s theorem of $\exp(\mathrm{i}t\xi H_{2})$ up to
first-order, then

$\begin{split}H_{1}\exp(\mathrm{i}t\xi
H_{2})\vec{v}&=H_{1}\vec{v}+\int_{0}^{t}\mathrm{i}\xi
H_{1}H_{2}\exp\left(\mathrm{i}\alpha\xi H_{2}\right)\vec{v}d\alpha\\\
&=H_{1}\vec{v}+\int_{0}^{t}\mathrm{i}\xi[H_{1},H_{2}]\exp\left(\mathrm{i}\alpha\xi
H_{2}\right)\vec{v}d\alpha+\int_{0}^{t}\mathrm{i}\xi
H_{2}H_{1}\exp\left(\mathrm{i}\alpha\xi
H_{2}\right)\vec{v}d\alpha.\end{split}$

The norm can be estimated using Eq. 41 as

$\begin{split}\left\|H_{1}\exp(\mathrm{i}t\xi
H_{2})\vec{v}\right\|&\leq\|H_{1}\vec{v}\|+\int_{0}^{t}|\xi|\|[H_{1},H_{2}]\exp\left(\mathrm{i}\alpha\xi
H_{2}\right)\vec{v}\|d\alpha+\int_{0}^{t}|\xi|\|H_{2}\|\left\|H_{1}\exp\left(\mathrm{i}\alpha\xi
H_{2}\right)\vec{v}\right\|d\alpha\\\
&\leq\|H_{1}\vec{v}\|+\widetilde{C}_{1}\|\xi\|\|\vec{v}\|t+\int_{0}^{t}\widetilde{C}_{1}|\xi|\|D_{1}\exp\left(\mathrm{i}\alpha\xi
H_{2}\right)\vec{v}\|d\alpha\\\
&\quad+\int_{0}^{t}|\xi|\|H_{2}\|\left\|H_{1}\exp\left(\mathrm{i}\alpha\xi
H_{2}\right)\vec{v}\right\|d\alpha.\end{split}$ (48)

Define $M(t):=\left\|H_{1}\exp(\mathrm{i}t\xi H_{2})\vec{v}\right\|$, and it
follows from Eq. 44 that

$\begin{split}\|D_{1}\exp\left(\mathrm{i}t\xi H_{2}\right)\vec{v}\|\leq
M(t)+\|\vec{v}\|.\end{split}$

Thus Eq. 48 can be rewritten as

$M(t)\leq\|H_{1}\vec{v}\|+2\widetilde{C}_{1}|\xi|\|\vec{v}\|t+\int_{0}^{t}|\xi|\left(\widetilde{C}_{1}+\|H_{2}\|\right)M(\alpha)\,d\alpha.$

Since $\|H_{1}\vec{v}\|+2\widetilde{C}_{1}|\xi|\|\vec{v}\|t$ is non-decreasing
with respect to $t$, applying Gronwall’s inequality yields the bound

$M(t)\leq\left(\|H_{1}\vec{v}\|+2\widetilde{C}_{1}|\xi|\|\vec{v}\|t\right)\exp\left(|\xi|(\widetilde{C}_{1}+\|H_{2}\|)t\right).$

Finally taking $t=1$ and applying the condition on $\xi$, the desired result
is achieved

$\left\|H_{1}\exp(\mathrm{i}\xi
H_{2})\vec{v}\right\|\leq\left(\|H_{1}\vec{v}\|+2\widetilde{C}_{1}|\xi|\|\vec{v}\|\right)\exp\left(|\xi|(\widetilde{C}_{1}+\|H_{2}\|)\right)\leq
2(\|H_{1}\vec{v}\|+\|\vec{v}\|).$

3\. We first show that it suffices to only prove the first inequality in Eq.
46. Let
$\vec{w}=\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}$.
Since $H_{2}$ is bounded, $\|H_{2}\vec{w}\|$ is directly bounded by
$\widetilde{C}\|\vec{v}\|$ . By Eqs. 41 and 47, and the fact
$\left\lVert\vec{w}\right\rVert=\left\lVert\vec{v}\right\rVert$, we have

$\|[H_{1},H_{2}]\vec{w}\|\leq\widetilde{C}_{1}(\|D_{1}\vec{w}\|+\|\vec{w}\|)\leq\widetilde{C}\left(\sqrt{\|\vec{w}\|\|H_{1}\vec{w}\|}+\|\vec{w}\|\right)=\widetilde{C}\left(\sqrt{\|\vec{v}\|\|H_{1}\vec{w}\|}+\|\vec{v}\|\right).$
(49)

Similarly Eq. 42 gives

$\|[H_{1},[H_{1},H_{2}]]\vec{w}\|\leq\widetilde{C}_{2}(\|H_{1}\vec{w}\|+\|\vec{w}\|)=\widetilde{C}_{2}(\|H_{1}\vec{w}\|+\|\vec{v}\|).$
(50)

Furthermore,

$\begin{split}\|[H_{2},[H_{2},H_{1}]]\vec{w}\|&\leq\|H_{2}[H_{2},H_{1}]\vec{w}\|+\|[H_{2},H_{1}]H_{2}\vec{w}\|\\\
&\leq\widetilde{C}\|[H_{2},H_{1}]\vec{w}\|+\widetilde{C}_{1}(\|D_{1}H_{2}\vec{w}\|+\|H_{2}\vec{w}\|)\\\
&\leq\widetilde{C}(\|D_{1}\vec{w}\|+\|\vec{w}\|)+\widetilde{C}(\|H_{1}H_{2}\vec{w}\|+\|H_{2}\vec{w}\|)\\\
&\leq\widetilde{C}(\|H_{1}\vec{w}\|+\|\vec{w}\|)+\widetilde{C}(\|[H_{1},H_{2}]\vec{w}\|+\|H_{2}H_{1}\vec{w}\|+\|\vec{w}\|)\\\
&\leq\widetilde{C}(\|H_{1}\vec{w}\|+\|\vec{w}\|)+\widetilde{C}(\|D_{1}\vec{w}\|+\|\vec{w}\|)\\\
&\leq\widetilde{C}(\|H_{1}\vec{w}\|+\|\vec{w}\|)=\widetilde{C}(\|H_{1}\vec{w}\|+\|\vec{v}\|).\end{split}$
(51)

Therefore we only need to bound $\|H_{1}\vec{w}\|$ further by
$\widetilde{C}(\|H_{1}\vec{v}\|+\|\vec{v}\|)$.

Notice that $\|H_{1}\exp(\mathrm{i}\xi H_{1})\vec{v}\|=\|\exp(\mathrm{i}\xi
H_{1})H_{1}\vec{v}\|=\|H_{1}\vec{v}\|$, together with Eq. 45,

$\|H_{1}\exp(\mathrm{i}\xi H_{l})\vec{v}\|\leq
2(\|H_{1}\vec{v}\|+\|\vec{v}\|)$ (52)

for $H_{l}$ being either $H_{1}$ or $H_{2}$. Then we have the recursive
relation

$\begin{split}\left\|H_{1}\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|&\leq
2\left\|H_{1}\left[\prod_{k=1}^{K-1}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|+2\left\|\left[\prod_{k=1}^{K-1}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|\\\
&=2\left\|H_{1}\left[\prod_{k=1}^{K-1}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|+2\left\|\vec{v}\right\|.\end{split}$
(53)

By applying this estimation $K$ times, we obtain the desired result

$\left\|H_{1}\left[\prod_{k=1}^{K}\exp\left(\mathrm{i}h\xi_{k}H_{l_{k}}\right)\right]\vec{v}\right\|\leq\widetilde{C}\left(\|H_{1}\vec{v}\|+\|\vec{v}\|\right).$
(54)

∎

Now we are ready to state our main theorems for the vector norm estimates.

######  Theorem 3.

For any vector $v$ and time step size
$h\leq(\|f_{1}\|_{\infty}+\|f_{2}\|_{\infty})^{-1}(\widetilde{C}_{1}+\|H_{2}\|)^{-1}/2$,
there exists a constant $\widetilde{C}$ such that

$\displaystyle\|U_{s,1}(h,0)\vec{v}-U(h,0)\vec{v}\|$
$\displaystyle\leq\widetilde{C}h^{2}(\|H_{1}\vec{v}\|+\|\vec{v}\|),$
$\displaystyle\|U_{g,1}(h,0)\vec{v}-U(h,0)\vec{v}\|$
$\displaystyle\leq\widetilde{C}h^{2}(\sqrt{\|\vec{v}\|\|H_{1}\vec{v}\|}+\|\vec{v}\|),$
$\displaystyle\|U_{s,2}(h,0)\vec{v}-U(h,0)\vec{v}\|$
$\displaystyle\leq\widetilde{C}h^{3}(\|H_{1}\vec{v}\|+\|\vec{v}\|),$
$\displaystyle\|U_{g,2}(h,0)\vec{v}-U(h,0)\vec{v}\|$
$\displaystyle\leq\widetilde{C}h^{3}(\|H_{1}\vec{v}\|+\|\vec{v}\|).$

###### Proof.

We start with the error representations (Lemmas 3, 4, 5 and 6). By multiplying
a vector $\vec{v}$ on the right of the error representations, bounding all the
$f_{j}^{(k)}$ 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

$\displaystyle\|U_{s,1}(h,0)\vec{v}-U(h,0)\vec{v}\|$
$\displaystyle\leq\widetilde{C}\theta_{s,1}h^{2},$
$\displaystyle\|U_{g,1}(h,0)\vec{v}-U(h,0)\vec{v}\|$
$\displaystyle\leq\widetilde{C}\theta_{g,1}h^{2},$
$\displaystyle\|U_{s,2}(h,0)\vec{v}-U(h,0)\vec{v}\|$
$\displaystyle\leq\widetilde{C}\theta_{s,2}h^{3},$
$\displaystyle\|U_{g,2}(h,0)\vec{v}-U(h,0)\vec{v}\|$
$\displaystyle\leq\widetilde{C}\theta_{g,2}h^{3},$

where $\theta_{s}$ and $\theta_{g}$ are linear combinations of the terms in
the form of Eq. 46, with an exception that $\theta_{g,1}$ only consists terms
like $\|[H_{1},H_{2}]\vec{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.

######  Theorem 4.

Let $T>0$ be the evolution time, $\vec{\psi}(t)$ be the exact solution of the
dynamics Eq. 1, and the dynamics Eq. 1 is discretized via standard and
generalized Trotter formulae with $L$ time steps such that the time step size
$h=T/L$ is bounded by
$(\|f_{1}\|_{\infty}+\|f_{2}\|_{\infty})^{-1}(\widetilde{C}_{1}+\|H_{2}\|)^{-1}/2$.
Then there exists a constant $\widetilde{C}$ such that

$\displaystyle\left\|\left(\prod_{l=1}^{L}U_{s,1}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)-U(T,0)\vec{\psi}(0)\right\|\leq\widetilde{C}\frac{T^{2}}{L}\left(\sup_{t\in[0,T]}\|H_{1}\vec{\psi}(t)\|+\|\vec{\psi}(0)\|\right),$
(55)
$\displaystyle\left\|\left(\prod_{l=1}^{L}U_{g,1}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)-U(T,0)\vec{\psi}(0)\right\|\leq\widetilde{C}\frac{T^{2}}{L}\left(\sup_{t\in[0,T]}\sqrt{\|\vec{\psi}(0)\|\|H_{1}\vec{\psi}(t)\|}+\|\vec{\psi}(0)\|\right),$
(56)
$\displaystyle\left\|\left(\prod_{l=1}^{L}U_{s,2}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)-U(T,0)\vec{\psi}(0)\right\|\leq\widetilde{C}\frac{T^{3}}{L^{2}}\left(\sup_{t\in[0,T]}\|H_{1}\vec{\psi}(t)\|+\|\vec{\psi}(0)\|\right),$
(57)
$\displaystyle\left\|\left(\prod_{l=1}^{L}U_{g,2}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)-U(T,0)\vec{\psi}(0)\right\|\leq\widetilde{C}\frac{T^{3}}{L^{2}}\left(\sup_{t\in[0,T]}\|H_{1}\vec{\psi}(t)\|+\|\vec{\psi}(0)\|\right).$
(58)

#### Remark

Although the error for the first-order generalized Trotter formula can also be
bounded by $\sup_{t\in[0,T]}\|H_{1}\vec{\psi}(t)\|+\|\vec{\psi}(0)\|$ as the
other schemes by a direct application of the Cauchy-Schwarz inequality, we
keep it as
$\sup_{t\in[0,T]}\sqrt{\|\vec{\psi}(0)\|\|H_{1}\vec{\psi}(t)\|}+\|\vec{\psi}(0)\|$
which promotes a better dependence on $\|H_{1}\vec{\psi}(t)\|$. This
improvement is achievable only for the first-order generalized Trotter formula
since its error only contains terms depending on $[H_{1},H_{2}]$ which process
a better bound as in Eq. 46, while the other schemes also contains terms
depending on $H_{1}$ and/or the nested commutators $[H_{1},[H_{1},H_{2}]]$ and
$[H_{1},[H_{2},H_{1}]]$.

###### 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 that $\|\vec{\psi}(t)\|=\|\vec{\psi}(0)\|$ for all $t\in[0,T]$, we
obtain

$\begin{split}&\quad\left\|\left(\prod_{l=1}^{L}U_{g,2}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)-U(T,0)\vec{\psi}(0)\right\|\\\
&=\left\|\left(\prod_{l=1}^{L}U_{g,2}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)-\left(\prod_{l=1}^{L}U\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)\right\|\\\
&\leq\sum_{k=1}^{L}\Bigg{\|}\left(\prod_{l=k+1}^{L}U_{g,2}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\left(\prod_{l=1}^{k}U\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)\\\
&\quad\quad\quad\quad-\left(\prod_{l=k}^{L}U_{g,2}\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\left(\prod_{l=1}^{k-1}U\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)\Bigg{\|}\\\
&\leq\sum_{k=1}^{L}\left\|\left(\prod_{l=1}^{k}U\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)-U_{g,2}\left(\frac{kT}{L},\frac{(k-1)T}{L}\right)\left(\prod_{l=1}^{k-1}U\left(\frac{lT}{L},\frac{(l-1)T}{L}\right)\right)\vec{\psi}(0)\right\|\\\
&=\sum_{k=1}^{L}\left\|\left(U\left(\frac{kT}{L},\frac{(k-1)T}{L}\right)-U_{g,2}\left(\frac{kT}{L},\frac{(k-1)T}{L}\right)\right)\vec{\psi}((k-1)T/L)\right\|\\\
&\leq\widetilde{C}\frac{T^{3}}{L^{3}}\sum_{k=1}^{L}\left(\left\|H_{1}\vec{\psi}((k-1)T/L)\right\|+\left\|\vec{\psi}((k-1)T/L)\right\|\right)\\\
&\leq\widetilde{C}\frac{T^{3}}{L^{2}}\left(\sup_{t\in[0,T]}\left\|H_{1}\vec{\psi}(t)\right\|+\left\|\vec{\psi}(0)\right\|\right).\end{split}$
(59)

∎

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 $\left\lVert
H_{2}\right\rVert$, and the dependence on $H_{1}$ only appears in the form of
$\|H_{1}\vec{\psi}\|$. Such a difference implies that the vector norm bounds
can be much sharper than the operator norm bounds when $\left\lVert
H_{1}\right\rVert$ is very large, but $\left\lVert
H_{1}\vec{\psi}\right\rVert$ and $\left\lVert H_{2}\right\rVert$ 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 $\epsilon$.

##  6 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 $\epsilon$-approximation of
the wavefunction at time $T\sim\mathcal{O}(1)$, where $D=[0,1]$ with periodic
boundary conditions. Throughout the section we make the following assumptions.

  1. 1.

$M_{\text{eff}}(t)$ is positive function, and is uniformly bounded from below.

  2. 2.

$M_{\text{eff}}(t),\omega(t)$ are second-order continuously differentiable
functions in $t$ with uniformly bounded function and derivative values up to
second order.

  3. 3.

$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 $[-\Delta,[-\Delta,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 $\widetilde{C}$ or the big-O notation $\mathcal{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\leq k\leq n-1$. Then the semi-discretized dynamics becomes
$\mathrm{i}\partial_{t}\vec{\psi}(t)=H(t)\vec{\psi}(t)$. Here the $k$-th entry
of $\vec{\psi}(t)$ will be an approximation of the exact wavefunction
evaluated at $t$ and $x=(k-1)/n$. $H(t)=f_{1}(t)H_{1}+f_{2}(t)H_{2}$ with

$H_{1}=n^{2}\left(\begin{array}[]{ccccc}2&-1&&&-1\\\ -1&2&-1&&\\\
&\ddots&\ddots&\ddots&\\\ &&-1&2&-1\\\ -1&&&-1&2\\\ \end{array}\right),$ (60)

and

$H_{2}=\text{diag}(V(0),V(1/n),\cdots,V((n-1)/n)).$ (61)

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 1 with $H_{1}=D_{1}^{\dagger}D_{1}$,

$D_{1}=n\left(\begin{array}[]{ccccc}-1&&&&1\\\ 1&-1&&&\\\ &\ddots&\ddots&&\\\
&&1&-1&\\\ &&&1&-1\end{array}\right),$ (62)

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.

######  Lemma 8.

Consider $H_{1}$ and $H_{2}$ defined in Eqs. 60 and 61, then 1 is satisfied
under the rescaled 2-norm $\|\cdot\|_{\star}$ with $D_{1}$ defined in Eq. 62.

Lemma 8 can be directly used to provide the scaling of the Hamiltonians and
their commutators in terms of $n$. Using the fact that the matrix 2-norms can
be estimated by considering its 1-norm (the maximum absolute column sum) and
$\infty$-norm (the maximum absolute row sum) by virtue of the fact that
$\|M\|_{2}\leq\sqrt{\|M\|_{1}\|M\|_{\infty}}$, we obtain the following Lemma.
Notice that this result is consistent with the continuous picture that their
continuous analogs $[-\Delta,V]$ and $[V,[-\Delta,V]]$ are differential
operators of the first order and $[-\Delta,[-\Delta,V]]$ of second order.

######  Lemma 9.

Consider $H_{1}$, $H_{2}$ and $D_{1}$ defined in Eqs. 60, 61 and 62, then

$\|H_{1}\|=\mathcal{O}(n^{2}),\quad\|H_{2}\|=\mathcal{O}(1),\quad\|D_{1}\|=\mathcal{O}(n),$
(63)

and

$\|[H_{1},H_{2}]\|=\mathcal{O}(n),\quad\|[H_{2},[H_{1},H_{2}]]\|\leq
2\|[H_{2}\|\|[H_{1},H_{2}]]\|=\mathcal{O}(n),\quad\|[H_{1},[H_{1},H_{2}]]\|=\mathcal{O}(n^{2}).$
(64)

######  Remark 3.

Although $\|H_{1}\|$ depends quadratically in $n$, the time-independent
simulations $\exp(-\mathrm{i}H_{1})$ and $\exp(-\mathrm{i}H_{2})$ can still be
performed efficiently. For $H_{1}$, it can be diagonalized under Fourier
transform. Specifically, let $F=(\omega^{jk}/\sqrt{n})_{j,k=0}^{n-1}$ be the
Fourier transform unitary matrix with $\omega=\exp(2\pi\mathrm{i}/n)$, then
$H_{1}=F\Lambda F^{\dagger}$ where
$\Lambda=\text{diag}(2n^{2}(1-\cos(2\pi\mathrm{i}j/n)))_{j=0}^{n-1}$.
Therefore $\exp(-\mathrm{i}H_{1})=F\exp(-\mathrm{i}\Lambda)F^{\dagger}$ can be
simulated efficiently, by first applying inverse quantum Fourier transform
$F^{\dagger}$, then applying fast-forwarding techniques [13, 1] for
$\exp(-\mathrm{i}\Lambda)$, 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(-\mathrm{i}H_{1})$ and $\exp(-\mathrm{i}H_{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._
$\|\vec{\widetilde{\psi}}(t)-(\phi(t,k/n))_{k=0}^{n-1}\|_{\star}$, where
$\vec{\widetilde{\psi}}(t)$ is the numerical solution after spatial and time
discretization at time $t$, and $\phi(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
$(\phi(t,k/n))_{k=0}^{n-1}$ is a discrete representation of the function
$\phi$, which is normalized under continuous $L^{2}$ norm rather than discrete
2-norm. Furthermore, if we encode the spatial discretized solution
$\vec{\psi}$ into a quantum state, then the normalized condition requires
$\ket{\psi}\sim\frac{1}{\sqrt{n}}\vec{\psi}$, thus
$\|\ket{\psi}\|\sim\frac{1}{\sqrt{n}}\|\vec{\psi}\|=\|\vec{\psi}\|_{\star}$,
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.

######  Lemma 10.

Let the exact solution of the Schrödinger equation with Hamiltonian Eq. 3 be
$\phi(t,x)$. Then

  1. 1.

for any $0\leq t\leq T,x\in[0,1]$,

$\begin{split}&\quad|n^{2}(\phi(t,x+1/n)-2\phi(t,x)+\phi(t,x-1/n))-\Delta\phi(t,x)|\\\
&\leq\frac{1}{3n^{2}}\sup_{y\in[0,1]}\left|\frac{\partial^{4}}{\partial
x^{4}}\phi(t,y)\right|.\end{split}$ (65)

  2. 2.

Let $\vec{\psi}(t)$ denote the solution of the dynamics

$\mathrm{i}\partial_{t}\vec{\psi}(t)=(f_{1}(t)H_{1}+f_{2}(t)H_{2})\vec{\psi}(t)$
(66)

where $H_{1}$ and $H_{2}$ are the discretized Hamiltonian defined in Eq. 60
and Eq. 61, then for any $0\leq t\in T$,

$\|(\phi(t,k/n))_{k=0}^{n-1}-\vec{\psi}(t)\|_{\star}\leq\frac{t}{3n^{2}}\|f_{1}\|_{\infty}\sup_{s\in[0,t],y\in[0,1]}\left|\frac{\partial^{4}}{\partial
x^{4}}\phi(s,y)\right|.$ (67)

This lemma shows that in order to make the vector error induced by the spatial
discretization bounded by $\epsilon$, it suffices to choose

$n=\mathcal{O}\left(\frac{T^{1/2}}{\epsilon^{1/2}}\left(\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|\right)^{1/2}\right).$ (68)

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 $\vec{\psi}(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 $\vec{\psi}_{s,p}(t)$ and $\vec{\psi}_{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
$\widetilde{C}>0$, independent of $T,L,H_{1},n,\vec{\psi},\phi$, such that

  1. 1.

$\begin{split}\|U_{s,1}(T,0)-U(T,0)\|&\leq\widetilde{C}\frac{n^{2}T^{2}}{L},\\\
\|U_{g,1}(T,0)-U(T,0)\|&\leq\widetilde{C}\frac{nT^{2}}{L},\\\
\|U_{s,2}(T,0)-U(T,0)\|&\leq\widetilde{C}\frac{n^{2}T^{3}}{L^{2}},\\\
\|U_{g,2}(T,0)-U(T,0)\|&\leq\widetilde{C}\frac{n^{2}T^{3}}{L^{2}}.\end{split}$
(69)

  2. 2.

$\begin{split}\|\vec{\psi}_{s,1}(T)-\vec{\psi}(T)\|_{\star}&\leq\widetilde{C}\frac{T^{2}}{L}\left((T+1)\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|+\sup\left|\frac{\partial^{2}\phi}{\partial
x^{2}}\right|+\sup_{y\in[0,1]}\left|\phi(0,y)\right|\right),\\\
\|\vec{\psi}_{g,1}(T)-\vec{\psi}(T)\|_{\star}&\leq\widetilde{C}\frac{T^{2}}{L}\left[\left((T+1)\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|+\sup\left|\frac{\partial^{2}\phi}{\partial
x^{2}}\right|\right)^{1/2}\left(\sup_{y\in[0,1]}\left|\phi(0,y)\right|\right)^{1/2}+\sup_{y\in[0,1]}\left|\phi(0,y)\right|\right],\\\
\|\vec{\psi}_{s,2}(T)-\vec{\psi}(T)\|_{\star}&\leq\widetilde{C}\frac{T^{3}}{L^{2}}\left((T+1)\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|+\sup\left|\frac{\partial^{2}\phi}{\partial
x^{2}}\right|+\sup_{y\in[0,1]}\left|\phi(0,y)\right|\right),\\\
\|\vec{\psi}_{g,2}(T)-\vec{\psi}(T)\|_{\star}&\leq\widetilde{C}\frac{T^{3}}{L^{2}}\left((T+1)\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|+\sup\left|\frac{\partial^{2}\phi}{\partial
x^{2}}\right|+\sup_{y\in[0,1]}\left|\phi(0,y)\right|\right),\end{split}$ (70)

where the notation $\sup$ without any subscript should be interpreted as
$\sup_{t\in[0,T],x\in[0,1]}$.

######  Remark 4.

The condition that $L$ should be sufficiently large is to ensure that the
lowest order term in the error bounds are dominant and to allow us to discard
the higher order terms. This can be guaranteed by requiring the desired level
of error $\epsilon$ to be sufficiently small in the complexity estimate later.

###### Proof.

1\. The result follows by combining Theorem 2 and the scaling of the matrix
norms provided in Lemma 9.

2\. According to Theorem 4 and the fact that
$\|\vec{\psi}(0)\|_{\star}\leq\sup_{y\in[0,1]}|\psi(0,y)|$, we only need to
bound $\|H_{1}\vec{\psi}(t)\|_{\star}$ for any $t\in[0,T]$. Let
$r(t,x)=n^{2}(\phi(t,x+1/n)-2\phi(t,x)+\phi(t,x-1/n))-\Delta\phi(t,x)$ where
$\phi(t,x)$ is the exact solution before any discretization. By Lemma 10,

$\begin{split}\|H_{1}\vec{\psi}(t)\|_{\star}&\leq\|H_{1}(\vec{\psi}(t)-(\phi(t,k/n))_{k=0}^{n-1})\|_{\star}+\|H_{1}(\phi(t,k/n))_{k=0}^{n-1}\|_{\star}\\\
&\leq\|H_{1}(\vec{\psi}(t)-(\phi(t,k/n))_{k=0}^{n-1}\|_{\star}+\|((\Delta\phi(t,k/n))_{k=0}^{n-1}\|_{\star}+\|(r(t,k/n))_{k=0}^{n-1}\|_{\star}\\\
&\leq\widetilde{C}\left(t\sup_{s\in[0,T],y\in[0,1]}\left|\frac{\partial^{4}}{\partial
x^{4}}\phi(s,y)\right|+\sup_{y\in[0,1]}\left|\frac{\partial^{2}}{\partial
x^{2}}\phi(t,y)\right|+\frac{1}{n^{2}}\sup_{s\in[0,T],y\in[0,1]}\left|\frac{\partial^{4}}{\partial
x^{4}}\phi(s,y)\right|\right)\\\
&\leq\widetilde{C}\left((T+1)\sup_{s\in[0,T],y\in[0,1]}\left|\frac{\partial^{4}}{\partial
x^{4}}\phi(s,y)\right|+\sup_{s\in[0,T],y\in[0,1]}\left|\frac{\partial^{2}}{\partial
x^{2}}\phi(s,y)\right|\right)\end{split}$ (71)

∎

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 $\epsilon$. In particular, the
operator norm error bound involves the derivatives of the exact solution of
interest $\phi$. 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 $\epsilon$-approximation in
rescaled 2-norm of the solution $\phi(t,x)$. Let $L_{\text{ope},s,1}$ and
$L_{\text{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_{\text{vec},s,1}$ and $L_{\text{vec},g,1}$ denote the
estimates from vector norm error bounds, and $L_{\text{ope},s,2}$,
$L_{\text{ope},g,2}$, $L_{\text{vec},s,2}$, $L_{\text{vec},g,2}$ are the
corresponding estimates for second-order schemes. Then for sufficiently small
$\epsilon$,

$\begin{split}&L_{\text{ope},s,1}=\mathcal{O}\left(\frac{T^{3}}{\epsilon^{2}}\left(\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|\right)\right),\\\
&L_{\text{ope},g,1}=\mathcal{O}\left(\frac{T^{5/2}}{\epsilon^{3/2}}\left(\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|\right)^{1/2}\right),\\\
&L_{\text{vec},s,1}=\mathcal{O}\left(\frac{T^{2}}{\epsilon}\left((T+1)\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|+\sup\left|\frac{\partial^{2}\phi}{\partial
x^{2}}\right|+\sup_{y\in[0,1]}\left|\phi(0,y)\right|\right)\right),\\\
&L_{\text{vec},g,1}=\mathcal{O}\left(\frac{T^{2}}{\epsilon}\left[\left((T+1)\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|+\sup\left|\frac{\partial^{2}\phi}{\partial
x^{2}}\right|\right)^{1/2}\left(\sup_{y\in[0,1]}\left|\phi(0,y)\right|\right)^{1/2}+\sup_{y\in[0,1]}\left|\phi(0,y)\right|\right]\right),\end{split}$
(72)

and

$\begin{split}&L_{\text{ope},s,2}=L_{\text{ope},g,2}=\mathcal{O}\left(\frac{T^{2}}{\epsilon}\left(\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|\right)^{1/2}\right),\\\
&L_{\text{vec},s,2}=L_{\text{vec},g,2}=\mathcal{O}\left(\frac{T^{3/2}}{\epsilon^{1/2}}\left((T+1)\sup\left|\frac{\partial^{4}\phi}{\partial
x^{4}}\right|+\sup\left|\frac{\partial^{2}\phi}{\partial
x^{2}}\right|+\sup_{y\in[0,1]}\left|\phi(0,y)\right|\right)^{1/2}\right).\end{split}$
(73)

###### Proof.

The complexity can be estimated by requiring both error bounds in Lemma 10 and
Lemma 11 to be smaller than $\epsilon$. First, by requiring the right hand
side of Eq. 67 to be bounded by $\epsilon$, 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 $\epsilon$, we obtain the complexity estimates. ∎

Theorem 5 clearly illustrates the advantage of vector norm error bounds in
terms of the desired level of error $\epsilon$. More precisely, the total
number of required Trotter steps estimated from vector norm bounds only scales
$\mathcal{O}(1/\epsilon^{1/p})$ for $p$-th order schemes. This is because the
operator norm error bounds depend on the spatial discretization $n$, where
$n=\mathcal{O}(1/\epsilon^{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 $\epsilon$ in
Table 2, where the simulation time $T$ is $\mathcal{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
$\operatorname{poly}(1/\epsilon)$. 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 $\epsilon$ for all the schemes
we consider.

######  Remark 5 (_a priori_ estimates of the solution $\phi$).

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}\equiv 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 $\epsilon$ between the standard and generalized Trotter formulae.
As discussed before, this is due to the fact that $\left\lVert
H_{1}\right\rVert$ and $\left\lVert[H_{1},[H_{1},H_{2}]]\right\rVert$ 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.

##  7 Numerical Results

To illustrate the difference between the operator norm and vector norm, we
consider the following Hamiltonian

$H(t)=-\frac{1}{2}(2+\sin(at+0.5))\Delta+(1+\cos(t))V(x),\quad
V(x)=1-\cos(x),\quad x\in[-\pi,\pi]$ (74)

with periodic boundary conditions. Here $a>0$ controls the magnitude of the
derivatives $\|f_{1}^{\prime}\|_{\infty}$ and
$\|f_{1}^{\prime\prime}\|_{\infty}$. 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 $-\Delta$ 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 $\vec{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 $\|[H_{1},[H_{1},H_{2}]]\|$ grows
quadratically with respect the the number of spatial grids while
$\|[H_{1},H_{2}]\|$ scales linearly, which agrees with Lemma 9. However, the
vector norms $\|[H_{1},[H_{1},H_{2}]]\vec{v}\|_{\star}$,
$\|[H_{1},H_{2}]\vec{v}\|_{\star}$ remain of the same scale. This behavior is
not restricted to the specific spatial discretization. Moreover, the scalings
of $\|D_{1}\vec{v}\|_{\star}$ and $\|H_{1}\vec{v}\|_{\star}$ are found to be
the same as those of $\|[H_{1},[H_{1},H_{2}]]\vec{v}\|_{\star}$ and
$\|[H_{1},H_{2}]\vec{v}\|_{\star}$. This verifies the in assumptions Eqs. 41
and 42, which is also proved for the finite difference scheme in Lemma 8.

(a) Finite difference discretization (b) Fourier discretization

Figure 1: Operator and vector norms of a smooth vector for various numbers of
spatial grids $n$. $\|[H_{1},[H_{1},H_{2}]]\|$ and $\|[H_{1},H_{2}]\|$ scales
quadratically and linearly with respect to $n$, but the vector norms do not
grow as $n$ gets larger.

We then verify the scaling of the errors with respect to $n$. The initial
wavefunction is $\phi(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 $\left\lVert
f_{1}\right\rVert_{\infty},\left\lVert
f_{1}^{\prime}\right\rVert_{\infty},\left\lVert
f_{1}^{\prime\prime}\right\rVert_{\infty}$ are comparable, and there is no
significant difference in the preconstants. However, when $a=10$, $\left\lVert
f_{1}^{\prime\prime}\right\rVert_{\infty}$ is one order of magnitude larger
than $\left\lVert f_{1}^{\prime}\right\rVert_{\infty},\left\lVert
f_{1}\right\rVert_{\infty}$. 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.

Moreover, we compare the scaling of the number of Trotter steps for various
precision $\epsilon$, measuring the relative error via the vector norm. We fix
$a=10$, $T=0.16$, and consider the precision $\epsilon$ as $2^{-10}$,
$2^{-12}$, $2^{-14}$, $2^{-16}$, $2^{-18}$, $2^{-20}$ and take
$n\propto\epsilon^{-0.5}$ as $2^{5}$, $2^{6}$, $2^{7}$, $2^{8}$, $2^{9}$ and
$2^{10}$. As is presented in Fig. 3, both second-order Trotter formulae
requires the number of Trotter steps $L=\mathcal{O}(\epsilon^{-0.5})$ while it
requires $L=\mathcal{O}(\epsilon^{-1})$ for both first-order Trotter formulae.
These results agree with Theorem 5.

(a) $a=1$. Finite difference discretization (b) $a=1$. Fourier discretization

(c) $a=10$. Finite difference discretization (d) $a=10$. Fourier
discretization

Figure 2: Relative Errors in the operator and vector norms. In the legend, “g”
stands for the generalized Trotter formula and “s” for the standard Trotter
formula. The error in operator norm is labeled as “ope” while the one in
vector norm as “vec”. First Row: $a=1$ with slowly varying control functions.
Second Row: $a=10$ with fast varying control functions. Figure 3: The number
of Trotter steps required to achieve various precision for the relative error
in the vector norm. The spatial discretization is finite difference. Both
second-order Trotter formulae scales proportionally to $\epsilon^{-0.5}$ while
the first-order formulae scales as $\epsilon^{-1}$, which agrees with the
theoretical bounds.

##  8 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
$\partial_{t}^{k}H_{j}(s)$ and $\partial_{t}^{k^{\prime}}H_{j}(s^{\prime})$
commute for any $j=1,2$ and $k,k^{\prime},s,s^{\prime}$ (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 $\mathcal{C}_{k}$ denote the set of the norms of all
possible $k$-th order nested commutators of $H_{1}$ and $H_{2}$, for example
$\mathcal{C}_{0}=\\{\|H_{1}\|,\|H_{2}\|\\}$,
$\mathcal{C}_{1}=\\{\|[H_{1},H_{2}]\|\\}$, and
$\mathcal{C}_{2}=\\{\|[H_{1},[H_{1},H_{2}]]\|,\|[H_{2},[H_{2},H_{1}]]\|\\}$.
For $p$-th order schemes, we expect that the one-step operator norm error
bounds for the standard and generalized Trotter formula scales as
$\mathcal{O}(\alpha_{s,p}h^{p+1})$, $\mathcal{O}(\alpha_{g,p}h^{p+1})$,
respectively. Here $\alpha_{s,p}$ is a linear combination of terms in the set
$\bigcup_{k=0}^{p}\mathcal{C}_{k}$, while $\alpha_{g,p}$ is expressed as a
linear combination of terms in the set $\bigcup_{k=1}^{p}\mathcal{C}_{k}$.
Hence the difference lies in whether $\mathcal{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 1 in this work) becomes

$\|[\underbrace{H_{1},[H_{1},\cdots,[H_{1}}_{k\text{
repeats}},H_{2}]]\cdots]\vec{v}\|\leq\mathcal{O}\left(\|H_{1}^{k/2}\vec{v}\|+\|\vec{v}\|\right)$
(75)

for any $1\leq k\leq 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 $\epsilon$ may become much weaker by employing higher order
discretization schemes. In such a scenario, the spectral norms $\left\lVert
H_{1}\right\rVert$ and $\left\lVert H_{2}\right\rVert$ 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 $\eta$-electron
sub-manifold. Although much smaller than that on the entire space, the
operator norm on the $\eta$-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.

## Acknowledgments:

This work was partially supported by the Department of Energy under Grant No.
DE-SC0017867 (D.A.,L.L.), under Grant No. DE-AC02-05CH11231, No. FWP-
NQISCCAWL, and by the National Science Foundation under the QLCI program
through grant number OMA-2016245 (L.L.). We thank Andrew Childs, Jinpeng Liu,
Yuan Su, Yu Tong for helpful discussions.

## References

  * [1] G. R. Ahokas.  Improved algorithms for approximate quantum Fourier transforms and sparse Hmailtonian simulations.  University of Calgary, 2004.  doi:10.11575/PRISM/22839.
  * [2] T. Albash and D. A. Lidar.  Adiabatic quantum computation.  Rev. Mod. Phys., 90:015002, 2018.  doi:10.1103/RevModPhys.90.015002.
  * [3] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders.  Efficient quantum algorithms for simulating sparse Hamiltonians.  Commun. Math. Phys., 270(2):359–371, 2007.  doi:10.1007/s00220-006-0150-x.
  * [4] D. W. Berry and A. M. Childs.  Black-box Hamiltonian simulation and unitary implementation.  Quantum Information & Computation, 12(1-2):29–62, 2012.  doi:10.26421/QIC12.1-2.
  * [5] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma.  Exponential improvement in precision for simulating sparse Hamiltonians.  In Proceedings of the forty-sixth annual ACM symposium on Theory of computing, pages 283–292, 2014.  doi:10.1145/2591796.2591854.
  * [6] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma.  Simulating Hamiltonian dynamics with a truncated Taylor series.  Phys. Rev. Lett., 114:090502, 2015.  doi:10.1103/PhysRevLett.114.090502.
  * [7] D. W. Berry, A. M. Childs, and R. Kothari.  Hamiltonian simulation with nearly optimal dependence on all parameters.  Proceedings of the 56th IEEE Symposium on Foundations of Computer Science, pages 792–809, 2015.  doi:10.1109/FOCS.2015.54.
  * [8] D. W. Berry, A. M. Childs, Y. Su, X. Wang, and N. Wiebe.  Time-dependent Hamiltonian simulation with $l^{1}$-norm scaling.  Quantum, 4:254, 2020.  doi:10.22331/q-2020-04-20-254.
  * [9] D. W. Berry, R. Cleve, and S. Gharibian.  Gate-efficient discrete simulations of continuous-time quantum query algorithms.  Quantum Information and Computation, 14(1-2):1–30, 2014.  doi:10.26421/QIC14.1-2-1.
  * [10] J. Bourgain.  On growth of sobolev norms in linear schrödinger equations with smooth time dependent potential.  Journal d’Analyse Mathématique, 77(1):315–348, 1999.  doi:10.1007/BF02791265.
  * [11] E. Campbell.  Random compiler for fast Hamiltonian simulation.  Phys. Rev. Lett., 123(7):070503, 2019.  doi:10.1103/PhysRevLett.123.070503.
  * [12] C.-F. Chen, H.-Y. Huang, R. Kueng, and J. A. Tropp.  Quantum simulation via randomized product formulas: Low gate complexity with accuracy guarantees.  2020\.  arXiv:2008.11751.
  * [13] A. M. Childs, R. Cleve, E. Deotto, E. Farhi, S. Gutmann, and D. A. Spielman.  Exponential algorithmic speedup by a quantum walk.  In Proceedings of the thirty-fifth annual ACM symposium on Theory of computing, pages 59–68, 2003.  doi:10.1145/780542.780552.
  * [14] A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su.  Toward the first quantum simulation with quantum speedup.  Proc. Nat. Acad. Sci., 115:9456–9461, 2018.  doi:10.1073/pnas.1801723115.
  * [15] A. M. Childs, A. Ostrander, and Y. Su.  Faster quantum simulation by randomization.  Quantum, 3:182, 2019.  doi:10.22331/q-2019-09-02-182.
  * [16] A. M. Childs and Y. Su.  Nearly optimal lattice simulation by product formulas.  Phys. Rev. Lett., 123(5):050503, 2019.  doi:10.1103/PhysRevLett.123.050503.
  * [17] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu.  Theory of trotter error with commutator scaling.  Phys. Rev. X, 11:011020, 2021.  doi:10.1103/PhysRevX.11.011020.
  * [18] B. F. Curchod and T. J. Martínez.  Ab initio nonadiabatic quantum molecular dynamics.  Chemical reviews, 118(7):3305–3336, 2018.  doi:10.1021/acs.chemrev.7b00423.
  * [19] C. M. A. Dantas, I. A. Pedrosa, and B. Baseia.  Harmonic oscillator with time-dependent mass and frequency and a perturbative potential.  Physical Review A, 45(3):1320–1324, 1992.  doi:10.1103/PhysRevA.45.1320.
  * [20] S. Descombes and M. Thalhammer.  An exact local error representation of exponential operator splitting methods for evolutionary problems and applications to linear Schrödinger equations in the semi-classical regime.  BIT Numer. Math., 50(4):729–749, 2010.  doi:10.1007/s10543-010-0282-4.
  * [21] D. Dong and I. R. Petersen.  Quantum control theory and applications: a survey.  IET Control Theory & Applications, 4(12):2651–2671, 2010.  doi:10.1049/iet-cta.2009.0508.
  * [22] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser.  Quantum computation by adiabatic evolution.  2000\.  arXiv:quant-ph/0001106.
  * [23] M. Feng.  Complete solution of the Schrödinger equation for the time-dependent linear potential.  Physical Review A, 64(3):034101 EP –, 2001.  doi:10.1103/PhysRevA.64.034101.
  * [24] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204, 2019.  doi:10.1145/3313276.3316366.
  * [25] E. Hairer, C. Lubich, and G. Wanner.  Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, volume 31.  Springer, 2006.  doi:10.1007/3-540-30666-8.
  * [26] E. Hairer, S. P. Nørsett, and G. Wanner.  Solving ordinary differential equation I: nonstiff problems, volume 8.  Springer, 1987.  doi:10.1007/978-3-540-78862-1.
  * [27] M. Hochbruck and C. Lubich.  On Magnus integrators for time-dependent Schrödinger equations.  SIAM J. Numer. Anal., 41(3):945–963, 2003.  doi:10.1137/S0036142902403875.
  * [28] J. Huyghebaert and H. De Raedt.  Product formula methods for time-dependent Schrödinger problems.  J. Phys. A, 23(24):5777–5793, 1990.  doi:10.1088/0305-4470/23/24/019.
  * [29] T. Jahnke and C. Lubich.  Error bounds for exponential operator splittings.  BIT, 40(4):735–744, 2000.  doi:10.1023/A:1022396519656.
  * [30] J.-Y. Ji, J. K. Kim, S. P. Kim, and K.-S. Soh.  Exact wave functions and nonadiabatic Berry phases of a time-dependent harmonic oscillator.  Physical Review A, 52(4):3352–3355, 1995.  doi:10.1103/PhysRevA.52.3352.
  * [31] I. D. Kivlichan, N. Wiebe, R. Babbush, and A. Aspuru-Guzik.  Bounding the costs of quantum simulation of many-body physics in real space.  J. Phys. A Math. Theor., 50:305301, 2017.  doi:10.1088/1751-8121/aa77b8.
  * [32] A. W. Knapp.  Basic Real Analysis.  Springer Science & Business Media, 2005.  doi:10.1007/0-8176-4441-5.
  * [33] R. J. LeVeque.  Finite difference methods for ordinary and partial differential equations.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2007.  Steady-state and time-dependent problems.  doi:10.1137/1.9780898717839.
  * [34] S. Lloyd.  Universal quantum simulators.  Science, pages 1073–1078, 1996.  doi:10.1126/science.273.5278.1073.
  * [35] G. H. Low.  Hamiltonian simulation with nearly optimal dependence on spectral norm.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 491–502, 2019.  doi:10.1145/3313276.3316386.
  * [36] G. H. Low and I. L. Chuang.  Optimal Hamiltonian simulation by quantum signal processing.  Phys. Rev. Lett., 118:010501, 2017.  doi:10.1103/PhysRevLett.118.010501.
  * [37] G. H. Low and N. Wiebe.  Hamiltonian simulation in the interaction picture.  2019\.  arXiv:1805.00675.
  * [38] Y. Maday, J. Salomon, and G. Turinici.  Monotonic time-discretized schemes in quantum control.  Numer. Math., 103(2):323–338, 2006.  doi:10.1007/s00211-006-0678-x.
  * [39] A. B. Magann, M. D. Grace, H. A. Rabitz, and M. Sarovar.  Digital quantum simulation of molecular dynamics and control.  2020\.  arXiv:2002.12497.
  * [40] R. Montalto.  On the growth of Sobolev norms for a class of linear Schrödinger equations on the torus with superlinear dispersion.  Asymptotic Analysis, 108:85–114, 2018.  doi:10.3233/ASY-181470.
  * [41] M. A. Nielsen and I. Chuang.  Quantum computation and quantum information, 2000.  doi:10.1017/CBO9780511976667.
  * [42] M. A. Nielsen, M. R. Dowling, M. Gu, and A. C. Doherty.  Optimal control, geometry, and quantum computing.  Phys. Rev. A, 73(6):062323, 2006.  doi:10.1103/PhysRevA.73.062323.
  * [43] S. Pang and A. N. Jordan.  Optimal adaptive control for quantum metrology with time-dependent hamiltonians.  Nature communications, 8(1):1–9, 2017.  doi:10.1038/ncomms14695.
  * [44] I. A. Pedrosa.  Exact wave functions of a harmonic oscillator with time-dependent mass and frequency.  Physical Review A, 55(4):3219–3221, 1997.  doi:10.1103/PhysRevA.55.3219.
  * [45] I. A. Pedrosa, G. P. Serra, and I. Guedes.  Wave functions of a time-dependent harmonic oscillator with and without a singular perturbation.  Physical Review A, 56(5):4300–4303, 1997.  doi:10.1103/PhysRevA.56.4300.
  * [46] D. Poulin, A. Qarry, R. Somma, and F. Verstraete.  Quantum simulation of time-dependent Hamiltonians and the convenient illusion of Hilbert space.  Phys. Rev. Lett., 106(17):170501, 2011.  doi:10.1103/PhysRevLett.106.170501.
  * [47] J. Roland and N. J. Cerf.  Quantum search by local adiabatic evolution.  Phys. Rev. A, 65(4):042308, 2002.  doi:10.1103/PhysRevA.65.042308.
  * [48] E. Runge and E. K. U. Gross.  Density-functional theory for time-dependent systems.  Phys. Rev. Lett., 52:997, 1984.  doi:10.1103/PhysRevLett.52.997.
  * [49] B. Şahinoğlu and R. D. Somma.  Hamiltonian simulation in the low energy subspace.  2020\.  arXiv:2006.02660.
  * [50] A. Schulze-Halberg.  Form-Preserving Transformations of Time-Dependent Schrödinger Equation with Time- and Position-Dependent Mass.  Communications in Theoretical Physics, 43(4):657–665, 2005.  doi:10.1088/0253-6102/43/4/017.
  * [51] Y. Su, H.-Y. Huang, and E. T. Campbell.  Nearly tight trotterization of interacting electrons.  2020\.  arXiv:2012.09194.
  * [52] M. Suzuki.  General decomposition theory of ordered exponentials.  Proc. Japan Acad., 69:161–166, 1993.  doi:10.2183/pjab.69.161.
  * [53] M. Thalhammer.  High-order exponential operator splitting methods for time-dependent Schrödinger equations.  SIAM J. Numer. Anal., 46(4):2022–2038, 2008.  doi:10.1137/060674636.
  * [54] J. W. Thomas.  Numerical partial differential equations: finite difference methods, volume 22 of Texts in Applied Mathematics.  Springer-Verlag, New York, 1995.  doi:10.1007/978-1-4899-7278-1.
  * [55] D. Wecker, M. B. Hastings, N. Wiebe, B. K. Clark, C. Nayak, and M. Troyer.  Solving strongly correlated electron models on a quantum computer.  Phys. Rev. A, 92:062318, 2015.  doi:10.1103/PhysRevA.92.062318.
  * [56] N. Wiebe, D. Berry, P. Høyer, and B. C. Sanders.  Higher order decompositions of ordered operator exponentials.  J. Phys. A, 43(6):065203, 2010.  doi:10.1088/1751-8113/43/6/065203.
  * [57] W. Zhu and H. Rabitz.  A rapid monotonically convergent iteration algorithm for quantum optimal control over the expectation value of a positive definite operator.  J. Chem. Phys., 109(2):385–391, 1998.  doi:10.1063/1.476575.

##  Appendix 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=\mathcal{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 $\epsilon$
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

$\|H_{1}\|=\mathcal{O}(n^{2}),\quad\|H_{2}\|=\mathcal{O}(1),\quad\|[H_{1},H_{2}]\|=\mathcal{O}(n),$
$\|[H_{1},[H_{1},H_{2}]]\|=\mathcal{O}(n^{2}),\quad\|[H_{2},[H_{2},H_{1}]]\|=\mathcal{O}(n),$

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
$\mathcal{O}(\epsilon^{-1/2})$. Plugging this back into the complexity
estimates leads to the overall scaling in terms of $\epsilon$, 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

$\frac{h^{3}}{12}\|[H_{2},[H_{2},H_{1}]]\|+\frac{h^{3}}{24}\|[H_{1},[H_{1},H_{2}]]\|,$

thus the global Trotter error is bounded by

$\left(\frac{1}{12}\|[H_{2},[H_{2},H_{1}]]\|+\frac{1}{24}\|[H_{1},[H_{1},H_{2}]]\|\right)\frac{T^{3}}{L^{2}}=\mathcal{O}\left(\frac{n^{2}}{L^{2}}\right).$

To bound this by $\epsilon$, it suffices to choose

$L=\mathcal{O}\left(\frac{n}{\epsilon^{1/2}}\right)=\mathcal{O}\left(\frac{1}{\epsilon}\right).$

[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

$\mathcal{O}\left(h^{2}\left(\|H_{1}\vec{v}\|_{\star}+\|D_{1}\vec{v}\|_{\star}+\|\vec{v}\|_{\star}\right)\right)=\mathcal{O}\left(\frac{1}{L^{2}}\left(\|H_{1}\vec{v}\|_{\star}+\|\vec{v}\|_{\star}\right)\right).$

We remark that [29] does not track explicitly the dependence on $T$. Noticing
that $\|H_{1}\vec{v}\|_{\star}$ and $\|\vec{v}\|_{\star}$ are independent of
$\epsilon$ and $n$ (shown in the proof of Lemma 11), the number of required
Trotter steps scales

$L=\mathcal{O}\left(\frac{1}{\epsilon^{1/2}}\right).$

#### 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
$\epsilon$-approximation of the exact evolution operator is

$L=\mathcal{O}\left(\frac{\widetilde{\alpha}_{\text{comm}}^{1/p}T^{1+1/p}}{\epsilon^{1/p}}\right),$

where

$\widetilde{\alpha}_{\text{comm}}=\sum_{\gamma_{1},\cdots,\gamma_{p+1}=1}^{2}\|[H_{\gamma_{p+1}},\cdots[H_{\gamma_{2}},H_{\gamma_{1}}]]\|.$

Straightforward bounds for these $p$-th nested commutators are that

$\|[H_{\gamma_{p+1}},\cdots[H_{\gamma_{2}},H_{\gamma_{1}}]]\|=\mathcal{O}(\|H_{1}\|^{p-2}\|[H_{1},[H_{1},H_{2}]]\|)=\mathcal{O}(n^{2p-2}),$

which results in

$L=\mathcal{O}\left(\frac{n^{2-2/p}}{\epsilon^{1/p}}\right)=\mathcal{O}\left(\frac{1}{\epsilon}\right).$

Notice that the scaling of $\epsilon$ 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. If the potential $V(x)$ has
higher regularity, we expect better bounds to exist for general nested
commutators, just like the case of $[H_{1},H_{2}]$ and
$[H_{1},[H_{1},H_{2}]]$. In particular, although we do not present complete
proof in this paper, a continuous analog as well as discretization under
Fourier basis suggests that the norm of $p$-th order nested commutator
$\|[H_{1},\cdots,[H_{1},H_{2}]]\|$ is bounded by $\mathcal{O}(\|D_{1}^{p}\|)$
if $V(x)$ is $(2p)$-th order continuously differentiable. In that case the
complexity can be improved to $L=\mathcal{O}(n/\epsilon^{1/p})$, although
there is still a linear dependence on $n$.

#### Truncated Taylor series

[5, Theorem 1] shows that to obtain an $\epsilon$-approximation of the exact
evolution operator using truncated Taylor series, the query complexity is

$\mathcal{O}\left(d^{2}\|H\|_{\max}\frac{\log(d\|H\|_{\max}/\epsilon)}{\log\log(d\|H\|_{\max}/\epsilon)}\right).$

Here $d$ is the sparsity of the Hamiltonian, $\|H\|_{\max}$ denotes the
largest matrix element of $H$ in absolute value. Notice that
$\|H_{1}\|_{\max}=\mathcal{O}(n^{2})$ since every non-zero entry of $H_{1}$ is
either $n^{2}$ or $(-2n^{2})$, and $\|H_{2}\|_{\max}=\mathcal{O}(1)$, we have
$\|H\|_{\max}=\mathcal{O}(n^{2})$. Therefore the query complexity becomes

$\mathcal{O}\left(n^{2}\frac{\log(n^{2}/\epsilon)}{\log\log(n^{2}/\epsilon)}\right)=\widetilde{\mathcal{O}}\left(n^{2}\right)=\widetilde{\mathcal{O}}\left(\frac{1}{\epsilon}\right).$

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$ than
$\mathcal{O}(1/\epsilon)$ [31, Theorem 4], which results in a total complexity
$\widetilde{\mathcal{O}}(1/\epsilon^{2})$ [31, Theorem 3 & 4]. The scaling can
be improved if $V(x)$ becomes smoother.

#### Quantum signal processing

[36, Theorem 3] proposes a quantum signal processing approach for time-
independent Hamiltonian simulation with optimal query complexity in all
parameters, which is

$\mathcal{O}\left(d\|H\|_{\max}+\frac{\log(1/\epsilon)}{\log\log(1/\epsilon)}\right).$

Here $d$ is the sparsity of the Hamiltonian, $\|H\|_{\max}$ denotes the
largest matrix element of $H$ in absolute value. Notice that
$\|H_{1}\|_{\max}=\mathcal{O}(n^{2})$ since every non-zero entry of $H_{1}$ is
either $n^{2}$ or $(-2n^{2})$, and $\|H_{2}\|_{\max}=\mathcal{O}(1)$, we have
$\|H\|_{\max}=\mathcal{O}(n^{2})$. Therefore the query complexity becomes

$\mathcal{O}\left(n^{2}+\frac{\log(n^{2}/\epsilon)}{\log\log(n^{2}/\epsilon)}\right)={\mathcal{O}}\left(n^{2}+\log(n^{2}/\epsilon)\right)=\mathcal{O}\left(\frac{1}{\epsilon}+\log(1/\epsilon)\right)=\mathcal{O}\left(\frac{1}{\epsilon}\right).$

#### Interaction picture

[37, Theorem 7] shows that by applying truncated Dyson series to simulate
time-independent Hamiltonian $H_{1}+H_{2}$ in the interaction picture rather
than the Schrödinger picture, it requires

$\mathcal{O}\left(\|H_{2}\|\frac{\log(\|H_{2}\|/\epsilon)}{\log\log(\|H_{2}\|/\epsilon)}\right)$

queries to $H_{2}$ and

$\mathcal{O}\left(\|H_{2}\|\frac{\log(\|H_{2}\|/\epsilon)}{\log\log(\|H_{2}\|/\epsilon)}\log\left(\frac{\|H_{1}\|+\|H_{2}\|}{\epsilon}\right)\right)$

queries to the unitary time evolution $e^{-\mathrm{i}sH_{1}}$. Therefore the
query complexity is logarithmic in $n$ and thus the scaling in terms of
$\epsilon$ is still poly-logarithmic. Note that the number of time steps is
included in the oracle HAM-T and scales as $\mathcal{O}(\left\lVert
H_{1}\right\rVert)$ [37, Lemma 6].

###  A.2 Time-dependent schemes

#### Time-dependent second order Trotter formulae

[28, Eq. (A12-A14)] show that for generalized second-order Trotter formula
applied to the model Eq. 3 with time-independent mass and time-dependent
frequency (in particular, $f_{2}(t)H_{2}$ and $f_{2}(s)H_{2}$ commute for any
$t$ and $s$), the one-step local Trotter error scales as

$\mathcal{O}\left(h^{3}\left(\|[H_{1},H_{2}]\|+\|[H_{1},[H_{1},H_{2}]]\|+\|[H_{2},[H_{2},H_{1}]]\|\right)\right),$

thus the global error scales

$\mathcal{O}\left(h^{2}\left(\|[H_{1},H_{2}]\|+\|[H_{1},[H_{1},H_{2}]]\|+\|[H_{2},[H_{2},H_{1}]]\|\right)\right)=\mathcal{O}\left(\frac{n^{2}}{L^{2}}\right).$

To bound this by $\epsilon$, it suffices to choose

$L=\mathcal{O}\left(\frac{n}{\epsilon^{1/2}}\right)=\mathcal{O}\left(\frac{1}{\epsilon}\right).$

The second order complexity estimate from [56] is a special case of their
general high order result. We will show the general case later.

[55, Appendix A] proves an improved operator norm error bound for the second
order standard Trotter formula. The one-step local Trotter error is bounded by

$\left(\frac{1}{24}\sup\|H^{\prime\prime}(s)\|+\frac{1}{12}\sup\|[H^{\prime}(s),H(s)]\|+\sup\|[H_{1}(s),[H_{1}(s),H_{2}(s)]\|+\|[H_{2}(s),[H_{2}(s),H_{1}(s)]\|\right)h^{3},$

thus the global error scales

$\mathcal{O}\left(h^{2}\left(\|H_{1}\|+\|[H_{1},H_{2}]\|+\|[H_{1},[H_{1},H_{2}]]\|+\|[H_{2},[H_{2},H_{1}]]\|\right)\right)=\mathcal{O}\left(\frac{n^{2}}{L^{2}}\right).$

To bound this by $\epsilon$, it suffices to choose

$L=\mathcal{O}\left(\frac{n}{\epsilon^{1/2}}\right)=\mathcal{O}\left(\frac{1}{\epsilon}\right).$

#### Time-dependent high order Trotter formula

[56, Theorem 1] proves that, to simulate a system with Hamiltonian
$H(t)=\sum_{j=1}^{m}H_{j}(t)$ within operator spectral norm error $\epsilon$
using a $2k$-th order standard Trotter formula, the total number of
exponentials is

$2m5^{k-1}\left\lceil 5k\Lambda
T\left(\frac{5}{3}\right)^{k}\left(\frac{\Lambda
T}{\epsilon}\right)^{1/(2k)}\right\rceil$

where

$\Lambda=\sup_{p=0,1,\cdots,2k}\left(\sup_{t}\left(\sum_{j=1}^{m}\|\partial_{t}^{p}H_{j}(t)\|\right)^{1/(p+1)}\right).$

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 $\epsilon$ into the big-$\mathcal{O}$
notation, in the case of the Schrödinger equation with a time-dependent
effective mass, the total number of Trotter steps becomes

$\mathcal{O}\left(\Lambda\left(\frac{\Lambda}{\epsilon}\right)^{1/2k}\right).$

It remains to estimate the scaling of $\Lambda$. By noticing
$\partial_{t}^{p}H_{j}(t)=f_{j}^{(p)}(t)H_{j}$, we obtain that
$\left(\sum_{j=1}^{m}\|\partial_{t}^{p}H_{j}(t)\|\right)$ is dominated by
$H_{1}=\mathcal{O}(n^{2})$, and

$\Lambda=\mathcal{O}\left(\sup_{p=0,1,\cdots,2k}(n^{2})^{1/(p+1)}\right)=\mathcal{O}\left(n^{2}\right).$

Therefore the total number of Trotter steps becomes

$\mathcal{O}\left(\frac{n^{2+1/k}}{\epsilon^{1/(2k)}}\right)=\mathcal{O}\left(\frac{1}{\epsilon^{1+1/k}}\right).$

#### Truncated Dyson series

[37, Theorem 9] shows that to obtain an $\epsilon$ approximation of the exact
evolution operator with success probability at least $(1-\epsilon)$ using
truncated Dyson series method, the query complexity is

$\mathcal{O}\left(d\|H\|_{\max,\infty}T\frac{\log(d\|H\|_{\max,\infty}T/\epsilon)}{\log\log(d\|H\|_{\max,\infty}T/\epsilon)}\right).$

Here $d$ is the sparsity of the Hamiltonian, and
$\|H\|_{\max,\infty}=\sup_{t\in[0,T]}\|H(t)\|_{\max}$, 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}=\mathcal{O}(n^{2})$ because
every non-zero entry of $H_{1}$ is either $n^{2}$ or $(-2n^{2})$, we have
$\|H(t)\|_{\max,\infty}=\mathcal{O}(\|H_{1}\|_{\max})=\mathcal{O}(n^{2})$,
then the query complexity becomes

$\mathcal{O}\left(n^{2}\frac{\log(n^{2}/\epsilon)}{\log\log(n^{2}/\epsilon)}\right)=\widetilde{\mathcal{O}}\left(n^{2}\right)=\widetilde{\mathcal{O}}\left(\frac{1}{\epsilon}\right).$

#### Rescaled Dyson series

[8, Theorem 10] shows that to obtain an $\epsilon$ approximation of the exact
evolution operator using rescaled Dyson series method, the query complexity is

$\mathcal{O}\left(d\|H\|_{\max,1}\frac{\log(d\|H\|_{\max,1}/\epsilon)}{\log\log(d\|H\|_{\max,1}/\epsilon)}\right).$

Here $d$ is the sparsity of the Hamiltonian,
$\|H\|_{\max,1}=\int_{0}^{T}\|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}=\mathcal{O}(n^{2})$ because every non-
zero entry of $H_{1}$ is either $n^{2}$ or $(-2n^{2})$, we have
$\|H\|_{\max,1}=\mathcal{O}(n^{2})$. Therefore the query complexity becomes

$\mathcal{O}\left(n^{2}\frac{\log(n^{2}/\epsilon)}{\log\log(n^{2}/\epsilon)}\right)=\widetilde{\mathcal{O}}\left(n^{2}\right)=\widetilde{\mathcal{O}}\left(\frac{1}{\epsilon}\right).$

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.

##  Appendix 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

$\partial_{h}U(h,0)=(-\mathrm{i}f_{1}(h)H_{1}-\mathrm{i}f_{2}(h)H_{2})U(h,0),$
(76)

and its numerical analogs $\partial_{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

$\begin{split}\partial_{h}U_{g,1}(h,0)&=-\mathrm{i}f_{2}(h)H_{2}U_{g,1}(h,0)+\exp\left(-\mathrm{i}\int_{0}^{h}f_{2}(s)dsH_{2}\right)\left(-\mathrm{i}f_{1}(h)H_{1}\right)\exp\left(-\mathrm{i}\int_{0}^{h}f_{1}(s)dsH_{1}\right)\\\
&=\left(-\mathrm{i}f_{1}(h)H_{1}-\mathrm{i}f_{2}(h)H_{2}\right)U_{g,1}(h,0)\\\
&\quad+\exp\left(-\mathrm{i}\int_{0}^{h}f_{2}(s)dsH_{2}\right)E_{g,1}(h)\exp\left(-\mathrm{i}\int_{0}^{h}f_{1}(s)dsH_{1}\right)\end{split}$
(77)

where $E_{g,1}(h)$ is defined as

$E_{g,1}(h)=\mathrm{i}f_{1}(h)\left[\exp\left(\operatorname{ad}_{\mathrm{i}\int_{0}^{h}f_{2}(s^{\prime})ds^{\prime}H_{2}}\right)H_{1}-H_{1}\right].$
(78)

By applying Lemma 2 to Eq. 76 and Eq. 77, one obtains

$U_{g,1}(h,0)=U(h,0)+\int_{0}^{h}U(h,s)\exp\left(-\mathrm{i}\int_{0}^{s}f_{2}(s^{\prime})ds^{\prime}H_{2}\right)E_{g,1}(s)\exp\left(-\mathrm{i}\int_{0}^{s}f_{1}(s^{\prime})ds^{\prime}H_{1}\right)ds.$
(79)

The representation of $E_{g,1}$, by Taylor’s theorem (Lemma 1), reads

$\begin{split}E_{g,1}(h)&=\int_{0}^{h}f_{1}(h)f_{2}(s)\left(\exp\left(\operatorname{ad}_{\mathrm{i}\int_{0}^{h}f_{2}(s^{\prime})ds^{\prime}H_{2}}\right)([H_{1},H_{2}])\right)ds.\end{split}$
(80)

∎

###### Proof of Lemma 3.

One starts by taking derivative of $U_{s,1}(h,0)$ with respect to $h$, which
reads

$\begin{split}\partial_{h}U_{s,1}(h,0)&=\left(-\mathrm{i}f_{2}(h)H_{2}-\mathrm{i}hf_{2}^{\prime}(h)H_{2}\right)\exp\left(-\mathrm{i}hf_{2}(h)H_{2}\right)\exp\left(-\mathrm{i}hf_{1}(h)H_{1}\right)\\\
&\quad+\exp\left(-\mathrm{i}hf_{2}(h)H_{2}\right)\left(-\mathrm{i}f_{1}(h)H_{1}-\mathrm{i}hf_{1}^{\prime}(h)H_{1}\right)\exp\left(-\mathrm{i}hf_{1}(h)H_{1}\right)\\\
&=\left(-\mathrm{i}f_{2}(h)H_{2}-\mathrm{i}f_{1}(h)H_{1}\right)U_{s,1}(h,0)\\\
&\quad+\exp\left(-\mathrm{i}hf_{2}(h)H_{2}\right)E_{s,1}(h)\exp\left(-\mathrm{i}hf_{1}(h)H_{1}\right),\end{split}$
(81)

where $E_{s,1}(h)$ is defined as

$\begin{split}E_{s,1}(h)=\mathrm{i}f_{1}(h)\left[\exp\left(\operatorname{ad}_{\mathrm{i}hf_{2}(h)H_{2}}\right)H_{1}-H_{1}\right]-\mathrm{i}hf_{1}^{\prime}(h)H_{1}-\mathrm{i}hf_{2}^{\prime}(h)H_{2}.\end{split}$
(82)

By applying Lemma 2 to Eq. 76 and Eq. 81, one has

$U_{s,1}(h,0)=U(h,0)+\int_{0}^{h}U(h,s)\exp\left(-\mathrm{i}sf_{2}(s)H_{2}\right)E_{s,1}(s)\exp\left(-\mathrm{i}sf_{1}(s)H_{1}\right)ds.$
(83)

It remains to derive the representation of $E_{s,1}$. The representation of
$E_{s,1}$ can be derived from Taylor’s theorem up to first-order. By Lemma 1

$\begin{split}&\quad\exp\left(\operatorname{ad}_{\mathrm{i}hf_{2}(h)H_{2}}\right)H_{1}-H_{1}\\\
&=\int_{0}^{h}\mathrm{i}f_{2}(s)\left(\exp\left(\operatorname{ad}_{\mathrm{i}sf_{2}(s)H_{2}}\right)([H_{2},H_{1}])\right)ds+\int_{0}^{h}\mathrm{i}sf_{2}^{\prime}(s)\left(\exp\left(\operatorname{ad}_{\mathrm{i}sf_{2}(s)H_{2}}\right)([H_{2},H_{1}])\right)ds.\end{split}$
(84)

Therefore, one has

$\begin{split}E_{s,1}(h)&=\int_{0}^{h}f_{1}(h)f_{2}(s)\left(\exp\left(\operatorname{ad}_{\mathrm{i}sf_{2}(s)H_{2}}\right)([H_{1},H_{2}])\right)ds-\mathrm{i}hf_{1}^{\prime}(h)H_{1}-\mathrm{i}hf_{2}^{\prime}(h)H_{2}\\\
&\quad+\int_{0}^{h}sf_{1}(h)f_{2}^{\prime}(s)\left(\exp\left(\operatorname{ad}_{\mathrm{i}sf_{2}(s)H_{2}}\right)([H_{1},H_{2}])\right)ds.\end{split}$
(85)

∎

Before proceeding, we first define the following quantities needed in the
error representations of the second order standard and generalized Trotter
formulae

$\begin{split}E_{s,2}(h)&=\mathrm{i}\int_{0}^{h}f_{1}^{\prime\prime}(s)(h-s)H_{1}ds-\frac{\mathrm{i}}{8}\int_{0}^{h}f_{1}^{\prime\prime}(s/2)(2h-s)H_{1}ds-\frac{\mathrm{i}}{4}\int_{0}^{h}f_{2}^{\prime\prime}(s/2)(2h-s)H_{2}ds\\\
&\quad-\frac{\mathrm{i}}{8}\int_{0}^{h}\left[f_{1}^{\prime\prime}(s/2)\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)H_{1}\right](2h-s)ds\\\
&\quad+\frac{1}{4}\int_{0}^{h}\left[f_{1}^{\prime}(s/2)f_{2}(s/2)\left(\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)H_{1}\right)\right]hds\\\
&\quad+\mathrm{i}\int_{0}^{h}\left[f_{2}^{\prime\prime}(s)\exp\left(\operatorname{ad}_{i\frac{s}{2}f_{1}(s/2)H_{1}}\right)H_{2}\right](h-s)ds\\\
&\quad+\frac{1}{2}\int_{0}^{h}(f_{1}^{\prime}(s/2)f_{2}(s/2)+f_{1}(s/2)f_{2}^{\prime}(s/2))\left(\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)[H_{1},H_{2}]\right)(h-s)ds\\\
&\quad-\frac{1}{2}\int_{0}^{h}\left(f_{2}^{\prime}(s)f_{1}(s/2)+f_{2}(s)f_{1}^{\prime}(s/2)\right)\left(\exp\left(\operatorname{ad}_{i\frac{s}{2}f_{1}(s/2)H_{1}}\right)[H_{1},H_{2}]\right)(h-s)ds\\\
&\quad+\frac{\mathrm{i}}{2}\int_{0}^{h}\left[f_{1}(s/2)f_{2}^{2}(s/2)\left(\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)[H_{2},[H_{1},H_{2}]]\right)\right](h-s)ds\\\
&\quad-\frac{\mathrm{i}}{4}\int_{0}^{h}\left[f_{2}(s)f_{1}^{2}(s/2)\left(\exp\left(\operatorname{ad}_{i\frac{s}{2}f_{1}(s/2)H_{1}}\right)[H_{1},[H_{1},H_{2}]]\right)\right](h-s)ds\\\
&\quad+\frac{1}{8}\int_{0}^{h}\left[f_{1}^{\prime}(s/2)f_{2}^{\prime}(s/2)\left(\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)H_{1}\right)\right]shds\\\
&\quad+\frac{1}{4}\int_{0}^{h}\left(f_{1}^{\prime}(s/2)f_{2}^{\prime}(s/2)+\frac{1}{2}f_{1}(s/2)f_{2}^{\prime\prime}(s/2)\right)\left(\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)[H_{1},H_{2}]\right)s(h-s)ds\\\
&\quad-\frac{1}{4}\int_{0}^{h}\left(f_{2}^{\prime}(s)f_{1}^{\prime}(s/2)+\frac{1}{2}f_{2}(s)f_{1}^{\prime\prime}(s/2)\right)\exp\left(\operatorname{ad}_{i\frac{s}{2}f_{1}(s/2)H_{1}}\right)[H_{1},H_{2}]s(h-s)ds\\\
&\quad+\frac{\mathrm{i}}{2}\int_{0}^{h}\left[f_{1}(s/2)f_{2}(s/2)f_{2}^{\prime}(s/2)\left(\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)[H_{2},[H_{1},H_{2}]]\right)\right]s(h-s)ds\\\
&\quad-\frac{\mathrm{i}}{4}\int_{0}^{h}\left[f_{2}(s)f_{1}(s/2)f_{1}^{\prime}(s/2)\left(\exp\left(\operatorname{ad}_{i\frac{s}{2}f_{1}(s/2)H_{1}}\right)[H_{1},[H_{1},H_{2}]]\right)\right]s(h-s)ds\\\
&\quad+\frac{\mathrm{i}}{8}\int_{0}^{h}\left[f_{1}(s/2)f_{2}^{\prime
2}(s/2)\left(\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)[H_{2},[H_{1},H_{2}]]\right)\right]s^{2}(h-s)ds\\\
&\quad-\frac{\mathrm{i}}{16}\int_{0}^{h}\left[f_{2}(s)f_{1}^{\prime
2}(s/2)\left(\exp\left(\operatorname{ad}_{i\frac{s}{2}f_{1}(s/2)H_{1}}\right)[H_{1},[H_{1},H_{2}]]\right)\right]s^{2}(h-s)ds,\end{split}$
(86)

and

$\begin{split}E_{g,2}(h)&=-\frac{h}{2}f_{1}(0)\int_{0}^{h}f_{2}^{\prime}(s)ds[H_{1},H_{2}]+\frac{h}{4}f_{2}(0)\int_{0}^{h}f_{1}^{\prime}(s/2)ds[H_{1},H_{2}]\\\
&\quad-
f_{2}(h)\int_{0}^{h}(f_{1}^{\prime}(s)-\frac{1}{4}f_{1}^{\prime}(s/2))\left(\exp\left(\operatorname{ad}_{\mathrm{i}\int_{s/2}^{s}f_{1}(s^{\prime})ds^{\prime}H_{1}}\right)[H_{1},H_{2}]\right)(h-s)ds\\\
&\quad+\frac{1}{2}f_{1}(h/2)\int_{0}^{h}f_{2}^{\prime}(s)\left(\exp\left(\operatorname{ad}_{-\mathrm{i}\int_{0}^{s}f_{2}(s^{\prime})ds^{\prime}H_{2}}\right)[H_{1},H_{2}]\right)(h-s)ds\\\
&\quad-\mathrm{i}f_{2}(h)\int_{0}^{h}\left(f_{1}(s)-\frac{1}{2}f_{1}(s/2)\right)^{2}\left(\exp\left(\operatorname{ad}_{\mathrm{i}\int_{s/2}^{s}f_{1}(s^{\prime})ds^{\prime}H_{1}}\right)[H_{1},[H_{1},H_{2}]]\right)(h-s)ds\\\
&\quad+\frac{\mathrm{i}}{2}f_{1}(h/2)\int_{0}^{h}f_{2}^{2}(s)\left(\exp\left(\operatorname{ad}_{-\mathrm{i}\int_{0}^{s}f_{2}(s^{\prime})ds^{\prime}H_{2}}\right)[H_{2},[H_{2},H_{1}]]\right)(h-s)ds.\end{split}$
(87)

###### Proof of Lemma 5.

One first compute the derivative with respect to $h$ of $U_{s,2}$

$\begin{split}\partial_{h}U_{s,2}&=\left(-\mathrm{i}\frac{1}{2}f_{1}(h/2)-\mathrm{i}\frac{h}{4}f_{1}^{\prime}(h/2)\right)H_{1}\exp\left(-\mathrm{i}\frac{h}{2}f_{1}(h/2)H_{1}\right)\exp\left(-\mathrm{i}hf_{2}(h/2)H_{2}\right)\exp\left(-\mathrm{i}\frac{h}{2}f_{1}(h/2)H_{1}\right)\\\
&\quad+\exp\left(-\mathrm{i}\frac{h}{2}f_{1}(h/2)H_{1}\right)\left(-\mathrm{i}f_{2}(h/2)-\mathrm{i}\frac{h}{2}f_{2}^{\prime}(h/2)\right)H_{2}\exp\left(-\mathrm{i}hf_{2}(h/2)H_{2}\right)\exp\left(-\mathrm{i}\frac{h}{2}f_{1}(h/2)H_{1}\right)\\\
&\quad+\exp\left(-\mathrm{i}\frac{h}{2}f_{1}(h/2)H_{1}\right)\exp\left(-\mathrm{i}hf_{2}(h/2)H_{2}\right)\left(-\mathrm{i}\frac{1}{2}f_{1}(h/2)-\mathrm{i}\frac{h}{4}f_{1}^{\prime}(h/2)\right)H_{1}\exp\left(-\mathrm{i}\frac{h}{2}f_{1}(h/2)H_{1}\right)\\\
&=\left(-\mathrm{i}f_{1}(h)H_{1}-\mathrm{i}f_{2}(h)H_{2}\right)U_{s,2}\\\
&\quad+\exp\left(-\mathrm{i}\frac{h}{2}f_{1}(h/2)H_{1}\right)E_{s,2}(h)\exp\left(-\mathrm{i}hf_{2}(h/2)H_{2}\right)\exp\left(-\mathrm{i}\frac{h}{2}f_{1}(h/2)H_{1}\right),\end{split}$
(88)

where $E_{s,2}(h)$ is defined as

$\begin{split}E_{s,2}(h)&=\mathrm{i}f_{1}(h)H_{1}-\mathrm{i}\frac{1}{2}f_{1}(h/2)H_{1}-\mathrm{i}\frac{h}{4}f_{1}^{\prime}(h/2)H_{1}-\left(\mathrm{i}\frac{1}{2}f_{1}(h/2)+\mathrm{i}\frac{h}{4}f_{1}^{\prime}(h/2)\right)\exp\left(\operatorname{ad}_{-\mathrm{i}hf_{2}(h/2)H_{2}}\right)H_{1}\\\
&\quad+\mathrm{i}f_{2}(h)\exp\left(\operatorname{ad}_{i\frac{h}{2}f_{1}(h/2)H_{1}}\right)H_{2}-\left(\mathrm{i}f_{2}(h/2)+\frac{\mathrm{i}h}{2}f_{2}^{\prime}(h/2)\right)H_{2}.\end{split}$
(89)

Similar as the proofs for first-order formulae, applying Lemma 2 gives

$U_{s,2}(h,0)=U(h,0)+\int_{0}^{h}U(h,s)\exp\left(-\mathrm{i}\frac{s}{2}f_{1}(s/2)H_{1}\right)E_{s,2}(s)\exp\left(-\mathrm{i}sf_{2}(s/2)H_{2}\right)\exp\left(-\mathrm{i}\frac{s}{2}f_{1}(s/2)H_{1}\right)ds.$
(90)

The rest of the proof follows straightforward calculations. To be exact, one
then applies the Taylor’s theorem (Lemma 1) to expand each term in $E_{s,2}$
to second-order in terms of $h$ with respect to 0. The first three terms can
be expressed as

$\mathrm{i}f_{1}(h)H_{1}=\mathrm{i}f_{1}(0)H_{1}+\mathrm{i}hf_{1}^{\prime}(0)H_{1}+\mathrm{i}\int_{0}^{h}f_{1}^{\prime\prime}(s)(h-s)H_{1}ds,$
(91)
$-\mathrm{i}\frac{1}{2}f_{1}(h/2)H_{1}=-\mathrm{i}\frac{1}{2}f_{1}(0)H_{1}-\mathrm{i}\frac{h}{4}f_{1}^{\prime}(0)H_{1}-\mathrm{i}\frac{1}{8}\int_{0}^{h}f_{1}^{\prime\prime}(s/2)(h-s)H_{1}ds,$
(92)
$-\mathrm{i}\frac{h}{4}f_{1}^{\prime}(h/2)H_{1}=-\mathrm{i}\frac{h}{4}f_{1}^{\prime}(0)H_{1}-\mathrm{i}\frac{h}{8}\int_{0}^{h}f_{1}^{\prime\prime}(s/2)H_{1}ds,$
(93)

Similarly, let us apply the Taylor theorem to the fourth term in $E_{s,2}$,
which is the sum of

$\begin{split}&-\mathrm{i}\frac{1}{2}f_{1}(h/2)\exp\left(\operatorname{ad}_{-\mathrm{i}hf_{2}(h/2)H_{2}}\right)H_{1}\\\
=&-\frac{\mathrm{i}}{2}f_{1}(0)H_{1}-\frac{\mathrm{i}h}{4}f_{1}^{\prime}(0)H_{1}+\frac{h}{2}f_{1}(0)f_{2}(0)[H_{1},H_{2}]-\frac{\mathrm{i}}{2}\int_{0}^{h}ds\frac{1}{4}f_{1}^{\prime\prime}(s/2)\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)H_{1}(h-s)\\\
&-\frac{\mathrm{i}}{2}\int_{0}^{h}dsf_{1}(s/2)\left(\mathrm{i}f_{2}(s/2)+\frac{\mathrm{i}s}{2}f_{2}^{\prime}(s/2)\right)^{2}\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)[H_{2},[H_{1},H_{2}]](h-s)\\\
&-\frac{\mathrm{i}}{2}\int_{0}^{h}ds\left(\mathrm{i}(f_{1}^{\prime}(s/2)f_{2}(s/2)+f_{1}(s/2)f_{2}^{\prime}(s/2))+\frac{\mathrm{i}s}{2}f_{1}^{\prime}(s/2)f_{2}^{\prime}(s/2)+\frac{\mathrm{i}s}{4}f_{1}(s/2)f_{2}^{\prime\prime}(s/2)\right)\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\times\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)[H_{1},H_{2}](h-s),\end{split}$
(94)

and

$\begin{split}&-\mathrm{i}\frac{h}{4}f_{1}^{\prime}(h/2)\exp\left(\operatorname{ad}_{-\mathrm{i}hf_{2}(h/2)H_{2}}\right)H_{1}\\\
=&-\frac{\mathrm{i}h}{4}f_{1}^{\prime}(0)H_{1}-\frac{\mathrm{i}h}{4}\int_{0}^{h}ds\frac{1}{2}f_{1}^{\prime\prime}(s/2)\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)H_{1}\\\
&-\frac{\mathrm{i}h}{4}\int_{0}^{h}ds\mathrm{i}f_{1}^{\prime}(s/2)\left(f_{2}(s/2)+\frac{s}{2}f_{2}^{\prime}(s/2)\right)\exp\left(\operatorname{ad}_{-\mathrm{i}sf_{2}(s/2)H_{2}}\right)H_{1},\end{split}$
(95)