Strong Error Bounds for Trotter & Strang-Splittings and Their Implications for Quantum Chemistry
Daniel Burgarth
Department Physik, Friedrich-Alexander-Universität Erlangen-Nürnberg, Staudtstraße 7, 91058 Erlangen, Germany
Center for Engineered Quantum Systems, Macquarie University, 2109 NSW, Australia Paolo Facchi
Dipartimento di Fisica, Università di Bari, I-70126 Bari, Italy
INFN, Sezione di Bari, I-70126 Bari, Italy
Alexander Hahn
Center for Engineered Quantum Systems, Macquarie University, 2109 NSW, Australia
Mattias Johnsson
Center for Engineered Quantum Systems, Macquarie University, 2109 NSW, Australia
Kazuya Yuasa
Department of Physics, Waseda University, Tokyo 169-8555, Japan

Efficient error estimates for the Trotter product formula are central in quantum computing, mathematical physics, and numerical simulations. However, the Trotter error's dependency on the input state and its application to unbounded operators remains unclear. Here, we present a general theory for error estimation, including higher-order product formulas, with explicit input state dependency. Our approach overcomes two limitations of the existing operator-norm estimates in the literature. First, previous bounds are too pessimistic as they quantify the worst-case scenario. Second, previous bounds become trivial for unbounded operators and cannot be applied to a wide class of Trotter scenarios, including atomic and molecular Hamiltonians.  Our method enables analytical treatment of Trotter errors in chemistry simulations, illustrated through a case study on the hydrogen atom. Our findings reveal: (i) for states with fat-tailed energy distribution, such as low-angular-momentum states of the hydrogen atom, the Trotter error scales worse than expected (sublinearly) in the number of Trotter steps; (ii) certain states do not admit an advantage in the scaling from higher-order Trotterization, and thus, the higher-order Trotter hierarchy breaks down for these states, including the hydrogen atom's ground state; (iii) the scaling of higher-order Trotter bounds might depend on the order of the Hamiltonians in the Trotter product for states with fat-tailed energy distribution. Physically, the enlarged Trotter error is caused by the atom's ionization due to the Trotter dynamics. Mathematically, we find that certain domain conditions are not satisfied by some states so higher moments of the potential and kinetic energies diverge. Our analytical error analysis agrees with numerical simulations, indicating that we can estimate the state-dependent Trotter error scaling genuinely.

§ INTRODUCTION

Solving the Schrödinger equation is the most fundamental endeavor in quantum physics.
However, depending on the model, this task is extremely difficult or even impossible to perform analytically.
Therefore, methods to find approximate solutions are of paramount importance.
Besides perturbation theory, the Trotter product formula serves as a powerful tool to approximate the time-evolution of a quantum system.
By subdividing the total time into small steps and by switching between simpler components, it allows the implementation of the dynamics under the full Hamiltonian.
The advantage of this approach is that the approximation becomes more accurate as the smaller time steps are chosen.
Thus, one can tune the accuracy of the Trotter approximation according to the problem at hand.
Due to this, the Trotter product formula has become the basic building block in many areas of physics.
To name three of them:
(i) Quantum computing — One of the most promising and powerful potential applications of quantum computers is simulating the dynamics of quantum systems.
Already in 1982, Feynman conjectured an exponential quantum advantage from this approach [1].
Lloyd then delivered an explicit algorithm in 1996 to implement Feynman's idea [2] on the basis of the Trotter product formula.
Ever since, it has been an active area of research and the list of papers referring to the exponential quantum advantage is rapidly growing [3].
A particularly important subclass of quantum simulation is that of quantum chemistry — the simulation of chemical materials [4, 5, 6, 7, 8, 9, 10].
It has been proposed as a potential pathway to developing new pharmaceuticals [11], catalysts [12], or fertilizers [13, 14].
All these simulations can be understood in terms of the Trotter product formula.
(ii) Mathematical physics — Here, the study of the convergence of the Trotter product formula has a long history.
See, for instance, Refs. [15, 16] for some early works.
Apart from elementary interest in the Trotter product formula due to its fundamental importance for quantum physics, it also has practical implications for mathematical physics.
For instance, it forms the basis for the path-integral formulation of quantum mechanics [17, 18, 19].
As such, the Trotter product formula can be the foundation for quantum field theory [20, 21, 22, 23].
(iii) Numerical simulations — The standard algorithms to solve the Schrödinger equation numerically are based on the Trotter product formula.
Most prominently, the split-step method (also known as strang-splitting) is used whenever the potential energy is reasonably well understood.
This is the case for instance in the simulation of atoms and molecules [24, 25, 26, 27, 28] or quantum optics [29].
Furthermore, numerical algorithms for studying quantum many-body systems are based on Trotterization [30, 31, 32, 33].

In all practical simulations of quantum systems, it is important to understand the Trotter error.
That is, how close is the Trotter evolution to the actual target evolution?
On one hand, in quantum simulation, one uses the Quantum Phase Estimation algorithm to extract the spectral information of a Hamiltonian of interest.
The controlled unitary gates in the quantum circuit are implemented by Trotterizing between noncommuting elementary gates corresponding to a representation of the target Hamiltonian.
The Trotter error then determines the necessary resources (gate count) of the algorithm.
Due to the lack of large fault-tolerant quantum computers existing to date, providing such resource estimates is a difficult problem.
For this reason, the current approach is to use heuristics, so that resource estimates can be made based on a set of assumptions on the Trotter error [34, 13, 14].
On the other hand, for numerical simulations, one directly Trotterizes between noncommuting parts of the target Hamiltonian or a truncated (discretized) version of it.
Here, the Trotter error determines the convergence of the algorithm to a reasonable solution.
Therefore, efficient implementations of the Trotter product formula and a realistic estimate for the Trotter error are central in the development of both quantum technology and numerical physics.
The current state-of-the-art approaches are bounds, which quantify the Trotter error in terms of the operator norm of the involved commutators [35, 36, 37].
However, for particular input states, these results are far from tight.
This is because the operator norm quantifies the worst-case scenario corresponding to the input state with the largest error.
For this reason, a state-dependent analysis leads to a more realistic and tighter estimate of the actual Trotter error.
This has also been noticed e.g. in Ref. [35], where the authors conclude: “The spectral-norm error bound […] would be overly pessimistic if we simulate with a low-energy initial state. It would then be beneficial to change to the error metric of the Euclidean distance to avoid the worst-case propagation.”
A simple example where this can be observed is given in Sec. <ref>:
For two Hamiltonians $A$ and $B$, we can have $[A, B]$ large, but $A\varphi\approx0$ and $B\varphi\approx0$ for the input state $\varphi$ of interest.
Even more convincing is the example of atomic and molecular Hamiltonians.
Here, the commutator norm bounds diverge and a state-dependent analysis is actually mandatory.
See for instance Refs. [21, 38, 39, 40] or the discussion in Sec. <ref>.

These considerations lead to the immediate desire to be able to quantify state-dependent Trotter convergence speeds in a general scenario.
In this paper, we solve this issue.
We derive state-dependent Trotter bounds that apply to arbitrary Hamiltonians with eigenstates.
These bounds might have wide-ranging applications in quantum technology, mathematical physics, and numerical simulations.
They refine existing bounds in the literature, where it is usually assumed that the Trotter error scales as $N^{-1}$, i.e. with the inverse number of Trotter steps $N$.
This is justified by comparing the Taylor expansions of Trotterized and target dynamics, which coincide up to the first order.
Due to the importance of efficient Trotterization, a lot of effort has been put into improving this method and achieving agreement of Trotterized and target dynamics in higher Taylor orders as well.
This led to what is called higher-order Trotter-Suzuki decompositions [35].
In practice, the second-order Trotterization (also known as strang-splitting or split-step method) is the most commonly used higher-order Trotter scheme.
This is because it improves the asymptotic scaling to $N^{-2}$, while still being similarly cost-efficient as the first-order Trotterization.
As for the standard first-order Trotter product formula, state-dependent bounds for higher-order schemes are not available in the literature yet.
We also derive general state-dependent Trotter bounds for higher-order Trotter schemes, which are valid for arbitrary Hamiltonians with eigenstates.
This combines both advantages of faster error scaling due to higher-order Trotterization and smaller error prefactor due to an explicit dependency on the input state.

Of course, upper bounds are not a mathematical proof of the asymptotic convergence speed.
Therefore, we bolster our general theory with a detailed case study of chemistry simulation.
This numerical analysis convincingly indicates the tightness of the asymptotic scaling of our bounds.
As an application of our results, we study the hydrogen atom, which is the paradigmatic example to look at.
Historically, it has been the first quantum system, for which the Schrödinger equation was solved analytically [41, 42].
This makes the hydrogen atom interesting from a fundamental perspective.
But also for practical purposes, the hydrogen atom is of outstanding importance.
It is simple enough to gain analytical insights while still reflecting many features of more complicated systems.
Therefore, its wave functions form the basis for our understanding of multi-electron atoms [43].
For example, many molecular Hamiltonians are modeled using linear combinations of atomic orbitals (LCAO), which are derived from the hydrogen eigenfunctions [44].
By studying certain eigenstates of the Hamiltonian of the hydrogen atom, we find that the asymptotic Trotter error scaling with $N$ differs from the predictions in the literature.
More precisely, the scaling with $N$ is determined by the quantum number $\ell$ of the orbital angular momentum.
In particular, we find that all eigenstates with $\ell\geq2$ (d-orbitals and higher) admit the desired $N^{-1}$ scaling.
However, for s-orbitals ($\ell=0$), the Trotter error only scales as $N^{-1/4}$.
This is remarkable as it implies that the Trotter error for the ground state is larger than expected.
Notice that the ground state plays a particularly important role in chemistry simulations as it encodes many interesting properties of an atom or a molecule.

We then ask the question whether higher-order Trotter methods provide relief and admit faster scalings.
This leads to two more revelations.
First, we find an asymmetry in the scaling of the Trotter bound.
This shows that the ordering of the Trotter product could be crucial.
Second, for s-orbitals including the ground state, we do not find any improvement from higher-order schemes.
In fact, the Trotter error scaling remains $N^{-1/4}$.
This shows that our observation on the slower ground-state convergence carries over to $p$th-order product formulas and puts higher-order schemes into question for certain ground-state chemistry simulations.

§.§ Relation to previous results

Historically, the study of product formulas started in the field of pure mathematics.
However, applications of product formulas are now ubiquitous in many fields, such as applied mathematics (see strang-splitting or symmetric integrators), physics (see Trotterization), and computer science (see split-step method).
In recent years, all these applied fields have contributed to advancing our understanding of product formulas.
Nevertheless, there has been a lack of transfer of the results among these research areas.
A prominent example where this can be observed is that of chemistry simulations, which have only been studied from the perspective of physics and computer science.
In this paper, we employ methods from different fields, shed light on their connections, and discuss the implications for chemistry simulations from different angles.
We hope that this will help to start a discourse among the different areas employing product formulas.

The main contribution of physics and computer science to product formulas has been the development of more and more efficient Trotter error bounds to improve digital simulations.
This led to the seminal and influential work of Childs et al. [35].
These results have been further improved very recently by Schubert and Mendl [36] and Zhuk et al. [37].
Furthermore, the physical origin of Trotter errors in hopping models is discussed in Refs. <cit.>.
The error bounds provided in these references quantify the Trotter error in some matrix norm, usually the operator norm, and give rise to a commutator scaling.
Such bounds are particularly well-suited for simulations of many-body systems, local observables, or quantum Monte-Carlo [35].
If we wanted to apply these norm bounds in the context of chemistry simulations, we would first truncate (discretize) the atomic or molecular Hamiltonian of interest at some finite level.
Afterward, we would either simulate the truncated Hamiltonian directly by Trotterization or go to a second-quantized picture.
In any case, the commutator of the noncommuting components of the resulting operator would determine the Trotter error.
The hope is that for a large enough truncation dimension, the actual target dynamics gets approximated well.
With the same logic, the so-computed Trotter error would be close to the true Trotter error of the full system.
However, this logic can sometimes be wrong and additional care has to be taken.
In fact, the norm Trotter errors simply diverge for unbounded operators such as the ones considered in chemistry simulations.
See Eq. (<ref>) for an example.
Intuitively, this happens because norm error bounds increase with extending truncation dimension, which ultimately leads to a diverging error.
This shows that there is a very intricate interplay between the error of truncation and the error of Trotterization [40].
The circumstance that norm error bounds are not the right quantities to look at for product formulas is well known in pure [15, 16, 21] and applied mathematics [45, 46, 47, 48, 49] as well as mathematical physics [39, 21].
For this reason, it is necessary to study state-dependent (also known as strong in the mathematical literature) error bounds for chemistry simulations.
So far, Trotter methods considering low-energy subspaces have been developed [50, 51] and some effort on state-dependent Trotter bounds has been made [45, 52, 53, 40, 54, 46, 47, 49, 48].
However, these results are not general enough to apply to atomic and molecular Hamiltonians in quantum chemistry.
In particular, the bounds in Refs. [50, 51] still diverge, whereas the bounds in Refs. [45, 52, 53, 40, 54, 46] still require at least one of the two operators to be bounded, which is not the case in chemistry simulations.
Furthermore, the assumptions on the multi-commutators in Refs. [47, 48, 49] are not satisfied for chemistry simulations.
On the contrary, our bounds are directly applicable to chemistry simulations.
They apply to general $p$th-order product formulas as well and provide conditions under which the Trotter error scales as $N^{-p}$.
For the case where these conditions are violated, we derive error bounds with slower fractional scalings.
Notice that fractional Trotter error bounds are known in the mathematical physics literature for operator-norm convergence of semigroup dynamics [55, 56].
However, these only apply to imaginary time-evolution and thus cannot be used to analyze the dynamics of atoms and molecules.
In addition, a particular state-dependent bound with a fractional scaling has been found in Ref. <cit.> under strong assumptions, which are not satisfied by chemistry Hamiltonians.

Due to the generality, our results can explain some open questions in the literature.
For example, the authors of Ref. [28] numerically found larger absolute ground-state Trotter errors that imply increased resource estimates for quantum chemistry.
So far, an explanation of this circumstance is lacking.
Furthermore, the problem of slower Trotter convergence for the ground state of the hydrogen atom has already been observed numerically in 2006 [57].
Nevertheless, the explanation for this effect remained open.
Both effects are aligned with our investigation and can be explained by our results.
Notice that similar issues are known for the Numerov method in numerical simulation [58], which might be due to related reasons.
In addition, a recent analysis by Lee et al. [59] suggests that a proof for an exponential quantum advantage is yet to be found.
In their conclusion, the authors state [59]:
“Numerical calculations are neither mathematical proof of asymptotics with respect to size and error, nor can we exclude exponential quantum advantage in specific problems.
However, our results suggest that without new and fundamental insights, there may be a lack of generic exponential quantum advantage in this task.”
Our paper aims to kick off such a discussion and provides ideas on how to tackle problems in quantum chemistry analytically.
Notice that our results only imply a polynomially increased runtime and do not affect a potential exponential quantum advantage of quantum simulations.

§.§ Structure of the paper

The remainder of the paper is structured as follows.
In Sec. <ref>, we present state-dependent error bounds for the standard (first-order) Trotter product formula (Thm. <ref>).
These bounds are valid even for unbounded operators such as atomic and molecular Hamiltonians under certain constraints on the input state.
Furthermore, we show how to obtain state-dependent error bounds when the input state does not satisfy the conditions (Thm. <ref>).
In this case, the scaling of the Trotter error gets slowed down.
We apply these results to the hydrogen atom in Sec. <ref> and find slower Trotter scalings for the ground state and other states with low angular momenta.
Afterwards, these findings are verified numerically.
Section <ref> is dedicated to the second-order Trotterization.
We develop state-dependent bounds here, which are also applicable to unbounded Hamiltonians (Thm. <ref>).
Again, these bounds rely on assumptions on the input state of interest, and error bounds with slower scalings are obtained when the assumptions are not satisfied (Thm. <ref>).
These general results are complemented by the application to the hydrogen atom in Sec. <ref>.
We find that there is no improvement from the second-order strategy for the ground state.
However, excited states can admit a provably advantageous scaling given the ordering of the operators in the Trotter product is correct.
In Sec. <ref>, we present a method to obtain state-dependent bounds for higher-order Trotter-Suzuki product formulas.
We explicitly compute the bounds for a fourth-order Trotterization (Thm. <ref>).
Furthermore, we provide a loose estimate for arbitrary-order Trotterization, that shows the desired asymptotic scaling if the input state satisfies certain conditions (Thm. <ref>).
These conditions are not satisfied for the ground state of the hydrogen atom, indicating the breakdown of the higher-order Trotter hierarchy.
Section <ref> then gives more details on the methods and proof ideas used to obtain the results.
Possible generalizations and an outlook on how our methodology might apply to other areas of research are given in Sec. <ref>.
Finally, we conclude in Sec. <ref>.
The technical details of the derivations are collected in the appendices.
Appendix <ref> provides the proof of the first-order Trotter bounds presented in Sec. <ref> (Thm. <ref> and Thm. <ref>).
How these results explicitly carry over to the hydrogen atom is discussed in Appendix <ref>.
The method to obtain the higher-order Trotter bounds is then outlined in detail in Appendix <ref>.
This enables us to prove the second-order Trotter bounds presented in Sec. <ref> (Thm. <ref> and Thm. <ref>), and the loose bound on the general $p$th-order Trotter error (Thm. <ref>).
Afterward, we provide a Mathematica script to explicitly calculate higher-order bounds with small prefactors in Appendix <ref>.
This script is used to compute the fourth-order Trotter bounds (Thm. <ref>) as well as the bounds for a particular sixth-order product formula.

§ FIRST-ORDER TROTTER BOUNDS

For two Hamiltonians $H_1$ and $H_2$, the commonly used bound in the literature to estimate the Trotter error is [35, 60]

\begin{align}
\left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}\right)^{N}
-\rme^{-\rmi t(H_{1}+H_{2})}
\right\|\nonumber \\
\le\frac{t^{2}}{2N}\Vert [H_{1},H_{2}]\Vert,
\label{eq:operator_norm_bound}
\end{align}

where $\Vert A\Vert=\sup_{\Vert\psi\Vert=1}\Vert A\psi\Vert$ is the operator norm of an operator $A$ with $\Vert\psi\Vert=\sqrt{\langle\psi\vert\psi\rangle}$ the (Euclidean) norm of a vector $\psi$.
If $A$ is a matrix, $\Vert A\Vert$ computes its largest singular value.
Thus, the standard error measure in Eq. (<ref>) considers a worst-case scenario of the input state with the largest Trotter error.
Clearly, by incorporating the explicit dependency on the input state, a tighter and more accurate error analysis is possible.
In particular, consider the case mentioned in the introduction, where $[H_1, H_2]$ is large but $H_1\varphi\approx 0$ and $H_2\varphi\approx 0$ for an input state $\varphi$ of interest.
A simple example of this scenario is given by the two block matrices

\begin{equation}
\medskip
0 & \epsilon\\
\epsilon & A
\end{pmatrix},
\qquad
\medskip
0 & \epsilon\\
\epsilon & B
\end{pmatrix},
\label{eq:example}
\end{equation}

$\Vert[A, B]\Vert\gg 1$ and $\Vert\epsilon\Vert\ll 1$.
Here, $\Vert[H_1,H_2]\Vert\gg 1$ due to the noncommuting blocks $A$ and $B$.
Thus, the norm error (<ref>) is large in this case.
However, consider a state

\begin{equation}
\varphi=\begin{pmatrix}
\medskip
\varphi_0 \\
\end{pmatrix},
\label{eq:state}
\end{equation}

which is only supported in the zero blocks of $H_1$ and $H_2$.
For such a state, one has $\Vert H_1\varphi\Vert=\Vert H_2\varphi\Vert=\Vert\epsilon\Vert\ll 1$, which would lead to a small state-dependent Trotter error.
This example shows that norm bounds can indeed be very loose in situations, where state-dependent bounds are actually small.
Intuitively, one would expect this to happen for low-energy input states [50, 51].
Therefore, it is natural to consider the Trotter error evaluated on eigenstates of the target Hamiltonian.

Furthermore, we will see in Sec. <ref> that state-dependent bounds are necessary to describe problems in quantum chemistry.
This is because in general we have $\Vert[H_1, H_2]\Vert=\infty$ since the involved operators are unbounded.
Our state-dependent bounds for the Trotter product formula are valid for both bounded and unbounded Hamiltonians.
For the case of unbounded Hamiltonians, we need to impose a condition on the input state.
This is because an unbounded operator $A$ is only defined on a certain (dense) subspace $\mathcal{D}(A)\subset\mathcal{H}$ of the Hilbert space $\mathcal{H}$.
The space $\mathcal{D}(A)$ is called the domain of $A$ and incorporates many important physical properties of the operator, such as boundary conditions.

As a simple example, consider the standard one-dimensional position operator
$X:\psi(x)\mapsto x \psi(x)$ on wave functions $\psi\in \mathcal{H}=L^2(\mathbb{R})$.
The operator $X$ has to map a wave function into a wave function.
That is, we have to require that  $X \psi$ is again square-integrable, i.e. $x\psi(x)\in L^2(\mathbb{R})$.
This imposes a condition on the possible input states of $X$, which restricts its domain $\mathcal{D}(X)=\{\psi \in L^2(\mathbb{R}) : x\psi(x)\in L^2(\mathbb{R})\}$.
Obviously, $\mathcal{D}(X)$ is a strict subset of the entire $L^2(\mathbb{R})$.
For example, the square-integrable function $\psi(x)= 1/(1+|x|)$ does not belong to $\mathcal{D}(X)$, because $x\psi(x)$ is not a wave function ($\|X\psi\|=\infty$).

Notice that all Hamiltonians on finite-dimensional Hilbert spaces are automatically bounded.
Bounded operators are defined on the entire Hilbert space so that we do not have to worry about their domains.
Therefore, the domain conditions in our results are trivially satisfied and can be dropped for all bounded Trotter pairs, in particular for finite-dimensional systems.
Conversely, for unbounded self-adjoint operators, the domain is always a strict subspace of the Hilbert space and the domain conditions are mandatory.

We now present our state-dependent bound for the first-order Trotter product formula.
For simplicity, we initially consider eigenstates and lift these bounds to general states later.
To this end, we define the eigenstate-dependent error

\begin{equation}
\xi_N(t;\varphi)\equiv
\left\Vert \left[\left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}\right)^{N}-\rme^{-\rmi th}\right]\varphi\right\Vert,\label{eq:definition_standard_error}
\end{equation}

where $h$ is an eigenstate of $H_1+H_2$, i.e. $(H_1+H_2)\varphi=h\varphi$.
Theorem <ref> gives a bound on $\xi_N(t;\varphi)$ and is a generalization of Ref. <cit.> to unbounded Hamiltonians.

Let $H_{1}$ be self-adjoint on $\mathcal{D}(H_{1})$ and $H_{2}$
be self-adjoint on $\mathcal{D}(H_{2})$.
Let $\varphi$ be an eigenstate of $H_{1}+H_{2}$ with eigenvalue $h$, i.e. $(H_{1}+H_{2})\varphi=h\varphi$.
If $\varphi\in\mathcal{D}(H_{1}^{2})\cap\mathcal{D}(H_{2}^{2})$,

\begin{equation*}
\xi_N(t;\varphi)\le\frac{t^{2}}{2N}\,\Bigl(\Vert(H_{1}-g)^{2}\varphi\Vert+\Vert(H_{2}-h+g)^{2}\varphi\Vert\Bigr),
\end{equation*}

for all $t,g\in\mathbb{R}$, and the Trotter product formula converges on $\varphi$.

This is proved in Appendix <ref>.
The idea is to shift the energy of $H_1+H_2$ in $\varphi$ to zero.
By this approach, the target evolution applied to $\varphi$ becomes the identity.
The Trotter error is bounded by $N$ times the error in a single step, which can be bounded under the above domain conditions on $\varphi$ by using a second-order Taylor expansion for the short-time evolution

\begin{equation}
\rme^{-\rmi \frac{t}{N} H_j}\varphi= \varphi -\rmi \frac{t}{N} H_j \varphi + \mathcal{O}\!\left(\frac{t^2}{N^2}\|H_j^2\varphi\|\right),
\label{eq:Taylor2t}
\end{equation}

for $j=1,2$.

One might wonder why Thm. <ref> does not give rise to a commutator scaling.
However, one can show that general state-dependent Trotter bounds with  a commutator scaling cannot exist:

There exist self-adjoint operators $H_{1}$ and $H_{2}$, and vectors $\varphi\in\mathcal{D}(H_1H_2)\cap\mathcal{D}(H_2H_1)$, such that a  Trotter error bound of the form

\begin{equation}
\xi_N(t;\varphi)\leq \omega\Vert [H_1,H_2]\varphi\Vert,
\label{eq:no-go_commutator_scaling}
\end{equation}

with $\omega>0$, does not hold.

Essentially, the reason why a general bound (<ref>) cannot exist is because $[H_1, H_2]\varphi$ does not carry enough information.
Therefore, we can easily construct counter-examples, where such a bound fails to apply.
To this end, we have to find an example with $[H_1, H_2]\varphi=0$ but $[\rme^{-\rmi t H_1},\rme^{-\rmi t H_2}]\varphi\neq0$.
Then, the right hand side in Eq. (<ref>) is zero, even though the Trotter formula is not exact.
We now give two examples of this situation.
Example 1 — Consider the two $3\times3$ Hamiltonians

\begin{equation*}
H_1 = \begin{pmatrix}
0 & 1 & 0 \\
1 & 0 & 1 \\
0 & 1 & 1
\end{pmatrix},
\quad
H_2 = \begin{pmatrix}
0 & -1 & 0 \\
-1 & 0 & -1 \\
0 & -1 & 0
\end{pmatrix},
\end{equation*}

and the input state

\begin{equation*}
\varphi = \begin{pmatrix}
1 \\
0 \\
\end{pmatrix}.
\end{equation*}

Notice that $\varphi$ is an eigenstate of $H_1+H_2$ with eigenvalue zero.
Furthermore, we have $\Vert [H_1,H_2]\varphi\Vert=0$, but $\xi_N(t;\varphi)>0$ if $N<\infty$ and $t\neq0$.
Example 2 — While Example 1 only goes wrong for certain input states, it is even possible to construct examples where $\Vert [H_1,H_2]\varphi\Vert=0$ but $\xi_N(t;\varphi)>0$ on a dense set of states.
For this, consider the model from Ref. [61]: $H_1=-P^2$ and $H_2=P^2+P$, where $P$ is the one-dimensional momentum operator defined on the real half-line $L^2(\mathbb{R}_+)$.
The operators $H_1$ and $H_2$ are self-adjoint on a common dense domain $D$ <cit.>, and  $[H_1,H_2]\varphi=0$ for all $\varphi\in \mathcal{D}(0,+\infty)\subset\mathcal{D}(H_1H_2)\cap\mathcal{D}(H_2H_1)$, which is dense in $L^2(\mathbb{R}_+)$.
Thus, the right hand side of Eq. (<ref>) is zero.
However, as shown in Ref. [61], the sum $H_1+H_2=P$ is not essentially self-adjoint on $D$ and Trotter does not converge strongly on $D$.
Therefore, the actual Trotter error is larger than zero even in the limit $N\rightarrow\infty$.

Notice that there might not always be eigenstates, i.e. if $H_1+H_2$ only admits a purely continuous spectrum.
However, this does not happen in chemistry problems, where the states of bounded electrons (“orbitals”) constitute a set of eigenstates for the atomic or molecular target Hamiltonian.
Note also that Thm. <ref> naturally extends to Trotter products with more than two operators.
See the discussion around Eq. (<ref>).

To treat generic input states, we can apply Thm. <ref> to superpositions of eigenstates.

For states that are superpositions of eigenstates of $H_1+H_2$,

\begin{equation*}
\psi=\sum_{\ell=1}^L c_{\ell}\varphi_{\ell},
\end{equation*}

with $L$ finite or $L=\infty$,
we can use the triangle and Cauchy-Schwarz inequalities to obtain for all $g\in\mathbb{R}$

\begin{align}
\xi_N(t;\psi)
\le\left(\sum_{\ell=1}^L  \xi_N(t;\varphi_\ell)^2\right)^{1/2}
\nonumber\\
\le\frac{t^{2}}{2N}\left[\sum_{\ell=1}^L\Bigl(
\Vert H_{1}(g)^{2}\varphi_{\ell}\Vert
+\Vert H_{2}(g)^{2}\varphi_{\ell}\Vert\Bigr)^2\right]^{1/2},
\label{eq:superposition_states}
\end{align}

where $H_1(g)=H_1-g$, $H_2(g)=H_2-h+g$, and the right-hand side is defined as $+\infty$ if it diverges [62].

Theorem <ref> requires the input state to be in the domain intersection $\varphi\in\mathcal{D}(H_{1}^{2})\cap\mathcal{D}(H_{2}^{2})$, so that
$\|H_j^2\varphi\|<\infty$ ($j=1,2$) and the bound does not diverge.
However, this might not always be the case.
In particular, we will see in Sec. <ref> that this domain condition is violated by the ground state of the hydrogen atom.
For such a case, we establish the following theorem, which works under relaxed domain conditions.
It comes at the price of a slower Trotter convergence speed.
The scaling is ruled by the high-energy distribution of the input state $\varphi$,

\begin{equation*}
\operatorname{\mathrm{Prob}}_{\varphi}(|H_j|\geq \Lambda)
\end{equation*}

where $\theta$ is the Heaviside step function, and
with $\Omega\subset\mathbb{R}$ are
the spectral measures of the Hamiltonians $H_j$ ($j=1,2$) at the input state $\varphi$.

Let $H_1$ be self-adjoint on $\mathcal{D}(H_1)$ and $H_2$ be self-adjoint on $\mathcal{D}(H_2)$.
Let $\varphi$ be an eigenstate of $H_1+H_2$ with eigenvalue $h$, i.e. $(H_1+H_2)\varphi=h\varphi$.
Furthermore, let the tails of the spectral measures $\mu_{j,\varphi}$ of $H_j$ at $\varphi$ decay as

\begin{equation}
\mu_{j,\varphi}(\{|\lambda|\geq\Lambda\})
=\mathcal{O}\!\left(\frac{1}{\Lambda^{2\delta}}\right),\quad j=1,2,
\label{eq:tail_bound_1}
\end{equation}

for $\Lambda>0$ and some $\delta>1$.

\begin{equation*}
\xi_N(t;\varphi)
\medskip
\displaystyle
\mathcal{O}\!\left(\frac{t^\delta}{N^{\delta-1}}\right),&1<\delta<2,\\
\medskip
\displaystyle
\mathcal{O}\!\left(\frac{t^2}{N}\sqrt{\log (N/t)}\right),&\delta=2,\\
\displaystyle
\mathcal{O}\!\left(\frac{t^2}{N}\right),&\delta>2.
\end{cases}
\end{equation*}

The proof for $\delta\in(1,2]$ is provided in Appendix <ref>.
We use the same strategy of Thm. <ref>.
However, the bound on the single-step error will now require a finer analysis than the second-order Taylor expansion (<ref>), since $\|H_j^2\varphi\|$ is no longer assured to be bounded.
By using the spectral theorem, one can show that the high-energy decay of the spectral measures implies instead the following error on the short-time evolution:

\begin{equation*}
\rme^{-\rmi \frac{t}{N} H_j}\varphi
=\varphi -\rmi \frac{t}{N} H_j \varphi + \mathcal{O}\!\left(\frac{t^\delta}{N^\delta}\right),
\end{equation*}

for $\delta<2$, and an analogous error with an additional logarithmic factor for $\delta=2$, whence the theorem follows.
Finally, the case $\delta>2$ is covered by Thm. <ref>, since in this case the decay (<ref>) implies that $\varphi\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2^2)$.

Note that the conditions on the tails of the spectral distributions in Eq. (<ref>) imply that $\varphi\in\mathcal{D}(|H_1|^{\gamma_1})\cap\mathcal{D}(|H_2|^{\gamma_2})$ for all $\gamma_j\in[0,\delta)$.
If $\delta\le2$, it is not guaranteed that $\varphi\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2^2)$, and one obtains a convergence slower (by a logarithmic factor) than the convergence $\mathcal{O}(1/N)$ given by Thm. <ref>.

Notice that the assumption that $\varphi$ is an eigenstate of $H_1+H_2$ and thus $\varphi\in\mathcal{D}(H_1)\cap\mathcal{D}(H_2)$ implies the decay (<ref>) of the spectral distributions with $\delta=1$.
Therefore, the condition (<ref>) contains no additional information and is immaterial for $\delta\leq 1$.
As a matter of fact, one can prove the convergence of the Trotter product formula without assuming (<ref>) (or equivalently for $\delta\leq 1$).
However, in such a case the rate of convergence $\xi_N(t;\varphi)\to 0$ as $N\to\infty$ can be arbitrarily slow.

If Eq. (<ref>) converges and $\psi$ is fully supported in the part of the Hilbert space corresponding to the point spectrum of $H_1+H_2$, then the Trotter error scaling is determined by the domain conditions in Thm. <ref> and Thm. <ref>.
Notice that this is always the case for any bounded $H_1+H_2$, for which the condition (<ref>) holds for any arbitrary large $\delta$.
Moreover, this is true for chemistry simulations, where the input state of quantum phase estimation or the Hartree state aims to approximate the target eigenstate well.
In particular, one is usually interested in the space of bounded electron states of a chemical Hamiltonian, which corresponds to its point spectrum.

Together, Thm. <ref> and Thm. <ref> allow us to study the Trotter error for arbitrary eigenstates of the target Hamiltonian.
The bounds can be lifted to generic states employing Cor. <ref>.
Our bound from Thm. <ref> is tight in the sense that it is saturated by the example provided in Eq. (<ref>):
If $\epsilon=0$, then the state $\varphi$ in Eq. (<ref>) has eigenvalue $h=0$, i.e. $(H_1+H_2)\varphi=0$.
Furthermore, $H_1^2\varphi=H_2^2\varphi=0$.
Therefore, choosing $g=0$ leads to the bound $\xi_N(t;\varphi)=0$, which is exact.
Thus, our bounds allow for more efficient quantum simulation due to an explicit and tight state dependency.
In the following, we apply them to the chemistry simulation of the hydrogen atom.

§ APPLICATION OF THE FIRST-ORDER BOUNDS TO THE HYDROGEN ATOM

The Hamiltonian of the hydrogen atom is given by

\begin{equation*}
H_\mathrm{hydrogen} = H_1 + H_2,
\end{equation*}

with (temporarily re-establishing the Planck constant $\hbar$)

\begin{equation}
\label{eq:hamiltonian_kin}
\end{equation}

\begin{equation}
\label{eq:hamiltonian_pot}
\end{equation}

$m_\rme$ is the electron mass (or equivalently the reduced electron mass), $e$ is the elementary (electron) charge, and

\begin{equation*}
\end{equation*}

is the Bohr radius, with $\varepsilon_0$ being the vacuum permittivity.
Furthermore, $\Delta=\nabla^2$ is the Laplace operator, and $r=\vert \mathbf{r} \vert$ denotes the distance from the nucleus, which is placed at $\mathbf{r}=(0, 0, 0)$.
The eigenfunctions $\Psi_{n\ell m}(\mathbf{r})$
of the hydrogen atom have eigenenergies

\begin{equation*}
\end{equation*}

and are labeled by the principal quantum number $n=1,2,\ldots$, as well as $\ell=0,1,\ldots,n-1$ and $m=-\ell,-\ell+1,\ldots,\ell-1,\ell$, which determine the orbital angular momentum.

To implement the dynamics under the Hamiltonian $H_\mathrm{hydrogen}=H_1+H_2$ of the hydrogen atom, one Trotterizes between the kinetic energy $H_1$ and the potential energy $H_2$.
This approach is called split-step method and is commonly used in the numerical analysis of quantum chemistry, see e.g. Refs. [28, 24, 26].
For the advantages of this approach over other methods, see Refs. [24, 26, 27, 25].
If we want to know the speed of Trotter convergence, we can compute the distance between the Trotterized dynamics and the actual dynamics.
For finite-dimensional systems, this can be bounded by Eq. (<ref>).
However, it is easy to see that for the hydrogen atom, we have

\begin{equation}
\Vert [H_{1},H_{2}]\Vert
\propto\left\Vert \left[\frac{\Delta}{2},\frac{1}{r}\right]\right\Vert
=\left\Vert \frac{1}{r^2}\frac{\rmd}{\rmd r}\right\Vert
\label{eq:commutator_diverge}
\end{equation}

Physically, this comes from the fact that the commutator is a differential operator (involving radial momentum) and has a singularity at $r=0$:
There are wave functions in the Hilbert space with arbitrary large radial momentum and/or large support (high electron density) close to this singularity at the nucleus.
From a mathematical perspective, this is not a surprise, since both $H_1$ and $H_2$ are unbounded and we need state-dependent bounds.
See, for instance, Refs. [21, 39, 40].
Since the task of quantum chemistry simulation is to find the spectrum of a Hamiltonian, it is perfectly reasonable to look at the convergence speed
on its eigenstates.
That is, the figure of merit for the hydrogen atom is

\begin{equation*}
\xi_N(t;\Psi_{n\ell m}) = \left\Vert
\left[
\left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}\right)^{N} -\rme^{-\rmi t E_n}
\right]\Psi_{n\ell m}
\right\Vert.
\end{equation*}

It is well known that the eigenfunctions of the hydrogen atom are given by the product of the radial wave function $R_{n\ell}$ and the spherical harmonics $Y_{\ell m}$.
The spherical part $Y_{\ell m}$ only accounts for the degeneracies in the spectrum of the hydrogen atom, so different $m$ values for fixed $n$ and $\ell$ lead to the same Trotter dynamics.
Therefore, we can restrict our attention to the radial part $R_{n\ell}$ when computing $\xi_N$.
Now, we apply our bound from Thm. <ref> to this Trotter scenario.
For example, we could consider the excited state $\Psi_{320}$ of the hydrogen atom, for which we can bound

\begin{equation*}
\xi_N(t;\Psi_{320})\leq \frac{1}{40 t_0^2}\frac{t^2}{N}, \qquad t_0=\frac{m_\mathrm{e}a_0^2}{\hbar}.
\end{equation*}

As opposed to the standard bound from Eq. (<ref>), which just gives the trivial bound for the hydrogen atom due to $\|[H_1, H_2]\|=\infty$, this bound is nontrivial and has the expected scaling of $\mathcal{O}(t^2/N)$.

The ground state is of particular importance for quantum chemistry applications.
However, our bound from Thm. <ref> diverges in this case as the domain conditions are not satisfied by the ground state $\Psi_{100}$ of the hydrogen atom.
Therefore, we need to perform a more refined analysis through Thm. <ref>.
This indicates that the Trotter product formula converges slower than $\mathcal{O}(N^{-1})$ for the ground state $\Psi_{100}$ of the hydrogen atom.
In fact, we find a slower scaling of the Trotter error for input states with low angular momenta.
That is, all s-orbitals only admit the $N^{-1/4}$ scaling and all p-orbitals lead to the $N^{-3/4}$ scaling.
Importantly, we emphasize that the quantum number $\ell$ of the orbital angular momentum alone determines the state-dependent Trotter scaling.
We summarize our findings in Table <ref>.

Scalings of the analytical first-order Trotter error bounds for the eigenfunctions $\Psi_{n\ell m}$ of the hydrogen atom.
See Sec. <ref> and specifically Eqs. (<ref>)–(<ref>).
We find different scalings depending on the input states, which are solely determined by the orbital angular momentum $\ell$.
In particular, the scalings for the s- and p-orbitals are slower than the expected $N^{-1}$ scaling.

c]@c@  Orbital angular momen-
tum quantum number $\bm{\ell}$       c]@c@ Scaling of bound on
Trotter error $\bm{\xi_N(t;\Psi_{n\ell m})}$

$\ell = 0$ (s-orbitals)     $\mathcal{O}(N^{-1/4})$

$\ell = 1$ (p-orbitals)     $\mathcal{O}(N^{-3/4})$

$\ell \geq 2$ (d-orbitals and higher)    $\mathcal{O}(N^{-1})$

We also provide explicit error bounds in Sec. <ref> of the Methods section.
See Eqs. (<ref>)–(<ref>).
These bounds allow us to compute the necessary resources to achieve a chemical accuracy for an algorithm for a chemistry simulation.
For instance, assume one would like to implement the dynamics under the hydrogen atom Hamiltonian for the ground state $\Psi_{100}$ for a total evolution time $t$ of the order of microseconds.
To reach a chemical accuracy [63],
one would need Trotter steps $N$ of the order $10^{3}$–$10^{4}$.

Since the results in Table <ref> are just upper bounds, it makes sense to convince ourselves numerically that we indeed obtain a slower scaling for the ground state than for higher excited states.
To compute the Trotter error, we numerically integrate the Schrödinger equation using XMDS2 [64], an open-source package for solving multidimensional partial differential equations.
For these simulations, we use a grid-based method and work in the Hartree atomic units, so that $\hbar=m_\rme=a_0=1$ and the Hamiltonians become $H_1=-\Delta/2$ and $H_2=-1/r$.
We diagonalize the kinetic energy $H_2=-\Delta/2$ by representing it in its spectral basis of spherical Bessel functions.
The number of Bessel functions used corresponds to the number of grid points in position space as well as the number of modes in the spectral space, resulting in a discretized system suitable for numerical integration.
More details on the numerical simulation methods can be found in Sec. <ref> of the Methods section.

We simulated the system multiple times, each using a different number of modes while keeping the radial cutoff constant.
The results are shown in Fig. <ref>.
The ground-state Trotter error initially scales as $N^{-1/4}$, but asymptotes to $N^{-1}$.
This transition occurs later and later (at a larger number of the Trotter steps $N$) as the number of Bessel modes (numerical cutoff dimension) increases.
Notice that the eigenstate Trotter error always asymptotically scales as $N^{-1}$ for finite-dimensional systems [40].
Therefore, for a finite number of Bessel modes, we eventually see the $N^{-1}$ scaling for large enough $N$.
For this reason, a later transition point to the $N^{-1}$ scaling with increased system dimension is the only way to numerically observe slower convergence rates for the full unbounded Hamiltonians.

[mark size=1.5, scale=1]
xlabel=Trotter steps $N$,
ylabel=First-order Trotter error $\xi_N(t;\Psi_{100})$,
label style=font=,
tick label style=font=,
x post scale=1,
y post scale=1,
legend style=at=(0.4,-0.2),anchor=north,
legend columns=3,
legend cell align=left,
xmin=2, xmax=3000,
[color=1, thick, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_groundstate_100.csv;
[color=2, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_groundstate_200.csv;
[color=3, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_groundstate_300.csv;
[color=4, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_groundstate_400.csv;
[color=5, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_groundstate_800.csv;
[color=0, ultra thick, domain=2:3000, densely dashed]0.3*x^(-1);
[color=6, ultra thick, domain=2:3000]2.04665*x^(-1/2) + 0.832107/x + 1.34365/x^(1/4);
[color=0, ultra thick, domain=2:3000]0.31/x + 0.1/x^(1/2) + 0.1/x^(1/4);

$100$ modes,
$200$ modes,
$300$ modes,
$400$ modes,
$800$ modes,
analytic bound,
$\infty$ mode scaling,

[color=6,rotate=-8] at (4,4.75) upper bound;
[color=0,rotate=-8] at (4,3.5) scaling for $\infty$ modes;
[-stealth, ultra thick] (3,2) node[below, align=left] increasing
# of modes – (4.2,3);
[color=0,rotate=-25] at (0.7,3.15) $\mathcal{O}(N^{-1})$;

The Trotter error at time $t=1$ in the Hartree atomic units $\hbar=m_\mathrm{e}=a_0=1$, as a function of the Trotter steps for the ground state $\Psi_{100}$ of the hydrogen atom.
The radial cutoff in the simulations is $R=30$.
We show five different levels of discretization characterized by the number of Bessel modes.
For reference, we show the slopes of $N^{-1}$ (grey dashed line) and our bound (brown solid line), which scales as $N^{-1/4}$.
See Eqs. (<ref>)–(<ref>).
We see that the asymptotics for any finite discretization initially start in a consistent way with our bound, and eventually go as $N^{-1}$, but $N$ at which this transition happens becomes larger with increasing number of modes.
This provides an evidence that the scaling in the infinite-mode limit is indeed slower than $N^{-1}$.
To indicate a potential curve for the infinite-mode limit, we show a slope (gray solid line) with the same scaling as our analytic bound.

[mark size=1.5, scale=1]
xlabel=Trotter steps $N$,
ylabel=First-order Trotter error $\xi_N(t;\Psi_{210})$,
label style=font=,
tick label style=font=,
x post scale=1,
y post scale=1,
legend style=at=(0.5,-0.2),anchor=north,
legend columns=3,
legend cell align=left,
xmin=2, xmax=3000,
[color=1, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_n2_l1_10.csv;
[color=2, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_n2_l1_20.csv;
[color=3, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_n2_l1_50.csv;
[color=4, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_n2_l1_200.csv;
[color=5, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_1st_order_n2_l1_800.csv;
[color=0, ultra thick, domain=2:3000, densely dashed]0.003*x^(-1);
[color=6, ultra thick, domain=2:3000]0.145959/x + 0.133831/x^(3/4);

$10$ modes,
$20$ modes,
$50$ modes,
$200$ modes,
$800$ modes,
analytic bound,

The Trotter error at time $t=1$ in the Hartree atomic units $\hbar=m_\mathrm{e}=a_0=1$, as a function of the Trotter steps for the state $\Psi_{210}$ of the hydrogen atom.
The radial cutoff in the simulations is $R=30$.
We show five different levels of discretization characterized by the number of Bessel modes; however, beyond 200 modes (red) the results are indistinguishable.
The grey dashed line shows the slope of the $N^{-1}$ scaling, and we see that the Trotter error scales as $N^{-1}$ in all cases, even though the state $\Psi_{210}$ does not satisfy the domain conditions for Thm. <ref>.
Our analytic bound, which scales as $N^{-3/4}$, is depicted in brown.

[mark size=1.5, scale=1]
xlabel=Trotter steps $N$,
ylabel=First-order Trotter error $\xi_N(t;\Psi_{320})$,
label style=font=,
tick label style=font=,
x post scale=1,
y post scale=1,
legend style=at=(0.5,-0.2),anchor=north,
legend columns=3,
legend cell align=left,
xmin=2, xmax=300,
[color=1, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_n3_l2_10.csv;
[color=2, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_n3_l2_20.csv;
[color=3, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_n3_l2_50.csv;
[color=4, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_n3_l2_800.csv;
[color=0, ultra thick, domain=2:300, densely dashed]0.0006*x^(-1);
[color=6, ultra thick, domain=2:300]0.023591/x;

$10$ modes,
$20$ modes,
$50$ modes,
$800$ modes,
analytic bound,

The Trotter error at time $t=1$ in the Hartree atomic units $\hbar=m_\mathrm{e}=a_0=1$, as a function of the Trotter steps for the state $\Psi_{320}$ of the hydrogen atom.
The radial cutoff in the simulations is $R=40$.
We use four different levels of discretization characterized by the number of Bessel modes; however, beyond 50 modes (green) the results are indistinguishable.
The grey dashed line shows the slope of the $N^{-1}$ scaling, and we see that the Trotter error scales as $N^{-1}$ in all cases as expected.
For reference, we plot our analytic bound in brown.
It also shows an $N^{-1}$ behavior.

We also investigate the Trotter error scaling for the excited state $\Psi_{210}$ of the hydrogen atom with $n=2$, $\ell=1$, and $m=0$ in Fig. <ref>, as well as for the excited state $\Psi_{320}$ with $n=3$, $\ell=2$, and $m=0$ in Fig. <ref>.
In both cases, we find a $N^{-1}$ scaling independent of the number of Bessel modes.
This behavior shows that the Trotter errors for the ground and excited states scale differently in the full system limit of infinitely many Bessel modes.
More specifically, excited states admit $N^{-1}$ scaling, whereas the ground state only scales as $N^{-1/4}$.
These results are rather remarkable as the ground state is usually the target of quantum chemistry simulations.
Its Trotter error scales slower than expected, indicating that more resources (quantum gates) are needed for accurate quantum chemistry simulation.
Ultimately, this leads to a polynomially increased runtime of ground-state quantum chemistry algorithms.
Furthermore, the increased gate count can intensify the errors of the quantum computation due to environmental decoherence or gate imperfections.

Physically, the slower scaling of the Trotter error in the ground state, and in all s-orbitals, can be explained by the fact that the electron is close to the $1/r$ singularity.
Due to this, it has very high fluctuations in both potential and kinetic energies.
In particular, we will see in Appendix <ref> that the fourth moments, i.e. $\|H_j^2\varphi\|$, are diverging.
When performing the Trotterization, in each period the state evolves under the bare potential and kinetic Hamiltonians.
Due to its huge energy fluctuations, there is a significant probability for the electron to be kicked out of the space of the bound states by this evolution.
In turn, the hydrogen atom gets ionized.
As we increase the number of Trotter steps $N$, each Trotter cycle becomes shorter and we approximate the target dynamics under the Hamiltonian of the hydrogen atom better.
Therefore, the ionization effect decreases in $N$.
If we truncate the system to make it amenable to numerical simulations, we will always cut off certain modes corresponding to free electron states.
As a result, these modes are not attained from the Trotter evolution and a smaller ionization effect is observed numerically.
Going to higher and higher truncation dimensions allows for more room in the space of free electron states.
In turn, we can numerically see a stronger ionization effect, which leads to a worse approximation of the target dynamics.
This explains the scaling behavior we see numerically in Fig. <ref>.
Ultimately, in the limit of infinitely many modes, there is always a significant probability of leaking out of the bound space no matter how large the Trotter number $N$ is.
Therefore, the full Trotter problem converges slower than $N^{-1}$.
We numerically observe this ionization behavior and show the results in Fig. <ref>.
It is noticeable that the s-orbitals admit particularly strong ionization behaviors.
This is consistent with our analytical bounds as well as the numerical results in Figs. <ref>–<ref>.

[mark size=1.5, scale=1]
xlabel=Trotter steps $N$,
label style=font=,
tick label style=font=,
x post scale=1,
y post scale=1,
legend style=at=(0.5,-0.2),anchor=north,
legend columns=3,
legend cell align=left,
xmin=2, xmax=30,
[color=1, mark=*, ultra thin] table[x=N, y=error, col sep=comma] Ionization_n1_l0.csv;
[color=2, mark=*, ultra thin] table[x=N, y=error, col sep=comma] Ionization_n2_l0.csv;
[color=3, mark=*, ultra thin] table[x=N, y=error, col sep=comma] Ionization_n3_l0.csv;
[color=4, mark=*, ultra thin] table[x=N, y=error, col sep=comma] Ionization_n4_l0.csv;
[color=5, mark=*, ultra thin] table[x=N, y=error, col sep=comma] Ionization_n3_l2.csv;
[color=6, mark=*, ultra thin] table[x=N, y=error, col sep=comma] Ionization_n4_l2.csv;

$n=1$, $\ell=0$,
$n=2$, $\ell=0$,
$n=3$, $\ell=0$,
$n=4$, $\ell=0$,
$n=3$, $\ell=2$,
$n=4$, $\ell=2$,

Numerical simulations of the ionization probability as a function of the number of Trotter steps.
The initial state $\varphi$ is prepared in various energy eigenstates $\Psi_{n\ell m}$ of the hydrogen atom, and then we let it evolve to $\varphi(t)=U_N(t)\varphi$ by the Trotter evolution $U_N(t)$, with the Hamiltonians $H_1$ and $H_2$ given by Eqs. (<ref>) and (<ref>), for a total time $t=1$ in the Hartree atomic units $\hbar=m_\mathrm{e}=a_0=1$.
The ionization probability is given by $1-\|P_{\mathrm{bd}}\varphi(t)\|^2$, where $P_{\mathrm{bd}}=\sum_{n\ell m} |\Psi_{n\ell m}\rangle \langle\Psi_{n\ell m}|$ is the projection on the space of bound states of the hydrogen atom.
For the evolution under the non-Trotterized Hamiltonian $H = H_1 + H_2$, the ionization fraction was zero up to numerical errors $(<10^{-10})$ as expected.
For the Trotterized evolution, $\varphi(t)$ acquires a nonzero component out of the space spanned by the bound states, resulting in a nonzero ionization probability, as shown in the figure.
The $\ell=0$ eigenstates ionize much more heavily than the $\ell=2$ eigenstates, but in both cases, the ionization rate decreases with the number of Trotter steps as the Trotter approximation approaches the true evolution.

Mathematically, the slower convergence of the ground state can be explained by the fact that both $H_1$ and $H_2$ are unbounded operators.
It is well known in mathematical physics that, for unbounded operators, the Trotter product formula can potentially converge slower (see e.g. Ref. [39, 45]) or even not at all (see e.g. Refs. [65, 66]).
In our context, slower convergence arises from the circumstance that the ground state of the hydrogen atom is not in the domain of $H_1^2$ and $H_2^2$, see Appendix <ref>.
On the contrary, the state $\Psi_{320}$ with $n=3$, $\ell=2$, and $m=0$ considered in Fig. <ref> satisfies this domain condition.
This allowed us to derive analytical bounds that scale as $N^{-1}$, which is consistent with the numerical simulation.

§ SECOND-ORDER TROTTER BOUNDS

A more refined method to simulate quantum systems that is commonly used in both classical and quantum simulations is the second-order Trotterization.
It is obtained by symmetrizing the Trotter product as

\begin{equation}
\mathcal{S}_N^{(2)}(t)=\left(\rme^{-\rmi \frac{t}{2N}H_1}\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{2N}H_1}\right)^N,
\label{eq:second_order}
\end{equation}

or alternatively,

\begin{equation}
\tilde{\mathcal{S}}_N^{(2)}(t)=\left(\rme^{-\rmi \frac{t}{2N}H_2}\rme^{-\rmi \frac{t}{N}H_1}\rme^{-\rmi \frac{t}{2N}H_2}\right)^N.\label{eq:second_order_reversed}
\end{equation}

Through this approach, one obtains a tighter error scaling of $\mathcal{O}(t^3/N^2)$ in the case of bounded operators.
In particular <cit.>,

\begin{align}
\Vert\mathcal{S}_N^{(2)}(t) - \rme^{-\rmi t(H_1+H_2)}\Vert \le{}& \frac{t^3}{12N^2}\Vert [H_2, [H_2, H_1]]\Vert\nonumber\\
&{}+\frac{t^3}{24N^2}\Vert [H_1, [H_1, H_2]]\Vert,
\label{eq:norm_bound_2nd_order}
\end{align}

and analogously for $\tilde{\mathcal{S}}_N^{(2)}(t)$.
In the following, we will focus on the $H_1H_2H_1$ scheme given by the product formula $\mathcal{S}_N^{(2)}(t)$.
Nevertheless, we can always derive equivalent results for the $H_2H_1H_2$ scheme $\tilde{\mathcal{S}}_N^{(2)}(t)$ by exchanging $H_1\leftrightarrow H_2$.
In fact, both product formulas only differ by boundary terms.
As in the case of the first-order Trotterization, the bound in Eq. (<ref>) diverges for unbounded operators such as the kinetic and potential energies of the hydrogen atom.
See also Ref. [39].
Therefore, we develop explicitly state-dependent bounds for the second-order Trotterization, which enable us to study such problems.
For this, we define

\begin{equation}
\xi_N^{(2)}(t;\varphi)\equiv \Vert [\mathcal{S}_N^{(2)}(t)-\rme^{-\rmi ht}]\varphi\Vert,
\label{eq:second_order_error}
\end{equation}

where again $(H_1+H_2)\varphi=h\varphi$.
We then find the following state-dependent bound on the second-order Trotter product.

Let $H_1$ and $H_2$ be self-adjoint on their respective domains $\mathcal{D}(H_1)$ and $\mathcal{D}(H_2)$, and let $(H_1+H_2)\varphi=h\varphi$.
For $g\in\mathbb{R}$, define ${H}_1(g)=H_1-g$ and ${H}_2(g)=H_2-h+g$.
Then, the state-dependent Trotter error of $\mathcal{S}_N^{(2)}(t)$ can be bounded for all $g\in\mathbb{R}$ and $t\geq 0$ by

\begin{align*}
\xi_N^{(2)}(t;\varphi)\leq \frac{t^3}{N^2}\,\biggl(&
\frac{1}{24}\Vert {H}_1(g)^3\varphi\Vert + \frac{1}{8}\Vert {H}_2(g){H}_1(g)^2\varphi\Vert\nonumber\\
&\qquad\qquad\qquad\quad\ \
{}+\frac{1}{12}\Vert {H}_2(g)^3\varphi\Vert\biggr),
\end{align*}

provided that $\varphi\in\mathcal{D}(H_1^3)\cap\mathcal{D}(H_2H_1^2)\cap\mathcal{D}(H_2^3)$.

This is proved in Appendix <ref>.
For the proof, we again shift the energy of $H_1+H_2$ with respect to $\varphi$ to zero to turn the target evolution on $\varphi$ to the identity $I$.
Then, we use integration by parts to rewrite the error operator $\mathcal{S}_N^{(2)}(t)-I$.
By noticing that the boundary terms disappear, we are only left with a remainder, that admits the desired scaling of $\mathcal{O}(t^3/N^2)$.
Explicitly bounding this remainder gives the theorem.

As in the case of the first-order Trotterization, it is possible to state a no-go theorem for the existence of the second-order Trotter bounds with commutator scaling by making the same argument as in Thm. <ref>.
In fact, this reasoning can be extended to arbitrary higher-order product formulas.
Furthermore, the bound from Thm. <ref> can be extended to generic input states by considering superpositions of eigenstates.
See Cor. <ref>.

Notice that the domain conditions on the input state $\varphi$ in Thm. <ref> are even stronger than the ones for the first-order Trotterization in Thm. <ref>.
For this reason, the ground state $\Psi_{100}$ of the hydrogen atom does not satisfy them and we expect a slower scaling behavior also for the second-order Trotterization.
As in the case of the first-order Trotterization (Thm. <ref>), we derive domain conditions and bounds for a slower fractional scaling.

Let $H_1$ be self-adjoint on $\mathcal{D}(H_1)$ and $H_2$ be self-adjoint on $\mathcal{D}(H_2)$.
Let $\varphi$ be an eigenstate of $H_1+H_2$ with eigenvalue $h$, i.e. $(H_1+H_2)\varphi=h\varphi$.
Let the spectral measures $\mu_{j,\varphi}$ of $H_j$ at $\varphi$ decay as

\begin{equation}
\mu_{j,\varphi}(\{|\lambda|\geq \Lambda\}) = \mathcal{O}\!\left(\frac{1}{\Lambda^{2\delta}}\right),\quad  j=1,2,
\label{eq:tail_bound_2}
\end{equation}

for $\Lambda>0$ and some $\delta\in(1,2]$.

\begin{equation*}
\xi_N(t;\varphi)
\medskip
\displaystyle
\mathcal{O}\!\left(\frac{t^\delta}{N^{\delta-1}}\right),
\displaystyle
\mathcal{O}\!\left(\frac{t^2}{N}\sqrt{\log (N/t)}\right), &\delta=2.
\end{cases}
\end{equation*}

Furthermore, if the decay (<ref>) holds for $\delta\geq 2$, $\varphi\in \mathcal{D}(H_1^2)$, and
the spectral measure $\mu_{2,H_1^2\varphi}$ of  $H_2$ at the vector $H_1^2\varphi$ decays as

\begin{equation*}
\mu_{2,H_1^2\varphi}(\{|\lambda|\geq \Lambda\})=\mathcal{O}\!\left(\frac{1}{\Lambda^{2(\delta-2)}}\right),
\end{equation*}

for  $\Lambda>0$, then,

\begin{equation*}
\xi_N(t;\varphi)=\begin{cases}
\medskip
\displaystyle
\mathcal{O}\!\left(\frac{t^\delta}{N^{\delta-1}}\right),
\medskip
\displaystyle
\mathcal{O}\!\left(\frac{t^3}{N^2}\sqrt{\log (N/t)}\right), &\delta=3,\\
\displaystyle
\mathcal{O}\!\left(\frac{t^3}{N^2}\right), &\delta>3.
\end{cases}
\end{equation*}

The proof is similar to the proof of Thm. <ref>.
The case $\delta>3$ is covered by Thm. <ref>, since in this case we are sure that $\varphi\in\mathcal{D}(H_1^3)\cap\mathcal{D}(H_2H_1^2)\cap\mathcal{D}(H_2^3)$ and the domain condition for Thm. <ref> is fulfilled.
The remaining nontrivial parts of the theorem are proved by the method described in  Appendix <ref>.

We now apply our second-order Trotter bounds to the problem of simulating the hydrogen atom.

§ APPLICATION OF THE SECOND-ORDER BOUNDS TO THE HYDROGEN ATOM

Applied to the eigenfunctions $\Psi_{n\ell m}$ of the hydrogen atom, Thm. <ref> results in scalings slower than $\mathcal{O}(N^{-2})$ for all $\ell\leq 3$.
As in the case of the first-order Trotterization, the scaling is again entirely determined by the quantum number $\ell$ of the orbital angular momentum.
We summarize our findings in Table <ref>.

Scalings of the analytical second-order Trotter error bounds for the  eigenfunctions $\Psi_{n\ell m}$ of the hydrogen atom.
The fractional scalings are derived from Thm. <ref>.
We find different scalings depending on the input states, which are solely determined by the quantum number $\ell$ of the orbital angular momentum.
In particular, the scalings for all s- to f-orbitals with $\ell\le3$ are slower than the expected $N^{-2}$ scaling.
Interestingly, the scalings for the f-orbitals slightly differ depending on whether we use the $H_1H_2H_1$ product formula $\mathcal{S}_N^{(2)}(t)$ [see the $\xi_N^{(2)}(t;\Psi_{n\ell m})$ column] or the $H_2H_1H_2$ product formula $\tilde{\mathcal{S}}_N^{(2)}(t)$ [see the $\tilde{\xi}_N^{(2)}(t;\Psi_{n\ell m})$ column].

c]@c@  Orbital angular momen-
tum quantum number $\bm{\ell}$
c]@c@ Scaling of
$\bm{\xi_N^{(2)}(t;\Psi_{n\ell m})}$
c]@c@ Scaling of
$\bm{\tilde{\xi}_N^{(2)}(t;\Psi_{n\ell m})}$

$\ell = 0$ (s-orbitals)     $\mathcal{O}(N^{-1/4})$     $\mathcal{O}(N^{-1/4})$

$\ell = 1$ (p-orbitals)     $\mathcal{O}(N^{-3/4})$     $\mathcal{O}(N^{-3/4})$

$\ell = 2$ (d-orbitals)    $\mathcal{O}(N^{-5/4})$     $\mathcal{O}(N^{-5/4})$

$\ell = 3$ (f-orbitals)    $\mathcal{O}(N^{-3/2})$     $\mathcal{O}(N^{-7/4})$

$\ell \geq 4$ (g-orbitals and higher)    $\mathcal{O}(N^{-2})$     $\mathcal{O}(N^{-2})$

We numerically investigate the ground-state Trotter error in Fig. <ref>.
Interestingly, we observe the same $N^{-1/4}$ scaling behavior as in the case of the first-order Trotterization.
This shows that the second-order Trotterization is not advantageous in this case.
See also Fig. <ref> below.

[mark size=1.5, scale=1]
xlabel=Trotter steps $N$,
ylabel=Second-order Trotter error $\xi_N^{(2)}(t;\Psi_{100})$,
label style=font=,
tick label style=font=,
x post scale=1,
y post scale=1,
legend style=at=(0.5,-0.2),anchor=north,
legend columns=3,
legend cell align=left,
xmin=2, xmax=3000,
[color=1, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_groundstate_100.csv;
[color=2, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_groundstate_200.csv;
[color=3, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_groundstate_300.csv;
[color=4, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_groundstate_400.csv;
[color=5, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_groundstate_800.csv;
[color=0, ultra thick, densely dashed, domain=2:3000]0.15*x^(-2);
[color=6, ultra thick, domain=2:3000]0.5*(2.04665*x^(-1/2) + 0.832107/x + 1.34365/x^(1/4));

$100$ modes,
$200$ modes,
$300$ modes,
$400$ modes,
$800$ modes,
analytic bound,

The error of the second-order Trotter product formula at time $t=1$ in the Hartree atomic units $\hbar=m_\mathrm{e}=a_0=1$, as a function of the Trotter steps for the ground state $\Psi_{100}$ of the hydrogen atom. See Eq. (<ref>). The radial cutoff in the simulations is $R=30$. We show five different levels of discretization characterized by the number of Bessel modes. For reference, we show the slopes of $N^{-2}$ (grey dashed line) and our analytic bound (brown solid line), which scales as $N^{-1/4}$. Notice that the second-order bound just computes to half the first-order bound. We see that the asymptotics for any finite discretization initially start in alignment with our bound, and eventually go as $N^{-2}$, but $N$ at which this transition happens becomes larger with increasing number of modes. This provides an evidence that the scaling in the infinite-mode limit is indeed slower than $N^{-2}$. In particular, we expect no improvement in the asymptotic scaling from the second-order Trotterization over the first-order Trotterization, see Fig. <ref>.

An interesting observation about the bounds in Thm. <ref> is that they are asymmetric under the exchange of $H_1$ and $H_2$.
This means that also the domain conditions on the state $\varphi$ change if we implement this change in the second-order Trotter product formula (<ref>).
Indeed, we find an asymmetry in the scaling for the eigenstate $\Psi_{430}$ of the hydrogen atom, see Table <ref>.
However, we do not see this effect numerically and both product formulas seem to admit the same scaling.
We simulated the Trotter error for $\Psi_{430}$ for the $H_1H_2H_1$ scheme $\mathcal{S}_N^{(2)}(t)$ and performed the same simulation for the $H_2H_1H_2$ scheme $\tilde{\mathcal{S}}_N^{(2)}(t)$.
The results are shown in Fig. <ref>.
Both curves are indistinguishable; nevertheless, $\Psi_{430}$ exhibits a scaling which is slower than $N^{-2}$.
However, we conjecture that there exist quantum systems, which indeed admit an asymmetric scaling for certain input states.

[mark size=1.5, scale=1]
xlabel=Trotter steps $N$,
ylabel=Second-order Trotter error $\mathord{\mathop{\xi}\limits^{\scriptscriptstyle(\sim)}}_N^{(2)}(t;\Psi_{430})$,
label style=font=,
tick label style=font=,
x post scale=1,
y post scale=1,
legend style=at=(0.5,-0.2),anchor=north,
legend columns=2,
legend cell align=left,
xmin=2, xmax=3000,
[color=1, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_ABA_n4_l3_200.csv;
[color=1, mark=x, only marks, mark size=4] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_BAB_n4_l3_200.csv;
[color=2, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_ABA_n4_l3_400.csv;
[color=2, mark=x, only marks, mark size=4] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_BAB_n4_l3_400.csv;
[color=3, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_ABA_n4_l3_800.csv;
[color=3, mark=x, only marks, mark size=4] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_BAB_n4_l3_800.csv;
[color=4, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_ABA_n4_l3_1200.csv;
[color=4, mark=x, only marks, mark size=4] table[x=N, y=error, col sep=comma] Trotter_Error_2nd_order_BAB_n4_l3_1200.csv;
[color=0, mark=*, ultra thick, mark options=line width=0.4pt] table[x=N, y=error, col sep=comma] Scaling_n4_l3_Three_Half.csv;
[color=0, mark=x, ultra thick, mark options=line width=0.4pt, mark size=4] table[x=N, y=error, col sep=comma] Scaling_n4_l3_Seven_Quarter.csv;
[color=0, ultra thick, domain=2:3000, densely dashed]0.01*x^(-2);

$200$ modes $\mathcal{S}_N^{(2)}$,
$200$ modes $\tilde{\mathcal{S}}_N^{(2)}$,
$400$ modes $\mathcal{S}_N^{(2)}$,
$400$ modes $\tilde{\mathcal{S}}_N^{(2)}$,
$800$ modes $\mathcal{S}_N^{(2)}$,
$800$ modes $\tilde{\mathcal{S}}_N^{(2)}$,
$1200$ modes $\mathcal{S}_N^{(2)}$,
$1200$ modes $\tilde{\mathcal{S}}_N^{(2)}$,

The error of the second-order Trotter product formula at time $t=1$ in the Hartree atomic units $\hbar=m_\mathrm{e}=a_0=1$, as a function of the Trotter steps for the eigenstate $\Psi_{430}$ of the hydrogen atom.
The radial cutoff in the simulations is $R=20$.
We show four different levels of discretization characterized by the number of Bessel modes.
For reference, we show the slope of $N^{-2}$ (grey dashed line).
The results of the $H_1H_2H_1$ scheme $\mathcal{S}_N^{(2)}(t)$ of the second-order Trotterization in Eq. (<ref>) are shown by dots $\bullet$.
The results of the $H_2H_1H_2$ scheme $\tilde{\mathcal{S}}_N^{(2)}(t)$ of the second-order Trotterization in Eq. (<ref>) are shown by crosses $\times$.
Both schemes show indistinguishable results.
However, we see again that the asymptotics for any finite discretization initially start slower, and eventually go as $N^{-2}$, but $N$ at which this transition happens becomes larger with increasing number of modes.
This provides an evidence that the scaling in the infinite-mode limit is indeed slower than $N^{-2}$.
For reference, we show the predicted scaling of $N^{-3/2}$ for $\mathcal{S}_N^{(2)}(t)$ (grey line with dots) and the predicted scaling of $N^{-7/4}$ for $\tilde{\mathcal{S}}_N^{(2)}(t)$ (grey line with crosses).

From the previous discussion, it seems like higher-order Trotter schemes are more restricted in terms of their domain conditions than the lower-order schemes.
To investigate this further, we develop a method to obtain state-dependent higher-order Trotter bounds for arbitrary order $p$.
This is done by generalizing the approach for the second-order results.
In the next section, we explain how to derive such bounds and explicitly study the example of a fourth-order Trotter product formula.

§ HIGHER-ORDER TROTTERIZATION

A common practice in quantum simulation is to employ higher-order Trotter schemes [35].
These are obtained by symmetrizing the Trotter product, which results in the cancellation of low-order terms in the Taylor expansion of the exponentials [67].
For two bounded Hamiltonians $H_1$ and $H_2$, the error of the $p$th-order Trotter product formula $\mathcal{S}^{(p)}(t)$ scales as $\mathcal{O}(t^{p+1}/N^p)$ [35].
Schemes with a higher order $p$ might have many different solutions for the short switching times between $H_1$ and $H_2$ [67].
An example of fourth-order scheme ($p=4$) is [68]

\begin{align}
\mathcal{S}_N^{(4)}(t)=\Bigl(&\rme^{-\rmi\frac{t}{N}\frac{\tau}{2}H_1}\rme^{-\rmi \frac{t}{N}\tau H_2}\rme^{-\rmi \frac{t}{N}\frac{1-\tau}{2}H_1}\rme^{-\rmi \frac{t}{N}(1-2\tau)H_2}\nonumber\\
{}\times \rme^{-\rmi\frac{t}{N}\frac{1-\tau}{2}H_1}\rme^{-\rmi\frac{t}{N}\tau H_2}\rme^{-\rmi\frac{t}{N}\frac{\tau}{2}H_1}\Bigr)^N,
\label{eq:fourth_order}
\end{align}

where $\tau=1/(2-2^{1/3})$.
If the two Hamiltonians $H_1$ and $H_2$ are bounded, the Trotter error scales as

\begin{equation*}
\Vert\mathcal{S}_N^{(4)}(t) - \rme^{-\rmi t(H_1+H_2)}\Vert = \mathcal{O}\!\left(\frac{t^5}{N^4}\right).
\end{equation*}

Notice that the better scalings of the higher-order schemes in the number of Trotter steps $N$ simultaneously go along with having to implement more unitaries per Trotter cycle.
A recent study [67] suggests that an 8th-order expansion could lead to optimal results requiring minimal resources.

Nevertheless, the results mentioned above are usually not applicable to unbounded operators, see e.g. Refs. [39, 40].
As before, we need to perform a state-dependent analysis in this case.
In Sec. <ref> of the Methods section, and in Appendix <ref> in more detail, we present a method to obtain such state-dependent bounds systematically for general $p$th-order product formulas through repeated integration by parts.
This method is inspired by Refs. [69, 40].
By generalizing these ideas, we are able to carry them over to unbounded operators provided the input state satisfies certain domain conditions.

In general, a $p$th-order Trotter product is given in the form

\begin{equation}
\mathcal{S}_N^{(p)}(t)=\left(\rme^{-\rmi\frac{t}{N}\tau_MH_{\sigma(M)}}\cdots \rme^{-\rmi\frac{t}{N}\tau_2H_{\sigma(2)}}\rme^{-\rmi\frac{t}{N}\tau_1H_{\sigma(1)}}\right)^N,
\label{eq:higher_order_evolution}
\end{equation}

\begin{equation}
\sigma(j)=\begin{cases}
\medskip
\displaystyle
1,&\text{if } j \text{ is odd},\\
\displaystyle
2,&\text{if } j \text{ is even},
\end{cases}
\label{eq:parity}
\end{equation}

and $\sum_{i=1}^{\lceil M/2\rceil}\tau_{2i-1}=\sum_{i=1}^{\lfloor M/2\rfloor}\tau_{2i}=1$, where $\lceil x\rceil$ is the ceiling and $\lfloor x\rfloor$ is the floor function.
This product is to be compared with $\rme^{-\rmi t(H_1+H_2)}$.
We then define the $p$th-order state-dependent Trotter error by

\begin{equation}
\xi_N^{(p)}(t;\varphi)=\Vert[\mathcal{S}_N^{(p)}(t)-\rme^{-\rmi th}]\varphi\Vert,
\label{eq:higher_order_error}
\end{equation}

where again $(H_1+H_2)\varphi=h\varphi$.

Note that some $\tau_j$ in the product formula (<ref>) can be negative.
Indeed, in the fourth-order Trotter product in Eq. (<ref>), $\tau_3=\tau_5=1-\tau$ and $\tau_4=1-2\tau$ are negative.
In fact, it has been shown by Suzuki that product formulas with $p>2$ necessarily have to involve some negative $\tau_j$ <cit.>.
Many examples of $p$th-order product formulas can be found in the literature, see e.g. Refs. [70, 71, 68, 67].
For instance, in the case of Suzuki's first fractal method for generating higher-order product formulas [70, 68, 67], the number of exponentials involved in the product formula of an even order $p$ is given by $M=2\cdot3^{p/2-1}+1$ [67].
In any case, our method for bounding $\xi_N^{(p)}(t;\varphi)$ works for any $p$th-order product formula.
Although our method allows us to set up and solve the equations for the switching times between $H_1$ and $H_2$ in a $p$th-order Trotter product, we assume that these are known here.

A Mathematica script that explicitly computes the Trotter bounds according to our procedure for given switching times is outlined in Appendix <ref>.
This method produces bounds with small prefactors that can be used for practical state-dependent Hamiltonian-simulation purposes.
Nevertheless, we here provide a loose state-dependent bound, which only aims to show the scaling.

Let $H_1$ be self-adjoint on $\mathcal{D}(H_1)$ and $H_2$ be self-adjoint on $\mathcal{D}(H_2)$.
Let $\varphi$ be an eigenstate of $H_1+H_2$ with eigenvalue $h$, i.e. $(H_1+H_2)\varphi=h\varphi$.
For $g\in\mathbb{R}$, define ${H}_1(g)=H_1-g$ and ${H}_2(g)=H_2-h+g$.
Then, the state-dependent error of a $p$th-order Trotter product $\mathcal{S}_N^{(p)}(t)$, consisting of $M$ exponentials of $H_1$ and $H_2$ per Trotter cycle,
can be bounded for all $g\in\mathbb{R}$ and $t\geq 0$ by

\begin{equation*}
\xi_N^{(p)}(t;\varphi)\leq\frac{1}{(p+1)!}\frac{(\tau_{*}t)^{p+1}}{N^p}K_{p+1}(\varphi),
\end{equation*}

where $\tau_*=\sum_{j=1}^M|\tau_j|$ and

{K}_p(\varphi)=\max_{1\leq j_1\leq\dots\leq j_p\leq M}\|H_{\sigma(j_p)}(g)\cdots H_{\sigma(j_1)}(p)\varphi\|,

provided that $\varphi$ is in the domains of all the involved operator products.

This is proved in Appendix <ref>.

The bound in Thm. <ref> admits the desired scaling of

\begin{equation*}
\xi_N^{(p)}(t;\varphi)=\mathcal{O}\!\left(\frac{t^{p+1}}{N^p}\right),
\end{equation*}

and is directly applicable to any $p$th-order Trotter scheme.
However, as pointed out before, it can be very loose.
To obtain tighter bounds, we follow the method outlined in Sec. <ref>, which is a generalization of the strategy we developed to obtain the second-order bounds in Sec. <ref>.
For the first-order Trotterization, this method leads to the same result as Thm. <ref>.
Applied to the fourth-order Trotter product presented in Eq. (<ref>), it results in the following theorem.

Let $H_1$ and $H_2$ be self-adjoint on their respective domains $\mathcal{D}(H_1)$ and $\mathcal{D}(H_2)$, and let $\varphi$ be an eigenstate of $H_1+H_2$ with eigenvalue $h$, i.e. $(H_1+H_2)\varphi=h\varphi$.
For $g\in\mathbb{R}$, define ${H}_1(g)=H_1-g$ and ${H}_2(g)=H_2-h+g$.
Then, the fourth-order state-dependent Trotter product formula $\mathcal{S}_N^{(4)}(t)$ in Eq. (<ref>) is bounded by

\begin{align}
\xi_N^{(4)}(t;\varphi)\leq\frac{t^5}{N^4}\,\biggl(&
\frac{a_0}{{17\,280}}\Vert{H}_1(g)^5 \varphi\Vert
\nonumber\\
&{}+ \frac{a_1}{1\,152}\Vert {H}_1(g)^2{H}_2(g){H}_1(g)^2\varphi\Vert
\nonumber\\
&{}+ \frac{a_1}{1\,728}\Vert {H}_1^2(g){H}_2(g)^3\varphi\Vert
\nonumber\\
&{}+ \frac{a_2}{1\,728} \Vert {H}_1(g){H}_2(g){H}_1(g)^3\varphi\Vert
\nonumber\\
&{}+ \frac{a_2}{576} \Vert {H}_1(g){H}_2(g)^2{H}_1(g)^2\varphi\Vert
\nonumber\\
&{}+ \frac{a_2}{864} \Vert {H}_1(g){H}_2(g)^4\varphi\Vert
\nonumber\\
&{}+ \frac{a_3}{6\,912} \Vert {H}_2(g){H}_1(g)^4\varphi\Vert
\nonumber\\
&{}+ \frac{a_4}{576} \Vert {H}_2(g){H}_1(g){H}_2(g){H}_1(g)^2\varphi\Vert
\nonumber\\
&{}+ \frac{a_4}{864} \Vert {H}_2(g){H}_1(g){H}_2(g)^3\varphi\Vert
\nonumber\\
&{}+ \frac{a_5}{48} \Vert {H}_2(g)^2{H}_1(g)^3\varphi \Vert
\nonumber\\
&{}+ \frac{a_6}{1\,728} \Vert {H}_2(g)^3{H}_1(g)^2\varphi \Vert
\nonumber\\
&{}+ \frac{a_7}{4\,320} \Vert {H}_2(g)^5\varphi\Vert
\biggr),
\label{eq:bound_4th_order}
\end{align}

\begin{align*}
a_0&=23+19\cdot2^{1/3}+17\cdot 2^{2/3},\\
a_1&=4+3\cdot2^{1/3}+2\cdot 2^{2/3},\\
a_2&=14+11\cdot2^{1/3}+9\cdot 2^{2/3},\\
a_3&=68+55\cdot2^{1/3}+44\cdot 2^{2/3},\\
a_4&=18+14\cdot2^{1/3}+11\cdot 2^{2/3},\\
a_6&= 226+180\cdot2^{1/3}+143\cdot 2^{2/3},\\
a_7&= 330+263\cdot2^{1/3}+209\cdot 2^{2/3},
\end{align*}

for all $t,g\in\mathbb{R}$, if $\varphi$ is in the domains of all the involved operator products in Eq. (<ref>).

This is computed by using the Mathematica script provided in Appendix <ref>.
It follows exactly the method outlined in Sec. <ref> of the Methods section.
For the numerical values of the coefficients, see Appendix <ref>.

We also computed error bounds for a $6$th-order product formula numerically with our Mathematica script.
See Appendix <ref>.

[mark size=1.5, scale=1]
xlabel=Trotter steps $N$,
ylabel=Fourth-order Trotter error $\xi_N^{(4)}(t;\Psi_{100})$,
label style=font=,
tick label style=font=,
x post scale=1,
y post scale=1,
legend pos=south west,
legend columns=1,
legend cell align=left,
xmin=2, xmax=2500,
[color=1, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_4th_order_groundstate_100.csv;
[color=2, thick, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_4th_order_groundstate_200.csv;
[color=3, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_4th_order_groundstate_300.csv;
[color=4, mark=*, only marks] table[x=N, y=error, col sep=comma] Trotter_Error_4th_order_groundstate_400.csv;
[color=5, mark=*, only marks] table[x=N, y=error, col sep=comma]Trotter_Error_4th_order_groundstate_800.csv;
[color=0, ultra thick, domain=2:2500]2*x^(-1/4);
[color=0, ultra thick, domain=2:2500, densely dashed]2*x^(-4);

$100$ modes,
$200$ modes,
$300$ modes,
$400$ modes,
$800$ modes,

The error of the fourth-order Trotter product formula at time $t=1$ in the Hartree atomic units $\hbar=m_\mathrm{e}=a_0=1$, as a function of the Trotter steps for the ground state $\Psi_{100}$ of the hydrogen atom.
See Eq. (<ref>).
The radial cutoff in the simulations is $R=30$.
We show five different levels of discretization characterized by the number of Bessel modes.
For reference, we show the slopes for $N^{-4}$ (grey dashed line) and $N^{-1/4}$ (grey solid line).
We see that the Trotter error for any finite discretization initially starts as $N^{-1/4}$, and eventually goes as $N^{-4}$, but $N$ at which this transition happens becomes larger with increasing number of modes.
This provides an evidence that the scaling in the infinite-mode limit is indeed slower than $N^{-4}$.
In particular, we expect no improvement in the asymptotic scaling from the fourth-order Trotterization over the first- or second-order Trotterization, see Figs. <ref> and <ref>.

As  expected, the domain conditions in Thm. <ref> on the input state $\varphi$ for the fourth-order Trotter product are even stronger than the conditions in Thm. <ref> for the second-order Trotterization.
Again, the ground state of the hydrogen atom does not satisfy these conditions and a slower Trotter scaling is expected.
We numerically investigate the fourth-order Trotterization in Fig. <ref>.
As before, we observe an error scaling of $\mathcal{O}(N^{-1/4})$ in the full system limit, which coincides with the first- and second-order Trotterizations.
In summary, we find

\begin{equation*}
\xi^{(p)}_N(t;\Psi_{100})
=\mathcal{O}(N^{-1/4}), \quad p=1,2,4.
\end{equation*}

This is a strong indication that the entire higher-order Trotter hierarchy breaks down for the ground state $\Psi_{100}$ of the hydrogen atom.
Indeed, we show in Appendix <ref> how to obtain error bounds that admit fractional scalings for any $p$th-order product formulas.
This is a generalization of Thm. <ref> and Thm. <ref>.
We summarize the implications of our results as follows.
Assume that we are given a $p$th-order product formula $\mathcal{S}_N^{(p)}(t)$ and an input state $\varphi$, which is an eigenstate of $H_1+H_2$.
On one hand, if $\varphi$ satisfies the domain conditions of the product formula $\mathcal{S}_N^{(q)}(t)$ with $q\leq p$, then the Trotter error scales at least as $\xi_N^{(p)}(t;\varphi)=\mathcal{O}(N^{-q})$.
On the other hand, if $\varphi$ does not satisfy the domain conditions of the product formula $\mathcal{S}_N^{(r)}(t)$ with $r\leq p$, the Trotter error $\xi_N^{(p)}(t;\varphi)$ scales slower than $\mathcal{O}(N^{-r})$.
In total,

\begin{equation*}
\xi_N^{(p)}(t;\varphi) = \mathcal{O}(N^{-\delta}),\quad q\leq\delta<r.
\end{equation*}

Appendix <ref> provides the proof for this statement including explicit error bounds for $q=1,2$, but the described procedure can be easily iterated to obtain error bounds for $q\ge3$.
Furthermore, the reasoning from Cor. <ref> to treat superpositions of eigenstates can be applied to $p$th-order product formulas as well.

Since the ground state $\Psi_{100}$ of the hydrogen atom does not even satisfy the domain conditions for the first-order Trotterization ($q=0$), any higher-order method only admits error bounds with $N^{-1/4}$ scaling.
From our numerical investigations, it seems like these error bounds reflect the true scaling.
See Fig. <ref>.
That is, the Trotter error for the ground state $\Psi_{100}$ of the hydrogen atom always scales the same and slower than $\mathcal{O}(N^{-1})$ irrespectively of the Trotter order $p$.
This makes the first-order Trotterization the most favorable one for quantum chemistry simulations, as it involves the smallest number of unitaries per Trotter cycle.
A comparison of the first-, second-, and fourth-order Trotter errors as functions of the total number of unitaries to implement can be found in Fig. <ref>.
Since the evolutions of the first- and second-order Trotterizations only differ by two boundary unitaries, they admit similar performances.
The fourth-order product formula has a worse error-to-resource ratio though.
See the white region of Fig. <ref>.

[mark size=1.5, scale=1]
xlabel=Total number of unitaries,
ylabel=Ground-state Trotter error $\xi_N^{(p)}(t;\Psi_{100})$,
label style=font=,
tick label style=font=,
x post scale=1,
y post scale=1,
legend pos=south west,
legend columns=1,
legend cell align=left,
xmin=4, xmax=18001,
enlarge y limits=false,
[color=1, mark=*, only marks] table[x=N, y=error, col sep=comma] Order_Comparison_1st_order_groundstate_800.csv;
[color=2, thick, mark=*, only marks] table[x=N, y=error, col sep=comma] Order_Comparison_2nd_order_groundstate_800.csv;
[color=3, mark=*, only marks] table[x=N, y=error, col sep=comma] Order_Comparison_4th_order_groundstate_800.csv;
[ultra thick, samples=50, smooth, color=0, dashed] coordinates (881,0.00000001) (881,1);
[fill=lightgray,draw=none] coordinates (881,0.00000001) (18001,0.00000001) (18001,1) (881,1) ;

First-order ($p=1$),
Second-order ($p=2$),
Fourth-order ($p=4$),

Comparison of the Trotter errors for different orders $p$ of Trotterization for the ground state $\Psi_{100}$ of the hydrogen atom.
For the simulation, we fixed the total evolution time to $t=1$ and worked in the Hartree atomic units $\hbar=m_\mathrm{e}=a_0=1$.
Furthermore, we truncated the potential and kinetic energies at $800$ Bessel modes and chose a radial cutoff of $R=30$.
We then compare the Trotter errors (<ref>) for the first-order ($p=1$, blue dots), second-order ($p=2$, orange dots), and fourth-order ($p=4$, green dots) Trotterizations as functions of the total number of unitaries in the Trotter product.
If boundary unitaries of two adjacent Trotter cycles are generated by the same Hamiltonian, they are combined and only counted as one unitary.
In the regime where all Trotter errors scale as $N^{-1/4}$ (white region), the first- and second-order Trotterizations perform almost the same and are both better than the fourth-order product formula.
Once the total number of unitaries exceeds a threshold (dashed line), the Trotter errors start to scale faster.
In this region (light gray), product formulas with higher $p$ will eventually outperform product formulas with lower $p$.
However, this behavior does not reflect the true error scaling, due to finite-size effects.

§ METHODS

In this section, we describe the methods we used to obtain our results.
This includes both details on the numerical simulations and the derivations of the analytical error bounds.
The technical details and lengthy calculations can be found in the appendices.

§.§ Numerical methods for the grid-based hydrogen simulations

To perform the grid-based simulation described in Sec. <ref>, we numerically solve the Schrödinger equation

\begin{equation}
\rmi  \frac{\rmd \psi({\mathbf{r}},t)} {\rmd t} = H(t) \psi({\mathbf{r}},t),  \quad \psi({\mathbf{r}}, 0) = \Psi_{n\ell m}(\mathbf{r})
\label{eqSEhydrogenNaturalUnits}
\end{equation}

over some time $t$, with $H(t)$ being alternately taken as $H_1$ and $H_2$ for the Trotter step evolution, or $H_1 + H_2$ for the actual evolution.
The numerical integration of Eq. (<ref>) was carried out using XMDS2 [64], an open-source package for solving multidimensional partial differential equations.

The time-dependent wave functions are given by $\psi_{\text{Trot}}({\mathbf{r}}, t) = U_N(t) \psi({\mathbf{r}}, 0) $ and $\psi({\mathbf{r}}, t) = V(t) \psi({\mathbf{r}}, 0) $, with the initial state chosen to be an eigenfunction $\psi({\mathbf{r}}, 0) = \Psi_{n\ell m}(\mathbf{r})$ of the hydrogen atom.
Here, $U_N(t)$ is the dynamics given by the Trotter product and $V(t)$ is the target dynamics under $H_1+H_2$.
In terms of these wave functions, the Trotter error (<ref>) at time $t$ becomes

\begin{equation}
\xi_N(t;\Psi_{n\ell m}) = \left( \int_{\mathbb{R}^3}  |\psi_{\text{Trot}}({\mathbf{r}}, t) - \psi({\mathbf{r}},t)|^2 \, \rmd^3{\mathbf{r}} \right)^{1/2}.
\label{eq:schroedinger_integral}
\end{equation}

As is well known, Eq. (<ref>) is separable for the hydrogen atom, and since all the terms in the Hamiltonian are spherically symmetric, the evolution conserves the angular momentum, meaning only the differential equation for the radial coordinate is relevant.
In spherical coordinates, the Laplacian (corresponding to $H_1$) acting on an eigenstate $\Psi_{n\ell m}$ is given by

\begin{equation*}
\Delta = \frac{\rmd^2}{\rmd r^2} + \frac{2}{r} \frac{\rmd}{\rmd r} - \frac{\ell(\ell+1)}{r^2}.
\end{equation*}

This form of the Laplacian allows the use of the spectral methods to numerically solve Eq. (<ref>), with the spherical Bessel functions as the spectral basis.
The number of Bessel functions used corresponds to the number of grid points in position space as well as the number of modes in the spectral space, resulting in a discretized system suitable for numerical integration.
The utility of the spherical Bessel modes is that, in the spectral basis, the action of the Laplacian is purely multiplicative.

The simulations were carried out on a radial domain $[0, R]$, where the radial cutoff $R$ was chosen to ensure that the wave function remained within the interval in both position and spectral space.
The integration of Eq. (<ref>) itself is performed exactly through quadrature integration and uses an adaptive fourth-fifth order Runge-Kutta integrator.
The numerical stability was tested by ensuring that the norm was stable, the wave function remained far from the boundary in both position and spectral bases, and the results were robust to changing the relative error tolerance of the stepper.

§.§ Bounds on the first-order Trotter convergence for the hydrogen atom

To obtain the scalings presented in Table <ref> in Sec. <ref>, we apply Thm. <ref>, which poses conditions on the input eigenstate.
Since not all the eigenstates of the hydrogen atom satisfy these conditions, some particular states need a special treatment.
An example of such a state is the ground state of the hydrogen atom.

In order to apply Thm. <ref> to the hydrogen atom, we first notice that $H_\mathrm{hydrogen}$ is self-adjoint on $\mathcal{D}(H_{1})$ and $\mathcal{D}(H_{1})\subset\mathcal{D}(H_{2})$ <cit.>.
By Trotter's theorem <cit.>, we thus know that the Trotter product formula converges, i.e. for all $\psi\in\mathcal{H}$,

\begin{equation*}
\lim_{N\rightarrow\infty}\left[\left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}\right)^{N}-\rme^{-\rmi t\left(H_{1}+H_{2}\right)}\right]\psi=0.
\end{equation*}

In addition, Thm. <ref> requires an eigenstate $\varphi$ to be in the domains of $H_1^2$ and $H_2^2$.
In the case of the eigenfunctions $\Psi_{n\ell m}$ of the hydrogen atom, we find this condition to be satisfied by all states with $\ell\geq2$.
However, for the kinetic energy $H_1$, we have that all the states with the quantum number $\ell=0,1$ of the orbital angular momentum do not satisfy this condition, i.e.  $\Psi_{n00},\Psi_{n1m}\not\in\mathcal{D}(H_1^2)$.
Furthermore, it is broken by the states with $\ell=0$ in the case of the potential energy, $\Psi_{n00}\not\in\mathcal{D}(H_2^2)$.
We show in Appendix <ref> how to apply Thm. <ref> to the Hamiltonian of the hydrogen atom and also derive bounds for $\Psi_{n00}$ and $\Psi_{n1m}$ separately.
This is done by applying Thm. <ref>, which only requires the state to be in the domains of $H_1$ and $H_2$ rather than those of $H_1^2$ and $H_2^2$, and is therefore suitable for treating $\Psi_{n00}$ and $\Psi_{n1m}$.

Succinctly, we obtain for the $\ell=0$ eigenstates (including the ground state $n=1$),

\begin{align}
\xi_N(t;\Psi_{n00})
\leq{}&
\frac{\tilde{t}^{\,5/4}}{N^{1/4}}
\frac{4}{\sqrt{5\sqrt{\pi}\,n^3}}
\frac{\tilde{t}^{\,3/2}}{N^{1/2}}
\sqrt{
\frac{4\pi}{3n^3}
\frac{\tilde{t}^{\,2}}{2N}
\left(
\sqrt{\frac{2}{n^7}}
\frac{1}{4n^4}
\right),\label{eq:bound_R_n0}
\end{align}

for the $\ell=1$ states,

\begin{align}
\xi_N(t;\Psi_{n1m})
\leq{}&
\frac{\tilde{t}^{\,7/4}}{N^{3/4}}
\sqrt{
\frac{64(n^2-1)}{189\sqrt{\pi}\,n^{5}}
\frac{\tilde{t}^{\,2}}{2N}
\left(
\sqrt{
\frac{6-4/n^2}{15n^3}
\sqrt{\frac{2}{3n^7}}
\frac{1}{4n^4}
\right),
\label{eq:bound_R_n1}
\end{align}

and finally for the $\ell\geq2$ states,

\begin{align}
&\xi_N(t;\Psi_{n\ell m})\nonumber\\
\leq
\frac{\tilde{t}^{\,2}}{2N}
\,\Biggl(
\frac{1}{8}
\sqrt{
\frac{(n+\ell)!\beta(\ell-3/2,\ell+11/2)}{n^7(n-\ell-1)![\Gamma(\ell+3/2)]^2}
\nonumber\\
\sqrt{
\frac{3-\ell(\ell+1)/n^2}{2n^3\ell(\ell+1)(\ell+1/2)(\ell+3/2)(\ell-1/2)}
\nonumber\\
\frac{1}{\sqrt{n^{7}(\ell+1/2)}}
\frac{1}{4n^4}
\Biggr),
\label{eq:bound_R_nl}
\end{align}

\begin{equation}
\tilde{t}=\frac{t}{t_0}, \quad \text{with} \quad t_0=\frac{m_\mathrm{e}a_0^2}{\hbar}
\end{equation}

is a reduced (dimensionless) time,
$\beta(x,y)=\Gamma(x)\Gamma(y)/\Gamma(x+y)$ is the Euler beta function and $\Gamma(x)$ is the gamma function.

§.§ Method for the higher-order Trotter bounds

Here, we present our method to obtain a bound on the state-dependent Trotter error $\xi_N^{(p)}(t;\varphi)$ defined in Eq. (<ref>) for an arbitrary order $p$.
The technical details and the proofs which form the basis of this section can be found in Appendix <ref>.
The idea of the method is to iterate the integration-by-part strategy developed in Ref. [69] to bound evolutions.

We first focus on a single Trotter cycle

\begin{equation}
\mathcal{S}_1^{(p)}(t/N)=\rme^{-\rmi\frac{t}{N}\tau_MH_{\sigma(M)}}\cdots \rme^{-\rmi\frac{t}{N}\tau_2H_{\sigma(2)}}\rme^{-\rmi\frac{t}{N}\tau_1H_{\sigma(1)}},
\label{eq:single_Trotter_cycle}
\end{equation}

and regard the sequence of the $M$ exponentials $\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}$ ($j=1,\ldots,M$) as the evolution under a piecewise-constant time-dependent Hamiltonian $H^{(p)}(s)$, defined by

\begin{equation}
=\tau_j H_{\sigma(j)}\ \ \mathrm{for}\ \ s\in[T_{j-1},T_j),
\label{eq:piecewise-constant}
\end{equation}

$T_j = j t/N$ for $j=1,\ldots,M$.
Without loss of generality, we shift the Hamiltonians $H_1$ and $H_2$ such that $(H_1+H_2)\varphi=0$, and the target evolution is the identity $I$ on the eigenstate $\varphi$.
Then, in Appendix <ref>, we show that the difference $\mathcal{S}_1^{(p)}(t/N)-I$
can be written, via iterating the integration by parts $q$ times, as

\begin{align}
\nonumber\\
\nonumber\\
+(-\rmi)^{q+1}\int_0^{T_M}\rmd s\,U^{(p)}(T_M,s)H^{(p)}(s)S_q(s),
\label{eq:iterated_formula}
\end{align}

\begin{gather}
S_k(s)=\int_0^s\rmd s_1\cdots \int_0^{s_{k-1}}\rmd s_k\,
H^{(p)}(s_1)\cdots H^{(p)}(s_k),
\label{eq:action}
\\
-\rmi\int_s^t\rmd u\,H^{(p)}(u)
\right),\nonumber
\end{gather}

and $\T$ denotes time-ordering.
Note that $\mathcal{S}_1^{(p)}(t/N)=U^{(p)}(T_M,0)$, and Eq. (<ref>) is essentially the Dyson expansion of $U^{(p)}(T_M,0)$.

A proper $p$th-order Trotter product formula should satisfy

\begin{equation}
S_k(T_M)=\frac{1}{k!}\frac{t^k}{N^k}(H_1+H_2)^k,\quad k=1,\ldots,p.
\label{eq:vanishing_actions}
\end{equation}

In other words, we can find a $p$th-order Trotter product formula, i.e. a set of times $\{\tau_j\}$, by imposing the conditions in Eq. (<ref>).
We hence have

\begin{equation*}
S_k(T_M)\varphi=0,\quad k=1,\ldots,p,
\end{equation*}

and the Trotter error $\xi_1^{(p)}(t/N;\varphi)=\|[\mathcal{S}_1^{(p)}(t/N)-I]\varphi\|$ on the input state $\varphi$ can be estimated by the last remainder term of the Dyson expansion (<ref>), with some $q\le p$.
If $\varphi\in\mathcal{D}\bm{(}H^{(p)}(s)S_p(s)\bm{)}$ for all $s\in[0,T_M)$ and the remainder operator with $q=p$ can be applied to the input state $\varphi$, the error of the single Trotter cycle is bounded by

\begin{equation}
\xi_1^{(p)}(t/N;\varphi)
\le\int_0^{T_M}\rmd s\,\Vert H^{(p)}(s) S_p(s) \varphi\Vert,
\label{eq:higher_order_bound_starting_point_single}
\end{equation}

and the overall error $\xi_N^{(p)}(t;\varphi)$ of the $p$th-order Trotter product $\mathcal{S}_N^{(p)}(t)$ is bounded by

\begin{equation}
\xi_N^{(p)}(t;\varphi)\le N\xi_1^{(p)}(t/N;\varphi),
\label{eq:higher_order_bound_starting_point}
\end{equation}

which follows from a standard telescoping sum identity, see e.g. Ref. <cit.>.

In order to bound Eq. (<ref>), we need to evaluate $S_p(s)$.
This can be done by explicitly calculating the integrals in Eq. (<ref>).
Recall that the single Trotter cycle $\mathcal{S}_1^{(p)}(t/N)$ consists of $M$ time slots, see Eq. (<ref>).
In the $j$th time slot, the Hamiltonian is constant, $H^{(p)}(s)=\tau_j H_{\sigma(j)}$, for time $t/N$, as defined in Eq. (<ref>).
For $s\in[T_{j-1},T_j)$ in the $j$th time slot, we get

\begin{align}
\mathop{\sum_{\bm{\beta}\in\mathbb{N}^j}}_{|\bm{\beta}|=k} \left(
\frac{t}{N}
\right)^{k-\beta_j}
\frac{1}{\bm{\beta}!}
\left[\left(s-\frac{(j-1)t}{N}\right) \tau_j H_{\sigma(j)}\right]^{\beta_j}
\nonumber\\[-2truemm]
&\qquad\qquad\quad\ \
(\tau_{j-1} H_{\sigma(j-1)})^{\beta_j}\cdots
(\tau_1 H_{\sigma(1)})^{\beta_1},
\label{eq:kth_action_explicit}
\end{align}

where $\bm{\beta}=(\beta_1,\dots,\beta_j)$ is a multi-index of nonnegative integers, $|\bm{\beta}|=\beta_1+\cdots+\beta_j$, and $\bm{\beta}!=\beta_1!\cdots\beta_j!$.
At the end of the Trotter cycle $s=T_M=Mt/N$, it reads

\begin{equation}
\left(
\frac{t}{N}
\right)^k
\mathop{\sum_{\bm{\beta}\in\mathbb{N}^M}}_{|\bm{\beta}|=k}
\frac{1}{\bm{\beta}!}
(\tau_M H_{\sigma(M)})^{\beta_M}
\cdots
\label{eq:kth_action_explicit_TM}
\end{equation}

Now, as we have an explicit expression of $S_p(s)$, we can estimate a bound on Eq. (<ref>).
By dividing the total time interval $[0, Mt/N)$ into $M$ time slots and by using the triangle inequality, we get

\begin{equation}
\xi_N^{(p)}(t;\varphi)\leq N \sum_{j=1}^M \int_{(j-1)t/N}^{j t/N}\rmd s\,\Vert \tau_j H_{\sigma(j)} S_p(s)\varphi\Vert.
\label{eq:bound_higher_order_integral}
\end{equation}

Equation (<ref>) gives a bound on a higher-order Trotterization.
Such a bound can be further simplified by using identities due to $(H_1+H_2)\varphi=0$, i.e.

\begin{equation*}
\end{equation*}

We provide a Mathematica script that performs this computation for given switching times $\{\tau_j\}$ in Appendix <ref>.
Furthermore, we give an example for $p=2$ in Appendix <ref>, where we perform the calculation analytically.
The last step is to shift back the energy of the target Hamiltonian so that indeed $(H_1+H_2)\varphi=h\varphi$.
This is trivial if the input state $\varphi$ satisfies the corresponding domain conditions.
It is done by distributing the energy $h$ to the individual Hamiltonians $H_1$ and $H_2$.
Since it is arbitrary in which way $h$ is distributed, we can replace $H_1\rightarrow {H}_1(g)=H_1-g$ and $H_2\rightarrow {H}_2(g)=H_2-h+g$ with any $g\in\mathbb{R}$.
If $\varphi$ does not satisfy the domain conditions and admits a fractional scaling $\mathcal{O}(N^{-\delta})$, the energy shift will add additional terms that scale faster than $N^{-\delta}$.
They are computed explicitly in Eq. (<ref>) of Appendix <ref>.
For example, for states $\varphi$ that yield scalings slower than $\mathcal{O}(N^{-1})$ in the case of the second-order Trotterization, combining Eq. (<ref>) with Eq. (<ref>) for $H_1$ gives

\begin{equation*}
\xi_N^{(2)}(t;\varphi)\leq \frac{1}{2}b_N(t;\varphi),
\end{equation*}

where $b_N(t;\varphi)$ is our bound for such states on the first-order Trotterization.
That is, $\xi_N(t;\varphi)\leq b_N(t;\varphi)$, see Eq. (<ref>).

§ GENERALIZATIONS AND OUTLOOK

While the methods derived in Sec. <ref> are used to compute bounds for higher-order Trotterization, they are more general than that.
For this reason, they might have further applications in various other fields.
In particular, Eq. (<ref>) provides a closed form for the Dyson expansion of any unitary $U$ generated by a piecewise-constant time-dependent Hamiltonian.
Since the Dyson expansion arises by repeatedly reinserting a differential equation (the Schrödinger equation) into itself, our method might be useful for solving general (piecewise-constant) linear differential equations.
Of course, in the case of non-self-adjoint generators, $U$ will not be unitary and some adjustments have to be made.
However, for applications in quantum mechanics with a piecewise-constant Hamiltonian, Eq. (<ref>) can be used directly.
For this, one only has to omit the map $\sigma$, which accounts for the periodicity of the Trotter Hamiltonian $H^{(p)}(s)$.
A potential application could be quantum control theory when considering time-independent drifts and piecewise-constant controls.

In the same manner, our methods from Sec. <ref> can be used to compute state-dependent bounds on Trotter errors for more than two operators.
For example, for a first-order Trotter product ($\tau_j=1$ for all $j=0,\ldots,L-1$) with a target Hamiltonian $H=\sum_{j=0}^{L-1} H_j$, we can compute the integral action at the $j$th time step as

\begin{equation}
S_1(s) = \left(s-\frac{(j-1)t}{N}\right) H_{j\ \text{mod}\ L} + \frac{t}{N}\sum_{i=1}^{j-1} H_{i\ \text{mod}\ L},
\label{eq:first_order_more_operators}
\end{equation}

for $s\in[T_{j-1},T_j)$.
Similar considerations can be carried out for higher-order schemes with $L$ summands in the target Hamiltonian as well.

For quantum simulations, there is a crucial consequence from the slower Trotter error scaling we observe for the ground state of the hydrogen atom.
One advantage of a quantum approach is the ability to use the large Hilbert space to store a big, hence very accurate, approximation of the Hamiltonian.
However, for states with low angular momenta, the higher the truncation dimension is, the more Trotter steps are necessary to arrive at the desired $N^{-1}$ scaling regime.
This shows that there is an intricate interplay between the error due to truncation and the Trotter error.
See also Ref. [40].
It would be interesting to see whether there is an optimal choice of parameters for the truncation dimension and the number of Trotter steps so that the total simulation error is minimal.
Some ideas from bosonic systems might potentially help tackle this problem [74, 75].

One might suspect that the slower scaling also occurs for other atoms and molecules than the hydrogen atom.
This is suggestive since the scaling arises from the unbounded nature of the involved Hamiltonians, such as the $1/r$ singularity of the Coulomb potential.
For this reason, one would expect that states with a large electron density close to Coulomb singularities admit slower Trotter scalings.
More specifically, for the eigenstates of many-electron systems, one usually makes an ansatz through a linear combination of atomic orbitals inspired by the eigenfunctions of the hydrogen atom [44, 43].
As we have already seen in Eq. (<ref>), such states can suffer from a slower scaling as well.
This gives rise to the hypothesis that molecular states overlapping the atomic orbitals with a high probability distribution close to one of the Coulomb singularities might reveal a reduced Trotter error scaling.
In addition, further electron-electron singularities can occur in more complicated atomic and molecular systems.
These could potentially slow the error scaling even further.
This conjecture is bolstered by the fact that some first signs of slower Trotter convergence for low-energy molecular orbitals have already been observed in Ref. [28].
The authors speculate that this issue could potentially be resolved by investigating higher-order Trotter schemes.
However, our results suggest that this is not the case.

Furthermore, one might expect similar effects of slower algorithmic convergence to appear in other quantum chemistry methods such as qubitization and linear combination of unitaries.
For example, a recent paper by Mukhopadhyay et al. [76] compares a qubitization algorithm with a Trotter algorithm.
In both cases, the approximation error diverges in the full-system limit of no truncation.
This indicates that the effects of unbounded operators play a role in this scheme as well.
It would be interesting to further analyze how exactly the truncation error enters here.

§ CONCLUSION

In summary, we have provided

* state-dependent bounds for the first-order Trotterization valid for unbounded operators,

* state-dependent bounds for the higher-order Trotterizations valid for unbounded operators, and

* bounds for the first- and higher-order Trotterizations, which admit slower fractional scalings when the input state does not satisfy the respective domain conditions.

Our methods naturally generalize to Trotter problems with more than two operators.
These results are very general and might have wide-ranging applications in different fields, where the Trotter product formula is used.
This includes quantum computing, quantum control, quantum field theory, open quantum systems, quantum optics, Floquet physics, as well as quantum many-body physics [40].
We discussed direct implications for quantum chemistry problems by studying the hydrogen atom.
For those applications, a state-dependent analysis is necessary as the involved Hamiltonians are unbounded.
Our results show three implications.

* There are input states that give slower scalings than the standard $N^{-1}$ scaling. For the hydrogen atom, the scaling is solely determined by the quantum number $\ell$ of the orbital angular momentum. Most importantly, the $\ell=0$ states including the ground state yield the $N^{-1/4}$ scaling. This behavior is due to the singularity of the Coulomb potential. It can be explained by the high electron density of low-angular-momentum states close to that singularity. This leads to diverging fourth moments in the kinetic and potential energies, due to which the Trotter dynamics can ionize the hydrogen atom. We conjecture that this behavior is generic to quantum chemistry as the Coulomb potential is ubiquitous there.

* Input states with slower Trotter convergence do not benefit from an enhanced scaling by the higher-order Trotter-Suzuki methods. For a $p$th-order product formula to give the scaling $\mathcal{O}(t^{p+1}/N^p)$, the moments of the kinetic and potential energies of order $2(p+1)$ must be finite. This condition does not hold for states with low $\ell$, such as the ground state.

* For states with finite higher moments of energy, both the first- and higher-order methods might lead to the desired $N^{-p}$ scaling. However, the order of the operators in the Trotter product can determine the provable scaling behavior. This is due to an asymmetry in the conditions for the higher moments.

Our findings imply that the first-order Trotterization is the most efficient scheme for the ground-state quantum chemistry, as it involves the least number of unitaries per Trotter cycle.
Furthermore, some previous resource and runtime estimates based on the $N^{-1}$ Trotter error scaling are too optimistic.
In practice, this means that the timescale towards useful chemistry simulations on quantum computers is longer than expected.
Nevertheless, our findings only imply a polynomially increased Trotter error and do not rule out a potential exponential quantum advantage.
In fact, the slower scaling also carries over to numerical simulations on classical computers.
Overall, we treat errors in digital simulations in an analytical way.
We hope that our analytical insights will inspire the development of more efficient (quantum) simulation algorithms and provide a methodology for the proof of exponential quantum advantage.

The authors thank Lauritz van Luijk for helpful feedback and for pointing out the argument for the non-existence of general state-dependent Trotter bounds with commutator scaling from Thm. <ref>.
DB and AH would like to thank Elliot Eklund and Ivan Kassal for their interesting and helpful discussions.
DB further acknowledges Reinhard F. Werner for a fruitful dialogue.
In addition, AH thanks Dominic Berry and Nathan K. Langford for stimulating conversations.

DB acknowledges funding from the Australian Research Council (project numbers FT190100106, DP210101367, CE170100009) and the Munich Quantum Valley project K8.
PF acknowledges financial
support from the PNRR MUR project
CN00000013 - National Centre for HPC, Big Data and
Quantum Computing, as well as partial support by Istituto Nazionale di Fisica Nucleare (INFN) through
the project “QUANTUM” and by the Italian National Group of Mathematical Physics (GNFM-INdAM).
AH was partially supported by the Sydney Quantum Academy.
KY acknowledges support by the Top Global University Project from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan, and by JSPS KAKENHI Grant Number JP18K03470 and JSPS KAKENHI Grant Number JP18KK0073 both from the Japan Society for the Promotion of Science (JSPS).

§ FIRST-ORDER TROTTER BOUNDS (THM. <REF> AND THM. <REF>)

In this appendix, we provide the proofs of Thm. <ref> and Thm. <ref> presented in Sec. <ref>.
We start by proving a lemma.
In the following, we denote by $I$ the identity operator.

Let $H_{1}$ be self-adjoint on $\mathcal{D}(H_{1})$ and $H_{2}$
be self-adjoint on $\mathcal{D}(H_{2})$. Let $\varphi\in\mathcal{D}(H_{1})\cap\mathcal{D}(H_{2})$
be an eigenstate of $H_{1}+H_{2}$ with eigenvalue $0$, i.e. $(H_{1}+H_{2})\varphi=0$.
Then, for all $t\in\mathbb{R}$,

\begin{equation}
\xi_N(t;\varphi)
\le N\left\Vert \left(\rme^{-\rmi \frac{t}{N}H_{1}}-I+\rmi\frac{t}{N}H_{1}\right)\varphi\right\Vert
+N \left\Vert \left(\rme^{\rmi\frac{t}{N}H_{2}}-I-\rmi\frac{t}{N}H_{2}\right)\varphi\right\Vert .
\label{eq:lemma_bound1}
\end{equation}

Moreover, if $\varphi\in\mathcal{D}(H_{1}^2)\cap\mathcal{D}(H_{2}^2)$, it is further bounded by

\begin{equation}
\xi_N(t;\varphi)
\le\frac{t^2}{2N}\,\Bigl(
\Vert H_1^2\varphi\Vert
\Vert H_2^2\varphi\Vert
\Bigr).
\label{eq:lemma_bound2}
\end{equation}

The definition of $\xi_N(t;\varphi)$ is given in Eq. (<ref>).
Under the assumption $(H_1+H_2)\varphi=0$, it reads

\begin{equation*}
\xi_N(t;\varphi)=\left\Vert\left[\left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}\right)^{N}-I\right]\varphi\right\Vert.
\end{equation*}

By telescoping the sum, we have

\begin{equation*}
D_N(t)\equiv\left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}\right)^{N}-I
=\sum_{\ell=1}^{N}\left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}\right)^{N-\ell}\left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}-I\right).
\end{equation*}

Then, using the invariance of the norm under unitaries and the assumption $(H_{1}+H_{2})\varphi=0$, we get

\begin{align}
\xi_N(t;\varphi)
\Vert D_N(t)\varphi\Vert
N\left\Vert \left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}-I\right)\varphi\right\Vert
\nonumber\\
&=N\left\Vert \left(\rme^{-\rmi \frac{t}{N}H_{1}}-\rme^{\rmi \frac{t}{N}H_{2}}\right)\varphi\right\Vert
\nonumber\\
&=N\left\Vert \left(\rme^{-\rmi \frac{t}{N}H_{1}}-I+\rmi\frac{t}{N}H_{1}-\rme^{\rmi\frac{t}{N}H_{2}}+I+\rmi\frac{t}{N}H_{2}\right)\varphi\right\Vert \nonumber\\
& \le N\left\Vert \left(\rme^{-\rmi \frac{t}{N}H_{1}}-I+\rmi\frac{t}{N}H_{1}\right)\varphi\right\Vert +N \left\Vert \left(\rme^{\rmi\frac{t}{N}H_{2}}-I-\rmi\frac{t}{N}H_{2}\right)\varphi\right\Vert,
\label{eq:first_order_bound}
\end{align}

which proves the first bound (<ref>).
Under the additional domain condition $\varphi\in\mathcal{D}(H_{1}^2)\cap\mathcal{D}(H_{2}^2)$, we can further bound it by Eq. (<ref>), using

\begin{equation}
(\rme^{-\rmi sH_j}-I+\rmi sH_j)\varphi
-\int_0^s\rmd u\,\rme^{-\rmi (s-u)H_j}uH_j^2 \varphi,
\label{eq:Taylor_first}
\end{equation}

and hence

\begin{equation}
\|(\rme^{-\rmi sH_j}-I+\rmi sH_j)\varphi\|
\le
\frac{1}{2}s^2\|H_j^2\varphi\|,
\label{eq:Taylor2}
\end{equation}

for all $s\in\mathbb{R}$ and $j=1,2$.

§.§ Proof of Thm. <ref>

We are now ready to prove Thm. <ref>.

We shift the eigenvalue of the target Hamiltonian $H_1+H_2$ with respect to $\varphi$ to zero and apply the second bound (<ref>) of Lemma <ref>:

\begin{align*}
\xi_N(t;\varphi)
={}&\left\Vert\left[\left(\rme^{-\rmi \frac{t}{N}H_2}\rme^{-\rmi \frac{t}{N}H_1}\right)^{N}-\rme^{-\rmi th}\right]\varphi\right\Vert
\nonumber\\
={}&\left\Vert\left[\left(\rme^{-\rmi \frac{t}{N}(H_2+g-h)}\rme^{-\rmi \frac{t}{N}(H_1-g)}\right)^{N}-I\right]\varphi\right\Vert
\nonumber\\
\le{}&
N\left\Vert \left(\rme^{-\rmi \frac{t}{N} (H_1-g)}-I+\rmi\frac{t}{N} (H_1-g)\right)\varphi\right\Vert
+N \left\Vert \left(\rme^{\rmi\frac{t}{N}(H_2+g-h)}-I-\rmi\frac{t}{N}(H_2+g-h)\right)\varphi\right\Vert
\nonumber\\
\le{}&\frac{t^{2}}{2N}\,\Bigl(\Vert (H_1-g)^{2}\varphi\Vert+\Vert (H_2+g-h)^{2}\varphi\Vert\Bigr).
\end{align*}

This proves Thm. <ref>.

§.§ Proof of Thm. <ref>

If the domain condition $\varphi\in\mathcal{D}(H_{1}^2)\cap\mathcal{D}(H_{2}^2)$ is not fulfilled, the bound (<ref>) is not useful and the second bound (<ref>) of Lemma <ref> does not hold.
In such a case, we need to perform a refined analysis to estimate the scaling of the first bound (<ref>) of Lemma <ref>.
Let us look at

\begin{equation}
\Delta_j(s)
=\|(\rme^{-\rmi sH_j}-I+\rmi sH_j)\varphi\|^2
=\int_{\mathbb{R}}|\rme^{-\rmi s\lambda}-1+\rmi s\lambda|^2
\,\rmd\mu_{j,\varphi}(\lambda),\quad j=1,2,
\label{eq:Delta1_def}
\end{equation}

for $s\in\mathbb{R}$, where

\begin{equation*}
\mu_{j,\varphi}(\Omega)= \| \mathbf{1}_{\Omega}(H_j)\,\varphi\|^2,
\end{equation*}

for a measurable $\Omega\subset \mathbb{R}$, is the spectral measure of $H_j$ at the state $\varphi$, with $\mathbf{1}_{\Omega}(x)$ being the characteristic function of $\Omega$, which equals 1 for $x\in\Omega$ and 0 for $x\notin\Omega$.
Note that the convergence of the integral in Eq. (<ref>) only requires $\varphi\in\mathcal{D}(H_j)$.
Since $f(x)= |\rme^{-\rmi x}-1+\rmi x|^2= f(-x)$, it suffices to study

\begin{equation*}
\Delta_j(s)
=\int_{\mathbb{R}_+} f(s\lambda)\,\rmd\nu_{j,\varphi}(\lambda),
\end{equation*}

where $\mathbb{R}_+=[0,+\infty)$, and

\begin{equation*}
\nu_{j,\varphi}(\Omega) = \mu_{j,\varphi}\bm{(}\Omega \cup (-\Omega) \bm{)},\quad \Omega \subset\mathbb{R}_+.
\end{equation*}

The asymptotic behavior of $\Delta_j(s)$ for small $s$ is ruled by the high-energy tail of the measure $\nu_{j,\varphi}$.
The following Lemma provides a bound on $\Delta_j(s)$, when we assume a suitable decay of $\nu_{j,\varphi}\bm{(}(\lambda,+\infty)\bm{)}$.

Let $\mu_\varphi(\Omega)$ be the spectral measure of a self-adjoint operator $H$ at the vector $\varphi$. Assume that

\begin{equation*}
\mu_{\varphi}(\{|x|\geq \lambda \}) = \nu_{\varphi}\bm{(}(\lambda,+\infty)\bm{)} \leq \left(\frac{\lambda_0}{\lambda}\right)^{2\delta},
\end{equation*}

for some $\lambda_0>0$ and $\delta>1$, and for all $\lambda>0$.

\begin{equation}
\Delta(s)
=\int_{\mathbb{R}_+}|\rme^{-\rmi s\lambda}-1+\rmi s\lambda|^2
\,\rmd\nu_{\varphi}(\lambda)
\label{eq:Delta(s)def}
\end{equation}

is bounded by

\begin{equation}
\Delta(s)
\le
\begin{cases}
\medskip
\displaystyle
\frac{\pi}{\Gamma(2\delta-1)\sin[(\delta-1)\pi]}(\lambda_0 s)^{2\delta},
\quad&1<\delta<2,\\
\medskip
\displaystyle
\frac{1}{4}\Lambda^2\|H\varphi\|^2s^4
g(\Lambda s),
\quad&\delta=2,\\
\displaystyle
\frac{1}{4}\|H^2\varphi\|^2s^4,&
\delta>2,
\end{cases}
\label{eqn:fractional_bound}
\end{equation}

for $s>0$, where $\Lambda>0$ is an arbitrary energy scale, and

\begin{equation*}
-x^4\Ci(x)+\left(1+\frac{1}{2}x^2\right)[(1-\cos x)^2+(x-\sin x)^2]-2x^3(x-\sin x)+\frac{3}{2}x^4,
\end{equation*}

$\Ci(x)=-\int_x^\infty\frac{\cos t}{t}\,\rmd t$.
Note that

\begin{equation*}
g(x)=-x^4\log x+\left(\frac{7}{4}-\gamma\right)x^4+\mathcal{O}(x^6),
\end{equation*}

as $x\downarrow 0$, with Euler's constant $\gamma=\Gamma'(1)\approx 0.577$.
In particular, one gets

\begin{equation}
g(x) \leq -x^4\log x+ 1.2003 x^4,
\label{eq:g(x)smallx}
\end{equation}

for $0\leq x\leq 1$.

For $1<\delta<2$, we use integration by parts.

\begin{equation*}
=|\rme^{-\rmi x}-1+\rmi x|^2
=(1-\cos x)^2+(x-\sin x)^2 .
\end{equation*}

Note that $f\in L^1(\mathbb{R}_+,\rmd\nu_{\varphi})$, $f'(\lambda)/\lambda^{2\delta}\in L^1(\mathrm{R}_+)$, $f(0)=0$, and $f(\lambda)=o(\lambda^{2\delta})$ as $\lambda\to+\infty$.
Since $\rmd\nu_{\varphi}(\lambda)= - \rmd[\nu_{\varphi}\bm{(}(\lambda,+\infty)\bm{)}]$, we get

\begin{align*}
\Delta(s)
-\int_{\mathbb{R}_+} \rmd[f(s\lambda)\nu_{\varphi}\bm{(}(\lambda,+\infty)\bm{)}]
f(0) \nu_{\varphi}\bm{(}(0,+\infty)\bm{)}
- \lim_{\lambda\rightarrow+\infty} f(s\lambda)\nu_{\varphi}\bm{(}(\lambda,+\infty)\bm{)}
+ \int_0^{+\infty}f'(x)\nu_{\varphi}\bm{(}(x/s,+\infty)\bm{)}\,\rmd x.
\end{align*}

\begin{align*}
\Delta(s)
\int_0^{+\infty}f'(x)\nu_{\varphi}\bm{(}(x/s,+\infty)\bm{)}\,\rmd x\\
\leq{}& (\lambda_0 s)^{2\delta} \int_0^{+\infty}\frac{f'(x)}{x^{2\delta}}\,\rmd x\\
\frac{\pi}{\Gamma(2\delta-1)\sin[(\delta-1)\pi]}(\lambda_0 s)^{2\delta}.
\end{align*}

For $\delta=2$, we fix $\Lambda>0$ and split the integration range as

\begin{equation*}
\Delta(s)
=\int_{[0,\Lambda]} f(s\lambda)\,\rmd\nu_{\varphi}(\lambda)
+\int_{(\Lambda,+\infty)} f(s\lambda)\,\rmd\nu_{\varphi}(\lambda)
\end{equation*}

Since $f(x)\le x^4/4$, the first term is estimated as

\begin{align*}
\Delta^{<}(s)
\int_{[0,\Lambda]} f(s\lambda)\,\rmd\nu_{\varphi}(\lambda)
\nonumber\\
\le{}&
\frac{1}{4}s^4
\int_{[0,\Lambda]} \lambda^4\,\rmd\nu_{\varphi}(\lambda)
\nonumber\\
\leq{}&
\frac{1}{4}\Lambda^2 s^4
\int_{[0,\Lambda]}\lambda^2\,\rmd\nu_{\varphi}(\lambda)
\nonumber\\
\le{}&
\frac{1}{4}\Lambda^2 s^4
\int_{\mathbb{R}_+}\lambda^2\,\rmd\nu_{\varphi}(\lambda)
\nonumber\\
\frac{1}{4}\Lambda^2\|H\varphi\|^2 s^4.
\intertext{Since
$f'(x)= 2x (1-\cos x) \geq 0$ for $x\geq0$,
the second term is bounded by}
\Delta^{>}(s)
\int_{(\Lambda,+\infty)} f(s\lambda)\,\rmd\nu_{\varphi}(\lambda)
\nonumber\\
={}& {-}\int_{(\Lambda,+\infty)}f(s\lambda)\,\rmd[\nu_{\varphi}\bm{(}(\lambda,+\infty)\bm{)}]
\nonumber\\
={}& f(s\Lambda)\nu_{\varphi}\bm{(}(\Lambda,+\infty)\bm{)} + s \int_{\Lambda}^{+\infty} f'(s\lambda)\nu_{\varphi}\bm{(}(\lambda,+\infty)\bm{)}\,\rmd\lambda
\nonumber\\
\leq{}& f(s\Lambda)\frac{\lambda_0^4}{\Lambda^4} + s \int_{\Lambda}^{+\infty} f'(s\lambda)\frac{\lambda_0^4}{\lambda^4}\,\rmd\lambda
\nonumber\\
={}& 4 \lambda_0^4 \int_{\Lambda}^{+\infty} \frac{f(s\lambda)}{\lambda^5}\,\rmd\lambda
\nonumber\\
4 (\lambda_0 s)^4\int_{s\Lambda}^{+\infty}\frac{f(x)}{x^5}\,\rmd x
\nonumber
\\
\frac{\lambda_0^4}{\Lambda^4}g(\Lambda s).
\end{align*}

Finally, if $\delta>2$ then $\varphi\in \mathcal{D}(H^2)$, and the bound is given by Eq. (<ref>) obtained in the proof of Lemma <ref>.

Let us collect some comments on this lemma in the following remark.

* $\lambda_0$ and $\Lambda$ have the dimension of energy.
$\Lambda$ is an arbitrary energy scale, which can be tuned to optimize the bound for $\delta=2$.
The optimal $\Lambda$ is the solution to the equation

\frac{f(s\Lambda)}{(s\Lambda)^6}=\frac{\|H\varphi\|^2}{8\lambda_0^4} \frac{1}{s^2}.

Since $f(x)\sim x^4/4$ as $x\to 0$, one gets that for $s\downarrow0$ the optimal value approaches $\Lambda \to \sqrt{2}\, \lambda_0^2 / \|H\varphi\|$, whence a quasi-optimal bound for small $s$ is given by

\Delta(s) \leq \frac{1}{2}(\lambda_0 s)^4 + \frac{ \|H\varphi\|^4}{4 \lambda_0^4} g\!\left(\frac{\sqrt{2}\,\lambda_0}{\|H\varphi\|} \lambda_0 s\right).

By using the upper bound (<ref>), one gets for  $s\leq \|H\varphi\|/ (\sqrt{2}\,\lambda_0^2)$

\begin{equation*}
\Delta(s) \leq (\lambda_0 s)^4 \left(-\log(\lambda_0 s) + 1.7003-\log\frac{\sqrt{2}\,\lambda_0}{\|H\varphi\|}\right) .
\end{equation*}

* The lemma can be extended to the range $\delta\leq 1$. Indeed, if one only assumes that $\varphi\in\mathcal{D}(H)$, so that $\Delta(s)$ in Eq. (<ref>) is finite, without any further assumptions on the decay of the spectral measure, one easily gets from Eq. (<ref>)

\begin{equation*}
\Delta(s)
=s^2 \int_{\mathbb{R}_+}|F(s\lambda)|^2
\lambda^2\,\rmd\nu_{\varphi}(\lambda),
\qquad F(x)= \rme^{\rmi x/2} \sinc(x/2)-1.
\end{equation*}

Since $F(x)$ is bounded and $F(x)\to 0$ as $x\to 0$, by dominated convergence one has that

\begin{equation*}
\Delta(s) =s^2 \|F(s H)H\varphi\|^2 = o(s^2).
\end{equation*}

Notice, however, that in such a case we have no control over the rate of convergence of $\Delta(s)/s^2$, which can be arbitrarily slow.

* If $\delta>2$, then $\varphi\in\mathcal{D}(H^2)$, and the bound from Thm. <ref> applies.
However, in this case, there might be a possibility of gaining a better scaling by a higher-order product formula.

Theorem <ref> is immediately proven, using the scaling analysis of Lemma <ref>.

By the same replacements as in the proof of Thm. <ref> and using the first bound (<ref>) of Lemma <ref>, we have

\begin{equation}
\xi_N(t;\varphi)
\le
\sqrt{\Delta_1(t/N)}
\sqrt{\Delta_2(t/N)}
\Bigr).
\label{eq:first_order_fractional_scaling}
\end{equation}

Under the conditions for Thm. <ref>, bounds on $\Delta_j(t/N)$ scale as Eq. (<ref>) of Lemma <ref>, i.e.,

\begin{equation*}
\Delta_j(t/N)
\begin{cases}
\medskip
\displaystyle
\mathcal{O}\bm{(}(t/N)^{2\delta}\bm{)},
\quad&1<\delta<2,\\
\displaystyle
\mathcal{O}\bm{(}(t/N)^4\log(N/t)\bm{)},
\quad&\delta=2,
\end{cases}
\end{equation*}

and we get the scalings of $\xi_N(t;\varphi)$ in Thm. <ref>.

§ ANALYTICAL TROTTER BOUNDS FOR THE HYDROGEN ATOM

In this section, we apply the Trotter error bounds from Thm. <ref> and Lemma <ref> (Thm. <ref>) to the Hamiltonian of the hydrogen atom, where we Trotterize between the kinetic energy (temporarily reinstating the Planck constant $\hbar$)

\begin{equation*}
\end{equation*}

and the potential energy

\begin{equation*}
H_2=-\frac{e^2}{4\pi\varepsilon_0 r}.
\end{equation*}

This results in the explicit bounds in Eqs. (<ref>)–(<ref>) presented in Sec. <ref> of the Methods section.
Recall that the Bohr radius $a_0$ is given by

\begin{equation*}
\end{equation*}

Furthermore, the hydrogen atom $H_1+H_2$ has eigenenergies

\begin{equation*}
E_n =-\frac{\hbar^2}{2m_\mathrm{e}a_0^2n^2},\quad n=1,2,\ldots,
\end{equation*}

corresponding to the radial hydrogen wave functions
$R_{n\ell}(r)=a_0^{-3/2}\mathcal{R}_{n\ell}(r/a_0)$ ($\ell=0,\ldots,n-1$) with

\begin{equation}
\mathcal{R}_{n\ell}(u)
\left(\frac{2}{n}\right)^{3/2}
\sqrt{
\frac{(n-\ell-1)!}{2n(n+\ell)!}
\rme^{-u/n}
\label{eq:hydrogen_eigenfunctions}
\end{equation}

\begin{equation*}
\end{equation*}

are the generalized Laguerre polynomials <cit.>.
The full eigenstates $\Psi_{n\ell m}$ of the hydrogen atom including the degeneracies (characterized by the magnetic quantum number $m=-\ell,-\ell+1,\ldots,\ell-1,\ell$) are given by

\begin{equation*}
\Psi_{n\ell m}(r, \theta, \phi) = R_{n\ell}(r) Y_{\ell m}(\theta, \phi),
\end{equation*}

where the spherical harmonics $Y_{\ell m}(\theta, \phi)$ are defined by

\begin{equation*}
Y_{\ell m}(\theta, \phi) = \sqrt{\frac{2\ell +1}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}} \rme^{\rmi m\phi} P_{\ell m}(\cos\theta),
\end{equation*}

with the associated Legendre polynomials

\begin{equation*}
P_{\ell m}(x) = \frac{(-1)^m}{2^\ell \ell!} (1-x^2)^{m/2} \frac{\rmd^{\ell+m}}{\rmd x^{\ell+m}}(x^2-1)^\ell.
\end{equation*}

However, the spherical part $Y_{\ell m}$ only accounts for the degeneracies in the spectrum of the hydrogen atom.
Different values of $m$ for fixed $(n,\ell)$ lead to the same Trotter dynamics.
For this reason, it suffices to bound the Trotter convergence speed for the radial part $R_{n\ell}(r)$.
The momentum distributions of the eigenstates of the hydrogen atom are given by $(a_0/\hbar)\Xi_{n\ell}(a_0p/\hbar)$ with [78]

\begin{equation}
\Xi_{n\ell}(u)
\frac{4^{2\ell+3}n^2(\ell!)^2(n-\ell-1)!}{2\pi(n+\ell)!}
\frac{(n^2u^2)^{\ell+1}}{(n^2u^2+1)^{2\ell+4}}
\left[
\frac{n^2u^2-1}{n^2u^2+1}
\right)
\right]^2,
\label{eq:momentum_distribution}
\end{equation}

where $C_n^{(\alpha)}(x)$ is the Gegenbauer polynomial <cit.>.
The radial wave functions and the momentum distributions are normalized as

\begin{equation*}
\int_0^\infty\rmd u\,u^2[\mathcal{R}_{n\ell}(u)]^2=1,\qquad
\int_0^\infty\rmd u\,\Xi_{n\ell}(u)=1.
\end{equation*}

We will use a bound on the generalized Laguerre polynomial [79]

\begin{equation}
\rme^{-\frac{1}{2}x}|L_n^{(\alpha)}(x)|
\le L_n^{(\alpha)}(0)
(x\ge0,\ \alpha\ge0),
\label{eqn:LaguerreBound}
\end{equation}

and a bound on the Gegenbauer polynomial <cit.>

\begin{equation}
\le C_n^{(\alpha)}(1)
(-1\le x\le1,\ \alpha>0).
\label{eqn:GegenbauerBound}
\end{equation}

In this section, we set $g=0$ in Thm. <ref> and define

\begin{equation*}
\xi_{N,j}(t;\Psi_{n\ell m}) = N\left\|\left(\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j} -I - \frac{\rmi}{\hbar}\frac{t}{N}H_j\right)\Psi_{n\ell m}\right\|
\le\frac{t^2}{2\hbar^2N}\|H_j^2 \Psi_{n\ell m}\|,\quad j=1,2,
\end{equation*}

for the kinetic and potential energies.
By computing these individually, we retrieve the Trotter bounds shown in Sec. <ref>.
Notice, however, that we shift the potential energy by the respective eigenvalue $h=E_n$ in $\xi_{N,2}(t;\Psi_{n\ell m})$ as we did in the proof of Thm. <ref> in Appendix <ref>.
We also define a characteristic energy scale and a time scale

\begin{equation}
\Lambda=\frac{\hbar^2}{m_\mathrm{e}a_0^2}, \qquad t_0=\frac{\hbar}{\Lambda}=\frac{m_\mathrm{e}a_0^2}{\hbar},
\label{eq:EnergyScale}
\end{equation}

and a reduced (dimensionless) time

\begin{equation*}
\tilde{t}=\frac{t}{t_0}.
\end{equation*}

§.§ Kinetic term for $\ell\ge 2$

For $\ell\ge2$, we have $\Psi_{n\ell m}\in\mathcal{D}(H_1^2)$ and we can just use Thm. <ref>.
That is, by setting $g=0$, we just have to bound $\|H_1^2 \Psi_{n\ell m}\|$ for

\begin{align*}
[\xi_{N,1}(t;\Psi_{n\ell m})]^2
\frac{t^4}{4\hbar^4N^2}\|H_1^2 \mathcal{R}_{n\ell}\|^2
\frac{\tilde{t}^{\,4}}{64N^2}
\int_0^\infty\rmd u\,
\Xi_{n\ell}(u),
\quad\ell\ge2.
\intertext{By bounding the Gegenbauer polynomial using Eq.~(\ref{eqn:GegenbauerBound}), we get}
[\xi_{N,1}(t;\Psi_{n\ell m})]^2
\frac{\tilde{t}^{\,4}}{64N^2}
\frac{4^{2\ell+3}n^2(\ell!)^2(n+\ell)!}{2\pi[(2\ell+1)!]^2(n-\ell-1)!}
\int_0^\infty\rmd u\,
\frac{(n^2u^2)^{\ell+1}}{(n^2u^2+1)^{2\ell+4}}
\nonumber\\
\frac{\tilde{t}^{\,4}}{16N^2}
\frac{(n+\ell)!\beta(\ell-3/2,\ell+11/2)}{n^7(n-\ell-1)![\Gamma(\ell+3/2)]^2},
\quad\ell\ge2,
\end{align*}

where $\beta(x,y)=\Gamma(x)\Gamma(y)/\Gamma(x+y)$ is the beta function.

§.§ Kinetic term for $\ell=0,1$

For $\ell=0,1$, we have diverging fourth moment $\|H_1^2 \Psi_{n\ell m}\|=\infty$, and $\Psi_{n\ell m}\not\in\mathcal{D}(H_1^2)$.
We hence need to estimate the finer bound, following Lemma <ref>.
The spectral measure $\nu_{1,\Psi_{n\ell m}}$ of the kinetic-energy operator $H_1$ at $\Psi_{n\ell m}$ is absolutely continuous,
$\rmd \nu_{1,\Psi_{n\ell m}}(\lambda)=\rho_{1,\Psi_{n\ell m}}(\lambda)\,\rmd\lambda$,
and is obtained from the momentum distribution of the eigenstate of the hydrogen atom (<ref>) as

\begin{equation*}
\rho_{1,\Psi_{n\ell m}}(\lambda)\,\rmd\lambda
\Xi_{n\ell}(u)\,\rmd u
\frac{4^{2\ell+3}n^2(\ell!)^2(n-\ell-1)!}{2\pi(n+\ell)!}
\frac{(2n^2\lambda/\Lambda)^{\ell+1}}{(2n^2\lambda/\Lambda+1)^{2\ell+4}}
\left[
\frac{2n^2\lambda/\Lambda-1}{2n^2\lambda/\Lambda+1}
\right)
\right]^2
\frac{\rmd\lambda}{\Lambda(2\lambda/\Lambda)^{1/2}},
\end{equation*}

where $\lambda=(\hbar^2/2m_\mathrm{e}a_0^2)u^2$ and
$\Lambda$ is the energy scale defined in Eq. (<ref>).
By Eq. (<ref>) and

we find that

\begin{equation}
\rho_{1,\Psi_{n\ell m}}(\lambda)
\le
\frac{2^{3\ell+3/2}(\ell!)^2(n+\ell)!}{\pi n^{2\ell+4}[(2\ell+1)!]^2(n-\ell-1)!}
\frac{\Lambda^{\ell+5/2}}{\lambda^{\ell+7/2}}.
\label{eqn:Tail1}
\end{equation}

This implies the decay condition of Lemma <ref> as

\begin{equation}
\nu_{1,\Psi_{n\ell m}}\bm{(}(\lambda,+\infty)\bm{)} = \int_\lambda^\infty \rho_{1,\Psi_{n\ell m}}(x)\,\rmd x \leq \left(\frac{\lambda_1}{\lambda}\right)^{2\delta_1},
\label{eqn:lambda1}
\end{equation}

\begin{equation*}
\delta_1
\medskip
\displaystyle
\frac{5}{4}, \quad&\ell=0,\\
\displaystyle
\frac{7}{4}, \quad&\ell=1.
\end{cases}
\end{equation*}

The bound (<ref>) from Lemma <ref> then  yields

\begin{equation}
\Delta_1(t/N)
\le
\frac{2^{3\ell+1}(\ell!)^2(n+\ell)!}{n^{2\ell+4}[(2\ell+1)!]^2(n-\ell-1)!}
\frac{2\ell+3}{\Gamma(\ell+7/2)}
\left(
\frac{\tilde{t}}{N}
\right)^{\ell+5/2},
\label{eqn:HydrogenBound1}
\end{equation}

so that

\begin{equation*}
\xi_{N,1}(t;\Psi_{n\ell m})
\leq
\frac{\tilde{t}^{\,\ell/2+5/4}}{N^{\ell/2+1/4}}
\sqrt{
\frac{2^{3\ell+1}(\ell!)^2(n+\ell)!}{n^{2\ell+4}[(2\ell+1)!]^2(n-\ell-1)!}
\frac{2\ell+3}{\Gamma(\ell+7/2)}
\quad\ell=0,1.
\end{equation*}

§.§ Potential term for $\ell\geq 1$

For $\ell\geq 1$, we have $\Psi_{n\ell m}\in\mathcal{D}(H_2^2)$ and can therefore apply Thm. <ref>.
We know from Ref. [81] that for $\ell\ge0$

\begin{align*}
\langle\Psi_{n\ell m}|\frac{1}{r}|\Psi_{n\ell m}\rangle & =\frac{1}{a_{0}n^{2}},\\
\langle\Psi_{n\ell m}|\frac{1}{r^{2}}|\Psi_{n\ell m}\rangle & =\frac{1}{a_{0}^{2}(\ell+1/2)n^{3}},
\intertext{and for $\ell\geq1$,}