\langle\Psi_{n\ell m}|\frac{1}{r^{3}}|\Psi_{n\ell m}\rangle & =\frac{1}{a_{0}^{3}\ell(\ell+1)(\ell+1/2)n^{3}},\\
\langle\Psi_{n\ell m}|\frac{1}{r^{4}}|\Psi_{n\ell m}\rangle& =\frac{3n^2-\ell(\ell+1)}{2a_{0}^{4}\ell(\ell+1)(\ell+1/2)(\ell-1/2)(\ell+3/2)n^5}.
\end{align*}

It is worth remarking that the last two moments diverge when $\ell=0$.
This is because s-orbitals have high probabilities for the electron
to be close to the center, where the potential diverges.
Using these moments, we can compute for $\ell\geq1$

\begin{align*}
\|(H_{2}-E_{n})^{2}\Psi_{n\ell m}\|^{2}
={}&\langle\Psi_{n\ell m}|(H_{2}-E_{n})^{4}|\Psi_{n\ell m}\rangle\\
={}&\langle\Psi_{n\ell m}|H_{2}^{4}|\Psi_{n\ell m}\rangle
-4E_{n}\langle\Psi_{n\ell m}|H_{2}^{3}|\Psi_{n\ell m}\rangle
\\
+6E_{n}^{2}\langle\Psi_{n\ell m}|H_{2}^{2}|\Psi_{n\ell m}\rangle
-4E_{n}^{3}\langle\Psi_{n\ell m}|H_{1}|\Psi_{n\ell m}\rangle+E_{n}^{4}
\end{align*}

to get

\begin{align*}
& \xi_{N,2}(t;\Psi_{n\ell m})\\
\frac{\tilde{t}^{\,2}}{2N}\sqrt{\frac{16n^5-32(\ell-1/2)(\ell+3/2)n^3+24\ell(\ell+1)(\ell-1/2)(\ell+3/2) n-7\ell(\ell+1)(\ell+1/2)(\ell-1/2)(\ell+3/2)}{16\ell(\ell+1)(\ell+1/2)(\ell-1/2)(\ell+3/2)n^8}}.
\end{align*}

§.§ Potential term for $\ell=0$

For $\ell=0$, we have a diverging fourth moment $\|H_2^2\Psi_{n00}\|=\infty$, and hence $\Psi_{n00}\not\in\mathcal{D}(H_2^2)$.
Therefore, we need to estimate the finer bound from Lemma <ref>.
The spectral measure $\nu_{2,\Psi_{n\ell m}}$ of the potential-energy operator $H_2$ at $\Psi_{n\ell m}$ is absolutely continuous,
$\rmd \nu_{2,\Psi_{n\ell m}}(\lambda)=\rho_{2,\Psi_{n\ell m}}(\lambda)\,\rmd\lambda$, and is obtained from the radial wave function (<ref>) as

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

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

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

This implies the decay condition of Lemma <ref> as

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

\begin{equation*}
\delta_2
\frac{3}{2},\quad\ell=0.
\end{equation*}

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

\begin{equation}
\Delta_2(t/N)
\le
\frac{4\pi}{3n^3}
\left(
\frac{\tilde{t}}{N}
\right)^3,\quad\ell=0.
\label{eqn:HydrogenBound2}
\end{equation}

Thus, $\xi_{N,2}(t,\Psi_{n00})$ is bounded by

\begin{align*}
\xi_{N,2}(t;\Psi_{n00})
\le{}&
\sqrt{\frac{4\pi}{3n^3}} \frac{\tilde{t}^{\,3/2}}{N^{1/2}},\quad\ell=0.
\end{align*}

§.§ Potential term with a shift

According to Thm. <ref>, we need to shift the potential Hamiltonian $H_2$ by the eigenvalue $E_n$ of the respective input eigenstate $\Psi_{n\ell m}$ to get the bound on the Trotter error.
This can be done as follows

\begin{align}
\xi_{N,j}(t;\Psi_{n\ell m})
N\left\|\left(\rme^{\frac{\rmi}{\hbar}\frac{t}{N}(H_j-g_j)} -I - \frac{\rmi}{\hbar}\frac{t}{N}(H_j-g_j)\right)\Psi_{n\ell m}\right\|
\nonumber\\
\left[
\left(
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j} -I
\right)
\left(
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}(H_j-g_j)}-\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j}
\frac{\rmi}{\hbar}\frac{t}{N}g_j
\right)
\right]
\Psi_{n\ell m}\right\|
\nonumber\\
\left[
\left(
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j} -I
\right)
\left(
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j}
\rme^{-\frac{\rmi}{\hbar}\frac{t}{N}g_j}
\frac{\rmi}{\hbar}\frac{t}{N}g_j
\right)
\right]
\Psi_{n\ell m}\right\|
\nonumber\\
\left(
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j} -I
\right)
\Psi_{n\ell m}
\left[
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j}
\left(
\rme^{-\frac{\rmi}{\hbar}\frac{t}{N}g_j}
\frac{\rmi}{\hbar}\frac{t}{N}g_j
\right)
\frac{\rmi}{\hbar}\frac{t}{N}g_j
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j}
\right]
\Psi_{n\ell m}
\biggr\|
\nonumber\\
\le{}&
\left(
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j} -I
\right)
\Psi_{n\ell m}
\right\|
\left|
\rme^{-\frac{\rmi}{\hbar}\frac{t}{N}g_j}
\frac{\rmi}{\hbar}\frac{t}{N}g_j
\right|
\frac{t}{\hbar}|g_j|
\|
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j}
\Psi_{n\ell m}
\|
\nonumber\\
\le{}&
\left(
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}H_j} -I
\right)
\Psi_{n\ell m}
\right\|
\frac{g_j^2t^2}{2N\hbar^2}
\frac{|g_j|t^2}{N\hbar^2}\|H_j\Psi_{n\ell m}\|.
\label{eq:first_order_with shift}
\end{align}

For instance, for the potential term $H_2$ with the shift $g_2=E_n=-\hbar^2/2m_\mathrm{e}a_0^2n^2$, we get

\begin{align*}
\xi_{N,2}(t;\Psi_{n\ell m})
\left(
\rme^{\frac{\rmi}{\hbar}\frac{t}{N}(H_2-E_n)} -I
\right)
\Psi_{n\ell m}
\right\|
\nonumber\\
\le{}&
\begin{cases}
\medskip
\displaystyle
\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/2}}}
\frac{1}{4n^4}
\right),
\displaystyle
\frac{\tilde{t}^{\,2}}{2N}
\left(
\sqrt{
\frac{3-\ell(\ell+1)/n^2}{2n^3\ell(\ell+1)(\ell+1/2)(\ell+3/2)(\ell-1/2)}
\frac{1}{n^{7/2}\sqrt{\ell+1/2}}
\frac{1}{4n^4}
\right),
\end{cases}
\end{align*}

These yield the bounds (<ref>)–(<ref>) presented in Sec. <ref> of the Methods section.

§.§ Hydrogen ground state

A particularly interesting case is the ground state of the hydrogen atom, $n=1$, $\ell=0$, for which we have

\begin{equation}
\xi_{N,1}(t;\Psi_{100})\le \frac{\tilde{t}^{\,5/4}}{N^{1/4}}
\frac{4}{\sqrt{5\sqrt{\pi}}},
\label{eqErrorBoundhydrogenGround1}
\end{equation}

\begin{equation}
\xi_{N,2}(t;\Psi_{100})\le
\frac{\tilde{t}^{\,3/2}}{N^{1/2}}
\sqrt{
\frac{4\pi}{3}
\frac{\tilde{t}^{\,2}}{2N}
\left(
\sqrt{2}
\frac{1}{4}
\right).
\label{eqErrorBoundhydrogenGround2}
\end{equation}

Notice that the total error is the sum of both. As we can see, the dominating part in the bound is of the order $\mathcal{O}(N^{-1/4})$.

§.§ Remark on the tightness of the bounds obtained from Lemma <ref>

For these hydrogen examples, the bounds (<ref>) and (<ref>) on $\Delta_j(s)$ are asymptotically tight for $s\to0$.
Indeed, the uniform bounds (<ref>) and (<ref>) on the spectral densities $\rho_{j,\Psi_{n\ell m}}(\lambda)$ allow us to apply Lebesgue's dominated convergence theorem to $\Delta_j(s)/s^{2\delta_j}$ as

\begin{align*}
\lim_{s\to0}\frac{1}{s^{2\delta_j}}\Delta_j(s)
\lim_{s\to0}
\int_0^\infty f(x)
\frac{\rho_{j,\varphi}(x/s)}{s^{2\delta_j+1}}
\,\rmd x
\nonumber\\
\int_0^\infty f(x)
\lim_{s\to0}
\frac{\rho_{j,\varphi}(x/s)}{s^{2\delta_j+1}}
\,\rmd x
\nonumber\\
2\delta_j (\lambda_j)^{2\delta_j}\int_0^\infty
\frac{f(x)}{x^{2\delta_j+1}}
\,\rmd x
\nonumber\\
\frac{\pi(\lambda_j)^{2\delta_j}}{\Gamma(2\delta_j-1)\sin[(\delta_j-1)\pi]},
\end{align*}

where $\lambda_j$ and $\delta_j$ ($j=1,2$) are introduced in Eqs. (<ref>) and (<ref>).
This shows that $\Delta_j(s)$ asymptotically approach the bounds.
Nevertheless, this does not necessarily lead to a tight bound on the convergence speed of the corresponding Trotter formula.

§ HIGHER-ORDER TROTTER BOUNDS AND PROOFS OF THM. <REF>, THM. <REF>, AND THM. <REF>

In this appendix, we show how to obtain state-dependent bounds for arbitrary higher-order Trotter product formulas.
We develop a general formalism, which we then apply to the first- and second-order Trotterizations (Thm. <ref> and Thm. <ref>).
Furthermore, we show how to obtain a loose state-dependent error bound for a general $p$th-order product formula (Thm. <ref>).
This follows the proofs of the error bounds with fractional scalings (Thm. <ref> and Thm. <ref>).
These proofs are based on the formalism presented before, and our proof method generalizes iteratively to the $p$th-order product formulas with errors scaling between $N^{-q+1}$ and $N^{-q}$, $q\in(1,p)$.
We explicitly carry out this computation for the cases $q=1$ and $q=2$.

§.§ General formalism

A general $p$th-order Trotter product, which aims to approximate the evolution $\rme^{-\rmi t(H_1+H_2)}$, 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,
\end{equation*}

with $\sum_{i=1}^{\lceil M/2\rceil}\tau_{2i-1}=\sum_{i=1}^{\lfloor M/2\rfloor}\tau_{2i}=1$ and $\sigma(j)$ defined in Eq. (<ref>).
We wish to estimate the error of this approximation on an input eigenstate $\varphi$ of $H_1+H_2$.
Without loss of generality, we shift the Hamiltonians $H_1$ and $H_2$ so that we have $(H_1+H_2)\varphi=0$.
Then, the Trotter product $\mathcal{S}_N^{(p)}(t)$ is to be compared with the identity $I$ on the input eigenstate $\varphi$.
As in the proofs of Thm. <ref> and Thm. <ref>, it suffices to study 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:1_trotter_cycle_appendix}
\end{equation}

instead of the full Trotter evolution.
This is because we can use a telescoping sum to obtain a bound on $\mathcal{S}_N^{(p)}(t)$ from a bound on $\mathcal{S}_1^{(p)}(t/N)$, that is

\begin{equation*}
\mathcal{S}_N^{(p)}(t) - I =
[\mathcal{S}_1^{(p)}(t/N)]^N - I = \sum_{k=0}^{N-1} [\mathcal{S}_1^{(p)}(t/N)]^k [\mathcal{S}_1^{(p)}(t/N) - I] ,
\end{equation*}

\begin{equation}
\xi_N^{(p)}(t)=\Vert[\mathcal{S}_N^{(p)}(t)-I] \varphi\Vert \le N\Vert[\mathcal{S}_1^{(p)}(t/N)-I] \varphi\Vert
\label{eq:telescope_higher_order}
\end{equation}

follows from the triangle inequality and the unitary equivalence of the Euclidean norm.

The idea of our approach is to treat the Trotterized evolution $\mathcal{S}_1^{(p)}(t/N)$ as a unitary evolution

\begin{equation*}
-\rmi\int_0^s\rmd u\,H^{(p)}(u)
\right)
\end{equation*}

generated by a piecewise-constant time-dependent Hamiltonian $H^{(p)}(s)$ [69].
For instance, in the case of the first-order Trotter product, we regard it as a unitary evolution generated by a time-dependent Hamiltonian $H^{(1)}(s)$ in which the Hamiltonian is switched as [40]

\begin{equation*}
\medskip
\displaystyle
\displaystyle
\end{cases}
\end{equation*}

In the case of the second-order Trotter product,

\begin{equation*}
\medskip
\displaystyle
\frac{1}{2} H_1,&s\in\left[0,\frac{t}{N}\right),\\
\medskip
\displaystyle
\displaystyle
\frac{1}{2} H_1,&s\in\left[\frac{2t}{N},\frac{3t}{N}\right).
\end{cases}
\end{equation*}

In the case of a general $p$th-order Trotter product, some of the $\tau_j$ can be negative.
In fact, Suzuki showed that every $p$th-order product formula of order $p\geq3$ has to have at least one negative $\tau_j$ <cit.>.
If $\tau_j$ is negative, the Hamiltonian $H^{(p)}(s)$ in the $j$th time slot is
$\tau_j H_{\sigma(j)}= -|\tau_j| H_{\sigma(j)}$ for time duration
$t/N$, and we have

\begin{equation}
\medskip
\displaystyle
\tau_1 H_{\sigma(1)},&s\in\left[0,\frac{t}{N}\right),\\
\medskip
\displaystyle
\cdots&\cdots\\
\medskip
\displaystyle
\tau_j H_{\sigma(j)},&s\in\left[\frac{(j-1)t}{N},\frac{jt}{N}\right),\\
\medskip
\displaystyle
\cdots&\cdots\\
\displaystyle
\tau_M H_{\sigma(M)},&s\in\left[\frac{(M-1)t}{N},\frac{Mt}{N}\right).
\end{cases}
\label{eq:piecesise-constant_Hamiltonian}
\end{equation}

Then, the evolution $U^{(p)}(s)$ generated by this piecewise-constant time-dependent Hamiltonian $H^{(p)}(s)$ yields
$U^{(p)}(M t/N)=\mathcal{S}_1^{(p)}(t/N)$ at time $s=T_M=M t/N$.
The Hamiltonian $H^{(p)}(s)$ is extended periodically from one Trotter cycle to the full Trotter product by $H^{(p)}(s+T_M)=H^{(p)}(s)$.

To bound the state-dependent Trotter error $\xi_1^{(p)}(t/N)$, we use the following Dyson expansion up to the $q$th order.

Let  $t \mapsto H(t)$ be a piecewise-constant Hamiltonian in $[0,s]$.
That is, there exists a partition $0=t_0<t_1<\dots<t_L=s$
of $[0,s]$, such that $H(t)=h_\ell$ for $t\in [t_{\ell-1},t_\ell)$, for $\ell=1,\ldots,L$ [and $H(s)=h_L$], with $h_\ell=h_\ell^\dag$. Let

\begin{equation*}
U(u)=\T\exp\!\left(-\rmi\int_0^u\rmd s_1\,H(s_1)\right),
\end{equation*}

for $u\in[0,s]$, be the unitary propagator generated by $H(t)$.

Let $q$ be a positive integer, and $\varphi\in\mathcal{D}(h_{\ell_1}h_{\ell_2} \cdots  h_{\ell_q} h_{\ell_{q+1}})$ for all $L\geq \ell_1\geq \ell_2 \geq \cdots \geq \ell_{q+1} \geq 1$.
Define the $k$th-order integral action on $\varphi$,

\begin{equation}
S_k(u) \varphi =\int_0^u\rmd s_1\cdots\int_0^{s_{k-1}}\rmd s_k\, H(s_1)\cdots H(s_k) \varphi,\quad k=1,\ldots,q,
\label{eq:action_appendix}
\end{equation}

for all $u\in[0,s]$.
Then, one has

\begin{align}
U(s)\varphi=\varphi+\sum_{k=1}^q(-\rmi)^kS_k(s)\varphi+(-\rmi)^{q+1}\int_0^s\rmd u\,U(s)U^\dag(u)H(u)S_q(u)\varphi.
\label{eq:lemma_error_operator}
\end{align}

For a constant Hamiltonian $H$, this lemma provides a remainder term of the Taylor expansion of the evolution operator, for $\varphi\in \mathcal{D}(H^{q+1})$:

\begin{equation}
\rme^{-\rmi sH}\varphi
=\varphi+\sum_{k=1}^q\frac{1}{k!}(-\rmi sH)^k\varphi
+(-\rmi)^{q+1}\int_0^s\rmd u\,\rme^{-\rmi(s-u)H}\frac{1}{q!}u^qH^{q+1}\varphi.
\label{eq:TaylorExpansion}
\end{equation}

It is already used in Eq. (<ref>) and will be used later.

We give two different derivations of the Dyson expansion (<ref>).
The first is a generalization of the integration-by-part strategy developed in Ref. [69].

Notice first that, by assumption, $H(s_1)\cdots H(s_k) \varphi$ is well defined and piecewise continuous in all variables for all $k=1,\ldots,q$, whence the integral actions $S_k(u)\varphi$ are piecewise differentiable. Moreover, $\varphi \in \mathcal{D}\bm{(}H(u) S_k(u)\bm{)}$ for all $u\in[0,s]$.
The difference between $U(s)$ and the identity $I$ on $\varphi$ can be estimated by [69]

\begin{equation*}
\biggr]_{u=0}^{u=s}
\varphi
=-\int_0^s\rmd u\,
\frac{\partial}{\partial u}[
\varphi
=-\rmi\int_0^s\rmd u\,
\end{equation*}

where we have used that $U^\dag(u)\varphi$ is piecewise continuosly differentiable with

\begin{equation*}
\frac{\rmd}{\rmd u}U^\dag(u)\varphi=\rmi U^\dag(u)H(u)\varphi.
\end{equation*}

Since $H(u)\varphi =\dot{S}_1(u)\varphi$ except at $u=t_k$, by integrating by parts one gets

\begin{align*}
&=-\rmi\int_0^s\rmd u\,
\nonumber\\
&=-\rmi S_1(s)\varphi
+(-\rmi)^2\int_0^s\rmd u\,
\nonumber
\intertext{By iterating the integration by parts $(q-1)$ more times, we get}
-\rmi S_1(s)\varphi
+(-\rmi)^2\int_0^s\rmd u\,
\nonumber\\
-\rmi S_1(s)\varphi
+(-\rmi)^3\int_0^s\rmd u\,
\nonumber\\
\nonumber\\
-\rmi S_1(s)\varphi
+(-\rmi)^{q+1}\int_0^s\rmd u\,
\nonumber
\end{align*}

This proves the lemma.

An alternate derivation of Lemma <ref> is to solve the backward Schrödinger equation

\begin{equation}
\frac{\partial}{\partial u}U(s,u) \varphi =\rmi U(s,u)H(u)\varphi,\qquad
\label{eq:BackwardSchrodingerEq}
\end{equation}

for $U(s,u)=U(s)U^\dag(u)$, iteratively.

The derivative of the unitary propagator on $\varphi$ in Eq. (<ref>) is piecewise continuous. By integrating the Schrödinger equation (<ref>), one has

\begin{equation*}
U(s,u)\varphi=\varphi-\rmi\int_u^s\rmd u_1\,U(s,u_1)H(u_1)\varphi.
\end{equation*}

By iteration,

\begin{align*}
-\rmi\int_u^s\rmd u_1\,H(u_1)\varphi
+(-\rmi)^2\int_u^s\rmd u_1\int_{u_1}^s\rmd u_2\,U(s,u_2)H(u_2)H(u_1)\varphi\\
(-\rmi)^k\int_u^s\rmd u_1\cdots\int_{u_{k-1}}^s\rmd u_k\,H(u_k)\cdots H(u_1)\varphi
\nonumber\\
{}+(-\rmi)^{q+1}\int_u^s\rmd u_1\cdots\int_{u_q}^s\rmd u_{q+1}\,U(s,u_{q+1})H(u_{q+1})\cdots H(u_1)\varphi.
\intertext{By reverting the order of the integrations,}
(-\rmi)^k\int_u^s\rmd u_k\int_u^{u_k}\rmd u_{k-1}\cdots\int_u^{u_2}\rmd u_1\,H(u_k)\cdots H(u_1)\varphi
\nonumber\\
{}+(-\rmi)^{q+1}\int_u^s\rmd u_{q+1}\int_u^{u_{q+1}}\rmd u_q\cdots\int_u^{u_2}\rmd u_1\,U(s,u_{q+1})H(u_{q+1})\cdots H(u_1)\varphi.
\end{align*}

By setting $u=0$, we get Eq. (<ref>).

Notice that the use of the backward Schrödinger equation is essential in our argument. As a matter of fact, even the forward Schrödinger equation itself

\begin{align*}
\frac{\rmd}{\rmd s}U(s)\varphi=-\rmi H(s)U(s)\varphi
\end{align*}

does not hold, unless one requires additional strong conditions on the operators $h_k$: that they have a common domain $D$, which is left invariant under their own unitary groups, and that $\varphi\in D$.
Moreover, in order to iterate it one needs even stronger assumptions. In any case, one would finally obtain by iteration

\begin{align*}
\varphi-\rmi\int_0^s\rmd s_1\,H(s_1)U(s_1)\varphi
\\
\varphi
+(-\rmi)^{q+1}\int_0^s\rmd s_1\cdots\int_0^{s_q}\rmd s_{q+1}\,H(s_1)\cdots H(s_{q+1})U(s_{q+1})\varphi.
\end{align*}

This leads to a different expression for the remainder term, which is not convenient for our purpose.
Indeed, the presence of the unitary $U(s_{q+1})$ on $\varphi$ before the action of the Hamiltonians makes it  difficult to bound the remainder term.
On the other hand, if the unitary comes last in the remainder term as in the expansion (<ref>), it is harmless in bounding the remainder term, as we will see below.

We use the Dyson expansion of Lemma <ref> for the Trotter evolution $U^{(p)}(s)$ generated by the piecewise-constant Hamiltonian $H^{(p)}(s)$ in Eq. (<ref>), and estimate the Trotter error $\xi_1^{(p)}(t/N)=\|[U^{(p)}(T_M)-I]\varphi\|$ for the single Trotter cycle.

Let $U^{(p)}(s)$ be the unitary propagator generated by the piecewise-constant Hamiltonian $H^{(p)}(s)$ in Eq. (<ref>) representing a $p$th-order Trotterized evolution with $M$ time slots in a single Trotter cycle.
The unitary reaches
$U^{(p)}(Mt/N)=\mathcal{S}_1^{(p)}(t/N)$ in Eq. (<ref>) at time $s=T_M=Mt/N$.

Let $q$ be a positive integer $q\le p$. Let the input state $\varphi$ satisfy $(H_1+H_2)\varphi=0$, and
$\varphi\in\mathcal{D}\bm{(}H^{(p)}(s_1)H^{(p)}(s_2)\cdots H^{(p)}(s_{q+1})\bm{)}$ for all $T_M> s_1 \geq s_2 \geq \cdots \geq s_{q+1}\geq 0$.
Then, one has

\begin{equation}
=R^{(p)}_q \varphi
=(-\rmi)^{q+1}\int_0^{T_M}\rmd s\,U^{(p)}(T_M)U^{(p)\dag}(s)H^{(p)}(s)S_q(s)\varphi,
\label{eqn:Dyson_error_q}
\end{equation}

where $S_q(s)\varphi$ is the $q$th-order integral action of the piecewise-constant Hamiltonian $H^{(p)}(s)$
defined in Eq. (<ref>).

If $q=p$, the error $\xi_N^{(p)}(t)$ of the overall $p$th-order Trotter product $\mathcal{S}_N^{(p)}(t)$ on the input state $\varphi$ is bounded by

\begin{equation}
\xi_N^{(p)}(t)
\le
\le
N\sum_{j=1}^M \int_{(j-1)t/N}^{jt/N}\rmd s\,\Vert \tau_j H_{\sigma(j)}S_p(s)\varphi\Vert,
\label{eqn:error_estimate_p}
\end{equation}

which is $\mathcal{O}(t^{p+1}/N^p)$.

Set $T_j=jt/N$. Notice first that the integral action $S_k(T_M)\varphi$ in Eq. (<ref>) for the piecewise-constant Hamiltonian $H^{(p)}(s)$ in Eq. (<ref>) is proportional to $t^k/N^k$.
Since a $p$th-order Trotter product $\mathcal{S}_1^{(p)}(t/N)=U^{(p)}(T_M)$ aims to approximate $\rme^{-\rmi\frac{t}{N}(H_1+H_2)}$ at $\mathcal{O}(t^{p+1}/N^{p+1})$, it formally requires that

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

We hence have

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

and the Dyson expansion of Lemma <ref> with $q\le p$ yields Eq. (<ref>).
If $q=p$, then $\varphi\in\mathcal{D}\bm{(}H^{(p)}(s)S_p(s)\bm{)}$ for all $s\in[0,T_M)$, and the single-cycle Trotter error $\xi_1^{(p)}(t/N)$ is bounded by

\begin{equation*}
\xi_1^{(p)}(t/N)
\|[U^{(p)}(T_M)-I]\varphi\|
\le\int_0^{T_M}\rmd s\,\|H^{(p)}(s)S_p(s)\varphi\|
\sum_{j=1}^M \int_{T_{j-1}}^{T_j}\rmd s\, \Vert \tau_j H_{\sigma(j)}S_p(s)\varphi\Vert.
\end{equation*}

Thus, Eq. (<ref>) gives the bound (<ref>).

In the case of unbounded operators,
Eq. (<ref>) might only be well defined for some $q<p$.
In such a case, we have to perform a more sophisticated estimation of the error, instead of simply bounding Eq. (<ref>) by triangle inequality as in Eq. (<ref>).
See the following subsections of this appendix.

In any case, to proceed, we need to explicitly compute the integral actions $S_k(s)$ defined in Eq. (<ref>) for the piecewise-constant Hamiltonian $H^{(p)}(s)$ in Eq. (<ref>) for a $p$th-order Trotter product. We proceed formally, by temporarily assuming that all operators are bounded. For unbounded operators the expressions obtained will be valid when applied to a vector belonging to a suitable domain.

For $s\in[(j-1)t/N,jt/N)$ in the $j$th time slot, we have

\begin{align*}
=\int_0^s\rmd s_1\cdots\int_0^{s_{k-1}}\rmd s_k\,
H^{(p)}(s_1)\cdots H^{(p)}(s_k)
\nonumber\\
\sum_{j\geq j_1 \geq \dots \geq j_k \geq 1 }
\int_{(j_1-1)t/N}^{\tilde{T}_{j_1}}\rmd s_1\cdots\int_{(j_k-1)t/N}^{\tilde{T}_{j_k}}\rmd s_k\,
\tau_{j_1}H_{\sigma(j_1)}\cdots
\tau_{j_k} H_{\sigma(j_k)},\quad
s\in\left[\frac{(j-1)t}{N},\frac{j t}{N}\right),
\intertext{where $\tilde{T}_{j_\ell}=j_\ell t/N$ if $j_{\ell-1}>j_{\ell}$,
or $\tilde{T}_{j_\ell}=s_{\ell-1}$ if $j_{\ell-1}=j_\ell$,
for $\ell=2,\ldots,k$, and $\tilde{T}_{j_1}=j_1 t/N$ if $j>j_1$, or $\tilde{T}_{j_1}=s$ if $j_1=j$.
By shifting the time integration in each slot, we get}
\sum_{j\geq j_1 \geq \dots \geq j_k \geq 1 }
\int_0^{s_{j_1}}\rmd s_1\cdots\int_0^{s_{j_k}}\rmd s_k\,
\tau_{j_1} H_{\sigma(j_1)}\cdots
\tau_{j_k} H_{\sigma(j_k)},\quad
s\in\left[\frac{(j-1)t}{N},\frac{j t}{N}\right),
\end{align*}

\begin{equation*}
\begin{cases}
\displaystyle
\medskip
s-\frac{(j-1)t}{N}, &\mathrm{if}\ \ell=1\ \mathrm{and}\ j_1=j,
\\
\displaystyle
\medskip
s_{\ell-1},&\mathrm{if}\ \ell\in\{2,\ldots,k\}\ \mathrm{and}\ j_\ell=j_{\ell-1},
\\
\displaystyle
\frac{t}{N}, &\mathrm{otherwise}.
\end{cases}
\end{equation*}

This yields

\begin{equation*}
= \left(\frac{t}{N}\right)^k
\sum_{j\geq j_1 \geq \dots \geq j_k \geq 1 } \frac{1}{m(j_1,\dots,j_k)}
(\tilde{\tau}_{j_k} H_{\sigma(j_k)}),\quad
s\in\left[\frac{(j-1)t}{N},\frac{j t}{N}\right),
\end{equation*}

where $\tilde{\tau}_{i}=\tau_i$ if $i< j$ and $\tilde{\tau}_{j}=\tau_j [s-(j-1)t/N]/{(t/N)}$, and $m(j_1,\dots,j_k)$ is a multiplicity factor which for each group of equal indices $j_\ell=j_{\ell+1}=\dots=j_{\ell+\beta-1}$ counts the factorial of its cardinality $\beta!$. By gathering together each group of equal indices, we finally get

\begin{equation}
S_k(s)= \left(
\frac{t}{N}
\right)^k
\mathop{\sum_{\bm{\beta}\in\mathbb{N}^j}}_{|\bm{\beta}|=k}
\frac{1}{\bm{\beta}!}
\left(\frac{s-(j-1)t/N}{t/N}\right)^{\beta_j}
\cdots
\quad
s\in\left[\frac{(j-1)t}{N},\frac{j t}{N}\right),
\label{eq:kth_action_explicit_expanded0}
\end{equation}

with $\bm{\beta}=(\beta_1,\dots,\beta_j)$ being a multi-index of nonnegative integers, $|\bm{\beta}|=\beta_1+\cdots+\beta_j$, and $\bm{\beta}!=\beta_1!\cdots\beta_j!$. This is Eq. (<ref>). For $s=T_M$, we obtain

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

which is Eq. (<ref>).

The usual strategy in the literature to obtain higher-order product formulas is by demanding that all terms up to order $p$ in the Taylor expansion of $U^{(p)}(T_M)$ vanish, see e.g. Refs. [67, 71].
That is, one requires

\begin{equation}
U^{(p)}(T_M) = I + \sum_{k=1}^p \frac{(-\rmi)^k}{k!}\frac{t^k}{N^k}(H_1+H_2)^k +\mathcal{R}_{p+1}(T_M),
\label{eq:Taylor_expansion}
\end{equation}

where the remainder $\mathcal{R}_{p+1}(T_M)$ is of order $\mathcal{O}(t^{p+1}/N^{p+1})$.
This requirement gives a system of equations for the switching coefficients $\tau_k$ to give rise to a $p$th-order product formula.
However, these equations are equivalent to the set of equations obtained by our approach through Eq. (<ref>).
This can be shown very easily.
It is obvious that the requirement in Eq. (<ref>) gives a $p$th-order product formula.
In fact, we showed that the overall Trotter error is of order $\mathcal{O}(t^{p+1}/N^p)$ under this assumption.
Thus, it remains to show that Eq. (<ref>) implies Eq. (<ref>).
Combining Eq. (<ref>) with Lemma <ref> gives

\begin{equation*}
\sum_{k=1}^p(-\rmi)^kS_k(T_M)+(-\rmi)^{p+1}\int_0^{T_M}\rmd u\,U(s)U^\dag(u)H(u)S_p(u)
\sum_{k=1}^p \frac{(-\rmi)^k}{k!}\frac{t^k}{N^k}(H_1+H_2)^k +\mathcal{R}_{p+1}(T_M).
\end{equation*}

However, on the left-hand side, each $S_k(T_M)=\mathcal{O}(t^{k}/N^k)$.
At the same time, each term $\frac{(-\rmi)^k}{k!}\frac{t^k}{N^k}(H_1+H_2)^k=\mathcal{O}(t^{k}/N^k)$ on the right-hand side.
Therefore, these terms have to agree with each other order by order and we retrieve Eq. (<ref>).
Nevertheless, as opposed to a Taylor expansion, our approach allows us to explicitly compute state-dependent error bounds.
This also allows us to compare how well different product formulas of the same order perform.

In Appendix <ref>, we provide a Mathematica script that automatically computes the bound in Eq. (<ref>), given the switching times $\tau_j$.
The fourth-order Trotter bound presented in Thm. <ref> is computed by this script.

§.§ First-order Trotter bound (Thm. <ref>)

Let us see how the general formalism presented above works for the first-order Trotterization.
Without loss of generality, we shift the Hamiltonians $H_1$ and $H_2$ so that $(H_1+H_2)\varphi=0$.
The switching times for the first-order Trotterization are given by $\tau_1=1$ and $\tau_2=1$, and hence, $T_1=t/N$ and $T_2=2t/N$.
In this case, the formula (<ref>) for the integral action gives

\begin{equation}
S_1(s) \varphi
\medskip
\displaystyle
sH_1 \varphi,&s\in[0,\frac{t}{N}),\\
\displaystyle
\left[\left(s-\frac{t}{N}\right)H_2
+\frac{t}{N}H_1 \right]\varphi ,&
\end{cases}
\label{eq:S1_first_order}
\end{equation}

It is easy to verify that we have

\begin{equation*}
S_1(2t/N) \varphi=\frac{t}{N}(H_1+H_2) \varphi =0.
\end{equation*}

assuming that $\varphi \in\mathcal{D}(H_1^2)\cap \mathcal{D}(H_2^2)$,
the error of the first-order Trotterization is bounded by Eq. (<ref>),

\begin{equation*}
\xi_N^{(1)}(t;\varphi)\leq N (I_1+I_2),
\end{equation*}

with $I_j$ ($j=1,2$) the integral over the $j$th time slot.
For the first time slot $s\in[0,\frac{t}{N})$, we have

\begin{equation*}
I_1=\int_0^{t/N}\rmd s\,s\Vert H_1^2\varphi\Vert =\frac{t^2}{2N^2}\Vert H_1^2\varphi\Vert.
\end{equation*}

For the second time slot $s\in[\frac{t}{N},\frac{2t}{N})$,

\begin{equation*}
= \int_{t/N}^{2t/N}\rmd s\left\Vert\left[
\left(s-\frac{t}{N}\right)H_2^2
\right]\varphi\right\Vert
\end{equation*}

where the second step follows by the fact that $H_2H_1\varphi = - H_2^2\varphi$, which also shows that for a zero-energy eigenstate of $H_1+H_2$, the domain condition $\varphi \in\mathcal{D}(H_2 H_1)$ is equivalent to $\varphi\in \mathcal{D}(H_2^2)$.
Summing up these two terms reproduces the bound (<ref>) of Lemma <ref>.

§.§ Second-order Trotter bound (Thm. <ref>)

The error bound on the second-order Trotterization in Thm. <ref> is obtained as an application of Lemma <ref>.

We shift the Hamiltonians $H_1$ and $H_2$ so that $(H_1+H_2)\varphi=0$, and apply the general formalism for bounding the error of a higher-order Trotterization described above.
The switching times for the second-order Trotterization are given by $\tau_1=1/2$, $\tau_2=1$, and $\tau_3=1/2$.
In this case, the formula (<ref>) for the integral action gives

\begin{gather}
\medskip
\displaystyle
\frac{1}{2} s H_1,&s\in[0,\frac{t}{N}),\\
\medskip
\displaystyle
\left(s-\frac{t}{N}\right)H_2+\frac{t}{2N}H_1,&
\displaystyle
\frac{1}{2}\left(s-\frac{t}{N}\right)H_1
\end{cases}
\label{eq:S1_second_order}
\\
\medskip
\displaystyle
\frac{1}{8}s^2H_1^2,&s\in[0,\frac{t}{N}),\\
\medskip
\displaystyle
\frac{1}{2}\left(s-\frac{t}{N}\right)^2H_2^2+\left(s-\frac{t}{N}\right)\frac{t}{2N}H_2H_1+\frac{t^2}{8N^2}H_1^2,&
\displaystyle
\frac{1}{8}\left(s-\frac{t}{N}\right)^2H_1^2
\end{cases}
\label{eq:S2_second_order}
\end{gather}

It is easy to verify that we have

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

and thus $S_1(3t/N)\varphi=0$ and $S_2(3t/N)\varphi=0$ for the eigenstate satisfying $(H_1+H_2)\varphi=0$.

Therefore, by assuming that

\begin{equation*}
\varphi\in\mathcal{D}(H_1^3)\cap\mathcal{D}(H_2 H_1^2)\cap\mathcal{D}(H_2^3),
\end{equation*}

the error of the second-order Trotterization is bounded by Eq. (<ref>),

\begin{equation*}
\xi_N^{(2)}(t;\varphi)\leq N (I_1+I_2+I_3),
\end{equation*}

with $I_j$ ($j=1,2,3$) the integral over the $j$th time slot.
For the first time slot $s\in[0,\frac{t}{N})$, we have

\begin{equation*}
I_1=\int_0^{t/N}\rmd s\,\frac{1}{16}s^2\Vert H_1^3\varphi\Vert  =\frac{t^3}{48N^3} \Vert H_1^3\varphi\Vert.
\end{equation*}

For the second time slot $s\in[\frac{t}{N},\frac{2t}{N})$,

\begin{align*}
I_2&= \int_{t/N}^{2t/N}\rmd s\left\Vert\left[
\frac{1}{2}\left(s-\frac{t}{N}\right)^2 H_2^3
\right]\varphi\right\Vert\\
&\leq \frac{t^3}{12N^3}\Vert H_2^3\varphi\Vert + \frac{t^3}{8N^3}\Vert H_2H_1^2\varphi\Vert,
\end{align*}

where the second step follows by the fact that $H_2^2H_1\varphi = - H_2^3\varphi$ [so that $\varphi\in\mathcal{D}(H_2^3)$ is equivalent to $\varphi\in\mathcal{D}(H_2^2 H_1)$] and the triangle inequality.
For the last time slot $s\in[\frac{2t}{N}, \frac{3t}{N})$, we get

\begin{align*}
I_3 &= \int_{2t/N}^{3t/N}\rmd s \left\Vert\left[
\frac{1}{16}\left(s-\frac{t}{N}\right)^2H_1^3
\right] \varphi\right\Vert\\
&=\frac{t^3}{48N^3}\Vert H_1^3\varphi\Vert,
\end{align*}

where we have used $H_1^2H_2\varphi=-H_1^3\varphi$, so that $\varphi\in\mathcal{D}(H_1^3)$ is equivalent to $\varphi\in\mathcal{D}(H_1^2 H_2)$, and $\varphi\in\mathcal{D}\bm{(}H_1H_2 (H_1+H_2)\bm{)}$ with $H_1H_2(H_1+H_2)\varphi=0$.
Summing up these three terms results in the second-order Trotter bound presented in Thm. <ref>.

§.§ A general bound on the $p$th-order Trotter error (Thm. <ref>)

By bounding each term appearing in Eq. (<ref>), we get Thm. <ref>.

We shift the Hamiltonians $H_1$ and $H_2$ so that $(H_1+H_2)\varphi=0$, and apply the general formalism for bonding the error of a higher-order Trotterization described above.
Bounding Eq. (<ref>) with Eq. (<ref>) by triangle inequality, we have

\begin{align*}
\xi_N^{(p)}(t)
\le{}& N\sum_{j=1}^M
\mathop{\sum_{\bm{\beta}\in\mathbb{N}^j}}_{|\bm{\beta}|=p} \frac{1}{\bm{\beta}!}
\left(\frac{t}{N}\right)^{p-\beta_j}
\int_0^{\frac{t}{N}}\rmd s\,
s^{\beta_j} |\tau_{j}|^{\beta_{j}+1}
\|H_{\sigma(j)}^{\beta_j+1}
\cdots
\nonumber\\
\mathop{\sum_{\bm{\beta}\in\mathbb{N}^j}}_{|\bm{\beta}|=p} \frac{1}{\bm{\beta}!}
\left(\frac{t}{N}\right)^{p+1}
\frac{|\tau_{j}|^{\beta_{j}+1}|\tau_{j-1}|^{\beta_{j-1}}\cdots|\tau_1|^{\beta_1}}{\beta_j+1}
\|H_{\sigma(j)}^{\beta_j+1}
\cdots
\nonumber\\
\mathop{\sum_{\bm{\beta}\in\mathbb{N}^M}}_{|\bm{\beta}|=p+1} \frac{1}{\bm{\beta}!}
\left(\frac{t}{N}\right)^{p+1}
\|H_{\sigma(M)}^{\beta_M}
\cdots
\nonumber\\
\le{}&N
\left(
\frac{t}{N}
\right)^{p+1}
\frac{1}{(p+1)!}
\mathop{\sum_{\bm{\beta}\in\mathbb{N}^M}}_{|\bm{\beta}|=p+1} \frac{|\bm{\beta}|!}{\bm{\beta}!} {|\tau_1|^{\beta_1}\cdots|\tau_M|^{\beta_M}}
\nonumber\\
\left(
\frac{t}{N}
\right)^{p+1}
\frac{(|\tau_1|+\cdots+|\tau_M|)^{p+1}}{(p+1)!}
\nonumber\\
\left(
\frac{\tau_*t}{N}
\right)^{p+1}
\frac{1}{(p+1)!}
\end{align*}

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)}\cdots H_{\sigma(j_1)}\varphi\|$.
The second last equality is the multinomial theorem, with ${|\bm{\beta}|!}/{\bm{\beta}!}$ being the multinomial coefficient.

§.§ Fractional scaling of the first-order Trotter bound (Thm. <ref>)

Theorem <ref> requires $\varphi\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2^2)$ and Thm. <ref> requires $\varphi\in\mathcal{D}(H_1^3)\cap\mathcal{D}(H_2H_1^2)\cap\mathcal{D}(H_2^3)$.
If these conditions are not met, we lose the $N^{-1}$ scaling of the first-order Trotter bound and the $N^{-2}$ scaling of the second-order Trotter bound.
To estimate the scalings in such cases, we perform refined analyses as done in Thm. <ref> for the first-order Trotterization.
The general formalism described above enables us to do it systematically for any general $p$th-order Trotterization.

Let us first look at the first-order Trotterization. In this case $\tau_1=\tau_2=1$,
and we have to compute the right-hand side of Eq. (<ref>), $R^{(p)}_q \varphi$, for $p=q=1$, where

\begin{equation*}
R^{(1)}_1= (-\rmi)^{2}\int_0^{2t/N}\rmd s\,U^{(1)}(2t/N)U^{(1)\dag}(s)H^{(1)}(s)S_1(s).
\end{equation*}

Here, the integral action $S_1(s)$ is given by Eq. (<ref>), and
the first-order Trotter evolution reads

\begin{equation*}
= \rme^{-\rmi\frac{t}{N}H_2}
\begin{cases}
\medskip
\displaystyle
\rme^{-\rmi(\frac{t}{N}-s)H_1},&s\in\left[0,\frac{t}{N}\right),\\
\displaystyle
\rme^{\rmi(s-\frac{t}{N})H_2},&
\end{cases}
\end{equation*}

Thus, one has

\begin{equation*}
{-}\int_0^{t/N}\rmd s\,\rme^{-\rmi\frac{t}{N}H_2}\rme^{-\rmi(\frac{t}{N}-s)H_1}sH_1^2
-\int_{t/N}^{2t/N}\rmd s\,\rme^{-\rmi(\frac{2t}{N}-s)H_2}\left[
\left(s-\frac{t}{N}\right)H_2^2+\frac{t}{N}H_2H_1
\right],
\end{equation*}

which can be integrated and gives

\begin{equation*}
R^{(1)}_1 =
\rme^{-\rmi\frac{t}{N}H_2}\left(
\rme^{-\rmi\frac{t}{N}H_1}-I+\rmi\frac{t}{N}H_1
\right)
\rme^{-\rmi\frac{t}{N}H_2}
\left(
\rme^{\rmi\frac{t}{N}H_2}-I-\rmi\frac{t}{N}H_2
\right)
\end{equation*}

Therefore, for $\varphi\in\mathcal{D}(H_1)\cap\mathcal{D}(H_2)$ with $(H_1+H_2)\varphi=0$,
one finally gets

\begin{equation*}
=R^{(1)}_1\varphi =
\rme^{-\rmi\frac{t}{N}H_2}\left(
\rme^{-\rmi\frac{t}{N}H_1}-I+\rmi\frac{t}{N}H_1
\right)
\varphi
\rme^{-\rmi\frac{t}{N}H_2}
\left(
\rme^{\rmi\frac{t}{N}H_2}-I-\rmi\frac{t}{N}H_2
\right)
\varphi.
\end{equation*}

This reproduces the bound (<ref>) on the first-order Trotter error, which is estimated as done in Appendix <ref>, even when the domain condition $\varphi\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2^2)$ for Thm. <ref> is not fulfilled.

§.§ Fractional scalings of the second-order Trotter bounds (Thm. <ref>)

Let us now estimate the second-order Trotter error when the domain condition $\varphi\in\mathcal{D}(H_1^3)\cap\mathcal{D}(H_2H_1^2)\cap\mathcal{D}(H_2^3)$ for Thm. <ref> is not fulfilled.
For the second-order Trotterization, the integral actions $S_1(s)$ and $S_2(s)$ are given by Eqs. (<ref>) and (<ref>), respectively.
Notice here that, for the second-order Trotter evolution,

\begin{equation*}
\medskip
\displaystyle
\rme^{\rmi \frac{s}{2} H_1},&s\in\left[0,\frac{t}{N}\right),\\
\medskip
\displaystyle
\rme^{\rmi\frac{t}{2N}H_1}
\rme^{\rmi(s-\frac{t}{N})H_2},
\displaystyle
\rme^{\rmi\frac{t}{2N}H_1}\rme^{\rmi\frac{t}{N}H_2}\rme^{\rmi \frac{1}{2}(s-\frac{2t}{N})H_1},&
\end{cases}
\end{equation*}

Then, we have the following lemma, which replaces the bound (<ref>) of Lemma <ref>, and can be used to estimate the second-order Trotter bound for $\varphi\not\in\mathcal{D}(H_1^3)\cap\mathcal{D}(H_2H_1^2)\cap\mathcal{D}(H_2^3)$.

The second-order Trotter bound can be estimated by

\begin{align}
\xi_N^{(2)}(t;\varphi)
\le{}&
\left(
\rme^{-\rmi\frac{t}{2N}H_1}
\right)
\varphi
\right\|
\nonumber\\
\left[
\rme^{-\rmi\frac{t}{2N}H_2}
\left(
\right)
\rme^{\rmi\frac{t}{2N}H_2}
\left(
\right)
\right]
\varphi
\right\|
\frac{t^2}{8N}
\left\|
\left(
\rme^{-\rmi\frac{t}{N}H_2}
\right)
\varphi
\right\|,
\label{eq:fractional_second_2}
\intertext{for $\varphi\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2)$, and by}
\xi_N^{(2)}(t;\varphi)
\le{}&
\left(
\rme^{-\rmi\frac{t}{2N}H_1}
\right)
\varphi
\right\|
\left[
\rme^{-\rmi\frac{t}{2N}H_2}
\left(
\right)
\rme^{\rmi\frac{t}{2N}H_2}
\left(
\right)
\right]
\varphi
\right\|,
\label{eq:fractional_second_1}
\end{align}

for $\varphi\in\mathcal{D}(H_1)\cap\mathcal{D}(H_2)$.

The second-order integral action $S_2(s)$ is given by Eq. (<ref>), and the right-hand side of Eq. (<ref>), $R^{(p)}_q \varphi$, for $p=q=2$ reads

\begin{align*}
R^{(2)}_2 ={}&
\rmi\int_0^{3t/N}\rmd s\,U^{(2)}(3t/N)U^{(2)\dag}(s)H^{(2)}(s)S_2(s)\nonumber\\
\rmi\int_0^{t/N}\rmd s\,\rme^{-\rmi\frac{t}{2N}H_1}\rme^{-\rmi\frac{t}{N}H_2}\rme^{-\rmi \frac{1}{2}(\frac{t}{N}-s)H_1}\frac{1}{16}s^2H_1^3
\nonumber\\
+\rmi\int_{t/N}^{2t/N}\rmd s\,\rme^{-\rmi\frac{t}{2N}H_1}\rme^{-\rmi(\frac{2t}{N}-s)H_2}\left[
\frac{1}{2}\left(s-\frac{t}{N}\right)^2H_2^3
\right]
\nonumber\\
+\rmi\int_{2t/N}^{3t/N}\rmd s\,\rme^{-\rmi \frac{1}{2}(\frac{3t}{N}-s)H_1}\left[
\frac{1}{16}\left(s-\frac{t}{N}\right)^2H_1^3
\right].
\end{align*}

By integrating, one gets

\begin{align*}
\rme^{-\rmi\frac{t}{2N}H_1}\rme^{-\rmi\frac{t}{N}H_2}
\left(
\rme^{-\rmi\frac{t}{2N}H_1}
\right)
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_2}
\left[
\rme^{-\rmi\frac{t}{2N}H_2}
\left(
\rmi\frac{t}{2N}
\right)
\rme^{\rmi\frac{t}{2N}H_2}
\left(
\rmi\frac{t}{2N}H_2
\right)
\right]
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_1}
\left(
\rme^{-\rmi\frac{t}{N}H_2}
\right)
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_1}
\left(
\rme^{-\rmi\frac{t}{N}H_2}
\right)
\nonumber\\
\left(
\rme^{\rmi\frac{t}{2N}H_1}
\rmi\frac{t}{2N}H_1
\frac{t^2}{8N^2}H_1^2
\right)
\nonumber\\
\left(
\rme^{-\rmi\frac{t}{2N}H_1}
\rmi\frac{t}{2N}H_1
\right)
\nonumber\\
\frac{t^2}{2N^2}
\left(
\rme^{-\rmi\frac{t}{2N}H_1}
\right)
\end{align*}

Therefore, for $\varphi\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2)$ with $(H_1+H_2)\varphi=0$, one automatically has that $\varphi\in \mathcal{D}\bm{(}X(H_1+H_2)\bm{)}$ for every operator $X$, and

\begin{align}
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_1}\rme^{-\rmi\frac{t}{N}H_2}
\left(
\rme^{-\rmi\frac{t}{2N}H_1}
\right)\varphi
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_2}
\left[
\rme^{-\rmi\frac{t}{2N}H_2}
\left(
\rmi\frac{t}{2N}
\right)
\rme^{\rmi\frac{t}{2N}H_2}
\left(
\rmi\frac{t}{2N}H_2
\right)
\right]\varphi
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_1}
\left(
\rme^{-\rmi\frac{t}{N}H_2}
\right)
\nonumber\\
\left(
\rme^{\rmi\frac{t}{2N}H_1}
\rmi\frac{t}{2N}H_1
\frac{t^2}{8N^2}H_1^2
\right)
\varphi.
\end{align}

Bounding it by the triangle inequality yields Eq. (<ref>).
This works for $\varphi\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2)$.

If the input state $\varphi$ is not in $\mathcal{D}(H_1^2)$, to estimate the second-order Trotter bound we can use the remainder term $R^{(p)}_q\varphi$ in Eq. (<ref>) for $p=2$ and $q=1$.
The first-order integral action $S_1(s)$ is given by Eq. (<ref>), and the remainder term $R^{(p)}_q$ for $p=2$ and $q=1$ reads

\begin{align*}
={}&{-}\int_0^{3t/N}\rmd s\,U^{(2)}(3t/N)U^{(2)\dag}(s)H^{(2)}(s)S_1(s)
\nonumber\\
{-}\int_0^{t/N}\rmd s\,\rme^{-\rmi\frac{t}{2N}H_1}\rme^{-\rmi\frac{t}{N}H_2}\rme^{-\rmi \frac{1}{2}(\frac{t}{N}-s)H_1} \frac{1}{4} sH_1^2
\nonumber\\
-\int_{t/N}^{2t/N}\rmd s\,
\rme^{-\rmi\frac{t}{2N}H_1}\rme^{-\rmi(\frac{2t}{N}-s)H_2}
\left[
\left(s-\frac{t}{N}\right)H_2^2
\right]
\nonumber\\
-\int_{2t/N}^{3t/N}\rmd s\,
\rme^{-\rmi \frac{1}{2}(\frac{3t}{N}-s)H_1}
\left[
\frac{1}{4}\left(s-\frac{t}{N}\right)H_1^2
\right].
\end{align*}

By integrating it, one has

\begin{align*}
\rme^{-\rmi\frac{t}{2N}H_1}\rme^{-\rmi\frac{t}{N}H_2}
\left(
\rme^{-\rmi\frac{t}{2N}H_1}
\right)
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_1}
\rme^{-\rmi\frac{t}{2N}H_2}
\left[
\rme^{-\rmi\frac{t}{2N}H_2}
\left(
\right)
\rme^{\rmi\frac{t}{2N}H_2}
\left(
\right)
\right]
\nonumber\\
\rmi\frac{t}{2N}\rme^{-\rmi\frac{t}{2N}H_1}
\left(
\rme^{-\rmi\frac{t}{N}H_2} -I
\right)
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_1}
\left(
\rme^{\rmi\frac{t}{2N}H_1}
\right)
\nonumber\\
\left(
\rme^{-\rmi\frac{t}{2N}H_1} -I
\right)
\end{align*}

whence, for $\varphi\in\mathcal{D}(H_1)\cap\mathcal{D}(H_2)$ with $(H_1+H_2)\varphi=0$, one gets

\begin{align*}
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_1}\rme^{-\rmi\frac{t}{N}H_2}
\left(
\rme^{-\rmi\frac{t}{2N}H_1}
\right)
\varphi
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_1}
\rme^{-\rmi\frac{t}{2N}H_2}
\left[
\rme^{-\rmi\frac{t}{2N}H_2}
\left(
\right)
\rme^{\rmi\frac{t}{2N}H_2}
\left(
\right)
\right]
\varphi
\nonumber\\
\rme^{-\rmi\frac{t}{2N}H_1}
\left(
\rme^{\rmi\frac{t}{2N}H_1}
\right)
\varphi.
\end{align*}

Bounding it by the triangle inequality yields Eq. (<ref>).

In order to prove Thm. <ref>, we need to bound

\begin{gather*}
\Delta^{(0)}(s)
=\|(\rme^{-\rmi sH}-1)\varphi\|^2
=\int_{\mathbb{R}} f_0(s\lambda)\,\rmd\mu_\varphi(\lambda),
\\
\Delta^{(1)}(s)
=\|(\rme^{-\rmi sH}-1+\rmi sH)\varphi\|^2
=\int_{\mathbb{R}} f_1(s\lambda)\,\rmd\mu_\varphi(\lambda),
\\
\Delta^{(2)}(s)
=\left\|\left(\rme^{-\rmi sH}-1+\rmi sH+\frac{1}{2}s^2H^2\right)\varphi\right\|^2
=\int_{\mathbb{R}} f_2(s\lambda)\,\rmd\mu_\varphi(\lambda),
\\
\tilde{\Delta}^{(2)}(s)
=\|[\rme^{-\rmi sH}(1+\rmi sH)-\rme^{\rmi sH}(1-\rmi sH)]\varphi\|^2
\end{gather*}

\begin{gather*}
=|\rme^{-\rmi x}-1|^2
=2(1-\cos x)
\le x^2,
\vphantom{\frac{1}{4}}
\\
=|\rme^{-\rmi x}-1+\rmi x|^2
=(1-\cos x)^2+(x-\sin x)^2
\le\frac{1}{4}x^4,
\\
\rme^{-\rmi x}-1+\rmi x+\frac{1}{2}x^2
\right|^2
=\left(1-\frac{1}{2}x^2-\cos x\right)^2+(x-\sin x)^2
\le\frac{1}{36}x^6,
\\
\tilde{f}_2(x)
\rme^{-\rmi x}(1+\rmi x)
\rme^{\rmi x}(1-\rmi x)
\sin x
-x\cos x
\le\frac{4}{9}x^6.
\end{gather*}

As done in the proof of Lemma <ref>
in Appendix <ref>, if

\begin{equation*}
\mu_{\psi}(\{|x|\geq \lambda\})
\le\left(\frac{\lambda_0}{\lambda}\right)^{2\delta},
\end{equation*}

for some $\lambda_0>0$, then we have

\begin{align}
\Delta^{(0)}(s)
\begin{cases}
\medskip
\displaystyle
\frac{2\pi\delta}{\Gamma(2\delta+1)\sin(\delta\pi)}(\lambda_0 s)^{2\delta},&
\medskip
\displaystyle
\Lambda^2s^2
\frac{2\lambda_0^2}{\Lambda^2}
g_0(\Lambda s),&
\delta=1,\\
\displaystyle
\|H\psi\|^2s^2,
\vphantom{\frac{1}{4}}&
\delta>1,
\end{cases}
\label{eq:bound_Delta0}
\\
\Delta^{(1)}(s)
\begin{cases}
\medskip
\displaystyle
\frac{\pi}{\Gamma(2\delta-1)\sin[(\delta-1)\pi]}(\lambda_0 s)^{2\delta},&
\medskip
\displaystyle
\frac{1}{4}\Lambda^2\|H\psi\|^2s^4
\frac{\lambda_0^4}{\Lambda^4}
g_1(\Lambda s),&
\delta=2,\\
\displaystyle
\frac{1}{4}\|H^2\psi\|^2s^4,&
\delta>2,
\end{cases}
\label{eq:bound_Delta1}
\\
\Delta^{(2)}(s)
\begin{cases}
\medskip
\displaystyle
\frac{\pi}{2\Gamma(2\delta-2)\sin(\delta\pi)]}
(\lambda_0 s)^{2\delta},&
\medskip
\displaystyle
\frac{1}{36}\Lambda^2\|H^2\psi\|^2s^6
\frac{\lambda_0^6}{6\Lambda^6}
g_2(\Lambda s),&
\delta=3,\\
\displaystyle
\frac{1}{36}\|H^3\psi\|^2s^6,&
\delta>3,
\end{cases}
\\
\tilde{\Delta}_2(s)
\begin{cases}
\medskip
\displaystyle
\frac{2^{2\delta-1}(\delta-2)\pi}{\Gamma(2\delta-1)\sin[(\delta-2)\pi]} (\lambda_0 s)^{2\delta},&
\medskip
\displaystyle
\frac{4}{9}\Lambda^2\|H^2\psi\|^2s^6
\frac{8\lambda_0^6}{3\Lambda^6}
\tilde{g}_2(\Lambda s),&
\delta=3,\\
\displaystyle
\frac{4}{9}\|H^3\psi\|^2s^6,&
\delta>3,
\end{cases}
\end{align}

\begin{align*}
&=-x^2\Ci(x)+(1-\cos x)+x\sin x
\vphantom{\frac{1}{4}}
\\
&=-x^2\log x+\left(\frac{3}{2}-\gamma\right)x^2+O(x^4),
\\
&=-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\\
&=-x^4\log x+\left(\frac{7}{4}-\gamma\right)x^4+O(x^6),
\\
-(12-6x^2+x^4)\cos x
-x(12+2x^2-x^4)\sin x
\vphantom{\frac{1}{4}}
\\
&=-x^6\log x+(2-\gamma)x^6+O(x^8),
\vphantom{\frac{1}{4}}
\\
\tilde{g}_2(x)
\frac{1}{8}[
\vphantom{\frac{1}{4}}
\\
\end{align*}

Theorem <ref> follows by using these bounds for Eqs. (<ref>) and (<ref>) of Lemma <ref>.
We summarize the scalings of the first- and second-order Trotter bounds for the hydrogen atom in Table <ref>.

Scalings of the state-dependent bounds on the first- and second-order Trotter errors for the hydrogen atom with $H_1=-\frac{\hbar^2}{2m_\mathrm{e}}\Delta$ and $H_2=-\frac{\hbar^2}{m_\mathrm{e}a_0r}$.
The state $\Psi_{n\ell m}$ is an eigenfunction of the Hamiltonian $H_1+H_2$ of the hydrogen atom.
The parameter $\delta_j$ specifies the scaling of the tail of the spectral density $\rho_{j,\Psi_{n\ell m}}(\lambda)\sim|\lambda|^{-2\delta_j-1}$ of the Hamiltonian $H_j$ ($j=1,2$) at $\Psi_{n\ell m}$.
Similarly, $\delta_{12}$ specifies the scaling of the tail of the spectral density $\rho_{2,H_1^2\Psi_{n\ell m}}(\lambda)\sim|\lambda|^{-2\delta_{12}-1}$ of the Hamiltonian $H_2$ at the vector $H_1^2\Psi_{n\ell m}$, while $\delta_{21}$ specifies the scaling of the tail of the spectral density $\rho_{1,H_2^2\Psi_{n\ell m}}(\lambda)\sim|\lambda|^{-2\delta_{21}-1}$ of the Hamiltonian $H_1$ at the vector $H_2^2\Psi_{n\ell m}$.
“$H_2H_1$” indicates the scaling of the bound on the first-order Trotter product with cycle $\rme^{-\rmi\frac{t}{N}H_2}\rme^{-\rmi\frac{t}{N}H_1}$ from Thm. <ref> and Thm. <ref>.
“$H_1H_2H_1$” and “$H_2H_1H_2$” indicate the scalings of the bounds on the second-order Trotter products with cycles $\rme^{-\rmi\frac{t}{2N}H_1}\rme^{-\rmi\frac{t}{N}H_2}\rme^{-\rmi\frac{t}{2N}H_1}$ and $\rme^{-\rmi\frac{t}{2N}H_2}\rme^{-\rmi\frac{t}{N}H_1}\rme^{-\rmi\frac{t}{2N}H_2}$, respectively, from Thm. <ref> and Thm. <ref>.
Blank items in the table indicate that the corresponding entry is not applicable.

c]@c@ $\ell$
c]@c@ $\delta_1$
c]@c@ $\delta_2$
c]@c@ $\delta_{12}$
c]@c@ $\delta_{21}$
c]@c@ $H_2H_1$
c]@c@ $H_1H_2H_1$
c]@c@ $H_2H_1H_2$

$0$   ${5}/{4}$   ${3}/{2}$           $N^{-1/4}$   $N^{-1/4}$   $N^{-1/4}$

$1$   ${7}/{4}$   ${5}/{2}$       ${1}/{4}$   $N^{-3/4}$   $N^{-3/4}$   $N^{-3/4}$

$2$   ${9}/{4}$   ${7}/{2}$   ${1}/{2}$   ${3}/{4}$   $N^{-1}$   $N^{-5/4}$   $N^{-5/4}$

$3$   ${11}/{4}$   ${9}/{2}$   ${1}/{2}$   ${5}/{4}$   $N^{-1}$   $N^{-3/2}$   $N^{-7/4}$

$\ell\ge4$   $\frac{1}{2}\ell+\frac{5}{4}$   $\ell+\frac{3}{2}$   $\ge\frac{5}{2}$   $\frac{1}{2}\ell-\frac{1}{4}$   $N^{-1}$   $N^{-2}$   $N^{-2}$

§.§ Fractional scalings of general $p$th-order Trotter bounds

The above strategy for estimating the fractional scaling works for a general $p$th-order Trotter product when $\varphi\not\in\mathcal{D}\bm{(}H^{(p)}(s)S_p(s)\bm{)}$.
Inserting the expression (<ref>) for the $q$th-order integral action into the remainder term (<ref>) for a $p$th-order Trotter product, we get

\begin{align}
\nonumber\\
&\ \
(-\rmi)^{q+1}\int_0^{T_M}\rmd s\,U^{(p)}(T_M)U^{(p)\dag}(s)H^{(p)}(s)S_q(s)\varphi
\nonumber\\
&\ \
\sum_{j=1}^M
\mathop{\mathop{\sum}_{\bm{\beta}\in\mathbb{N}^j}}_{|\bm{\beta}|=q}
\rme^{-\rmi\frac{t}{N}\tau_MH_{\sigma(M)}}
\cdots
\rme^{-\rmi\frac{t}{N}\tau_{j+1}H_{\sigma(j+1)}}
\nonumber\\
\int_0^{\frac{t}{N}}\rmd s\,
\rme^{-\rmi(\frac{t}{N}-s) \tau_j H_{\sigma(j)}}
\frac{1}{\beta_j!}
s^{\beta_j} (\tau_j H_{\sigma(j)})^{\beta_j+1}
\left(
\frac{t}{N}
\right)^{q-\beta_j}
\frac{1}{\beta_{j-1}!}
\cdots
\frac{1}{\beta_1!}
\varphi
\nonumber
\displaybreak[0]
\\
&\ \
\sum_{j=1}^M
\mathop{\mathop{\sum}_{\bm{\beta}\in\mathbb{N}^j}}_{|\bm{\beta}|=q}
\rme^{-\rmi\frac{t}{N}\tau_MH_{\sigma(M)}}
\cdots
\rme^{-\rmi\frac{t}{N}\tau_{j+1}H_{\sigma(j+1)}}
\nonumber\\
\Biggl[
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\Biggr]
\left(
\right)^{q-\beta_j}
\frac{1}{\beta_{j-1}!}
\cdots
\frac{1}{\beta_1!}
\varphi,
\label{eq:remainder_q}
\end{align}

where we have used the Taylor expansion (<ref>) of Remark <ref>.
It can be bounded by triangle inequality to yield an error bound

\begin{align*}
\xi_N^{(p)}(t;\varphi)
\le
\mathop{\mathop{\sum}_{\bm{\beta}\in\mathbb{N}^j}}_{|\bm{\beta}|=q}
\frac{t}{N}
\right)^{q-\beta_j}
\frac{|\tau_{j-1}|^{\beta_{j-1}}\cdots|\tau_1|^{\beta_1}}{\beta_{j-1}!\cdots\beta_1!}
\nonumber\\[-4truemm]
\Biggl\|
\Biggl[
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\Biggr]
\cdots
\varphi
\Biggr\|.
\end{align*}

This is just a rough bound.
By taking care of cancellations among terms due to $(H_1+H_2)\varphi=0$ before taking the norm, one can get a better bound.
For instance, for $q=1$, the error operator (<ref>) simplifies to

\begin{align*}
\sum_{j=1}^M
\rme^{-\rmi\frac{t}{N}\tau_MH_{\sigma(M)}}
\cdots
\rme^{-\rmi\frac{t}{N}\tau_{j+1}H_{\sigma(j+1)}}
\nonumber\\
\left[
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\rmi\frac{t}{N}
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\sum_{k=1}^{j-1}
\tau_k
\right]
\varphi
\nonumber
\displaybreak[0]
\\
\sum_{j=1}^M
\rme^{-\rmi\frac{t}{N}\tau_MH_{\sigma(M)}}
\cdots
\rme^{-\rmi\frac{t}{N}\tau_{j+1}H_{\sigma(j+1)}}
\nonumber\\
\left[
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\rmi\frac{t}{N}
\tau_j
\right)
\rmi\frac{t}{N}
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\sum_{k=1}^{j-1}
\tau_k
\right]
\varphi,
\end{align*}

and the Trotter error is bounded by

\begin{equation*}
\xi_N^{(p)}(t;\varphi)
\le
\left\|
\left[
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\rmi\frac{t}{N}
\tau_j
\right)
\rmi\frac{t}{N}
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\sum_{k=1}^{j-1}
\tau_k
\right]
\varphi
\right\|.
\end{equation*}

This is valid for $\varphi\in\mathcal{D}(H_1)\cap\mathcal{D}(H_2)$, and reproduces Eq. (<ref>) of Lemma <ref> for the first-order Trotterization and Eq. (<ref>) of Lemma <ref> for the second-order Trotterization.
Using the bounds in Eqs. (<ref>) and (<ref>), we see that it scales as
$\xi_N^{(p)}(t;\varphi)=\mathcal{O}(N^{-\delta})$ with some $\delta\in(0,1)$, if $\varphi\in\mathcal{D}(H_1)\cap\mathcal{D}(H_2)$ but $\varphi\not\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2^2)$, which is the case for the hydrogen eigenstates with $\ell=0,1$.
See Table <ref>.

For $q=2$, the error operator (<ref>) yields

\begin{align*}
\sum_{j=1}^M
\rme^{-\rmi\frac{t}{N}\tau_MH_{\sigma(M)}}
\cdots
\rme^{-\rmi\frac{t}{N}\tau_{j+1}H_{\sigma(j+1)}}
\nonumber\\
\Biggl[
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\nonumber\\
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\sum_{k=1}^{j-1}
\tau_kH_{\sigma(k)}
\nonumber\\
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\left(
\frac{1}{2}
\sum_{k=1}^{j-1}
\tau_k^2H_{\sigma(k)}^2
\sum_{k_1=2}^{j-1}
\sum_{k_2=1}^{k_1-1}
\tau_{k_1}H_{\sigma(k_1)}
\tau_{k_2}H_{\sigma(k_2)}
\right)
\Biggr]\,
\varphi
\nonumber\\
\sum_{j=1}^M
\rme^{-\rmi\frac{t}{N}\tau_MH_{\sigma(M)}}
\cdots
\rme^{-\rmi\frac{t}{N}\tau_{j+1}H_2}
\nonumber\\
\Biggl[
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\nonumber\\
\,\Biggl(
\sum_{k=1}^{j-1}
\Biggr)
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\nonumber\\
\,\Biggl(
\frac{1}{2}
\sum_{k_1=1}^{\lfloor\frac{j}{2}\rfloor}
\sum_{k_2=1}^{\lfloor\frac{j}{2}\rfloor}
\tau_{2k_1-1}
\tau_{2k_2-1}
\sum_{k_1=2}^{\lfloor\frac{j}{2}\rfloor}
\sum_{k_2=1}^{k_1-1}
\tau_{2k_1-1}
\tau_{2k_2}
\Biggr)
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\nonumber\\
\,\Biggl(
\frac{1}{2}
\sum_{k_1=1}^{\lceil\frac{j}{2}\rceil-1}
\sum_{k_2=1}^{\lceil\frac{j}{2}\rceil-1}
\tau_{2k_1}
\tau_{2k_2}
\sum_{k_1=1}^{\lceil\frac{j}{2}\rceil-1}
\sum_{k_2=1}^{k_1}
\tau_{2k_1}
\tau_{2k_2-1}
\Biggr)
\left(
\rme^{-\rmi\frac{t}{N}\tau_jH_{\sigma(j)}}
\right)
\Biggr]\,
\varphi,
\end{align*}

and the Trotter error is bounded by

\begin{align*}
\xi_N^{(p)}(t;\varphi)
\le{}&
\,\Biggl\|
\Biggl[
\left(
\rme^{-\rmi\frac{t}{N}\tau_{2i-1}H_1}
\right)
\nonumber\\
\,\Biggl(
\sum_{k=1}^{2i-2}
\Biggr)
\left(
\rme^{-\rmi\frac{t}{N}\tau_{2i-1}H_1}
\right)
\nonumber\\
\,\Biggl(
\frac{1}{2}
\sum_{k_1=1}^{i-1}
\sum_{k_2=1}^{i-1}
\tau_{2k_1-1}
\tau_{2k_2-1}
\sum_{k_1=2}^{i-1}
\sum_{k_2=1}^{k_1-1}
\tau_{2k_1-1}
\tau_{2k_2}
\Biggr)
\left(
\rme^{-\rmi\frac{t}{N}\tau_{2i-1}H_1}
\right)
\Biggr]\,
\varphi
\Biggr\|
\nonumber\\
\,\Biggl\|
\Biggl[
\left(
\rme^{-\rmi\frac{t}{N}\tau_{2i}H_2}
\right)
\nonumber\\
\,\Biggl(
\sum_{k=1}^{2i-1}
\Biggr)
\left(
\rme^{-\rmi\frac{t}{N}\tau_{2i}H_2}
\right)
\nonumber\\
\,\Biggl(
\frac{1}{2}
\sum_{k_1=1}^{i-1}
\sum_{k_2=1}^{i-1}
\tau_{2k_1}
\tau_{2k_2}
\sum_{k_1=1}^{i-1}
\sum_{k_2=1}^{k_1}
\tau_{2k_1}
\tau_{2k_2-1}
\Biggr)
\left(
\rme^{-\rmi\frac{t}{N}\tau_{2i}H_2}
\right)
\Biggr]\,
\varphi
\Biggr\|
\nonumber\\
\frac{t^2}{N}
\sum_{i=1}^{\lceil\frac{M}{2}\rceil}
\,\Biggl(
\frac{1}{2}
\sum_{k_1=1}^{i-1}
\sum_{k_2=1}^{i-1}
\tau_{2k_1}
\tau_{2k_2}
\sum_{k_1=1}^{i-1}
\sum_{k_2=1}^{k_1}
\tau_{2k_1}
\tau_{2k_2-1}
\Biggr)
\left\|
\left(
\rme^{-\rmi\frac{t}{N}\tau_{2i-1}H_1}
\right)
\varphi
\right\|
\nonumber\\
\frac{t^2}{N}
\sum_{i=1}^{\lfloor\frac{M}{2}\rfloor}
\,\Biggl(
\frac{1}{2}
\sum_{k_1=1}^i
\sum_{k_2=1}^i
\tau_{2k_1-1}
\tau_{2k_2-1}
\sum_{k_1=2}^i
\sum_{k_2=1}^{k_1-1}
\tau_{2k_1-1}
\tau_{2k_2}
\Biggr)
\left\|
\left(
\rme^{-\rmi\frac{t}{N}\tau_{2i}H_2}
\right)
\varphi
\right\|.
\end{align*}

This is valid for $\varphi\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2^2)$, and reproduces Eq. (<ref>) of Lemma <ref> for the second-order Trotterization.
It scales as $\xi_N^{(p)}(t;\varphi)=\mathcal{O}(N^{-\delta})$ with some $\delta\in(1,2)$, if $\varphi\in\mathcal{D}(H_1^2)\cap\mathcal{D}(H_2^2)$ but $\varphi\not\in\mathcal{D}(H_1^3)\cap\mathcal{D}(H_1H_2^2)\cap\mathcal{D}(H_2H_1^2)\cap\mathcal{D}(H_2^3)$, which is the case for the hydrogen eigenstates with $\ell=2,3$.
See Table <ref> again.

The method presented here is iterative and can be extended to all $q\leq p$ for a $p$th-order product formula.
Thereby, error bounds with fractional scalings are obtained for any state $\varphi\not\in\mathcal{D}\bm{(}H^{(p)}(s)S_p(s)\bm{)}$.

§ MATHEMATICA SCRIPT FOR COMPUTING HIGHER-ORDER TROTTER BOUNDS

In this appendix, we provide a script written in Wolfram Mathematica to compute the higher-order Trotter error.
As an input, one has to give the order $p$ of the product formula as well as the switching coefficients $\{\tau_j\}$.

§.§ Compute $S_p(s)$ at time step $j$

First, we compute $S_p(s)$ by Eq. (<ref>).
For this purpose, we use the string $a$ to denote the Hamiltonian $H_1$ and the string $b$ to denote the Hamiltonian $H_2$.
[breakable,sharp corners,colframe=white,left=1mm]
ClearAll[T, KK, t, n, s];

Definition of a power of the Hamiltonian.
H = If[EvenQ[K], "b", "a"];
Return[StringRepeat[H, power]],
Return[Hold[Ham[K, power]]];

All combinations to distribute the integer p in lists of length L.
This function solves the non-commutative integer partition problem to distribute
the exponents in the j-th term of the integral action.
AllSums = IntegerPartitions[p, len];
If[Length[AllSums] > 1,
AllSums = MapAt[Permutations, AllSums, All];
AllSums = Flatten[AllSums, 1],
AllSums = Permutations @@ AllSums;

Use the fact that the input state is an eigenvector of eigenvalue 0 to simplify.
expressionVector_(*integral action:array*)
Module[result, bookKeeping,
(*If the position is divisible by 2 but not by 4: move the entry to the left.*)
(*The next entry after this: move it to the right.*)
result = expressionVector;
bookKeeping = False;(*helper variable*)
\[And] Mod[i, 4] != 0,
result[[i - 1]] -= result[[i]];
result[[i]] = 0;
bookKeeping = True,
result[[i + 1]] -= result[[i]];
result[[i]] = 0;
bookKeeping = False;
, {i, 1, Length[expressionVector]}](*Table*);

Compute the integral action at time step j.
The result is stored in a vector of length 2^p, where each entry corresponds to a
product of operators H_1="a" and H_2="b", sorted in alphabetical order.
tauList_(*List of the switching coefficients:array*),
boundaryQ_ : False(*whether or not to compute the integral action
at the boundary:bool*),
simplifyQ_ : False(*Whether or not to use the function `ZeroEigenstate` to
simplify the result:bool*)
Module[{SummationList, tuples, delete, result, ii, L, AllSums,
times, Hamiltonians, jj, power, prefactor, HamTerm, position},
tuples = Map[StringJoin, Tuples[Sort@{"a", "b"}, order], 2];
delete = 0;(*deletion parameter q*)
result = ConstantArray[0, 2^order];
While[delete < order,
L = order - delete;
AllSums = Combinations[order, L];
For[ii = 1, ii <= Length[AllSums], ii++,
prefactor = 1/Times @@ Map[Factorial, AllSums[[ii]]];
power = AllSums[[ii, 1]];
times = T[KK[1]]^power;
Hamiltonians = {Ham[KK[1], power]};
power = AllSums[[ii, jj]];
times *= T[KK[jj]]^power;
Hamiltonians = AppendTo[Hamiltonians, Ham[KK[jj], power]]
, {jj, 2, L}
HamTerm = Map[StringJoin, Evaluate @@@ Hamiltonians, All];
position = Position[tuples, HamTerm][[1]];
result[[position]] += prefactor*times
, Evaluate[
Table[{KK[index], 1, KK[index - 1] - 1}, {index, 2, L}]
,{KK[1], 1, mMax}](*Prepend*)
If[boundaryQ == False,
result =
result /.
T[mMax] -> (time - Sum[Abs[T[ll]], {ll, 1, mMax - 1}])*
](*For ii*);
result = result /. Table[T[j] -> tauList[[j]], {j, 1, mMax}];
\end{Verbatim}
\end{tcolorbox}

\subsection{Compute the error bound $\xi_N^{(p)}(t;\varphi)$}
Having defined $S_p(s)$ for time $s$ in the $j$th time slot, we can now compute the error bounds using Eq.~\eqref{eqn:error_estimate_p} of<EMAIL_ADDRESS>Here, we will use the convention that $(H_1+H_2)\varphi=0$.
\begin{tcolorbox}[breakable,sharp corners,colframe=white,left=1mm]
\begin{Verbatim}[numbers=left, firstnumber=129]
Multiply the vector representing the integral action at time step j by the Hamiltonian
expressionVector_(*integral action:array*),
Module[{zeros, result},
zeros = ConstantArray[0, Length[expressionVector]];
result = Join[zeros, expressionVector],
result = Join[expressionVector, zeros];

Compute the actual error bound.
order_(*p:int*), mMax_(*M:int*),
tauList_(*List of the switching coefficients:array*),
simplifyQ_ : True(*Whether or not to use the function `ZeroEigenstate` to
simplify the result:bool*)
Module[{Action, integralBounds, bound, operators},
integralBounds = Table[Sum[Abs[tauList[[j]]], {j, 1, i}], {i, 1, Length[tauList]}];
integralBounds = PrependTo[integralBounds, 0];
bound = (t^(order + 1)/n^order)*Sum[
Sp[order, i, s, tauList, False, simplifyQ]
, {s, integralBounds[[i]], integralBounds[[i + 1]]}](*Integrate*)
, i](*MultiplyBy*)
, {i, 1, mMax}](*Sum*);
operators =
Tuples[Sort@{"a", "b"}, order + 1]
// Map[Append["\[CurlyPhi]||"]] // Map[Prepend["||"]]
, 2](*Map*);
Return[Chop[bound, 10^-MachinePrecision].operators];
\end{Verbatim}
\end{tcolorbox}

\subsection{Examples}\label{sec:Mathematica_examples}
To show how the code works, we compute the examples of $p=1,2,4,6$.
For $p=1$, we enter
\begin{tcolorbox}[breakable,sharp corners,colframe=white,left=1mm]
\begin{Verbatim}[numbers=left, firstnumber=175]
ErrorBound[1, 2, {1, 1}, True]
\end{Verbatim}
\end{tcolorbox}
\noindent
and get the output
\begin{equation*}
\frac{\text{$||$aa$\varphi ||$} t^2}{2 n}+\frac{\text{$||$bb$\varphi ||$} t^2}{2 n}.
\end{equation*}
For $p=2$, we enter
\begin{tcolorbox}[breakable,sharp corners,colframe=white,left=1mm]
\begin{Verbatim}[numbers=left, firstnumber=176]
ErrorBound[2, 3, {1/2, 1, 1/2}, True]
\end{Verbatim}
\end{tcolorbox}
\noindent
to obtain
\begin{equation*}
\frac{\text{$||$aaa$\varphi ||$} t^3}{24 n^2}+\frac{\text{$||$baa$\varphi ||$} t^3}{8 n^2}+\frac{\text{$||$bbb$\varphi ||$} t^3}{12 n^2}.
\end{equation*}
The case $p=4$ can be computed with
\begin{tcolorbox}[breakable,sharp corners,colframe=white,left=1mm]
\begin{Verbatim}[numbers=left, firstnumber=177]
s4 = 1/(2 - 2^(1/3));
times4 = {s4/2, s4, (1-s4)/2, 1-2s4, (1-s4)/2, s4, s4/2}//N;
ErrorBound[4, 7, times4, True]
\end{Verbatim}
\end{tcolorbox}
\noindent
and results in
\begin{align*}
&\frac{0.00427803 \text{$||$aaaaa$\varphi ||$} t^5}{n^4}
+\frac{0.00950917 \text{$||$aabaa$\varphi ||$} t^5}{n^4}
+\frac{0.00633945 \text{$||$aabbb$\varphi ||$}
+\frac{0.0243899 \text{$||$abaaa$\varphi ||$} t^5}{n^4}
\\
&+\frac{0.0731697 \text{$||$abbaa$\varphi ||$} t^5}{n^4}
+\frac{0.0487798 \text{$||$abbbb$\varphi ||$}
+\frac{0.0299684 \text{$||$baaaa$\varphi ||$} t^5}{n^4}
+\frac{0.092188 \text{$||$babaa$\varphi ||$} t^5}{n^4}
\\
&+\frac{0.0614587 \text{$||$babbb$\varphi ||$}
+\frac{0.212061 \text{$||$bbaaa$\varphi ||$} t^5}{n^4}
+\frac{0.393394 \text{$||$bbbaa$\varphi ||$} t^5}{n^4}
+\frac{0.22989 \text{$||$bbbbb$\varphi ||$}
\end{align*}
Finally, for $p=6$, we input
\begin{tcolorbox}[breakable,sharp corners,colframe=white,left=1mm]
\begin{Verbatim}[numbers=left, firstnumber=180]
s4 = 1/(2 - 2^(1/3));
s6 = 1/(2 - 2^(1/5));
times6 = {
s4*s6/2, s4*s6, (1-s4)s6/2, (1-2 s4)s6, (1-s4)s6/2, s4*s6,
s4*s6/2+s4*(1-2 s6)/2, s4 (1-2s6), (1-s4)(1-2s6)/2, (1-2s4)(1-2s6),
(1-s4)(1-2s6)/2, s4(1-2s6), s4(1-2s6)/2+s4*s6/2, s4*s6,
(1-s4)s6/2, (1-2s4)s6, (1-s4)s6/2, s4*s6, s4*s6/2
ErrorBound[6, 19, times6, True]
\end{Verbatim}
\end{tcolorbox}
\noindent
and obtain
\begingroup
\allowdisplaybreaks
\begin{align*}
&\frac{0.000342245 \text{$||$aaaaaaa$\varphi ||$} t^7}{n^6}+\frac{0.00466073 \text{$||$aaaabaa$\varphi ||$} t^7}{n^6}+\frac{0.00310716
\text{$||$aaaabbb$\varphi ||$} t^7}{n^6}+\frac{0.0128515 \text{$||$aaabaaa$\varphi ||$} t^7}{n^6}\\
&+\frac{0.0385546 \text{$||$aaabbaa$\varphi ||$}
t^7}{n^6}+\frac{0.0257031 \text{$||$aaabbbb$\varphi ||$} t^7}{n^6}+\frac{0.0128622 \text{$||$aabaaaa$\varphi ||$} t^7}{n^6}+\frac{0.0417926
\text{$||$aababaa$\varphi ||$} t^7}{n^6}\\
&+\frac{0.0278617 \text{$||$aababbb$\varphi ||$} t^7}{n^6}+\frac{0.0164871 \text{$||$aabbaaa$\varphi ||$}
t^7}{n^6}+\frac{0.0382987 \text{$||$aabbbaa$\varphi ||$} t^7}{n^6}+\frac{0.0240441 \text{$||$aabbbbb$\varphi ||$} t^7}{n^6}\\
\text{$||$abaaaaa$\varphi ||$} t^7}{n^6}+\frac{0.0241779 \text{$||$abaabaa$\varphi ||$} t^7}{n^6}+\frac{0.0161186 \text{$||$abaabbb$\varphi ||$}
t^7}{n^6}+\frac{0.0620135 \text{$||$ababaaa$\varphi ||$} t^7}{n^6}\\
&+\frac{0.18604 \text{$||$ababbaa$\varphi ||$} t^7}{n^6}+\frac{0.124027
\text{$||$ababbbb$\varphi ||$} t^7}{n^6}+\frac{0.0213565 \text{$||$abbaaaa$\varphi ||$} t^7}{n^6}+\frac{0.234396 \text{$||$abbabaa$\varphi ||$}
&+\frac{0.156264 \text{$||$abbabbb$\varphi ||$} t^7}{n^6}+\frac{0.168267 \text{$||$abbbaaa$\varphi ||$} t^7}{n^6}+\frac{0.205275
\text{$||$abbbbaa$\varphi ||$} t^7}{n^6}+\frac{0.0969135 \text{$||$abbbbbb$\varphi ||$} t^7}{n^6}\\
&+\frac{0.00379873 \text{$||$baaaaaa$\varphi ||$}
t^7}{n^6}+\frac{0.0501259 \text{$||$baaabaa$\varphi ||$} t^7}{n^6}+\frac{0.0334173 \text{$||$baaabbb$\varphi ||$} t^7}{n^6}+\frac{0.137266
\text{$||$baabaaa$\varphi ||$} t^7}{n^6}\\
&+\frac{0.411798 \text{$||$baabbaa$\varphi ||$} t^7}{n^6}+\frac{0.274532 \text{$||$baabbbb$\varphi ||$}
t^7}{n^6}+\frac{0.127458 \text{$||$babaaaa$\varphi ||$} t^7}{n^6}+\frac{0.420193 \text{$||$bababaa$\varphi ||$} t^7}{n^6}\\
\text{$||$bababbb$\varphi ||$} t^7}{n^6}+\frac{0.290088 \text{$||$babbaaa$\varphi ||$} t^7}{n^6}+\frac{0.467689 \text{$||$babbbaa$\varphi ||$}
t^7}{n^6}+\frac{0.258116 \text{$||$babbbbb$\varphi ||$} t^7}{n^6}\\
&+\frac{0.0631099 \text{$||$bbaaaaa$\varphi ||$} t^7}{n^6}+\frac{0.141983
\text{$||$bbaabaa$\varphi ||$} t^7}{n^6}+\frac{0.0946554 \text{$||$bbaabbb$\varphi ||$} t^7}{n^6}+\frac{0.36417 \text{$||$bbabaaa$\varphi ||$}
&+\frac{1.09251 \text{$||$bbabbaa$\varphi ||$} t^7}{n^6}+\frac{0.72834 \text{$||$bbabbbb$\varphi ||$} t^7}{n^6}+\frac{0.0509254
\text{$||$bbbaaaa$\varphi ||$} t^7}{n^6}+\frac{1.05772 \text{$||$bbbabaa$\varphi ||$} t^7}{n^6}\\
&+\frac{0.705143 \text{$||$bbbabbb$\varphi ||$}
t^7}{n^6}+\frac{0.67482 \text{$||$bbbbaaa$\varphi ||$} t^7}{n^6}+\frac{0.679037 \text{$||$bbbbbaa$\varphi ||$} t^7}{n^6}+\frac{0.281524
\text{$||$bbbbbbb$\varphi ||$} t^7}{n^6}.
\end{align*}
\endgroup
\end{widetext}

\bibliography{hydrogen.bib}
\end{document}
\begin{thebibliography}{10}

[1]
R.~P. Feynman, Simulating physics with computers, \href{https://doi.org/10.1007/bf02650179}{Int. J. Theor. Phys. \textbf{21},  467  (1982)}.

[2]
S. Lloyd, Universal Quantum Simulators, \href{https://doi.org/10.1126/science.273.5278.1073}{Science \textbf{273},  1073  (1996)}.

[3]
S. Lee, J. Lee, and G.~K.-L. Chan, Collection of papers referring to exponential quantum advantage in quantum chemistry, GitHub, \href{https://github.com/seunghoonlee89/Refs_EQA_GSQC}{https://github.com/seunghoonlee89/Refs\textunderscore EQA\textunderscore GSQC}, 2023.

[4]
B.~P. Lanyon, J.~D. Whitfield, G.~G. Gillett, M.~E. Goggin, M.~P. Almeida, I. Kassal, J.~D. Biamonte, M. Mohseni, B.~J. Powell, M. Barbieri, A. Aspuru-Guzik, and A.~G. White, Towards quantum chemistry on a quantum computer, \href{https://doi.org/10.1038/nchem.483}{Nature Chem. \textbf{2},  106  (2010)}.

[5]
D. Wecker, B. Bauer, B.~K. Clark, M.~B. Hastings, and M. Troyer, Gate-count estimates for performing quantum chemistry on small quantum computers, \href{https://doi.org/10.1103/physreva.90.022305}{Phys. Rev. A \textbf{90},  022305  (2014)}.

[6]
R. Babbush, J. McClean, D. Wecker, A. Aspuru-Guzik, and N. Wiebe, Chemical basis of {Trotter-Suzuki} errors in quantum chemistry simulation, \href{https://doi.org/10.1103/physreva.91.022311}{Phys. Rev. A \textbf{91},  022311  (2015)}.

[7]
D. Poulin, M.~B. Hastings, D. Wecker, N. Wiebe, A.~C. Doberty, and M. Troyer, The {Trotter} step size required for accurate quantum simulation of quantum chemistry, \href{https://doi.org/10.26421/QIC15.5-6-1}{Quant. Inf. Comp. \textbf{15},  361  (2015)}.

[8]
R. Babbush, N. Wiebe, J. McClean, J. McClain, H. Neven, and G.~K.-L. Chan, Low-Depth Quantum Simulation of Materials, \href{https://doi.org/10.1103/PhysRevX.8.011044}{Phys. Rev. X \textbf{8},  011044  (2018)}.

[9]
Y. Cao, J. Romero, J.~P. Olson, M. Degroote, P.~D. Johnson, M. Kieferov{\'{a}}, I.~D. Kivlichan, T. Menke, B. Peropadre, N.~P.~D. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik, Quantum Chemistry in the Age of Quantum Computing, \href{https://doi.org/10.1021/acs.chemrev.8b00803}{Chem. Rev. \textbf{119},  10856  (2019)}.

[10]
S. McArdle, S. Endo, A. Aspuru-Guzik, S.~C. Benjamin, and X. Yuan, Quantum computational chemistry, \href{https://doi.org/10.1103/RevModPhys.92.015003}{Rev. Mod. Phys. \textbf{92},  015003  (2020)}.

[11]
Y.-h. Lam, Y. Abramov, R.~S. Ananthula, J.~M. Elward, L.~R. Hilden, S.~O. Nilsson~Lill, P.-O. Norrby, A. Ramirez, E.~C. Sherer, J. Mustakis, and G.~J. Tanoury, Applications of Quantum Chemistry in Pharmaceutical Process Development: Current State and Opportunities, \href{https://doi.org/10.1021/acs.oprd.0c00222}{Org. Process Res. Dev. \textbf{24},  1496  (2020)}.

[12]
S. Ahn, M. Hong, M. Sundararajan, D.~H. Ess, and M.-H. Baik, Design and Optimization of Catalysts Based on Mechanistic Insights Derived from Quantum Chemical Reaction Modeling, \href{https://doi.org/10.1021/acs.chemrev.9b00073}{Chem. Rev. \textbf{119},  6509  (2019)}.

[13]
D.~W. Berry, C. Gidney, M. Motta, J.~R. McClean, and R. Babbush, Qubitization of Arbitrary Basis Quantum Chemistry Leveraging Sparsity and Low Rank Factorization, \href{https://doi.org/10.22331/q-2019-12-02-208}{Quantum \textbf{3},  208  (2019)}.

[14]
J. Lee, D.~W. Berry, C. Gidney, W.~J. Huggins, J.~R. McClean, N. Wiebe, and R. Babbush, Even More Efficient Quantum Computations of Chemistry Through Tensor Hypercontraction, \href{https://doi.org/10.1103/prxquantum.2.030305}{{PRX} {Q}uantum \textbf{2},  030305  (2021)}.

[15]
H.~F. Trotter, On the product of semi-groups of operators, \href{https://doi.org/10.1090/S0002-9939-1959-0108732-6}{Proc. Am. Math. Soc. \textbf{10},  545  (1959)}.

[16]
T. Kato, Trotter's product formula for an arbitrary pair of self-adjoint contraction semigroups,  in \href{https://searchworks.stanford.edu/view/746600}{\textit{Topics in Functional Analysis: Essays Dedicated to {M. G. Krein} on the Occasion of His 70th Birthday}}, edited by I. Gohberg and M. Kac (Academic Press, New York, 1978), pp.\ 185--195.

[17]
E. Nelson, Feynman Integrals and the {Schr\"odinger} Equation, \href{https://doi.org/10.1063/1.1704124}{J. Math. Phys. \textbf{5},  332  (1964)}.

[18]
G.~W. Johnson and M.~L. Lapidus, \href{https://global.oup.com/academic/product/the-feynman-integral-and-feynmans-operational-calculus-9780198515722}{\textit{{T}he {Feynman} {I}ntegral and {F}eynman's {O}perational {C}alculus}} (Oxford University Press, New York, 2002).

[19]
G.~B. Folland, \href{https://bookstore.ams.org/surv-149}{\textit{Quantum Field Theory: A Tourist Guide for Mathematicians}} (American Mathematical Society, Providence, RI, 2008).

[20]
R.~D. Klauber, \href{http://www.quantumfieldtheory.info/}{\textit{Student Friendly Quantum Field Theory}, 2nd ed.} (Sandtrove Press, Fairfield, 2013).

[21]
B. Simon, \href{https://bookstore.ams.org/chel-351-h}{\textit{Functional Integration and Quantum Physics}, 2nd ed.} (American Mathematical Society, Providence, RI, 2005).

[22]
F. Nicola and S.~I. Trapasso, On the Pointwise Convergence of the Integral Kernels in the {Feynman-Trotter} Formula, \href{https://doi.org/10.1007/s00220-019-03524-2}{Commun. Math. Phys. \textbf{376},  2277  (2020)}.

[23]
B. Gaveau, E. Mih{\'{o}}kov{\'{a}}, M. Roncadelli, and L.~S. Schulman, Path integral in a magnetic field using the {Trotter} product formula, \href{https://doi.org/10.1119/1.1630334}{Am. J. Phys. \textbf{72},  385  (2004)}.

[24]
R. Kosloff, Time-dependent quantum-mechanical methods for molecular dynamics, \href{https://doi.org/10.1021/j100319a003}{J. Phys. Chem. \textbf{92},  2087  (1988)}.

[25]
D.~J. Tannor, \href{https://uscibooks.aip.org/books/introduction-to-quantum-mechanics-a-time-dependent-perspective/}{\textit{Introduction to Quantum Mechanics: A Time-Dependent Perspective}} (University Science Books, California, 2007).

[26]
S. Choi and J. Van{\'{\i}}{\v{c}}ek, Efficient geometric integrators for nonadiabatic quantum dynamics. I. The adiabatic representation, \href{https://doi.org/10.1063/1.5092611}{J. Chem. Phys. \textbf{150},  204112  (2019)}.

[27]
J. Roulet, S. Choi, and J. Van{\'{\i}}{\v{c}}ek, Efficient geometric integrators for nonadiabatic quantum dynamics. {II}. The diabatic representation, \href{https://doi.org/10.1063/1.5094046}{J. Chem. Phys. \textbf{150},  204113  (2019)}.

[28]
H.~H.~S. Chan, R. Meister, T. Jones, D.~P. Tew, and S.~C. Benjamin, Grid-based methods for chemistry simulations on a quantum computer, \href{https://doi.org/10.1126/sciadv.abo7484}{Sci. Adv. \textbf{9},    (2023)}.

[29]
F. Fitzek, J.-N. Siem\ss, S. Seckmeyer, H. Ahlers, E.~M. Rasel, K. Hammerer, and N. Gaaloul, Universal atom interferometer simulation of elastic scattering processes, \href{https://doi.org/10.1038/s41598-020-78859-1}{Sci. Rep. \textbf{10},  22120  (2020)}.

[30]
F. Verstraete, J.~J. Garc\'{\i}a-Ripoll, and J.~I. Cirac, Matrix Product Density Operators: Simulation of Finite-Temperature and Dissipative Systems, \href{https://doi.org/10.1103/PhysRevLett.93.207204}{Phys. Rev. Lett. \textbf{93},  207204  (2004)}.

[31]
M. Zwolak and G. Vidal, Mixed-State Dynamics in One-Dimensional Quantum Lattice Systems: A Time-Dependent Superoperator Renormalization Algorithm, \href{https://doi.org/10.1103/PhysRevLett.93.207205}{Phys. Rev. Lett. \textbf{93},  207205  (2004)}.

[32]
G. Vidal, Efficient Simulation of One-Dimensional Quantum Many-Body Systems, \href{https://doi.org/10.1103/PhysRevLett.93.040502}{Phys. Rev. Lett. \textbf{93},  040502  (2004)}.

[33]
L. Cincio, J. Dziarmaga, and M.~M. Rams, Multiscale Entanglement Renormalization Ansatz in Two Dimensions: Quantum Ising Model, \href{https://doi.org/10.1103/PhysRevLett.100.240603}{Phys. Rev. Lett. \textbf{100},  240603  (2008)}.

[34]
M. Reiher, N. Wiebe, K.~M. Svore, D. Wecker, and M. Troyer, Elucidating reaction mechanisms on quantum computers, \href{https://doi.org/10.1073/pnas.1619152114}{Proc. Natl. Acad. Sci. U.S.A. \textbf{114},  7555  (2017)}.

[35]
A.~M. Childs, Y. Su, M.~C. Tran, N. Wiebe, and S. Zhu, Theory of {Trotter} Error with Commutator Scaling, \href{https://doi.org/10.1103/PhysRevX.11.011020}{Phys. Rev. X \textbf{11},  011020  (2021)}.

[36]
A. Schubert and C.~B. Mendl, Trotter error with commutator scaling for the {Fermi-Hubbard} model, \href{https://doi.doi/10.1103/PhysRevB.108.195105}{Phys. Rev. B \textbf{108},  195105  (2023)}.

[37]
S. Zhuk, N. Robertson, and S. Bravyi, Trotter error bounds and dynamic multi-product formulas for {Hamiltonian} simulation, \href{https://doi.org/10.48550/arXiv.2306.12569}{arXiv:2306.12569  (2023)}.

[38]
K. Ito and F. Kappel, The {Trotter-Kato} theorem and approximation of {PDEs}, \href{https://doi.org/10.1090/s0025-5718-98-00915-6}{Math. Comp. \textbf{67},  21  (1998)}.

[39]
T. Ichinose and H. Tamura, Note on the Norm Convergence of the Unitary {Trotter} Product Formula, \href{https://doi.org/10.1007/s11005-004-3760-2}{Lett. Math. Phys. \textbf{70},  65  (2004)}.

[40]
D. Burgarth, N. Galke, A. Hahn, and L. van Luijk, State-dependent {Trotter} limits and their approximations, \href{https://doi.org/10.1103/PhysRevA.107.L040201}{Phys. Rev. A \textbf{107},  l040201  (2023)}.

[41]
E. Sch\"odinger, Quantisierung als Eigenwertproblem, \href{https://doi.org/10.1002/andp.19263840404}{Ann. Phys. \textbf{384},  361  (1926)}.

[42]
W. Pauli, \"Uber das Wasserstoffspektrum vom Standpunkt der neuen Quantenmechanik, \href{https://doi.org/10.1007/BF01450175}{Z. Phys. \textbf{36},  336  (1926)}.

[43]
T. Helgaker, P. J{\o}rgensen, and J. Olsen, \href{https://doi.org/10.1002/9781119019572}{\textit{Molecular Electronic-Structure Theory}} (Wiley, West Sussex, 2000).

[44]
J.~E. Huheey, Inorganic chemistry: Principles of structure and reactivity, \href{https://doi.org/10.1021/ed050pA379.1}{J. Chem. Educ. \textbf{50},  A379  (1973)}.

[45]
T. Jahnke and C. Lubich, Error Bounds for Exponential Operator Splittings, \href{https://doi.org/10.1023/A:1022396519656}{{BIT} Numer. Math. \textbf{40},  735  (2000)}.

[46]
M. Thalhammer, High-Order Exponential Operator Splitting Methods for Time-Dependent {Schr\"odinger} Equations, \href{https://doi.org/10.1137/060674636}{{SIAM} J. Numer. Anal. \textbf{46},  2022  (2008)}.

[47]
E. Hansen and A. Ostermann, Exponential splitting for unbounded operators, \href{https://doi.org/10.1090/s0025-5718-09-02213-3}{Math. Comp. \textbf{78},  1485  (2009)}.

[48]
C. Neuhauser and M. Thalhammer, On the convergence of splitting methods for linear evolutionary {Schr\"odinger} equations involving an unbounded potential, \href{https://doi.org/10.1007/s10543-009-0215-2}{{BIT} Numer. Math. \textbf{49},  199  (2009)}.

[49]
E. Kieri, Stiff convergence of force-gradient operator splitting methods, \href{https://doi.org/10.1016/j.apnum.2015.03.005}{Appl. Numer. Math. \textbf{94},  33  (2015)}.

[50]
B. {\c{S}}ahino{\u{g}}lu and R.~D. Somma, Hamiltonian simulation in the low-energy subspace, \href{https://doi.org/10.1038/s41534-021-00451-w}{npj Quantum Inf. \textbf{7},  119  (2021)}.

[51]
C. Yi and E. Crosson, Spectral analysis of product formulas for quantum simulation, \href{https://doi.org/10.1038/s41534-022-00548-w}{npj Quantum Inf. \textbf{8},  37  (2022)}.

[52]
D. An, D. Fang, and L. Lin, Time-dependent unbounded {Hamiltonian} simulation with vector norm scaling, \href{https://doi.org/10.22331/q-2021-05-26-459}{Quantum \textbf{5},  459  (2021)}.

[53]
T. Hatomura, State-dependent error bound for digital quantum simulation of driven systems, \href{https://doi.org/10.1103/PhysRevA.105.L050601}{Phys. Rev. A \textbf{105},  L050601  (2022)}.

[54]
T. Hatomura, The first-order Trotter decomposition in the dynamical-invariant basis, \href{https://doi.org/10.48550/arXiv.2308.15100}{arXiv:2308.15100  (2023)}.

[55]
H. Neidhardt and V.~A. Zagrebnov, Fractional powers of self-adjoint operators and Trotter-Kato product formula, \href{https://doi.org/10.1007/BF01196384}{Integr. Equ. Oper. Theory \textbf{35},  209  (1999)}.

[56]
T. Ichinose, H. Neidhardt, and V.~A. Zagrebnov, Trotter-Kato product formula and fractional powers of self-adjoint generators, \href{https://doi.org/10.1016/S0022-1236(03)00177-0}{J. Funct. Anal. \textbf{207},  33  (2004)}.

[57]
S.~A. Chin and P. Anisimov, Gradient symplectic algorithms for solving the radial {Schr\"odinger} equation, \href{https://doi.org/10.1063/1.2150831}{J. Chem. Phys. \textbf{124},  054106  (2006)}.

[58]
E. Buend{\'{\i}}a and R. Guardiola, The numerov method and singular potentials, \href{https://doi.org/10.1016/0021-9991(85)90038-5}{J. Comput. Phys. \textbf{60},  561  (1985)}.

[59]
S. Lee, J. Lee, H. Zhai, Y. Tong, A.~M. Dalzell, A. Kumar, P. Helms, J. Gray, Z.-H. Cui, W. Liu, M. Kastoryano, R. Babbush, J. Preskill, D.~R. Reichman, E.~T. Campbell, E.~F. Valeev, L. Lin, and G.~K.-L. Chan, Evaluating the evidence for exponential quantum advantage in ground-state quantum chemistry, \href{https://doi.org/10.1038/s41467-023-37587-6}{Nat. Commun. \textbf{14},    (2023)}.

[60]
M. Suzuki, Decomposition formulas of exponential operators and {Lie} exponentials with some applications to quantum mechanics and statistical physics, \href{https://doi.org/10.1063/1.526596}{J. Math. Phys. \textbf{26},  601  (1985)}.

[61]
D. Burgarth, P. Facchi, M. Fraas, and R. Hillier, {Non-Markovian noise that cannot be dynamically decoupled by periodic spin echo pulses}, \href{https://doi.org/10.21468/SciPostPhys.11.2.027}{SciPost Phys. \textbf{11},  027  (2021)}.

[62]
The bound could diverge for two different reasons: (i) one of $\varphi_\ell$ does not satisfy the corresponding domain condition; (ii) $L=\infty$ and, even if the domain conditions are satisfied for all $\varphi_\ell$, the sum of the series diverges.

[63]
Chemical accuracy means an error of the order $10^{-3}\,\mathrm{Hartee}$~\cite[p.~267]{Helgaker2014}, where $1\,\mathrm{Hartree} = \hbar^2/(m_\rme e^2 a_0^2)$.

[64]
G.~R. Dennis, J.~J. Hope, and M.~T. Johnsson, {XMDS}2: Fast, scalable simulation of coupled stochastic partial differential equations, \href{https://doi.org/10.1016/j.cpc.2012.08.016}{Comput. Phys. Comm. \textbf{184},  201  (2013)}.

[65]
C. Zhu and J.~R. Klauder, Classical symptoms of quantum illnesses, \href{https://doi.org/10.1119/1.17221}{Am. J. Phys. \textbf{61},  605  (1993)}.

[66]
C. Arenz, D. Burgarth, P. Facchi, and R. Hillier, Dynamical decoupling of unbounded {Hamiltonians}, \href{https://doi.org/10.1063/1.5016495}{J. Math. Phys. \textbf{59},  032203  (2018)}.

[67]
M.~E.~S. Morales, P.~C.~S. Costa, D.~K. Burgarth, Y.~R. Sanders, and D.~W. Berry, Greatly improved higher-order product formulae for quantum simulation, \href{https://doi.org/10.48550/arXiv.2210.15817}{arXiv:2210.15817  (2022)}.

[68]
M. Suzuki, General theory of fractal path integrals with applications to many-body theories and statistical physics, \href{https://doi.org/10.1063/1.529425}{J. Math. Phys. \textbf{32},  400  (1991)}.

[69]
D. Burgarth, P. Facchi, G. Gramegna, and K. Yuasa, One bound to rule them all: from {Adiabatic} to {Zeno}, \href{https://doi.org/10.22331/q-2022-06-14-737}{Quantum \textbf{6},  737  (2022)}.

[70]
M. Suzuki, Fractal decomposition of exponential operators with applications to many-body theories and {Monte Carlo} simulations, \href{https://doi.org/10.1016/0375-9601(90)90962-n}{Phys. Lett. A \textbf{146},  319  (1990)}.

[71]
H. Yoshida, Construction of higher order symplectic integrators, \href{https://doi.org/10.1016/0375-9601(90)90092-3}{Phys. Lett. A \textbf{150},  262  (1990)}.

[72]
B.~C. Hall, \href{https://doi.org/10.1007/978-1-4614-7116-5}{\textit{Quantum Theory for Mathematicians}} (Springer, New York, 2013).

[73]
A. Hahn, D. Burgarth, and K. Yuasa, Unification of random dynamical decoupling and the quantum {Zeno} effect, \href{https://doi.org/10.1088/1367-2630/ac6b4f}{New J. Phys. \textbf{24},  063027  (2022)}.

[74]
Y. Tong, V.~V. Albert, J.~R. McClean, J. Preskill, and Y. Su, Provably accurate simulation of gauge theories and bosonic systems, \href{https://doi.org/10.22331/q-2022-09-22-816}{Quantum \textbf{6},  816  (2022)}.

[75]
B. Peng, Y. Su, D. Claudino, K. Kowalski, G.~H. Low, and M. Roetteler, Quantum Simulation of {Boson}-Related {Hamiltonians}: Techniques, Effective {Hamiltonian} Construction, and Error Analysis, \href{https://doi.org/10.48550/arXiv.2307.06580}{arXiv:2307.06580  (2023)}.

[76]
P. Mukhopadhyay, T.~F. Stetina, and N. Wiebe, Quantum Simulation of the First-Quantized {Pauli-Fierz} {Hamiltonian}, \href{https://doi.org/10.48550/arXiv.2306.11198}{arXiv:2306.11198  (2023)}.

[77]
M. Abramowitz and I.~A. Stegun, \href{https://store.doverpublications.com/0486612724.html}{\textit{Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables}} (Dover Publications, Washington, D.C., 1972).

[78]
B. Podolsky and L. Pauling, The Momentum Distribution in Hydrogen-Like Atoms, \href{https://doi.org/10.1103/PhysRev.34.109}{Phys. Rev. \textbf{34},  109  (1929)}.

[79]
T. Koornwinder, The Addition Formula for {Laguerre} Polynomials, \href{https://doi.org/10.1137/0508041}{{SIAM} J. Math. Anal. \textbf{8},  535  (1977)}.

[80]
G. Szeg\H{o}, \href{https://doi.org/10.1090/coll/023}{\textit{Orthogonal Polynomials}, 4th ed.} (American Mathematical Society, Providence, RI, 1975).

[81]
L. Pauling and E. Bright~Wilson, Jr., \href{https://store.doverpublications.com/0486134938.html}{\textit{Introduction to Quantum Mechanics with Applications to Chemistry}} (Dover Publications, New York, 1985).

\end{thebibliography}
\end{document}

\end{document}