# A Theory of Digital Quantum Simulations in the
Low-Energy Subspace

Weiyuan Gong Shuo Zhou∗ Tongyang Li Equal contribution.John A. Paulson School
of Engineering and Applied Sciences, Harvard University; Institute for
Interdisciplinary Information Sciences, Tsinghua University. Email:
<EMAIL_ADDRESS>on Frontiers of Computing Studies, School of
Computer Science, Peking University; Yuanpei College, Peking University.
Email<EMAIL_ADDRESS>author. Center on Frontiers of
Computing Studies, School of Computer Science, Peking University. Email:
<EMAIL_ADDRESS>
###### Abstract

Digital quantum simulation has broad applications in approximating unitary
evolutions of Hamiltonians. In practice, many simulation tasks for quantum
systems focus on quantum states in the low-energy subspace instead of the
entire Hilbert space. In this paper, we systematically investigate the
complexity of digital quantum simulation based on product formulas in the low-
energy subspace. We show that the simulation error depends on the effective
low-energy norm of the Hamiltonian for a variety of digital quantum simulation
algorithms and quantum systems, allowing improvements over the previous
complexities for full unitary simulations even for imperfect state
preparations. In particular, for simulating spin models in the low-energy
subspace, we prove that randomized product formulas such as qDRIFT and random
permutation require smaller Trotter numbers. Such improvement also persists in
symmetry-protected digital quantum simulations. We prove a similar improvement
in simulating the dynamics of power-law quantum interactions. We also provide
a query lower bound for general digital quantum simulations in the low-energy
subspace.

##  1 Introduction

Quantum computers possess the potential to simulate the dynamics of given
quantum systems more efficiently compared to their classical counterparts [20,
18, 25]. A main approach for simulating Hamiltonian evolution is the digital
quantum simulation, which maps the quantum systems into qubits and
approximates the evolution using digitalized quantum gates. The most standard
digitalization approaches are product formulas, including the Suzuki-Trotter
formulas [36, 38, 37, 23, 13, 10, 4], which decompose the target Hamiltonian
into a sum of non-commuting terms and perform a product of exponentials of
each individual term. Alternative principles have also led to advanced digital
quantum simulation approaches such as linear combinations of unitaries [15,
43], qubitization [27], truncated Taylor series [6], quantum signal processing
[26, 30], etc. Nevertheless, product formulas are still the most popular and
promising approaches due to their relatively simple implementations,
especially on near-term devices [34, 11].

Recent works provided refined analysis for standard product formulas [14] and
proposed improved digital quantum simulation algorithms including randomized
compiling [7, 9, 8], random permutation [12], symmetry transformation [39],
and multiproduct formula [28, 43, 21] based on product formulas. While most of
these simulation algorithms are important and efficient for approximating
Hamiltonian evolutions, the resource requirement depends heavily on the norm
of (the terms of) the Hamiltonian. However, a significant class of simulation
tasks in quantum physics do not explore high-energy states, suggesting a
potential to improve the resource demand under physically relevant assumptions
on energy scales. As pointed out by Ref. [35], standard product formulas
require smaller complexity to simulate the dynamics of $k$-local spin models
when the states are restricted within a low-energy subspace. Hitherto, a
general and complete theoretical framework for digital quantum simulations in
the low-energy subspace still requires further exploration.

To this end, we systematically investigate the Hamiltonian simulation problem
in the low-energy subspace. We analyze the Trotter number and gate complexity
of randomized and symmetry-protected digital quantum simulation and prove the
significant improvements compared to the best-known results for simulating
states chosen from the full Hilbert space. We show that similar improvement
can be obtained for quantum models such as the power-law model. In general, we
provide a lower bound for Hamiltonian simulation in the low-energy subspace
with logarithmic dependence on $\Delta$ and linear dependence on $t$. We
further prove that the improvement above persists even when the quantum states
are imperfectly prepared. Our results pave the way to more practical and
efficient digital quantum simulations, which is an essential step toward the
implementation of quantum simulation.

###  1.1 Our results

We begin with some basic concepts for digital quantum simulations. Given a
target Hamiltonian $H=\sum_{l=1}^{L}H_{l}$ consisting of $L$ terms, we
consider simulating an $n$-qubit quantum state $\rho$ undergoing a unitary
evolution $V(t)=e^{-iHt}$ for some time $t$ using a unitary $U$. To achieve
this goal, we split $t$ into $r$ (referred to as the Trotter number or the
step complexity) identical Trotter steps such that the step size satisfies
$\delta=t/r\ll 1$, and perform a product formula
$W({\delta})=e^{-i\delta_{q}H_{l_{q}}}\cdots e^{-i\delta_{1}H_{l_{1}}}$ in
each time step, where $l_{1},\ldots,l_{q}$ are chosen from $1,\ldots,L$. A
standard instance of product formula is the Lie-Trotter formula

$\displaystyle
U(t)=S_{1}(\delta)^{r}=\left(\prod_{l=1}^{L}\exp(-iH_{l}\delta)\right)^{r}.$

Suzuki systematically extended the Lie-Trotter formula to a family of $p$-th
order Suzuki-Trotter formulas defined recursively by

$\displaystyle S_{2}(\delta)$
$\displaystyle=\prod_{l=1}^{L}\exp\Bigl(-iH_{l}\frac{\delta}{2}\Bigr{missing})\prod_{l=L}^{1}\exp\Bigl(-iH_{l}\frac{\delta}{2}\Bigr{missing}),$
$\displaystyle S_{p}(\delta)$
$\displaystyle=S_{p-2}(u_{p}\delta)^{2}S_{p-2}((1-4u_{p})\delta)S_{p-2}(u_{p}\delta)^{2},$

where $u_{p}=(4-4^{1/(p-1)})^{-1}$ [36, 37, 38]. For large enough $r$, we have
$V(t)\approx U(t)=S_{p}(\delta)^{r}$.

The main challenge for implementing product formulas is to decide a proper $r$
to guarantee that the simulation error is bounded below the error threshold
$\epsilon$. In the general case, the simulation error for a unitary simulation
$U$ is given by $U^{\dagger}\rho U-V^{\dagger}\rho V$. For technical
simplicity, we consider the spectral norm $\norm{\cdot}$ as the error metric
in this paper. As $\norm{U^{\dagger}\rho U-V^{\dagger}\rho V}\leq\norm{U-V}$,
we consider the spectral distance between the unitaries $U$ and $V$ in the
following analysis. Given a Trotter number $r$, the Lie-Trotter formula
provides a second-order simulation approximation error
$O(\lambda_{H}^{2}t^{2}/r)$, where $\lambda_{H}=\sum_{l=1}^{L}\norm{H_{l}}\leq
L\lambda_{h}$ and $\lambda_{h}=\max_{l}\norm{H_{l}}$ [14]. Therefore, choosing
the Trotter number $r_{1}=O(t^{2}\lambda_{H}^{2}/\epsilon)$ is sufficient to
ensure the error below $\epsilon$. Similarly, choosing the Trotter number
$r_{p}=O((\lambda_{H}t)^{1+1/p}/\epsilon^{1/p})$ is sufficient to ensure the
error below $\epsilon$ for $p$-th order Lie-Suzuki-Trotter formulas.

It is sometimes possible to improve the complexity of product formulas with
further information or assumptions in the form of Hamiltonians and states. In
this paper, we consider digital quantum simulation in the low-energy subspace.
Without loss of generality, we assume the Hamiltonian $H$ is composed of $L$
positive semi-definite terms $H_{l}\geq 0$.111Otherwise, we can consider
simple shifts $H_{l}\to H_{l}+a_{l}I$ for all $H_{l}$’s for some positive
$a_{l}$ such that $H_{l}+a_{l}I\geq 0$, which only adds an additional global
phase to the evolution $V$. For the simulation task, we only care about the
low-energy quantum states that are below spectrum $\Delta\ll\norm{H}$. We
denote the projector to this subspace as $\Pi_{\leq\Delta}$. The formal
definition for the low-energy simulation is given as follows:

######  Definition 1.

Given a Hamiltonian $H$ with the above assumptions and a quantum state $\rho$
with energy lower than $\Delta\leq\norm{H}$, our goal is to find a simulation
channel $\mathcal{U}(\cdot)$ such that the simulation error is below some
threshold $\epsilon$, i.e., $\norm{\mathcal{U}(\rho)-V^{\dagger}\rho
V}\leq\epsilon$, where $V=e^{-iHt}$ is the target evolution of $H$ with
simulation time $t$.

In the special case when $\mathcal{U}(\cdot)$ in Definition 1 is a unitary
channel for some unitary $U$, we only need to consider the norm
$\norm{(U-V)\Pi_{\leq\Delta}}$ when estimating the upper bound for the
simulation error between the simulation unitary $U$ and the target evolution
$V$. A specific well-studied model is the spin Hamiltonian $H$ which has
$k$-local interactions and strengths for each individual interaction term
bounded below by $J$. Assuming each $H_{l}$ contains at most $M$ interaction
terms and the number of interaction terms acting on any spin bounded by $d>0$,
Ref. [35] showed that choosing a Trotter number
$r^{\prime}=O((t\Delta)^{1+1/p}/\epsilon^{1/p}+(t\sqrt{n})^{1+1/(2p+1)}/\epsilon^{1/(2p+1)})$
for a $p$-th order Suzuki-Trotter formula suffices to achieve $\epsilon$
accuracy, where $n$ is the number of qubits. Compared to the simulation
complexity for Suzuki-Trotter formulas for the general case simulation, we can
obtain an improvement on $n$ dependence for $p=1$ as $\lambda_{H}=O(n)$ for
$k$-local Hamiltonians. Ref. [35] also showed that other improvements are
possible for $p\geq 2$ when choosing the parameters $d$, $M$, and $J$ of the
$H$ within some regimes. This result indicates that standard product formulas
applied to low-energy states have smaller errors than the worst-case error for
$k$-local spin Hamiltonians.

Model | Approach | Unitary simulation | Low-energy simulation
---|---|---|---
$k$-local | PF |  $O\left(\frac{t^{1+1/p}}{\epsilon^{1/p}}n^{1/p}\right)$ [14] |  20cm$\tilde{O}\left(\frac{(t\Delta)^{1+1/p}}{\epsilon^{1/p}}\right)+$
$O\left(\frac{t^{1+1/(2p+1)}}{\epsilon^{1/(2p+1)}}n^{\frac{1}{2}+\frac{1}{4p+2}}\right)$ [35] |  |  |
$k$-local | qDRIFT |  $O\left(\frac{t^{2}}{\epsilon}n^{2}\right)$ [7] | $\tilde{O}\left(\frac{t^{2}}{\epsilon}\Delta^{2}\right)+O\left(\frac{t^{4/3}}{\epsilon^{1/3}}n^{2/3}\right)$
$k$-local | rand. perm. |  $O\left(\frac{L(tMJ)^{1+1/p}}{\epsilon^{1/p}}\right)$ [12]222This paper considered bounds in the diamond norm instead of the spectral norm. However, we can derive the complexity for the spectral norm case directly through their results. |  20cm$\tilde{O}\left(\frac{L(t\Delta)^{1+1/p}}{\epsilon^{1/p}}\right)+$
$O\left(\frac{t^{1+1/(2p+1)}}{\epsilon^{1/(2p+1)}}(LnJ)^{\frac{1}{2}+\frac{1}{4p+2}}\right)$ |  |  |
$k$-local | sym. prot. |  $O\left(\left(\frac{nt^{2}}{\epsilon}\right)^{\frac{1}{1+\theta}}+\frac{nt^{3/2}}{\epsilon^{1/2}}\right)$ [39] |  20cm$\tilde{O}\left(\left(\frac{\Delta t^{2}}{\epsilon}\right)^{\frac{1}{1+\theta}}+\frac{\Delta t^{3/2}}{\epsilon^{1/2}}\right)$+
$O\left(\left(\frac{nt^{3}}{\epsilon}\right)^{\frac{1}{2+\theta}}+\frac{t^{5/4}n^{1/2}}{\epsilon^{1/4}}\right)$ |  |  |
power-law | PF |  $O\left(\frac{(gt)^{1+1/p}}{\epsilon^{1/p}}n^{1/p}\right)$ [14] |  20cm$\tilde{O}\left(\frac{(gt)^{1+1/p}}{\epsilon^{1/p}}\Delta^{1+1/p}\right)+$
$O\left(\frac{t^{1+1/(2p+1)}}{\epsilon^{1/(2p+1)}}(gn)^{\frac{1}{2}+\frac{1}{4p+2}}\right)$ |  |  |

Table 1: A summary of low-energy simulation complexity (Trotter number) for
product formulas (PF), qDRIFT, random permutation (rand. perm.), and symmetry
protection (sym. prot.) approaches, with comparisons to the state-of-the-art
results. The models include $k$-local Hamiltonians and power-law models. The
results in this work are marked in the bold font. Here, $\tilde{O}(\cdot)$
omits the polylogarithmic dependence on the parameters. For $k$-local
Hamiltonians, we assume that $d$, $L$, and $J$ to be $O(1)$ except for the
random permutation approach. For the symmetry protection approach,
$0\leq\theta\leq 1$ is a constant depending on the structure of the
Hamiltonian. In the last line, $\alpha$ is the decay factor for the power-law
interactions. Here, $g$ is the interaction strength and $g=O(1)$, $O(\log n)$,
and $O(n^{1-\alpha/D})$ for $\alpha>D$, $\alpha=D$, and $\alpha<D$, where $D$
is the dimension of the system. For random approaches, we only list the
results considering the simulation errors in expectation. We derive the
Trotter number required for ensuring the random fluctuation in the following
sections. In addition, we provide corresponding gate complexities for the
results in the above table.

In this work, we conduct a systematic study of quantum simulation under the
low-energy subspace assumption by proving the robustness of low-energy
simulation under imperfect state preparation, proving the lower bound for low-
energy simulations, analyzing the complexity required for implementing various
digital quantum simulation approaches [7, 12, 39], and simulating different
models such as $k$-local Hamiltonians and power-law models. We summarize our
results in the following Table 1 with comparisons to the previous results. We
also carry out extensive numerical experiments to benchmark the low-energy
simulation using various algorithms in physically relevant models.

####  1.1.1 Robustness against imperfect state preparations

Concerning experimental implementations of quantum simulations, it is hard to
guarantee that the support for the input states is restricted exactly in the
subspace below the threshold $\Delta$. In particular, the state preparation
procedures on these near-term devices are imperfect, which produces a Gaussian
tail on the higher energy states in the energy spectrum. Quantitatively, the
probability distribution for eigenstates of energy $E$ higher than $\Delta$
satisfies a Gaussian distribution $\propto e^{-(E-\Delta)^{2}/\sigma^{2}}$,
where $\sigma\ll\Delta$ is the variance. In this work, we prove that such
imperfect state preparation causes an $O((\sigma/\Delta)^{2})$ multiplicative
extra error in Section 2.1, which is ignorable compared to the original low-
energy simulation error. This result indicates that the low-energy improvement
persists even under imperfect state preparations.

####  1.1.2 Applications

##### Randomized quantum compiling (qDRIFT).

The first approach we analyze is the qDRIFT algorithm proposed by Ref. [7]. We
recall that the first-order Lie-Trotter algorithm deterministically cycles
through every term $H_{l}$ in the Hamiltonian $H$ in each step. The qDRIFT
algorithm [7], however, approximates the target evolution using a channel
constructed by averaging products $U=U_{r}U_{r-1}\cdots U_{1}$. Here, each
unitary $U_{i}$ is a short-time evolution on a single term $H_{l}$ sampled
according to the probability distribution
$\\{p_{l}=\norm{H_{l}}/\lambda_{H}\\}_{l}$. It is proved that only
$O(\lambda_{H}^{2}t^{2}/\epsilon)$ steps suffice to ensure the error within
$\epsilon$ in expectation, which is explicitly independent of $L$. In this
work, we consider the performance of qDRIFT in the low-energy subspace. We
obtain the following complexity upper bound.

######  Theorem 1.

Let $H=\sum_{l=1}^{L}H_{l}$ be an $n$-qubit $k$-local Hamiltonian with
parameters $M$, $J$, and $d$ defined as above ($LMJ=O(n)$ and $dkJ=O(n)$). By
choosing the Trotter number

$\displaystyle
r_{\text{exp}}=\tilde{O}\left(\frac{L^{2}(\Delta+dkJ)^{2}t^{2}}{\epsilon}\right)+O\left(\frac{(LM)^{2/3}J^{4/3}t^{4/3}}{\epsilon^{1/3}}\right),$
(1)

we can ensure that the expected simulation error in the low-energy subspace
below $\Delta$ for the qDRIFT algorithm is bounded by $\epsilon$, i.e.,
$\norm{(V-\mathbb{E}[U])\Pi_{\leq\Delta}}\leq\epsilon$. Moreover, if we choose
the Trotter number

$\displaystyle
r_{\text{prob}}=\tilde{O}\left(\frac{L^{2}(\Delta+dkJ)^{2}t^{2}}{\epsilon^{2}}\left(n+\log(\frac{1}{\chi})\right)\right)+O\left(\frac{(LM)^{2/3}J^{4/3}t^{4/3}}{\epsilon^{2/3}}\left(n+\log(\frac{1}{\chi})\right)^{1/3}\right),$
(2)

then with probability at least $1-\chi$, we can guarantee that the simulation
error for a single product $U$ is bounded by $\epsilon$.

We leave the full proof for Theorem 1 in Section 3.1. Theorem 1 provides the
step complexity (Trotter number) required for simulating low-energy states of
$k$-local spin Hamiltonians both in the expectation case and in the case when
one needs to ensure the performance of every single simulation unitary $U$
with high confidence. Compared to the step complexity for implementing the
qDRIFT algorithm in the full Hilbert space, the above theorem indicates that
we can expect an improvement concerning the dependence on system size $n$ as
long as the strengths for each individual interaction term being $J=o(n)$,
which can be satisfied by most of the spin models. Although the step
complexity (1) in expectation and the step complexity (2) to control random
fluctuations are no better than that for the low-energy simulation using
standard product formulas, we remark that this deficiency originates from the
disadvantage of qDRIFT in simulating $k$-local spins even in the full Hilbert
space [9].

##### Random permutation.

Except for the qDRIFT algorithm, we consider the random permutation approach
in simulating $k$-local Hamiltonians [12]. While the qDRIFT approach
implements random evolution on each term of the Hamiltonian, the random
permutation approach randomly permutes the sequence of the terms in higher-
order Suzuki-Trotter formulas. In particular, for any permutation in the
permutation group $\sigma\in\mathcal{S}_{L}$, we denote

$\displaystyle S_{2}^{\sigma}(\delta)$
$\displaystyle=\prod_{l=1}^{L}\exp\Bigl(-iH_{\sigma(l)}\frac{\delta}{2}\Bigr{missing})\prod_{l=L}^{1}\exp\Bigl(-iH_{\sigma(l)}\frac{\delta}{2}\Bigr{missing}),$
$\displaystyle S_{p}^{\sigma}(\delta)$
$\displaystyle=S_{p-2}^{\sigma}(u_{p}\delta)^{2}S_{p-2}^{\sigma}((1-4u_{p})\delta)S_{p-2}^{\sigma}(u_{p}\delta)^{2},$

where $u_{p}=(4-4^{1/(p-1)})^{-1}$. In each time step $\delta$, we consider
averaging over all possible $\sigma\in\mathcal{S}_{L}$. For the simulation in
low-energy subspace, we estimate the simulation error using

$\displaystyle\norm{\left[V(t)-\left(\frac{1}{L!}\sum_{\sigma\in\mathcal{S}_{L}}S_{p}^{\sigma}(\delta)\right)^{r}\right]\Pi_{\leq\Delta}}.$
(3)

In practice, it is difficult to implement all the permutations $\sigma$ and
average the results as the number of all permutations $L!$ increases
exponentially with $L$. We consider a random sampling procedure similar to the
qDRIFT algorithm. We sample by implementing a random evolution
$U=U_{r}U_{r-1}\cdots U_{1}$ with each $U_{i}=S_{p}^{\sigma_{i}}(\delta)$ a
randomly chosen permuted Suzuki-Trotter formula and
$\sigma_{i}\sim\mathcal{S}_{L}$ i.i.d. random permutations. In expectation,
the evolution in each step is
$\frac{1}{L!}\sum_{\sigma\in\mathcal{S}_{L}}S_{p}^{\sigma}(\delta)$, which is
the same as the standard random permutation approach. Concerning implementing
this approach in the low-energy subspace, we obtain the following complexity
result.

######  Theorem 2.

Let $H=\sum_{l=1}^{L}H_{l}$ be an $n$-qubit $k$-local Hamiltonian with
parameters $M$, $J$, and $d$ ($LMJ=O(n)$ and $dkJ=O(n)$). By choosing the
Trotter number

$\displaystyle
r_{\text{exp}}=\tilde{O}\left(\frac{Lt^{1+1/p}(\Delta+dkJ)^{1+1/p}}{\epsilon^{1/p}}\right)+O\left(\frac{Lt^{1+1/(2p+1)}(dkMJ^{2})^{\frac{1}{2}+\frac{1}{4p+2}}}{\epsilon^{1/(2p+1)}}\right),$
(4)

we can ensure that the expected simulation error in the low-energy subspace
below $\Delta$ for the random permutation approach is bounded by $\epsilon$.
For the random sampling implementation of the random permutation algorithm, if
we choose the Trotter number

$\displaystyle
r_{\text{prob}}=\tilde{O}\left(\frac{(Lt(\Delta+dkJ))^{\frac{2p+2}{2p+1}}\left(n+\log\left(\frac{1}{\chi}\right)\right)^{\frac{1}{2p+1}})}{\epsilon^{\frac{2}{2p+1}}}\right)+O\left(\frac{(L^{2}t^{2}dkMJ^{2})^{\frac{2p+2}{4p+3}}\left(n+\log\left(\frac{1}{\chi}\right)\right)^{\frac{1}{4p+3}})}{\epsilon^{\frac{2}{4p+3}}}\right),$
(5)

then with probability at least $1-\chi$, we can ensure that the simulation
error for a single random sequence is bounded by $\epsilon$.

We leave the technical details for Theorem 2 in Section 3.1. Theorem 2 shows
the improvement over the step complexity of random permutation in the full
Hilbert space [12] concerning the dependence on systems size $n$. Compared to
the step complexity in (4) with that for standard product formulas in the low-
energy subspace [35], we can obtain an $L^{1/2p}$ reduction in the step
complexity. This originates from the reduction from implementing the random
permutation in the Hilbert space [12].

Based on the random permutation approach, Ref. [16] proposed an approach to
double the order of approximation via the randomized product formulas. In
particular, given $r\geq(5^{p/2-1}+\frac{1}{2})L\lambda_{h}t$ and $p$-th order
contributing product formula, this algorithm provides a $(2p+1)$-th order
simulation for the ideal evolution. For this approach, we also prove the
following corollary on the Trotter numbers required for low-energy
simulations.

######  Corollary 1.

Let $H=\sum_{l=1}^{L}H_{l}$ be an $n$-qubit $k$-local Hamiltonian with
parameters $M$, $J$, and $d$ ($LMJ=O(n)$ and $dkJ=O(n)$). By choosing the
Trotter number

$\displaystyle
r_{\text{exp}}=\tilde{O}\left(\frac{(Lt)^{1+1/(2p+1)}(\Delta+dkJ)^{1+1/(2p+1)}}{\epsilon^{1/(2p+1)}}\right)+O\left(\frac{(Lt)^{1+1/(4p+3)}(dkMJ^{2})^{\frac{2p+2}{4p+3}}}{\epsilon^{1/(4p+3)}}\right),$
(6)

we can ensure that the expected simulation error in the low-energy subspace
below $\Delta$ for the doubling order approach is bounded by $\epsilon$. And
the Trotter complexity $r_{\text{prob}}$ to ensure the algorithm converges
remains the same as that in Theorem 2.

##### Symmetry protection.

Other than the randomized digital quantum simulation approaches, another
approach to improve the performance of digital quantum simulations is to
implement a symmetry transformation in each time step [39].

Given a Hamiltonian $H$, we assume that it is invariant under a group of
unitary transformations denoted by $\mathcal{C}$. For each unitary $C$ chosen
from $\mathcal{C}$, we explicitly have $[H,C]\equiv 0$. The group
$\mathcal{C}$ represents a symmetry group of the system. According to Ref.
[39], we “rotate” each implementation of the Lie-Trotter formula
$S_{1}(\delta)$ via a symmetry transformation $C_{\mu}\in\mathcal{C}$ in each
time step. The simulation for $V=e^{-iHt}$ reads
$V(t)\approx\prod_{\mu=1}^{r}C_{\mu}^{\dagger}S_{1}(\delta)C_{\mu}$. Assuming
the simulation error coherent, the digital quantum simulation of an evolution
$V(t)$ may end up with the time evolution of a different Hamiltonian
$H_{\text{eff}}$. We denote $\kappa=H_{\text{eff}}-H$ and rewrite the
simulation under the assumption $\norm{\kappa}$ is small

$\displaystyle
V(t)\approx\prod_{\mu=1}^{r}C_{\mu}^{\dagger}S_{1}(\delta)C_{\mu}=\prod_{\mu=1}^{r}e^{-iC_{\mu}^{\dagger}H_{\text{eff}}C_{\mu}\delta}=\prod_{\mu=1}^{r}e^{-i(H+C_{\mu}^{\dagger}\kappa
C_{k})\delta}\approx
e^{-i\left(H+\frac{1}{r}\sum_{\mu=1}^{r}C_{\mu}^{\dagger}\kappa
C_{\mu}\right)t}.$

In the last line, we use the Baker-Campbell-Hausdorff (BCH) formula. The
simulation error for this symmetry protection approach can thus be represented
as

$\displaystyle\norm{V(t)-\prod_{\mu=1}^{L}C_{\mu}^{\dagger}S_{1}(\delta)C_{\mu}}.$

For the symmetry protection approach, the Trotter number required for the
$k$-local spin model is provided by Ref. [39]

$\displaystyle
r=\max\left\\{O\left(\left(\frac{nt^{2}}{\epsilon}\right)^{\frac{1}{1+\theta}}\right),O\left(\frac{nt^{3/2}}{\sqrt{\epsilon}}\right)\right\\},$

where $\theta\in[0,1]$ is a constant that depends only on the structure of the
Hamiltonian $H$ and the properties of the symmetry group $\mathcal{C}$. There
are several examples provided with a specific value of $\theta$. For example,
when we draw symmetry transformations randomly, the behavior of the error
would be analogous to a random walk, which results in $\theta=0.5$ under
specific settings. The rigorous proof for this intuition for the localized
Heisenberg model was provided in Ref. [39]. Following a similar derivation in
Ref. [21], one can explicitly obtain a construction where $\theta$ achieves
the maximal value of $1$.

Now, we consider the performance of this symmetry protection approach in the
low-energy subspace, which is also an open question in Ref. [39]. We estimate
the simulation error by computing the following error term

$\displaystyle\norm{\left(V(t)-\prod_{\mu=1}^{r}C_{\mu}^{\dagger}S_{1}(\delta)C_{\mu}\right)\Pi_{\leq\Delta}}.$

We prove the following theorem concerning the low-energy simulations of
$k$-local Hamiltonians using the symmetry protection approach.

######  Theorem 3.

Let $H=\sum_{l=1}^{L}H_{l}$ be an $n$-qubit $k$-local Hamiltonian with
parameters $M$, $J$, and $d$ ($LMJ=O(n)$ and $dkJ=O(n)$). By choosing the
Trotter number

$\displaystyle
r=\tilde{O}\left(\left(\frac{L(\Delta+JdkL)t^{2}}{\epsilon}\right)^{\frac{1}{1+\theta}}+\left(\frac{L^{2}MJdkt^{3}}{\epsilon}\right)^{\frac{1}{2+\theta}}+\frac{L(\Delta+JdkL)t^{\frac{3}{2}}}{\epsilon^{\frac{1}{2}}}+\frac{t^{\frac{5}{4}}L(JdkM)^{\frac{1}{2}}}{\epsilon^{\frac{1}{4}}}\right),$
(7)

we can ensure that the symmetry protection approach provides an approximation
within error $\epsilon$ in the low-energy subspace below $\Delta$.

The proof for Theorem 3 is provided in Section 3.2. The above theorem provides
the step complexity for the symmetry protection approach in the low-energy
simulation. Compared with the complexity required for full Hilbert space
simulations [39], the $\lambda_{H}=O(n)$ dependence is replaced by $L\Delta$
in two terms and $O(n^{1/2})$ in the other two terms, which enables an
improvement concerning the $n$ dependence. When $\theta>0$, the required
Trotter number is smaller than that for implementing the standard first-order
Lie-Trotter product formula in the low-energy subspace.

##### Power-law model.

Theorem 1, Theorem 2, and Theorem 3 consider different digital quantum
simulation approaches in simulating the low-energy quantum states of $k$-local
spin models. We also consider power-law models, which can be regarded as a
specific extension of $k$-local spin models at $k=2$, unbounded $d$, and
additional constraints on the interaction strength $J$ for each term. We
consider a $D$-dimensional power-law interaction model consisting of $n$
qubits. A power-law interaction
$H=\sum_{{\bm{i}},{\bm{j}}\in\Sigma}H_{{\bm{i}},{\bm{j}}}$ with exponent
satisfies:

$\displaystyle\norm{H_{{\bm{i}},{\bm{j}}}}\leq\begin{cases}1&\quad\text{if
}{\bm{i}}={\bm{j}},\\\
\frac{1}{\norm{{\bm{i}}-{\bm{j}}}_{2}^{\alpha}}&\quad\text{if
}{\bm{i}}\neq{\bm{j}},\end{cases}$

where ${\bm{i}},{\bm{j}}\in\Sigma$ are the qubit sites,
$\Lambda\subseteq\mathbb{R}^{D}$ is a $D$-dimensional square lattice, and
$H_{{\bm{i}},{\bm{j}}}$ is an operator supported on two sites
${\bm{i}},{\bm{j}}$. Some notable examples include Coulomb interactions
($\alpha=1$), dipole-dipole interactions ($\alpha=3$), and van der Waals
interactions ($\alpha=6$). According to Lemma H.1 of Ref. [14], we have

$\displaystyle
g=\norm{\max_{{\bm{i}}}\sum_{{\bm{j}}\neq{\bm{i}}}H_{{\bm{i}},{\bm{j}}}}\leq\begin{cases}O(n^{1-\alpha/D})&\quad\text{for
}0\leq\alpha<D,\\\ O(\log(n))&\quad\text{for }\alpha=D,\\\ O(1)&\quad\text{for
}\alpha>D.\\\ \end{cases}$

Here, $g$ is an upper bound of the strengths of the interactions associated
with a single spin qubit. We can observe that every term of the power-law
Hamiltonian is $2$-local. In the following, we assume that each part $H_{l}$
has strength at most $M$. We obtain the following complexity result for
simulating such a power-law model in the low-energy subspace.

######  Theorem 4.

Consider simulating a $D$-dimensional power-law Hamiltonian
$H=\sum_{l=1}^{L}H_{l}$ with interaction strength $g$ in the low-energy
subspace $\Delta$. Assume that each $H_{l}$ has strength at most $M$. Then,
the Trotter number required to ensure the simulation error below $\epsilon$ is

$\displaystyle
r=\tilde{O}\left(\frac{t^{1+\frac{1}{p}}L^{\frac{1}{p}}}{\epsilon^{\frac{1}{p}}}(\Delta+gq\log
q)^{1+\frac{1}{p}}\right)+O\left(\frac{t^{1+\frac{1}{2p+1}}}{\epsilon^{\frac{1}{2p+1}}}(LMg)^{\frac{1}{2}+\frac{1}{4p+2}}\right).$
(8)

The corresponding gate complexity is given as:

$\displaystyle
G=\tilde{O}\left(\frac{t^{1+\frac{1}{p}}L^{\frac{1}{p}+1}}{\epsilon^{\frac{1}{p}}}(\Delta+gq\log
q)^{1+\frac{1}{p}}\right)+O\left(\frac{Lt^{1+\frac{1}{2p+1}}}{\epsilon^{\frac{1}{2p+1}}}(LMg)^{\frac{1}{2}+\frac{1}{4p+2}}\right).$

The proof for the above theorem is provided in Section 3.3. Compared with the
best previously known complexity result for simulating the power-law model
[14], the step complexity in (8) shows an improvement concerning the $n$
dependence at $p=1$. Notice that the best result for full Hilbert space
simulations is obtained by simulating each term $H_{{\bm{i}},{\bm{j}}}$
separately using one term $H_{l}$. For $p>1$, we can also obtain some
advantage for specific choices of parameters $M$ and $\Delta$ in this case as
we assume that $H_{l}$’s consist of more than one term of
$H_{{\bm{i}},{\bm{j}}}$.

####  1.1.3 Lower bound for low-energy simulations

We now consider the necessary number of queries to the Hamiltonian $H$ to
simulate any given state $\rho$ in the low-energy subspace $\leq\Delta$. In
Section 4, we show the following theorem as the lower bound on the quantum
query complexity of simulations in the low-energy subspace:

######  Theorem 5 (Informal, see Theorem 6 for the formal version).

There exists a Hamiltonian $H$ such that simulating some state with constant
error within the low-energy threshold $\leq\Delta$ for some chosen scaled time
$\tau$ requires at least
$\Omega\left(\max\\{\tau,\frac{\log(\Delta)}{\log\log(\Delta)}\\}\right)$
queries to $H$.

The above theorem indicates that, similar to full space simulations [4],
simulating Hamiltonians in the low-energy subspace also requires a linear
number of queries to $H$ in the simulation time $\tau$. In addition, we also
obtain a logarithmic dependence on the threshold $\Delta$. Although this is
far from the upper bound, where we require $\operatorname{poly}(\Delta)$
queries, we remark that this gap also exists in the full space simulation [5].
The main idea to prove this theorem is to show a Hamiltonian of which exact
low-energy subspace simulations for any time $\tau>0$ enable one to compute
the parity of a string, which requires a linear number of queries to the
length of the string [3, 17].

###  1.2 Open questions

Our paper leaves several open questions for future investigations:

  * •

Ref. [42] developed a theory of average error for Hamiltonian simulation with
random inputs and showed that improvement is possible when considering the
average-case instead of the worst-case error. It is interesting to explore if
we can obtain further improvement if we consider the simulation error for low-
energy random states.

  * •

In this paper, we study the simulation of evolution $V$ on a low-energy
quantum state $\rho$. It is also meaningful to study the complexity for
estimating the expected value $\operatorname{Tr}(OV^{\dagger}\rho V)$ for some
observable $O$ instead of the full state $V^{\dagger}\rho V$ for a low-energy
state $\rho$.

  * •

Can we propose some Hamiltonians that are widely considered in the field of
quantum information such that we can obtain improvement in the step complexity
of digital quantum simulation under the low-energy assumption?

##### Roadmap.

The rest of the paper is organized as follows. In Section 2, we provide the
general framework for low-energy simulations. In particular, we recap some of
the results in [35] and prove several important lemmas for our results. We
clarify the difference between our work and previous works. We further prove
the robustness of the improvement originates from the low-energy assumption in
imperfect state preparation. In Section 3, we consider the applications of
low-energy settings for different digital quantum simulation approaches and
quantum models. In Section 3.1, we consider randomized approaches including
qDRIFT and random permutation for simulating $k$-local Hamiltonians, and
provide the proofs for Theorem 1, Theorem 2, and Corollary 1. In Section 3.2,
we consider the symmetry protection approach and provide the proof for Theorem
3. We show that the low-energy setting can provide improvements in power-law
Hamiltonians in Section 3.3 and prove Theorem 4. Finally, we prove the query
lower bound in Theorem 5 for simulating dynamics in the low-energy subspace in
Section 4.

##  2 The General Framework

##### Notations.

We consider simulating the Hamiltonian $H=\sum_{l=1}^{L}H_{l}$ composed of $n$
particles. We assume that each term $H_{l}$ contains at most $M$ interaction
terms and each interaction term acts on at most $k$ qubits. For each qubit, we
assume that the strength of interaction between this qubit and the rest qubits
is bounded by $g$ measured by the spectral norm. Without loss of generality,
we assume that $H_{l}\geq 0$ for any $l$. We also assume the number of terms
acting on each qubit is bounded by $d$ and the strength of each term is
bounded by $J$. Thus, we have $g\leq dJ$.

Throughout this paper, we use the notation $\tilde{O}(\cdot)$ which omits the
polylogarithmic dependence on the parameters. Given positive constant
$\Lambda$, we denote $\Pi_{\leq\Lambda}$ the projector to the subspace spanned
by eigenstates of $H$ with energy smaller than $\Lambda$. We further denote
$\Pi_{>\Lambda}=I-\Pi_{\leq\Lambda}$ be the projector to the orthogonal
subspace.

##### Simulating $k$-local Hamiltonians in the low-energy subspace.

For $k$-local Hamiltonians with bounded interaction strengths, the following
lemma holds as the backbone of the framework for analyzing the simulation
performance in the low-energy subspace.

######  Lemma 1 (Theorem 2.1 of [1]).

Given $H=\sum_{l=1}^{L}H_{l}$ defined above with parameters $g$, $L$, and $k$,
and any operator $A$, the following inequality holds

$\displaystyle\norm{\Pi_{>\Lambda^{\prime}}A\Pi_{\leq\Lambda}}\leq\norm{A}\cdot
e^{-\lambda(\Lambda-\Lambda-2R_{A})},$

where $R_{A}$ is the strengths of interaction terms acting on $A$ and
$\lambda=(2gk)^{-1}$. Here, $\Lambda^{\prime}\geq\Lambda\geq 0$ are two
positive values.

Based on Lemma 1, we can obtain some corollaries when $A$ is an evolution of
some terms $H_{l}$. We list these lemmas in Appendix A. Using these lemmas,
Ref. [35] obtained the following result concerning the performance of
arbitrary product formulas of length $q$ in the low-energy subspace. Given
$H=\sum_{l=1}^{L}H_{l}$ an $n$-qubit $k$-local Hamiltonian with parameters
$M$, $J$, and $d$ ($LMJ=O(n)$ and $dkJ=O(n)$), the simulation error between
the target evolution $V(\delta)$ and a product formula $W(\delta)$ of $p$-th
order error is bounded by [35]:

$\displaystyle\norm{(V(\delta)-W(\delta))\Pi_{\leq\Delta}}=O((L\Delta^{\prime}\delta)^{p+1}),$
(9)

where

$\displaystyle\Delta^{\prime}=\Delta+\frac{\alpha
LM\delta}{\lambda}+\frac{q\log
q}{\lambda}+\frac{q}{\lambda}\log\left(\frac{1}{J\delta}\right),$ (10)

$\alpha$ is some constant, and $\lambda=(2Jdk)^{-1}$. We can further deduce
the following lemma regarding the Trotter number required for low-energy
simulations based on the above error decomposition.

######  Lemma 2 (Eq. (111) of [35]).

Let $H=\sum_{l=1}^{L}H_{l}$ be an $n$-qubit $k$-local Hamiltonian with
parameters $M$, $J$, and $d$ ($LMJ=O(n)$ and $dkJ=O(n)$). By choosing the
Trotter number

$\displaystyle
r=\frac{t}{\delta}=\tilde{O}\left(\frac{t^{1+\frac{1}{p}}}{\epsilon^{\frac{1}{p}}}(L\Delta+LdkJq\log
q)^{1+\frac{1}{p}}\right)+O\left(\frac{t^{1+\frac{1}{2p+1}}}{\epsilon^{\frac{1}{2p+1}}}(L^{2}dMJ^{2})^{\frac{1}{2}+\frac{1}{4p+2}}\right),$
(11)

we can ensure that the $p$-th order product formula $W(\delta)$ provides an
approximation within error $\epsilon$ in the low-energy subspace below
$\Delta$.

Furthermore, as the implementation of the $p$-th order Suzuki-Trotter formula
requires $O(2\cdot 5^{p/2-1}L)$ gates, the gate complexity required is

$\displaystyle
G=\tilde{O}\left(\frac{Lt^{1+\frac{1}{p}}}{\epsilon^{\frac{1}{p}}}(L\Delta+LdkJq\log
q)^{1+\frac{1}{p}}\right)+O\left(\frac{Lt^{1+\frac{1}{2p+1}}}{\epsilon^{\frac{1}{2p+1}}}(L^{2}dMJ^{2})^{\frac{1}{2}+\frac{1}{4p+2}}\right),$

where the $p$ dependence is hidden as we assume it to be a constant. Notice
that $LMJ=O(n)$, the step complexity has an explicit
$n^{\frac{1}{2}\left(1+\frac{1}{2p+1}\right)}$ term. For $p=1$, the step
complexity obtained by Lemma 2 achieves a reduction concerning the dependence
on $n$ compared to the performance for simulating high energy states. For
general $p\geq 1$ and $d=O(n^{\frac{1}{2(p+1)})}$, we can also obtained a
similar improvement.

As an illustration for the reduction of simulating dynamics of Hamiltonian
assuming input states from the low-energy subspace, we perform numerical
experiments to benchmark the product formulas in simulating $k$-local
Hamiltonians. In particular, we study the following homogeneous Heisenberg
model without external fields:

$H=-\sum_{<i,j>}X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j},$ (12)

where $X_{i},Y_{i}$, and $Z_{i}$ are Pauli operators acting on the $i$-th
spin, and the summation is over every adjacent spin pair. Since the ground
energy $E_{0l}<0$ for every $H_{l}$, we implement the shift $H_{l}\rightarrow
H_{l}-E_{0l}I$ for every $H_{l}$. We empirically pick the threshold $\Delta=4$
to make sure that the corresponding low-energy subspace is neither (close to)
the full Hilbert space nor an empty subspace.

Figure 1: Quantum simulation in full Hilbert space vs. low-energy subspace:
Trotter error induced by product formulas for a $\mathbf{2\times 6}$
homogeneous Heisenberg spin-$\mathbf{1/2}$ model. The solid lines denote
simulation errors in the full Hilbert space while the dashed lines denote
errors in the low-energy subspace. The plot distinguishes different Trotter
numbers $r$ by the color of the lines.

As shown in the left-hand side of Figure 1, the Hamiltonian in (12) is
decomposed into three interaction terms: $H=H_{1}+H_{2}+H_{3}$, respectively
represented by yellow, green, and red bonds. Here, a bond represents all the
interactions between a pair of spins. In the right-hand side of Figure 1, we
plot the simulation $\epsilon$ error as a function of single-step evolution
time $\delta$ with respect to different Trotter numbers $r$. The time
evolution operator $V(\delta)$ is approximated by the second-order product
formula $S_{2}(\delta)$. We can observe that the projection to the low-energy
subspace significantly reduces the error compared to the full Hilbert space
for the same parameter configuration.

###  2.1 Low-energy simulation with imperfect state preparations

In the previous analysis, we focused on simulating the dynamics of low-energy
input states for a given Hamiltonian. We assumed the quantum state is chosen
in the low-energy subspace under a threshold $\Delta$. However, the state
preparation on near-term devices is not perfect, resulting in a Gaussian tail
distribution over higher energy states on the energy spectrum. It is natural
to ask if the improvement obtained by low-energy simulations persists under
the imperfect state preparation scenario. In the following, we provide an
affirmative answer to this question.

For technical simplicity, we focus on pure states and suppose the perfect
input state $\ket{\psi_{p}}$ has energy $\Delta$ and denote the actual input
state to be $\ket{\psi}$. In general, the mean energy and variance of $H$ in
the state $\ket{\psi}$ is

$\displaystyle E_{\psi}$ $\displaystyle=\bra{\psi}H\ket{\psi},$
$\displaystyle\sigma_{\psi}^{2}$
$\displaystyle=\bra{\psi}(H-\Delta)^{2}\ket{\psi}.$

In the case when $H$ is local and $\ket{\psi}$ has finite correlation length,
$\sigma_{\psi}^{2}$ will scale as $O(n)$ [2]. Assuming imperfect state
preparation, we adapt $E_{\psi}=\Delta$ and $\sigma_{\psi}^{2}\ll\Delta$ from
current analog state preparation schemes [19, 29]. Quantitatively, the support
of the imperfect input state is a Gaussian distribution with parameter
$\sigma_{\psi}$ around $\Delta$. By using the error representation in (9) and
(10), we can compute the simulation error $\text{Err}_{\text{imp}}$ as:

$\displaystyle\text{Err}_{\text{imp}}$
$\displaystyle\propto\int_{-\infty}^{\infty}e^{-(x-\Delta)^{2}/\sigma_{\psi}^{2}}(L\Delta_{x}\delta)^{p+1}dx,$
(13)

where

$\displaystyle\Delta_{x}=x+\frac{\alpha LM\delta}{\lambda}+\frac{q\log
q}{\lambda}+\frac{q}{\lambda}\log\left(\frac{1}{J\delta}\right).$

Recall that for standard Suzuki-Trotter formulas, $p$ is an even number.
Therefore, we can compute $\text{Err}_{\text{imp}}$ as

$\displaystyle\text{Err}_{\text{imp}}\propto(L\delta\sigma_{\psi})^{p+1}\left[\left(\frac{\Delta^{\prime}}{\sigma_{\psi}}\right)^{(p+1)}+\sum_{i=1}^{p/2}\binom{p+1}{2i}\left(\frac{\Delta^{\prime}}{\sigma_{\psi}}\right)^{(p+1-2i)}\cdot\prod_{j=1}^{i}\frac{1}{2(2j-1)}\right].$

As $\sigma_{\psi}\ll\Delta\leq\Delta^{\prime}$ and $p$ a constant, we can
discard the lower-order terms for $\Delta^{\prime}$ and obtain

$\displaystyle\text{Err}_{\text{imp}}\simeq(L\delta\Delta^{\prime})^{p+1}\left(1+\frac{p(p+1)}{4}\left(\frac{\sigma_{\psi}}{\Delta^{\prime}}\right)^{2}\right).$

We compare the above error decomposition with (9) and observe that we only get
an additional $O(\Delta^{{}^{\prime}(p-1)}\sigma_{\psi}^{2})$ term when the
state preparation is imperfect. Therefore, we conclude that for imperfect
state preparation, the improvement obtained by the low-energy simulation
assumption persists. We further claim that this robustness persists for all
the analyses on any following low-energy performance for various digital
quantum simulation approaches. This is because the low-energy simulation error
for any following algorithms has a polynomial scaling on $\Delta^{\prime}$. We
can thus employ similar calculations as (13) and prove the robustness.

##  3 Applications

Intuitively, the resource requirement for different simulations reduces when
the evolved quantum system does not consider high-energy states. Yet, Lemma 2
only provides strict proof for standard product formulas and $k$-local
Hamiltonians. In this work, we explore other digital quantum simulation
algorithms and Hamiltonians. In this section, we provide the technical details
for our results in applying low-energy analysis for different digital quantum
simulation algorithms and Hamiltonians.

###  3.1 Randomized product formulas

In this subsection, we analyze the performance of randomized simulation
approaches in the low-energy subspace. In particular, we provide the proof for
Theorem 1, Theorem 2, and Corollary 1.

##### qDRIFT algorithm.

We begin with recapping the qDRIFT algorithm in the below Algorithm 1.

0: Hamiltonian $H=\sum_{l=1}^{L}H_{l}$, parameter
$\lambda_{H}=\sum_{l}\norm{H_{l}}$, evolution time $t$, and number of steps
$r$

1: At each time interval $t/r$: evolve a random term in Hamiltonian
$U_{i}=e^{-i(t/r)X_{i}}$, where $X_{i}$ is randomly chosen as
$\lambda_{H}H_{l}/\norm{H_{l}}$ with probability $\norm{H_{l}}/\lambda_{H}$

1: The unstructured (randomly generated) product formula $U=U_{r}\ldots U_{1}$

Algorithm 1 The qDRIFT algorithm

We now analyze the performance of the qDRIFT algorithm in the low-energy
subspace and prove Theorem 1. We still consider the $k$-local Hamiltonians and
parameters $d$, $M$, $J$, and $d$ defined above. We carefully pick
$\Delta^{\prime}\geq\Delta$ a value to be determined later. We denote
$\overline{H}=\Pi_{\leq\Delta^{\prime}}H\Pi_{\leq\Delta^{\prime}}$ and
$\overline{H}_{l}=\Pi_{\leq\Delta^{\prime}}H_{l}\Pi_{\leq\Delta^{\prime}}$.
Therefore, $\norm{\overline{H}_{l}}\leq\norm{\overline{H}}=\Delta^{\prime}$.

Given a random sequence $U_{r}U_{r-1}\ldots U_{1}$ generated by Algorithm 1,
our goal is to estimate the error

$\displaystyle\norm{(U_{r}\ldots U_{1}-V)\Pi_{\leq\Delta}},$

where $V=e^{-iHt}$ is the ideal evolution. To begin with, we decompose this
error into three terms:

$\displaystyle\quad\ \norm{(U_{r}\ldots U_{1}-V)\Pi_{\leq\Delta}}$
$\displaystyle\leq\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}+\norm{(\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r})\Pi_{\leq\Delta}}+\norm{(\mathbb{E}[\overline{U}_{i}]^{r}-V)\Pi_{\leq\Delta}}$
$\displaystyle\leq\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}+\norm{\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r}}+\norm{\mathbb{E}[\overline{U}_{i}]^{r}-\overline{V}},$

where $\overline{V}=e^{-i\overline{H}t}$ and $\overline{U}_{i}$ are obtained
by running Algorithm 1 on $\overline{H}$ with corresponding
$\overline{H}_{l}$’s, and choosing
$\overline{X}_{i}=\lambda_{\overline{H}}\overline{H}_{l}/\norm{\overline{H}_{l}}$
in the $i$-th step with probability
$\norm{\overline{H}_{l}}/\lambda_{\overline{H}}$. Here, the first term is the
projection error for the low-energy subspace estimation. The second term is
the random fluctuation in the low-energy subspace. The third term is the
deterministic bias in the low-energy subspace.

We first consider the deterministic bias error term. For the projected
Hamiltonian $\overline{H}$,
$\lambda_{\overline{H}}=\sum_{l=1}^{L}\norm{\overline{H}_{l}}=O(L\Delta^{\prime})$
for some $\Delta^{\prime}\geq\Delta$ to be fixed later. It is worthwhile to
mention that for $k$-local Hamiltonians, $\lambda_{\overline{H}}$ is
explicitly independent of system size $n$ while $\lambda_{H}=O(n)$ [14].
According to Lemma 3.5 of Ref. [9], for Hamiltonian
$\overline{H}=\sum_{l=1}^{L}$ and $\overline{U}_{i}$, the deterministic bias
term is bounded by

$\displaystyle\norm{\mathbb{E}[\overline{U}_{i}]^{r}-\overline{V}}\leq
r\norm{\mathbb{E}[\overline{U}_{i}]-\overline{V}^{1/r}}\leq
O\left(r\cdot\frac{\lambda^{2}_{\overline{H}}t^{2}}{r^{2}}\right)=O\left(\frac{\lambda^{2}_{\overline{H}}t^{2}}{r}\right).$

Here $\overline{V}^{1/r}=e^{-iHt/r}=e^{iH\delta}$.

The next step is to provide a bound for the random fluctuation term
$\norm{\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r}}$.
For the randomized sequence generated by the qDRIFT algorithm projected into
the low-energy subspace $\overline{U}_{r}\ldots\overline{U}_{1}$, we denote
$B_{k}=(\mathbb{E}[\overline{U}_{i}])^{r-k}\overline{U}_{k}\cdots\overline{U}_{1}$.
Therefore, $B_{r}=\overline{U}_{r}\ldots\overline{U}_{1}$ and
$B_{0}=(\mathbb{E}[\overline{U}_{i}])^{r}$. We can verify that this sequence
$\\{B_{k}:k=0,\ldots,r\\}$ satisfies the martingale properties as follows:

  * •

Causality:  every $B_{k}$ completely depends on the previous information of
$B_{k-1},\ldots,B_{1}$, i.e., the choice of
$\overline{U}_{k-1},\ldots,\overline{U}_{1}$.

  * •

Status quo:  the expectation value for $B_{k}$ equals to $B_{k-1}$ conditioned
on the previous history, i.e.,
$\mathbb{E}[B_{k}|B_{k-1},\ldots,B_{1}]=B_{k-1}$ for all $k=1,\ldots r$.

The property for such martingale features similar properties to an unbiased
random walk. Based on Freedman’s inequality and its application to martingales
[22, 32, 40, 33, 24, 41], we have Lemma 10 in Appendix A. Starting from this
lemma, we consider the martingale provided by the qDRIFT algorithm defined
above. We define $C_{k}=B_{k}-B_{k-1}$. We observe that

$\displaystyle\norm{C_{k}}$ $\displaystyle=\norm{B_{k}-B_{k-1}}$
$\displaystyle=\norm{\mathbb{E}[\overline{U}_{i}]^{r-k}(\overline{U}_{k}-\mathbb{E}[\overline{U}_{i}])\overline{U}_{k-1}\ldots\overline{U}_{1}}$
$\displaystyle\leq\norm{\mathbb{E}[\overline{U}_{i}]^{r-k}}\norm{(\overline{U}_{k}-\mathbb{E}[\overline{U}_{i}])}\norm{\overline{U}_{k-1}\ldots\overline{U}_{1}}$
$\displaystyle=\norm{(\overline{U}_{k}-\mathbb{E}[\overline{U}_{i}])}$
$\displaystyle=2\delta\lambda_{\overline{H}}.$

We invoke Lemma 10 with $v$ bounded by $r\norm{C_{k}}^{2}$ and obtain that

$\displaystyle\Pr[\norm{B_{0}-B_{r}}\geq\epsilon]\leq
2d\exp\left(-\frac{3r\epsilon^{2}}{24(t\lambda_{\overline{H}})^{2}+4(t\lambda_{\overline{H}})\epsilon}\right),$

where $d=2^{n}$ is the dimension of the Hilbert space. With probability at
least $1-\chi$, we have

$\displaystyle\norm{\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r}}\leq\tilde{O}\left(\sqrt{\frac{\lambda_{\overline{H}}^{2}t^{2}}{r}\left(n+\log\left(\frac{1}{\chi}\right)\right)}\right).$
(14)

Finally, we provide the bound for the projection error $\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}$. We decompose
the projection error as follows:

$\displaystyle\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}$
$\displaystyle=$ $\displaystyle\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}U_{r-1}\ldots U_{1}+\overline{U}_{r}U_{r-1}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}$
$\displaystyle\leq$ $\displaystyle\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}U_{r-1}\ldots
U_{1})\Pi_{\leq\Delta}}+\norm{(\overline{U}_{r}U_{r-1}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}$
$\displaystyle\leq$
$\displaystyle\norm{(U_{r}-\overline{U}_{r})\Pi_{\leq\Delta}}+\norm{U_{r-1}\ldots
U_{1}-\overline{U}_{r-1}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}$
$\displaystyle\leq$
$\displaystyle\sum_{i=1}^{r}\norm{(U_{i}-\overline{U}_{i})\Pi_{\leq\Delta}}$
$\displaystyle=$ $\displaystyle
r\norm{(U_{i}-\overline{U}_{i})\Pi_{\leq\Delta}}.$

We now provide a bound for the truncation error in each step
$\norm{(U_{i}-\overline{U}_{i})\Pi_{\leq\Delta}}$. We assume that $U_{i}$
evolves on $H_{i}$ for some $i=1,\ldots,L$ and denote
$\delta^{\prime}=\lambda_{\overline{H}}\delta/\norm{H_{i}}$, the error can be
calculated as follows:

$\displaystyle\norm{(e^{-iH_{i}\lambda_{\overline{H}}\delta/\norm{H_{i}}}-e^{-i\overline{H}_{i}\lambda_{\overline{H}}\delta/\norm{H_{i}}})\Pi_{\leq\Delta}}$
$\displaystyle=$
$\displaystyle\norm{(e^{-iH_{i}\delta^{\prime}}-\Pi_{\leq\Delta^{\prime}}e^{-iH_{i}\delta^{\prime}}+\Pi_{\leq\Delta^{\prime}}e^{-iH_{i}\delta^{\prime}}-\Pi_{\leq\Delta^{\prime}}e^{-i\overline{H}_{i}\delta^{\prime}}+\Pi_{\leq\Delta^{\prime}}e^{-i\overline{H}_{i}\delta^{\prime}}-e^{-i\overline{H}_{i}\delta^{\prime}})\Pi_{\leq\Delta}}$
$\displaystyle\leq$
$\displaystyle\norm{\Pi_{>\Delta^{\prime}}e^{-iH_{i}\delta^{\prime}}\Pi_{\leq\Delta}}+\norm{\Pi_{\leq\Delta^{\prime}}(e^{-iH_{i}\delta^{\prime}}-e^{-i\overline{H}_{i}\delta^{\prime}})\Pi_{\leq\Delta}}+\norm{\Pi_{>\Delta^{\prime}}e^{-i\overline{H}_{i}\delta^{\prime}}\Pi_{\leq\Delta}}$
$\displaystyle\leq$ $\displaystyle
5\exp\left(-[\lambda(\Delta^{\prime}-\Delta)-\alpha sM]\right),$

where $\lambda=(2Jdk)^{-1}$ and $\alpha=eJ$. The last line follows from Lemma
7 and Lemma 8.

Now, we are ready to prove Theorem 1. We first consider the expectation error
and prove (1). In this case, we do not need to guarantee the random
fluctuation term. We compute the step complexity $r$ such that the projection
error is of the same scale as the deterministic bias, i.e.,

$\displaystyle\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}\simeq\norm{\mathbb{E}[\overline{U}_{i}]^{r}-\overline{V}}.$

Based on the previous bound on these two terms of error, we conclude that

$\displaystyle r\exp\left(-[\lambda(\Delta^{\prime}-\Delta)-\alpha\delta
M]\right)=O\left(\frac{\lambda_{\overline{H}}^{2}t^{2}}{r}\right).$

To achieve this, we can fix the value of $\Delta^{\prime}$ as follows:

$\displaystyle\Delta^{\prime}=\Delta+\frac{1}{\lambda}\alpha\delta
M+\frac{1}{\lambda}\log(\frac{r^{2}}{\lambda_{\overline{H}}^{2}t^{2}}).$

The total error is thus
$O\left(\frac{\lambda_{\overline{H}}^{2}t^{2}}{r}\right)$ with
$\lambda_{\overline{H}}=O(L\Delta^{\prime})$. Now, we consider three cases
separately to derive the final step complexity.

In the first case, $\Delta$ is the dominant term of $\Delta^{\prime}$. In this
case, we require

$\displaystyle\frac{L^{2}\Delta^{2}t^{2}}{r}\leq\epsilon\to
r_{1}\geq\frac{L^{2}\Delta^{2}t^{2}}{\epsilon}.$

In the second case, $\frac{1}{\lambda}\alpha\delta M$ is the dominant term of
$\Delta^{\prime}$. In this case, we require

$\displaystyle\frac{\alpha^{2}M^{2}t^{4}}{\lambda^{2}r^{3}}\leq\epsilon\to
r_{2}\geq L^{2/3}M^{2/3}J^{4/3}t^{4/3}\epsilon^{-1/3}.$

In the last case,
$\frac{1}{\lambda}\log(\frac{r^{2}}{\lambda_{\overline{H}}^{2}t^{2}})$ is the
dominant term of $\Delta^{\prime}$. Except for a mild polylogarithmic
correction, this term is the same as $1/\lambda$. In this case, we can prove
that the corresponding step complexity is not the dominant term.

In general, the Trotter step complexity is

$\displaystyle
r=\tilde{O}\left(\frac{L^{2}(\Delta+dkJ)^{2}t^{2}}{\epsilon}\right)+O\left(\frac{(LM)^{2/3}J^{4/3}t^{4/3}}{\epsilon^{1/3}}\right),$

which finishes the proof for (1). The corresponding gate complexity is

$\displaystyle
G=\tilde{O}\left(\frac{L^{3}(\Delta+dkJ)^{2}t^{2}}{\epsilon}\right)+O\left(\frac{(L^{5/2}M)^{2/3}J^{4/3}t^{4/3}}{\epsilon^{1/3}}\right),$

given that $O(L)$ gates are required to implement a Suzuki-Trotter formula in
each time step.

Now, we consider the step complexity required to ensure qDRIFT converges and
prove (2). If we want to further control the random fluctuation, we require

$\displaystyle r\exp\left(-[\lambda(\Delta^{\prime}-\Delta)-\alpha\delta
M]\right)=O\left(\sqrt{\frac{\lambda_{\overline{H}}t^{2}}{r}\cdot\left(n+\log(\frac{1}{\chi})\right)}\right).$

To achieve this, we can fix the value of $\Delta^{\prime}$ as follows:

$\displaystyle\Delta^{\prime}=\Delta+\frac{1}{\lambda}\alpha\delta
M+\frac{1}{2\lambda}\log(\frac{r^{3}}{\lambda_{\overline{H}}^{2}t^{2}(n+\log(1/\chi))}).$

Following a similar procedure to the average-case performance, we derive the
Trotter number to ensure a $1-\chi$ success probability for qDRIFT as

$\displaystyle
r=\tilde{O}\left(\frac{L^{2}(\Delta+dkJ)^{2}t^{2}}{\epsilon^{2}}\left(n+\log(\frac{1}{\chi})\right)\right)+O\left(\frac{(LM)^{2/3}J^{4/3}t^{4/3}}{\epsilon^{2/3}}\left(n+\log(\frac{1}{\chi})\right)^{1/3}\right),$
(15)

which finishes the proof for Theorem 1. In addition, the corresponding gate
complexity is

$\displaystyle
G=\tilde{O}\left(\frac{L^{3}(\Delta+dkJ)^{2}t^{2}}{\epsilon^{2}}\left(n+\log(\frac{1}{\chi})\right)\right)+O\left(\frac{(L^{5/2}M)^{2/3}J^{4/3}t^{4/3}}{\epsilon^{2/3}}\left(n+\log(\frac{1}{\chi})\right)^{1/3}\right).$

##### Random permutation.

In the previous qDRIFT algorithm, we only consider applying a random evolution
in each time step, which is a random version of the first-order Lie-Trotter
formula. A randomized approach for higher-order formulas, known as random
permutation, was proposed in Ref. [12]. As shown in Algorithm 2, this
algorithm randomly permutes the order of Hamiltonian terms within each block
to $S_{p}^{\sigma_{i}}(\delta=t/r)$ for a permutation
$\sigma_{i}\in\mathcal{S}_{L}$ from the permutation group.

0: Hamiltonian $H=\sum_{l=1}^{L}H_{l}$, evolution time $t$, and number of
steps $r$

1: At each time interval $\delta=t/r$: evolve a random term in Hamiltonian
$U_{i}=S_{p}^{\sigma_{i}}(\delta)$ with a randomly chosen sequence with
probability $1/L!$

1: The unstructured (randomly generated) product formula $U=U_{r}\ldots U_{1}$

Algorithm 2 The random permutation algorithm

We now consider the performance of this algorithm in the low-energy subspace
and prove Theorem 2. Similar to the previous section, we can decompose the
error into three terms:

$\displaystyle\quad\ \norm{(U_{r}\ldots U_{1}-V)\Pi_{\leq\Delta}}$
$\displaystyle\leq\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}+\norm{(\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r})\Pi_{\leq\Delta}}+\norm{(\mathbb{E}[\overline{U}_{i}]^{r}-V)\Pi_{\leq\Delta}}$
$\displaystyle\leq\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}+\norm{\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r}}+\norm{\mathbb{E}[\overline{U}_{i}]^{r}-\overline{V}},$

where $\overline{U}_{i}=S_{p}^{\sigma_{i}}(\delta)$ with Hamiltonian
$\overline{H}$ and $\overline{U}=e^{-i\overline{H}t}$. The first term is the
projection error for the low-energy subspace estimation. The second term is
the random fluctuation in the low-energy subspace. The third term is the
deterministic bias in the low-energy subspace. We introduce two new parameters
$\lambda_{h}=\max_{l}\norm{H_{l}}=JML=O(n)$ and
$\lambda_{\overline{h}}=\max_{l}\norm{\overline{H}_{l}}=\Delta^{\prime}$ for
$n$-qubit $k$-local Hamiltonian.

We begin with the deterministic bias error term. According to the direct
calculation [12, 36], we can derive that (see, e.g. Eq. (C14) of Ref. [9])

$\displaystyle\norm{\mathbb{E}[\overline{U}_{i}]-V^{1/r}}\leq
O\left(\left(\frac{t\lambda_{\overline{h}}}{r}\right)^{p+1}L^{p}\right).$

According to Lemma 3.5 of Ref. [9], for Hamiltonian
$\overline{H}=\sum_{l=1}^{L}\overline{H}_{l}$ and $\overline{U}_{i}$, the
deterministic bias term is bounded by

$\displaystyle\norm{\mathbb{E}[\overline{U}_{i}]^{r}-\overline{V}}\leq
r\norm{\mathbb{E}[\overline{U}_{i}]-\overline{V}^{1/r}}\leq
O\left(\left(\frac{L}{r}\right)^{p}\cdot(t\lambda_{\overline{h}})^{p+1}\right).$

Here, $\overline{U}^{1/r}=e^{-iHt/r}=e^{iH\delta}$.

The next step is to provide a bound for the random fluctuation term
$\norm{\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r}}$.
We still use the Freedman inequality for martingales, and define
$B_{k}=(\mathbb{E}[\overline{U}_{i}])^{r-k}\overline{U}_{k}\cdots\overline{U}_{1}$,
$B_{0}=(\mathbb{E}[\overline{U}_{i}])^{r}$ and $C_{k}=B_{k}-B_{k-1}$. For the
random permutation approach, each $C_{k}$ is bounded by

$\displaystyle\norm{C_{k}}\leq\norm{\overline{U}_{i}-\mathbb{E}[\overline{U}_{i}]}=\norm{\overline{U}_{i}-\overline{V}^{1/r}}+\norm{\mathbb{E}[\overline{U}_{i}]-\overline{V}^{1/r}}=O\left(\left(\frac{tL\lambda_{\overline{h}}}{r}\right)^{p+1}\right).$

Here, the last equation follows from the fact [12, 38] that
$\norm{\overline{U}_{i}-\overline{V}^{1/r}}=O\left(\left(\frac{tL\lambda_{\overline{h}}}{r}\right)^{p+1}\right)$
and
$\norm{\mathbb{E}[\overline{U}_{i}]-\overline{V}^{1/r}}=O\left(\left(\frac{t\lambda_{\overline{h}}}{r}\right)^{p+1}L^{p}\right)$.
Therefore, we can derive from Lemma 10 that with probability at least $1-\chi$
the random fluctuation scales as

$\displaystyle\norm{\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r}}\leq
O\left(\sqrt{\frac{(tLh_{\lambda}^{\prime})^{2p+2}}{r^{2p+1}}\cdot\left(n+\log(\frac{1}{\chi})\right)}\right).$
(16)

Finally, we provide the bound for the projection error $\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}$. According to
Lemma 9, we decompose the projection error as follows:

$\displaystyle\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}\leq r\cdot
5\exp(-\frac{1}{q}(\lambda(\Delta^{\prime}-\Delta)-\alpha^{\prime}\delta
ML-q\log q)),$

where $\lambda=(2Jdk)^{-1}$, $\alpha=eJ$, and $q(p)$ is the length of the
product formula in each step depending only on $p$.

Now, we are ready to derive Theorem 2. Similar to the previous section, we
still consider two cases:

In the first case, we consider the average-case performance of the random
permutation approach in the low-energy subspace and prove (4). In this case,
we ignore the random fluctuation and compute the step complexity $r$ such that
the projection error is of the same scale as the deterministic bias, i.e.,

$\displaystyle\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}\simeq\norm{\mathbb{E}[\overline{U}_{i}]^{r}-\overline{V}}.$

We can derive the bound on $\Delta^{\prime}$ and $r$ as

$\displaystyle\Delta^{\prime}$
$\displaystyle=\Delta+\frac{1}{\lambda}\alpha\delta ML+\frac{q}{\lambda}\log
q+\frac{q}{\lambda}\log\left(\frac{r^{p+1}}{t^{p+1}\lambda_{\overline{h}}^{p+1}L^{p}}\right),$
$\displaystyle r$ $\displaystyle=O\left(\frac{Lt^{1+1/p}(\Delta+dkJq\log
q)^{1+1/p}}{\epsilon^{1/p}}\right)+O\left(\frac{Lt^{1+1/(2p+1)}(dkMJ^{2})^{\frac{1}{2}+\frac{1}{4p+2}}}{\epsilon^{1/(2p+1)}}\right).$

This finishes the proof for (4), which provides an $O(L^{1/p})$ or
$O(L^{1/(2p+1)})$ reduction on the step complexity compared to standard
product formulas [35]. This means that as $p$ increases, the random
permutation will eventually lose its reduction—a similar observation to that
in Ref. [12]. The gate complexity requires an $O(L)$ overhead from the Trotter
number $r$.

In the second case, we want to further control the random fluctuation. Thus we
require

$\displaystyle\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}\simeq\norm{\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r}}.$

We can derive the bound on $\Delta^{\prime}$ and $r$ to ensure that the
simulation error is bounded by $\epsilon$ with probability at least $1-\chi$
as

$\displaystyle\Delta^{\prime}$
$\displaystyle=\Delta+\frac{1}{\lambda}\alpha\delta ML+\frac{q}{\lambda}\log
q+\frac{q}{2\lambda}\log\left(\frac{r^{2p+3}(n+\log(1/\chi))}{t^{2p+2}\lambda_{\overline{h}}^{2p+2}L^{2p+2}}\right),$
(17) $\displaystyle r$
$\displaystyle=O\left(\frac{(Lt)^{1+1/(2p+1)}(\Delta+dkJq\log
q)^{1+1/(2p+1)}(n+\log(1/\chi)^{1/(2p+1)})}{\epsilon^{2/(2p+1)}}\right)$
$\displaystyle\quad+O\left(\frac{(Lt)^{1+1/(4p+3)}(dkMJ^{2})^{\frac{2p+2}{4p+3}}(n+\log(1/\chi)^{1/(4p+3)})}{\epsilon^{2/(4p+3)}}\right),$
(18)

which finishes the proof for (5). The overhead for controlling randomized
fluctuation will quickly disappear as the order $p$ increases. Yet, the
advantage of this randomized approach compared to deterministic product
formulas would also disappear correspondingly. Similarly, The gate complexity
requires an $O(L)$ overhead from the Trotter number $r$.

##### Doubling the order of product formula.

In Ref. [16], Cho, Berry, and Huang proposed an approach to double the order
of digital quantum simulations via the randomized product formulas. We also
consider the performance of this approach in the low-energy subspace and
obtain Corollary 1. We refer to Ref. [16] for a detailed description of this
algorithm. As the algorithm requires the implementation of some random well-
designed unitaries
$\left\\{U_{h}^{(l)}=\exp(\alpha_{h,l}H_{h}^{(l)})\right\\}$ according to some
well-designed probability $\\{p_{h,l}\\}$, we have to assume that
$U_{h}^{(l)}$’s also satisfy the properties of $\exp(-iH_{l}\delta)$’s. We
still consider the error decomposition as follows:

$\displaystyle\quad\ \norm{(U_{r}\ldots U_{1}-V)\Pi_{\leq\Delta}}$
$\displaystyle\leq\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}+\norm{\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r}}+\norm{\mathbb{E}[\overline{U}_{i}]^{r}-\overline{V}}.$

According to Theorem 2 of Ref. [16] and the assumption that
$r\geq(5^{p/2-1}+\frac{1}{2})L\lambda_{\overline{h}}t$, we can first compute
the deterministic bias as

$\displaystyle\norm{\mathbb{E}[\overline{U}_{i}]^{r}-\overline{V}}\leq
O\left(\frac{(L\lambda_{\overline{h}}t)^{2p+2}}{r^{2p+1}}\right),$

For the random fluctuation, we still define
$B_{k}=\overline{U}_{r}\ldots\overline{U}_{1}$,
$B_{0}=(\mathbb{E}[\overline{U}_{i}])^{r}$, and $C_{k}=B_{k}-B_{k-1}$. We can
compute the norm of $C_{k}$ as

$\displaystyle\norm{C_{k}}\leq\norm{\overline{U}-\mathbb{E}[\overline{U}_{i}]}\leq
O\left(\frac{(L\lambda_{\overline{h}}t)^{p+1}}{r^{p+1}}\right).$

Therefore, we can derive that with probability at least $1-\chi$, the random
fluctuation scales as

$\displaystyle\norm{\overline{U}_{r}\ldots\overline{U}_{1}-\mathbb{E}[\overline{U}_{i}]^{r}}\leq
O\left(\sqrt{\frac{(tL\lambda_{\overline{h}})^{2p+2}}{r^{2p+1}}\cdot\left(n+\log(\frac{1}{\chi})\right)}\right).$

For the projection error, we have

$\displaystyle\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}\leq r\cdot
5\exp(-\frac{1}{q}(\lambda(\Delta^{\prime}-\Delta)-\alpha^{\prime}\delta
ML-q\log q)),$

where $\lambda=(2Jdk)^{-1}$, $\alpha=eJ$, and $q(p)$ is the length in each
step depending only on $p$.

To prove Corollary 1, we consider the step complexity required to ensure the
average-case performance and the convergence of the algorithm, respectively.
In the first setting, we ignore the random fluctuation and compute the step
complexity $r$ such that the projection error is of the same scale as the
deterministic bias, i.e.,

$\displaystyle\norm{(U_{r}\ldots
U_{1}-\overline{U}_{r}\ldots\overline{U}_{1})\Pi_{\leq\Delta}}\simeq\norm{\mathbb{E}[\overline{U}_{i}]^{r}-\overline{V}}.$

We can derive the bound on $\Delta^{\prime}$ and $r$ as

$\displaystyle\Delta^{\prime}$
$\displaystyle=\Delta+\frac{1}{\lambda}\alpha\delta ML+\frac{q}{\lambda}\log
q+\frac{q}{\lambda}\log\left(\frac{r^{2p+2}}{t^{2p+2}\lambda_{\overline{h}}^{2p+2}L^{2p+2}}\right),$
$\displaystyle r$
$\displaystyle=O\left(\frac{(Lt)^{1+1/(2p+1)}(\Delta+dkJq\log
q)^{1+1/(2p+1)}}{\epsilon^{1/(2p+1)}}\right)+O\left(\frac{(Lt)^{1+1/(4p+3)}(dkMJ^{2})^{\frac{2p+2}{4p+3}}}{\epsilon^{1/(4p+3)}}\right).$

This finishes the proof for $r_{\text{exp}}$ in (6). For $r_{\text{prob}}$, we
follow a similar proof toTheorem 2 to obtain the same result. The gate
complexity required correspondingly is

$\displaystyle G$
$\displaystyle=O\left(\frac{L(Lt)^{1+1/(2p+1)}(\Delta+dkJq\log
q)^{1+1/(2p+1)}}{\epsilon^{1/(2p+1)}}\right)+O\left(\frac{L(Lt)^{1+1/(4p+3)}(dkMJ^{2})^{\frac{2p+2}{4p+3}}}{\epsilon^{1/(4p+3)}}\right).$

(a) qDRIFT for system size $n=8$

(b) Random permutation for system size $n=8$

(c) qDRIFT for step number $r=10^{6}$

(d) Random permutation for step number $r=10^{3}$

Figure 2: Full Hilbert space vs. low-energy subspace: simulation error induced
by the qDRIFT and random permutation algorithms for localized Heisenberg
spin-$\mathbf{1/2}$ models. The solid lines denote errors in the full Hilbert
space while the dashed lines denote errors in the low-energy subspace. We
distinguish the total evolution time $t$ or the order of the product formula
$p$ by the color of the lines. (a) We plot the worst-case and the low-energy
simulation errors of the qDRIFT algorithm as a function of step number $r$.
(b) The same as the upper left subfigure, but for the random permutation
approach at $p=2$ and $p=4$. (c) We plot the worst-case and the low-energy
simulation errors of the qDRIFT algorithm as a function of system size $n$.
(d) The same as the lower left subfigure, but for the random permutation
approach at $p=2$ and $p=4$.

To benchmark the performance of random Hamiltonian simulation approaches in
the low-energy subspace, we further carry out extensive numerical experiments
concerning simulating the one-dimensional localized homogeneous Heisenberg
models of different system sizes $n\in[4,8]$ without external fields as
follows:

$H=\underbrace{-\sum_{i=1}^{n-1}X_{i}X_{i+1}}_{H_{1}}\underbrace{-\sum_{i=1}^{n-1}Y_{i}Y_{i+1}}_{H_{2}}\underbrace{-\sum_{i=1}^{n-1}Z_{i}Z_{i+1}}_{H_{3}},$
(19)

where $X_{i}$, $Y_{i}$, and $Z_{i}$ are Pauli operators acting on the $i$-th
spin, and the summation is over every adjacent spin pair. The Hamiltonian in
(19) is decomposed into three interaction terms: $H=H_{1}+H_{2}+H_{3}$,
containing all $X$, $Y$, and $Z$ terms, respectively, as indicated by the
curly brackets. Since the ground energy $E_{0l}<0$ for every $H_{l}$, we
implement the shift $H_{l}\rightarrow H_{l}-E_{0l}I$ for every $H_{l}$. We
empirically pick $\Delta=14$ (the shifted ground energy of the largest system
$n=8$) such that for system sizes $n\in[4,8]$, the low-energy subspace covers
neither (nearly) the whole spectrum nor a vacant subspace.

In Figure 2, we plot the qDRIFT error and the random permutation error as a
function of step number $r$ and system size $n$ in four subfigures. We fix the
total evolution time $t$ and randomly sample $10$ sequences for each choice of
parameters. For the qDRIFT approach, in order to distinguish between closely
spaced error bars, we made a tiny shift to the abscissa of nearby data points.
We can observe that the projection to the low-energy subspace reduces the
error compared to the full Hilbert space for the same parameter configuration.
We plot the numerical results for the qDRIFT approach in the left two figures
of Figure 2. We estimate that $\epsilon=O(nr^{-1/2})$ for full Hilbert space
simulations and $\epsilon=O(n^{2/3}r^{-1/2})$ for low-energy simulations. The
numerical performances for both cases have better scaling concerning the $n$
dependence than theoretical bounds as the latter one gives
$\epsilon=O(n^{3/2}r^{-1/2})$ and
$\epsilon=O(\max\\{n^{1/2}r^{-1/2},n^{3/2}r^{-3}\\})$ for the two cases as
shown in Ref. [9] and (15). We remark that we take the random fluctuation into
consideration here as the variance is much larger than the average value for
simulation errors. In the right two figures, we plot the numerical results for
the random permutation approach. We estimate that
$\epsilon=O(n^{p/2+1}r^{-p})$ and $\epsilon=O(nr^{-p})$ for worst-case and
low-energy simulations, which are again better than the theoretical bounds of
$\epsilon=O(n^{p+1}r^{-p})$ and
$\epsilon=O(\max\\{n^{1/2}r^{-p-1/2},n^{p+1}r^{-2p-3/2}\\})$ in Ref. [12] and
(18). We can conclude that numerical error for the low-energy subspace has a
better performance than the theoretical bounds in numerics and provides a
steady reduction in the simulation error with various parameters and settings.

###  3.2 Symmetry protections

In this subsection, we consider an alternative approach, symmetry protection
[39]. Given a Hamiltonian $H$, we consider the group of unitary
transformations denoted by $\mathcal{C}$ such that $[H,C]\equiv 0$ for each
unitary $C$ chosen from $\mathcal{C}$. As mentioned in Section 1.1, we
implement the Lie-Trotter formula $S_{1}(\delta)$ in each time step with a
symmetry transformation $C\in\mathcal{C}$. Explicitly, the simulation for
$V=e^{-iHt}$ reads

$\displaystyle
U(t)=\prod_{\mu=1}^{r}C_{\mu}^{\dagger}S_{1}(\delta)C_{\mu}\approx V(t),$

where $V(t)=e^{-iHt}$ is the ideal evolution.

We briefly recap the results in Ref. [39] for the intuition of reducing
simulation error using symmetry protection. Intuitively, $U(t)$ can be
represented as an evolution on a slightly different Hamiltonian
$H_{\text{eff}}$ and

$\displaystyle
U(t)=\prod_{\mu=1}^{r}e^{-iC_{\mu}^{\dagger}H_{\text{eff}}C_{\mu}\delta}=\prod_{\mu=1}^{r}e^{-i(H+C_{\mu}^{\dagger}\kappa
C_{k})\delta}\approx
e^{-i\left(H+\frac{1}{r}\sum_{\mu=1}^{r}C_{\mu}^{\dagger}\kappa
C_{\mu}\right)t}\approx V(t),$

where $\kappa=H_{\text{eff}}-H$.

For $k$-local Hamiltonians, we define the following quantities for a $k$-local
Hamiltonian [14]:

$\displaystyle\alpha=\sum_{\gamma_{1}=1}^{L}\sum_{\gamma_{2}=\gamma_{1}+1}^{L}\norm{[H_{\gamma_{2}},H_{\gamma_{1}}]}=O(n),$
$\displaystyle\beta=\sum_{\gamma_{1}=1}^{L}\sum_{\gamma_{2}=\gamma_{1}+1}^{L}\sum_{\gamma_{3}=\gamma_{2}}^{L}\norm{[H_{\gamma_{3}},[H_{\gamma_{2}},H_{\gamma_{1}}]]}=O(n),$
$\displaystyle\gamma=\sum_{\gamma_{1}=1}^{L}\sum_{\gamma_{2}=\gamma_{1}+1}^{L}\norm{[H,[H_{\gamma_{2}},H_{\gamma_{1}}]]}=O(n^{2}).$

Ref. [39] obtains the following result concerning the simulation error in each
time step:

######  Lemma 3 (Lemma 7 of [39]).

For all $\delta$ such that $\beta\delta\leq 2\alpha$ and
$\alpha^{2}\delta\leq\beta+\gamma$, the Trotter error in each step is
represented by

$\displaystyle
V(\delta)-S_{1}(\delta)=V(\delta)v_{0}\frac{\delta^{2}}{2}+\mathcal{V}(\delta),$

where
$v_{0}=\sum_{\gamma_{1}=1}^{L}\sum_{\gamma_{2}=\gamma_{1}+1}^{L}[H_{\gamma_{2}},H_{\gamma_{1}}]$
and $\mathcal{V}(\delta)\leq\iota\delta^{3}$ with $\iota=5(\gamma+\beta)/6$.

Based on the above result, the Trotterization error for the full Hilbert space
simulation is:

######  Lemma 4 (Theorem 1 of [39]).

The simulation error for the whole evolution using symmetry protection for
$k$-local Hamiltonians is bounded by:

$\displaystyle\norm{V(t)-\prod_{\mu=1}^{L}C_{\mu}^{\dagger}S_{1}(\delta)C_{\mu}}\leq\overline{v}_{0}\frac{t^{2}}{2r}+\iota\frac{t^{3}}{r^{2}}+2\left(\frac{1}{2}\norm{v_{0}}+\iota\frac{t}{r}\right)^{2}\frac{t^{4}}{r^{2}},$

where

$\displaystyle\overline{v}_{0}=\norm{\frac{1}{r}\sum_{\mu=1}^{r}C_{\mu}^{\dagger}V(\mu
s)^{\dagger}v_{0}V(\mu\delta)C_{\mu}}$

is the averaged first-order error term.

We then focus on $\overline{v}_{0}$. In general, we can expect that
$\norm{\overline{v}_{0}}\propto\norm{v_{0}}/r^{\theta}$ for $\theta\in[0,1]$.
Hence, we require a step complexity of

$\displaystyle
r=\max\left\\{\Theta\left(\left(\frac{nt^{2}}{\epsilon}\right)^{\frac{1}{1+\theta}}\right),\Theta\left(\frac{nt^{3/2}}{\sqrt{\epsilon}}\right)\right\\}.$

to bound the simulation error below $\epsilon$. There are several examples
given for the value this $\theta$. For example, when we draw symmetry
transformations randomly, the behavior of the error would be analogous to a
random walk, which results in $\theta=0.5$. The rigorous proof for this
intuition for localized Heisenberg model was provided in [39]. There are also
some instances when $\theta=1$ [39, 21].

Now, we consider the performance of this symmetry protection approach in the
low-energy subspace, which is also an open question in [39]. Our goal is to
compute the following error term

$\displaystyle\norm{\left(V(t)-\prod_{\mu=1}^{r}C_{\mu}^{\dagger}S_{1}(\delta)C_{\mu}\right)\Pi_{\leq\Delta}}.$

We start with the following observation:

######  Proposition 1.

Suppose $\mathcal{S}$ is the symmetry group for $H$, then $\mathcal{S}$ is a
subgroup of the symmetry group for $\overline{H}$.

We then decompose the error into three parts:

$\displaystyle\norm{\left(V(t)-\prod_{\mu=1}^{r}C_{\mu}^{\dagger}S_{1}(\delta)C_{\mu}\right)\Pi_{\leq\Delta}}$
$\displaystyle=\norm{(V(t)-\overline{V}(t))\Pi_{\leq\Delta}}+\norm{\left(\overline{V}(t)-\prod_{\mu=1}^{r}C_{\mu}^{\dagger}\overline{S}_{1}(\delta)C_{\mu}\right)\Pi_{\leq\Delta}}$
$\displaystyle+\norm{\left(\prod_{\mu=1}^{r}C_{\mu}^{\dagger}\overline{S}_{1}(\delta)C_{\mu}-\prod_{\mu=1}^{r}C_{\mu}^{\dagger}S_{1}(\delta)C_{\mu}\right)\Pi_{\leq\Delta}},$

where $\overline{S}_{1}(\delta)$ is the Lie-Trotter formula for
$\overline{H}=\sum_{l=1}^{L}\overline{H}_{l}$. The first part is $0$ according
to the property of the projector. The second part, according to the results of
the symmetry protection approach, gives:

$\displaystyle\epsilon_{2}=O\left(\frac{L\Delta^{\prime}t^{2}}{r^{1+\theta}}\right)+O\left(\frac{L^{2}\Delta^{{}^{\prime}2}t^{3}}{r^{2}}\right).$

According to Lemma 9 and the fact that $q=L$ for first-order Trotterizations,
we have

$\displaystyle\norm{(\overline{S}_{1}(\delta)-S_{1}(\delta))\Pi_{\leq\Delta}}=5\exp(-\frac{1}{L}(\lambda(\Delta^{\prime}-\Delta)-\alpha^{\prime}\delta
ML-L\log L)).$

Analogous to the intuition from symmetry transformation, we assume

$\displaystyle\quad\
\norm{\left(\prod_{\mu=1}^{r}S_{\mu}^{\dagger}\overline{S}_{1}(\delta)S_{\mu}-\prod_{\mu=1}^{r}S_{\mu}^{\dagger}S_{1}(\delta)S_{\mu}\right)\Pi_{\leq\Delta}}$
$\displaystyle=O\left(r^{1-\theta^{\prime}}\exp\left(-\frac{1}{L}(\lambda(\Delta^{\prime}-\Delta)-\alpha^{\prime}\delta
ML-L\log L)\right)\right)$

for some $\theta^{\prime}\in[0,1]$.

We can derive the step complexity here as $(\lambda=(2Jdk)^{-1})$:

$\displaystyle
r=\tilde{O}\left(\left(\frac{L(\Delta+JdkL)t^{2}}{\epsilon}\right)^{\frac{1}{1+\theta}}+\left(\frac{L^{2}MJdkt^{3}}{\epsilon}\right)^{\frac{1}{2+\theta}}+\frac{L(\Delta+JdkL)t^{\frac{3}{2}}}{\epsilon^{\frac{1}{2}}}+\frac{t^{\frac{5}{4}}L(JdkM)^{\frac{1}{2}}}{\epsilon^{\frac{1}{4}}}\right),$

which finishes the proof for Theorem 3. The gate complexity of this approach
depends on the property of each symmetry transformation $C_{\mu}$ and varies
from case to case. However, we can still observe that the step complexity
above is strictly better than that provided by either standard Trotterizations
or the symmetry-protected simulation in the full Hilbert space.

(a) Symmetry transformation for $n=8$

(b) Symmetry transformation for $r=10^{3}$

Figure 3: Full Hilbert space vs. low-energy subspace: Error induced by the
symmetry protection approach for Heisenberg spin-$\mathbf{1/2}$ models. The
solid lines denote errors in the full Hilbert space while the dashed lines
denote errors in the low-energy subspace. We distinguish the experiments with
different parameters and settings by the color of the lines. (a) We plot the
worst-case and the low-energy simulation errors as a function of step number
$r$. (b) We plot the worst-case and the low-energy simulation errors as a
function of system size $n$.

We conduct numerical experiments concerning simulating the one-dimensional
Heisenberg models of different system sizes $n\in[4,8]$ in (19). We consider
three different schemes of simulations: the standard first-order
Trotterization (“Standard”), the first-order Trotterization protected by a
random set symmetry transformation (“Random-ST”), and the first-order
Trotterization protected by the optimal set (“Optimal-SP”). For the “Random-
ST” approach, considering the $SU(2)$-symmetry of Heisenberg models, each step
we sample a $2\times 2$ matrix from norm distribution ensuring diagonal
conjugation, then scale the determinant to $1$ and finally tensorize the
normalized matrix to the $n$-th power. We randomly repeat this process 10
times and plot the average value as well as error bars. For the “Optimal-SP”
approach, we employ $C_{0}=U_{H}^{\otimes n},C_{\mu}=C_{0}^{\mu}$ as [39],
where $U_{H}$ denotes single-qubit Hadamard gate.

In Figure 3, we plot the symmetry protection error as a function of step
number $r$ and system size $n$ in two subfigures. In our numerical
simulations, we fix the total evolution time $t=1.0$. We can observe that the
projection to the low-energy subspace reduces the error compared to the full
Hilbert space for the same parameter configuration. We estimate that full
Hilbert space simulations and low-energy simulations respectively yield
$\epsilon=O(nr^{-1})$ and $\epsilon=O(n^{2/5}r^{-1})$ for the standard Trotter
formula, $\epsilon=O(nr^{-3/2})$ and $\epsilon=O(n^{1/5}r^{-3/2})$ for the
random set symmetry transformation, $\epsilon=O(nr^{-2})$ and
$\epsilon=O(n^{4/5}r^{-2})$ for the optimal set symmetry protection. The
numerical performances for all the three approaches conform to the scaling of
theoretical bounds, which give
$\epsilon=O(\max\\{nr^{-(1+\theta)},n^{2}r^{-2}\\})$ and
$\epsilon=O(\max\\{r^{-(1+\theta)},r^{-2},nr^{-(2+\theta)},n^{2}r^{-4}\\})$
for the two cases.

###  3.3 Simulating power-law models

Except for $k$-local Hamiltonians, we also show that simulating low-energy
states provides improvement for power-law interactions
$H=\sum_{{\bm{i}},{\bm{j}}\in\Sigma}H_{{\bm{i}},{\bm{j}}}$ with exponent
satisfying:

$\displaystyle\norm{H_{{\bm{i}},{\bm{j}}}}\leq\begin{cases}1&\quad\text{if
}{\bm{i}}={\bm{j}},\\\
\frac{1}{\norm{{\bm{i}}-{\bm{j}}}_{2}^{\alpha}}&\quad\text{if
}{\bm{i}}\neq{\bm{j}},\end{cases}$

where ${\bm{i}},{\bm{j}}\in\Sigma$ are the qubit sites,
$\Lambda\subseteq\mathbb{R}^{D}$ is a $D$-dimensional square lattice, and
$H_{{\bm{i}},{\bm{j}}}$ is an operator supported on two sites
${\bm{i}},{\bm{j}}$. Through direct computation [14], we have

$\displaystyle
g=\norm{\max_{{\bm{i}}}\sum_{{\bm{j}}\neq{\bm{i}}}H_{{\bm{i}},{\bm{j}}}}\leq\begin{cases}O(n^{1-\alpha/D})&\quad\text{for
}0\leq\alpha<D,\\\ O(\log(n))&\quad\text{for }\alpha=D,\\\ O(1)&\quad\text{for
}\alpha>D.\\\ \end{cases}$

Here, $g$ is an upper bound on the strengths of the interactions associated
with a single spin. We can observe that every term of the power-law
Hamiltonian is $2$-local. As assumed by the Theorem 4, each part $H_{l}$
contains at most $k$ qubits and strength at most $M$. We consider an operator
$A$, and denote $E_{A}$ as the subset of interaction terms in $H$ that do not
commute with $A$ and $R$ as the strengths of the terms in $E_{A}$. Then we can
bound $R$ by $R\leq 2g$ and $R\leq 2gn$ for any $H_{l}$ or $H_{l}^{n}$ for
arbitrary $n$. We consider the following term:

$\displaystyle\norm{e^{sH}Ae^{-sH}}.$

Using the Hadamard formula [31], we write

$\displaystyle\norm{e^{sH}Ae^{-sH}}=A+s[H,A]+\frac{s^{2}}{2!}[H,[H,A]]+\cdots=\sum_{\ell=0}^{\infty}\frac{s^{\ell}}{\ell!}K_{\ell}.$

For $\ell=0$, we have $\norm{K_{\ell}}=\norm{A}$, and for $\ell>0$, we have

$\displaystyle K_{\ell}\coloneqq\sum_{X_{1}\in
E_{A}}\sum_{X_{2}|X_{1}}\sum_{X_{3}|X_{2},X_{1}}\cdots\sum_{X_{\ell}|X_{\ell-1},\cdots,X_{1}}[h_{X_{\ell}},[h_{X_{\ell-1}},\cdots,[h_{X_{1}},A]\cdots]],$

where $h_{X}$ denotes a term acting on a two-qubit string $X$ and the sum
$\sum_{X_{j}|X_{j-1},\ldots,X_{1}}$ denotes a summation over the non-zero
terms in the commutator. We can observe that the norm of $\norm{K_{\ell}}$ can
be bounded by

$\displaystyle\norm{K_{\ell}}\leq\sum_{X_{1}\in
E_{A}}\norm{h_{X_{1}}}\cdots\sum_{X_{\ell}|X_{\ell-1}\cdots
X_{1}}\norm{h_{X_{\ell}}}\norm{A}.$

Now, we compute the value for each term. Notice that the sum of the norms of
$h_{X}$ that do not commute with $A$ is bounded by $R$. The sum of the norms
of $h_{X}$ that does not commute with another $h_{Y}$ is bounded by $2g$. We
have

$\displaystyle\sum_{X_{j}|X_{j-1},\ldots,X_{1}}\norm{h_{X_{j}}}\leq
R+2(j-1)g.$

Therefore, we have

$\displaystyle\norm{K_{\ell}}\leq\norm{A}R(R+2g)(R+4g)\cdots(R+2(\ell-1)g)=(2g)^{\ell}r(r+1)\cdots(r+\ell-1)\norm{A},$

where $r\coloneqq\frac{R}{2g}$. Thus, we have

$\displaystyle\norm{e^{sH}Ae^{-sH}}\leq\norm{A}\sum_{\ell=0}^{\infty}\frac{(2gs)^{\ell}}{\ell!}r(r+1)\cdots(r+\ell-1)=\frac{1}{(1-2gs)^{r}}.$

Formally, we have the following lemma:

######  Lemma 5.

Consider a Hamiltonian with $2$-local terms and suppose the interaction on
each qubit is bounded by strength $g$. For any $0\leq s\leq\frac{1}{2g}$, we
have $\norm{e^{sH}Ae^{-sH}}\leq\norm{A}(1-2gs)^{-R/2g}$, where $g$, $R$, and
$E_{A}$ are defined as above.

Following the relationship [1] that

$\displaystyle\norm{\Pi_{>\Lambda^{\prime}}A\Pi_{\leq\Lambda}}=\norm{\Pi_{>\Lambda^{\prime}}e^{-sH}e^{sH}Ae^{-sH}e^{sH}\Pi_{\leq\Lambda}}\leq\norm{e^{sH}Ae^{-sH}}\cdot
e^{-s(\Lambda^{\prime}-\Lambda)},$

we obtain the following lemma:

######  Lemma 6.

Consider a Hamiltonian with $2$-local terms and suppose the interaction on
each qubit is bounded by strength $g$. We have

$\displaystyle\norm{\Pi_{>\Lambda^{\prime}}A\Pi_{\leq\Lambda}}\leq\norm{A}\cdot
e^{-\lambda(\epsilon^{\prime}-\epsilon-2R)},$

where $\lambda\coloneqq(4g)^{-1}$.

According to the above lemma, we can deduce that

$\displaystyle\norm{\Pi_{>\Lambda^{\prime}}(H_{l})^{n}\Pi_{\leq\Lambda}}\leq((M)^{n}e^{-\frac{1}{4g}(\Lambda^{\prime}-\Lambda-2g)}\leq(e^{-1/2}M)^{n}e^{-\frac{1}{4g}(\Lambda^{\prime}-\Lambda)},$

where we assume that $\Lambda^{\prime}-\Lambda-2R\geq 0$.

Thus, we can deduce that

$\displaystyle\norm{\Pi_{>\Lambda^{\prime}}e^{-i\delta
H_{l}}\Pi_{\leq\Lambda}}$ $\displaystyle\leq
e^{-\lambda(\Lambda^{\prime}-\Lambda)}(e^{\alpha M\delta}-1),$
$\displaystyle\norm{\Pi_{>\Lambda^{\prime}}e^{-i\delta\overline{H}_{l}}\Pi_{\leq\Lambda}}$
$\displaystyle\leq 3e^{-\lambda(\Lambda^{\prime}-\Lambda)}(e^{\alpha
M\delta}-1).$

for $\alpha=e^{-1/2}$ a constant. Compared with Lemma 8 and Lemma 7 in
Appendix A, we find that $k$ is exactly replaced by a constant $2g$. Thus,
similar to the proof for Lemma 2, we can prove the step complexity in Theorem
4 as:

$\displaystyle
r=\tilde{O}\left(\frac{t^{1+\frac{1}{p}}L^{\frac{1}{p}}}{\epsilon^{\frac{1}{p}}}(\Delta+gq\log
q)^{1+\frac{1}{p}}\right)+O\left(\frac{t^{1+\frac{1}{2p+1}}}{\epsilon^{\frac{1}{2p+1}}}(LMg)^{\frac{1}{2}+\frac{1}{4p+2}}\right).$

It is straightforward to observe that the gate complexity required is

$\displaystyle
G=\tilde{O}\left(\frac{t^{1+\frac{1}{p}}L^{\frac{1}{p}+1}}{\epsilon^{\frac{1}{p}}}(\Delta+gq\log
q)^{1+\frac{1}{p}}\right)+O\left(\frac{Lt^{1+\frac{1}{2p+1}}}{\epsilon^{\frac{1}{2p+1}}}(LMg)^{\frac{1}{2}+\frac{1}{4p+2}}\right).$

Compared with the gate count required for full-energy simulation tasks in
[14], our result gets rid of the explicit overhead of $O(n^{2})$ gates to
implement each Trotter step. When we combine multiple terms in each $H_{l}$,
the overhead can be largely reduced.

Here, we perform numerical experiments to benchmark the product formulas in
simulating low-energy states for power-law Hamiltonians. In particular, we
study the following 2D power-law Heisenberg spin-$1/2$ model without external
fields (as shown in the left-hand side of Figure 4):

$H=-\sum_{<i,j>}\frac{1}{\|i-j\|_{2}^{\alpha}}(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}),$
(20)

where $X_{i},Y_{i}$, and $Z_{i}$ are Pauli operators acting on the $i$-th
spin, and the summation is over every spin pair. Since the ground energy
$E_{0l}<0$ for every $H_{l}$, we make a shift $H_{l}\rightarrow H_{l}-E_{0l}I$
for every term. We empirically pick $\Delta=15$ such that the low-energy
subspace contains eigenstates suitable.

We plot the Trotter error as a function of single-step evolution time at
different decay factors $\alpha$ in the right-hand side of Figure 4. The
Hamiltonian in (20) is decomposed into three interaction terms: the
horizontal, the vertical, and the remaining, respectively represented by
yellow, green, and red arrows as shown in the left-hand side of Figure 4. The
target evolution within each time step $V(\delta)$ is approximated by the
second-order product formula $S_{2}(\delta)$. We can observe that the
projection to the low-energy subspace significantly reduces the error compared
to the full Hilbert space for the same parameter configuration. Moreover, in
the full Hilbert space, the error increases as $\alpha$ grows, while in the
low-energy subspace, the opposite is true. This could be explained by the fact
that the norm of the Hamiltonian decreases as $\alpha$ increases. Therefore,
$\Delta=15$ becomes less stringent compared to the whole spectrum of the
Hamiltonian as $\alpha$ increases. For a fixed $\Delta$, the proportion of the
spectrum in the low-energy subspace increases for larger $\alpha$, which
results in a smaller improvement compared to the full Hilbert space
simulation.

Figure 4: Full Hilbert space vs. low-energy subspace: Trotter error induced by
the product formulas for a $\mathbf{3\times 3}$ power-law model. The solid
lines denote errors in the full Hilbert space while the dashed lines denote
errors in the low-energy subspace, and the plot distinguishes interaction
power $\alpha$ by the color of the lines.

##  4 A Lower Bound for Simulating in the Low-Energy Subspace

In the above discussion, we show that Hamiltonian simulation in the low-energy
subspace can provide improvement compared to full space simulation even with
imperfect state preparations. It is natural to ask what is the lower bound in
the resource requirement for low-energy simulations. Here, we provide a lower
bound on the query complexity to simulate any Hamiltonian $H$ for low-energy
simulations. In particular, we prove the following theorem:

######  Theorem 6 (Lower bound).

Given any positive number $D$, we can construct a $2$-sparse positive-
semidefinite Hamiltonian $H$ with $\norm{H}=\Theta(D)$ and zero minimal
eigenvalues such that simulating some state within the low-energy threshold
$\Delta$ and accuracy $\epsilon$ for some chosen scaled time $\tau$ requires
at least
$\Omega\left(\max\left\\{\tau,\frac{\log(\Delta/\epsilon)}{\log\log(\Delta/\epsilon)}\right\\}\right)$
queries to $H$.

###### Proof.

For the construction of the Hamiltonian, we adopt the idea of proving the
no–fast-forwarding theorem for computing parity using quantum computation [3,
17]. We start with the following Hamiltonian $H_{1}$ acting on
$(D+1)$-dimensional quantum states. We denote the basis of these quantum
states as $\ket{i}$ for $i\in\\{0,...,D\\}$. The nonzero entries for the
Hamiltonian are

$\displaystyle\bra{i+1}H_{1}\ket{i}=\bra{i}H_{1}\ket{i+1}=\frac{\sqrt{(D-i)(i+1)}}{2},\quad
i=0,...,D-1.$

We can observe that the Hamiltonian is actually an $J_{x}$ operator acting on
a spin-$D/2$ system. The operator norm of the Hamiltonian is exactly
$\frac{D}{2}$. There are two eigenstates corresponding to the eigenvalue
$-D/2$ and $D/2$, the former within which is the ground state. In addition, we
can find that $\bra{0}e^{iH_{1}t/2}\ket{D}=(\sin(t/2))^{D}$.

This Hamiltonian has been utilized as a component to construct lower bounds
for Hamiltonian simulation through the following construction $H_{2}$ [4, 5].
In particular, we consider the Hamiltonian $H_{2}$ generated by an $D$-bit
string $x=x_{1},...,x_{D}$ acting on the quantum state $\ket{i,j}$ with
$j\in\\{0,...,D\\}$ and $i\in\\{0,1\\}$. The nonzero entries of $H_{2}$ read:

$\displaystyle\bra{i,j}H_{2}\ket{i^{\prime},j+1}=\bra{i,j+1}H_{2}\ket{i^{\prime},j+1}=\frac{\sqrt{(D-i)(i+1)}}{2},\quad
j=0,...,D-1,$

where $i\oplus i^{\prime}=x_{i+1}$. This Hamiltonian, however, has a minimal
eigenvalue of $-D/2$. In order to make the Hamiltonian positive semidefinite,
we consider the following Hamiltonian

$\displaystyle H=H_{2}+\frac{D}{2}I,$

which adds an additional $e^{-i\Delta t}$ global phase to the evolution of any
input states. As we can always add back the inverse of a global phase to the
input state when we fix $t$, we simply ignore the effect of the term
$e^{-i\Delta t}$ in the following analysis. By definition, $\ket{0,0}$ is
connected to exactly one of $\ket{0,1}$ or $\ket{1,1}$. Moreover, it is
connected to the state $\ket{i,j}$ if and only if $x_{1}\oplus
x_{2}\oplus...\oplus x_{i}=j$. Therefore, $\ket{0,0}$ is exactly connected to
one of $\ket{D,0}$ or $\ket{D,1}$, determining the parity of the bit string
$x_{1}x_{2}...x_{D}$. The graph of this Hamiltonian contains two paths within
which the vertices are connected, one containing $\ket{0,0}$ and
$\ket{D,p(x)}$ for $p(x)$ the parity of $x$ and the other containing
$\ket{0,1}$ and $\ket{D,p(x)\oplus 1}$. Given a time $t$, the overlap between
$e^{iHt}\ket{0,0}$ and $\ket{D,p(x)}$ is
$\absolutevalue{\bra{D,p(x)}e^{-iHt}\ket{0,0}}=(\sin(t/2))^{D}$ while the
overlap between $e^{iHt}\ket{0,0}$ and $\ket{D,p(x)\oplus 1}$ is strictly
zero.

Notice that the Hamiltonian is $2$-sparse and we can only query at most two
values of $x_{i}$ by querying $H$ once. Therefore, simulating the evolution of
input state $\ket{0,0}$ for time $t=\pi$ yields an unbounded-error algorithm
for computing the parity of $x$, thus requiring $\Omega(D)$ queries when we
measure the output state in computational basis and consider the probability
of obtaining $\ket{D,p(x)}$ [4]. However, there are two issues with this
technique in our setting. The first issue is that the state $\ket{0,0}$ is not
in the low-energy subspace. It has energy $\bra{0,0}H\ket{0,0}=D/2$ due to the
second term $DI/2$. To address this issue, we consider the input state as a
mixed state of the minimal eigenstate $\ket{\psi_{\min}}$ and $\ket{0,0}$ as

$\displaystyle\rho_{\text{in}}=\left(1-\frac{2\Delta}{D}\right)\ket{\psi_{\min}}\bra{\psi_{\min}}+\frac{2\Delta}{D}\ket{0,0}\bra{0,0}.$

This state has energy exactly $\operatorname{Tr}(D\rho_{\text{in}})=\Delta$.
The second issue is that, in approximated simulation, we allow an $\epsilon$
error tolerance. Here, we use the techniques in [5] and consider the choice of
$D$ such that the overlap between the evolution of the second term
$\frac{2\Delta}{D}\ket{0,0}\bra{0,0}$ and $\ket{D,p(x)}$ the is larger than
$\epsilon$, which directly require
$\frac{2\Delta}{D}\cdot(\sin(t/2))^{D}\geq\epsilon$. Even if the simulation
error is allowed to be within $\epsilon$, we still have to figure out the
parity. In addition, the scaled time $\tau=\norm{H}t=Dt$. Thus we have

$\displaystyle\left(\sin\left(\frac{\tau}{2D}\right)\right)^{D}\leq\frac{\Delta\epsilon}{2D}.$

By taking the logarithm at two sides and using the approximation $\sin
x\approx x$, we have $D(\log
D-\log\tau)=O\left(\frac{1}{\Delta\epsilon}\right)$. Therefore, we have

$\displaystyle
D\geq\Omega\left(\max\left\\{\tau,\frac{\log(\Delta/\epsilon)}{\log\log(\Delta/\epsilon)}\right\\}\right).$

The lower bound for the number of queries to $H$ is thus
$\Omega(D)=\Omega\left(\max\left\\{\tau,\frac{\log(\Delta/\epsilon)}{\log\log(\Delta/\epsilon)}\right\\}\right)$,
which is what we expect. ∎

## Code Availability

Github repository: https://github.com/Qubit-Fernand/Digital-Quantum-
Simulation.

## Acknowledgements

We thank Andrew M. Childs, Alexey V. Gorshkov, Minh C. Tran, Chi-Fang Chen,
and Dong Yuan for helpful discussions. SZ and TL were supported by the
National Natural Science Foundation of China (Grant Numbers 62372006 and
92365117), and a startup fund from Peking University.

## References

  * [1] Itai Arad, Tomotaka Kuwahara, and Zeph Landau, _Connecting global and local energy distributions in quantum spin models on a lattice_ , Journal of Statistical Mechanics: Theory and Experiment 2016 (2016), no. 3, 033301, arXiv:1406.3898.
  * [2] Mari Carmen Bañuls, David A. Huse, and J. Ignacio Cirac, _Entanglement and its relation to energy variance for local one-dimensional hamiltonians_ , Physical Review B 101 (2020), no. 14, 144305, arXiv:1912.07639.
  * [3] Robert Beals, Harry Buhrman, Richard Cleve, Michele Mosca, and Ronald de Wolf, _Quantum lower bounds by polynomials_ , Journal of the ACM (JACM) 48 (2001), no. 4, 778–797, quant-ph/9802049.
  * [4] Dominic W. Berry, Graeme Ahokas, Richard Cleve, and Barry C. Sanders, _Efficient quantum algorithms for simulating sparse hamiltonians_ , Communications in Mathematical Physics 270 (2007), 359–371, arXiv:quant-ph/0508139.
  * [5] Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma, _Exponential improvement in precision for simulating sparse hamiltonians_ , Proceedings of the Forty-sixth Annual ACM Symposium on Theory of Computing, pp. 283–292, 2014, 1312.1414.
  * [6] Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma, _Simulating hamiltonian dynamics with a truncated taylor series_ , Physical Review Letters 114 (2015), no. 9, 090502, arXiv:1412.4687.
  * [7] Earl Campbell, _Random compiler for fast hamiltonian simulation_ , Physical Review Letters 123 (2019), 070503, arXiv:1811.08017.
  * [8] Chi-Fang Chen and Fernando G.S.L. Brandão, _Concentration for Trotter error_ , 2021, arXiv:2111.05324.
  * [9] Chi-Fang Chen, Hsin-Yuan Huang, Richard Kueng, and Joel A. Tropp, _Concentration for random product formulas_ , PRX Quantum 2 (2021), 040305, arXiv:2008.11751.
  * [10] Andrew M. Childs, _Quantum information processing in continuous time_ , Ph.D. thesis, Massachusetts Institute of Technology, Cambridge, MA, June 2004.
  * [11] Andrew M. Childs, Dmitri Maslov, Yunseong Nam, Neil J. Ross, and Yuan Su, _Toward the first quantum simulation with quantum speedup_ , Proceedings of the National Academy of Sciences 115 (2018), no. 38, 9456–9461, arXiv:1711.10980.
  * [12] Andrew M. Childs, Aaron Ostrander, and Yuan Su, _Faster quantum simulation by randomization_ , Quantum 3 (2019), 182, arXiv:1805.08385.
  * [13] Andrew M. Childs and Yuan Su, _Nearly optimal lattice simulation by product formulas_ , Physical Review Letters 123 (2019), no. 5, 050503, arXiv:1901.00564.
  * [14] Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe, and Shuchen Zhu, _Theory of Trotter error with commutator scaling_ , Physical Review X 11 (2021), 011020, arXiv:1912.08854.
  * [15] Andrew M. Childs and Nathan Wiebe, _Hamiltonian simulation using linear combinations of unitary operations_ , Quantum Information & Computation 12 (2012), no. 11-12, 901–924, arXiv:1202.5822.
  * [16] Chien Hung Cho, Dominic W. Berry, and Min-Hsiu Hsieh, _Doubling the order of approximation via the randomized product formula_ , 2022, arXiv:2210.11281.
  * [17] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Michael Sipser, _Limit on the speed of quantum computation in determining parity_ , Physical Review Letters 81 (1998), no. 24, 5442, quant-ph/9802045.
  * [18] Richard P. Feynman, _Simulating physics with computers_ , International Journal of Theoretical Physics 21 (1982), no. 6/7, 467–488.
  * [19] Yimin Ge, Jordi Tura, and J. Ignacio Cirac, _Faster ground state preparation and high-precision ground energy estimation with fewer qubits_ , Journal of Mathematical Physics 60 (2019), no. 2, 022202, arXiv:1712.03193.
  * [20] Iulia M. Georgescu, Sahel Ashhab, and Franco Nori, _Quantum simulation_ , Reviews of Modern Physics 86 (2014), no. 1, 153, arXiv:1308.6253.
  * [21] W. Gong, Yaroslav Kharkov, Minh C. Tran, Przemyslaw Bienias, and Alexey V. Gorshkov, _Improved digital quantum simulation by non-unitary channels_ , 2023, arXiv:2307.13028.
  * [22] De Huang, Jonathan Niles-Weed, Joel A. Tropp, and Rachel Ward, _Matrix concentration for products_ , Foundations of Computational Mathematics 22 (2022), no. 6, 1767–1799, arXiv:2003.05437.
  * [23] Jacky Huyghebaert and Hans De Raedt, _Product formula methods for time-dependent Schrödinger problems_ , Journal of Physics A: Mathematical and General 23 (1990), no. 24, 5777.
  * [24] Tarun Kathuria, Satyaki Mukherjee, and Nikhil Srivastava, _On concentration inequalities for random matrix products_ , 2020, arXiv:2003.06319.
  * [25] Youngseok Kim, Andrew Eddins, Sajant Anand, Ken Xuan Wei, Ewout Van Den Berg, Sami Rosenblatt, Hasan Nayfeh, Yantao Wu, Michael Zaletel, Kristan Temme, and Abhinav Kandala, _Evidence for the utility of quantum computing before fault tolerance_ , Nature 618 (2023), no. 7965, 500–505.
  * [26] Guang Hao Low and Isaac L. Chuang, _Optimal hamiltonian simulation by quantum signal processing_ , Physical Review Letters 118 (2017), no. 1, 010501, arXiv:1606.02685.
  * [27] Guang Hao Low and Isaac L. Chuang, _Hamiltonian simulation by qubitization_ , Quantum 3 (2019), 163, arXiv:1610.06546.
  * [28] Guang Hao Low, Vadym Kliuchnikov, and Nathan Wiebe, _Well-conditioned multiproduct Hamiltonian simulation_ , 2019, arXiv:1907.11679.
  * [29] Sirui Lu, Mari Carmen Banuls, and J. Ignacio Cirac, _Algorithms for quantum simulation at finite energies_ , PRX Quantum 2 (2021), no. 2, 020321, arXiv:2006.03032.
  * [30] John M. Martyn, Zane M. Rossi, Andrew K. Tan, and Isaac L. Chuang, _Grand unification of quantum algorithms_ , PRX Quantum 2 (2021), no. 4, 040203, arXiv:2105.02859.
  * [31] Willard Miller, _Symmetry groups and their applications_ , Academic Press, 1973.
  * [32] Roberto Imbuzeiro Oliveira, _The spectrum of random k-lifts of large graphs (with possibly large k)_ , 2009, arXiv:0911.4741.
  * [33] Iosif Pinelis, _Optimum bounds for the distributions of martingales in banach spaces_ , The Annals of Probability (1994), 1679–1706.
  * [34] John Preskill, _Quantum computing in the NISQ era and beyond_ , Quantum 2 (2018), 79, arXiv:1801.00862.
  * [35] Burak Şahino$\mathbf{u}$glu and Rolando D. Somma, _Hamiltonian simulation in the low-energy subspace_ , npj Quantum Information 7 (2021), no. 1, 119, arXiv:2006.02660.
  * [36] Masuo Suzuki, _Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics and statistical physics_ , Journal of Mathematical Physics 26 (1985), no. 4, 601–612.
  * [37] Masuo Suzuki, _Fractal decomposition of exponential operators with applications to many-body theories and monte carlo simulations_ , Physics Letters A 146 (1990), no. 6, 319–323.
  * [38] Masuo Suzuki, _General theory of fractal path integrals with applications to many-body theories and statistical physics_ , Journal of Mathematical Physics 32 (1991), no. 2, 400–407.
  * [39] Minh C. Tran, Yuan Su, Daniel Carney, and Jacob M. Taylor, _Faster digital quantum simulation by symmetry protection_ , PRX Quantum 2 (2021), no. 1, 010323, arXiv:2006.16248.
  * [40] Joel A. Tropp, _Freedman’s inequality for matrix martingales_ , Electronic Communications in Probability 16 (2011), no. 25, 262–270, arXiv:1101.3039.
  * [41] Joel A. Tropp, _An introduction to matrix concentration inequalities_ , Foundations and Trends® in Machine Learning 8 (2015), no. 1-2, 1–230, arXiv:1501.01571.
  * [42] Qi Zhao, You Zhou, Alexander F. Shaw, Tongyang Li, and Andrew M. Childs, _Hamiltonian simulation with random inputs_ , Physical Review Letters 129 (2022), no. 27, 270502, arXiv:2111.04773.
  * [43] Sergiy Zhuk, Niall Robertson, and Sergey Bravyi, _Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation_ , 2023, arXiv:2306.12569.

##  Appendix A Auxiliary Lemmas

In the following, we carefully pick
$\Delta^{\prime}\geq\Lambda^{\prime}\geq\Lambda\geq 0$ as some value depends
on $\Delta$ and $\Delta^{\prime}\geq\Delta$. We denote
$\overline{H}=\Pi_{\leq\Delta^{\prime}}H\Pi_{\leq\Delta^{\prime}}$ and
$\overline{H}_{l}=\Pi_{\leq\Delta^{\prime}}H_{l}\Pi_{\leq\Delta^{\prime}}$.
Therefore,
$\norm{\overline{H}_{l}}\leq\norm{\overline{H}}\leq\Delta^{\prime}$.

######  Lemma 7 (Lemma 1 of [35]).

For $k$-local Hamiltonian $H=\sum_{l=1}^{L}H_{l}$ with $H_{l}\geq 0$ and
$\Lambda^{\prime}\geq\Lambda\geq 0$. Then, for all $s\in\mathbb{R}$ and
$l=1,\ldots,L$,

$\displaystyle\norm{\Pi_{>\Lambda^{\prime}}e^{-i\delta
H_{l}}\Pi_{\leq\Lambda}}\leq
e^{-\alpha_{1}(\Lambda^{\prime}-\Lambda)/J}(e^{\alpha_{2}J\absolutevalue{\delta}n}-1),$

where $\alpha_{1}>0$ and $\alpha_{2}>0$ are constants.

######  Lemma 8 (Lemma 2 of [35]).

For $k$-local Hamiltonian $H=\sum_{l=1}^{L}H_{l}$ with $H_{l}\geq 0$ and
$\Delta^{\prime}\geq\Lambda^{\prime}\geq\Lambda\geq 0$. Then, for all
$\delta\in\mathbb{R}$ and $l=1,\ldots,L$,

$\displaystyle\norm{\Pi_{\leq\Lambda^{\prime}}(e^{-i\delta\overline{H}_{l}}-e^{-i\delta
H_{l}})\Pi_{\leq\Lambda}}$ $\displaystyle\leq
e^{-\alpha_{1}(\Lambda^{\prime}-\Lambda)/J}(e^{\alpha_{2}J\absolutevalue{\delta}n}-1),$
$\displaystyle\norm{\Pi_{>\Lambda^{\prime}}e^{-i\delta\overline{H}_{l}}\Pi_{\leq\Lambda}}$
$\displaystyle\leq
3e^{-\alpha_{1}(\Lambda^{\prime}-\Lambda)/J}(e^{\alpha_{2}J\absolutevalue{\delta}n}-1),$

where $\alpha_{1}>0$ and $\alpha_{2}>0$ are constants.

Based on these two lemmas, we consider approximating errors for product
formulas in the low-energy subspaces. We first give the following generic
product formulas of $q>1$ terms:

$\displaystyle W(\bm{\delta})$ $\displaystyle=e^{-i\delta_{q}H_{l_{q}}}\cdots
e^{-i\delta_{1}H_{l_{1}}},$ $\displaystyle\overline{W}(\bm{\delta})$
$\displaystyle=e^{-i\delta_{q}\overline{H}_{l_{q}}}\cdots
e^{-i\delta_{1}\overline{H}_{l_{1}}},$

where $\bm{\delta}=\delta_{1},\ldots,\delta_{q},\delta_{j}\in\mathbb{R}$ and
$1\leq l_{j}\leq L$. Given
${\bm{\Lambda}}=(\Lambda_{1},...,\Lambda_{q}),\Lambda_{j}\geq 0$, we define

$\displaystyle W^{\bm{\Lambda}}(\bm{\delta})$
$\displaystyle=\Pi_{\leq\Lambda_{q}}e^{-i\delta_{q}H_{l_{q}}}\cdots\Pi_{\leq\Lambda_{1}}e^{-i\delta_{1}H_{l_{1}}},$
$\displaystyle\overline{W}^{\bm{\Lambda}}(\bm{\delta})$
$\displaystyle=\Pi_{\leq\Lambda_{q}}e^{-i\delta_{q}\overline{H}_{l_{q}}}\cdots\Pi_{\leq\Lambda_{1}}e^{-i\delta_{1}\overline{H}_{l_{1}}}.$

Regarding these product formulas, we have the following lemma:

######  Lemma 9 (Corollaries 1, 2, 3, and 4 of [35]).

Let $\tau>0$, $\Delta\geq 0$, $\lambda=(2Jdk)^{-1}$, and $\alpha=eJ$. Then, if
${\bm{\Lambda}}$ satisfies
$\Lambda_{j}-\Lambda_{j-1}\geq\frac{1}{\lambda}(\alpha\absolutevalue{\delta_{j}}M+\log(q/\tau))$,
$\Lambda=\Delta$, and $\Delta^{\prime}\geq\Lambda_{q}$, we have

$\displaystyle\norm{(W^{\bm{\Lambda}}(\bm{\delta})-W(\bm{\delta}))\Pi_{\leq\Delta}}\leq\tau,$
$\displaystyle\norm{(\overline{W}^{\bm{\Lambda}}(\bm{\delta})-W^{\bm{\Lambda}}(\bm{\delta}))\Pi_{\leq\Delta}}\leq\tau,$
$\displaystyle\norm{(\overline{W}(\bm{\delta})-\overline{W}^{\bm{\Lambda}}(\bm{\delta}))\Pi_{\leq\Delta}}\leq
3\tau.$

Combing the three inequalities here, we obtain the result. For $\tau>0$,
$\Delta\geq 0$, $\lambda=(2Jdk)^{-1}$, $\alpha=eJ$, and
$\absolutevalue{s}=\sum_{j=1}^{q}\absolutevalue{\delta_{j}}$. Then, if
$\Delta^{\prime}\geq\Delta+\frac{1}{\lambda}(\alpha\absolutevalue{\bm{\delta}}M+q\log(q/\tau))$,
we have

$\displaystyle\norm{(\overline{W}(\bm{\delta})-W(\bm{\delta}))\Pi_{\leq\Delta}}\leq
5\tau.$

######  Lemma 10 (Corollary 3.4 of [9]).

Given a martingale $\\{B_{k}:k=0,\ldots,r\\}\in\mathbb{C}_{d\times d}$ with
$d=2^{n}$ and assume that $C_{k}=B_{k}-B_{k-1}$ satisfies $\norm{C_{k}}\leq R$
and $\norm{\sum_{k-1}^{r}\mathbb{E}_{k-1}C_{k}C_{K}^{\dagger}}\leq v$. Here
$\mathbb{E}_{k-1}$ is the average conditioned on $B_{k-1},\ldots,B_{1}$. Then,
we have

$\displaystyle\Pr[\norm{B_{r}-B_{0}}\geq\epsilon]\leq
2d\exp\left(\frac{3\epsilon^{2}}{6v+2R\epsilon}\right).$