aff1]Department of Computer Science, University of Maryland aff2]Institute for
Advanced Computer Studies, University of Maryland aff3]Joint Center for
Quantum Information and Computer Science, University of Maryland
aff4]Department of Physics, University of Maryland

# Faster quantum simulation by randomization

Andrew M. Childs  Aaron Ostrander  Yuan Su  [  [  [  [

###### Abstract

Product formulas can be used to simulate Hamiltonian dynamics on a quantum
computer by approximating the exponential of a sum of operators by a product
of exponentials of the individual summands. This approach is both
straightforward and surprisingly efficient. We show that by simply randomizing
how the summands are ordered, one can prove stronger bounds on the quality of
approximation for product formulas of any given order, and thereby give more
efficient simulations. Indeed, we show that these bounds can be asymptotically
better than previous bounds that exploit commutation between the summands,
despite using much less information about the structure of the Hamiltonian.
Numerical evidence suggests that the randomized approach has better empirical
performance as well.

##  1 Introduction

Simulating quantum dynamics is one of the major potential applications of
quantum computers. The apparent intractability of simulating quantum dynamics
with a classical computer led Feynman [19] and others to propose the idea of
quantum computation. Lloyd gave the first explicit quantum algorithm for
simulating the dynamics of local Hamiltonians [25], and later work showed that
the more general class of sparse Hamiltonians can also be simulated
efficiently [1]. Quantum simulation can be applied to understand the behavior
of various physical systems—including many-body physics [31], quantum
chemistry [36, 30, 2], and quantum field theory [23]—and designing new quantum
algorithms [16, 18, 21, 9, 11].

The main ingredient in Lloyd’s algorithm is the Lie product formula, which
provides a first-order approximation to the exponential of a sum as a product
of exponentials of the summands. Given Hermitian operators
$H_{1},\ldots,H_{L}$ (which we refer to as the summands of the Hamiltonian
$H=\sum_{j=1}^{L}H_{j}$) and a complex number $\lambda$, the Lie product
formula

$S_{1}(\lambda):=\prod_{j=1}^{L}\exp(\lambda H_{j}),$ (1)

approximates the exponentiation

$V(\lambda):=\exp\biggl{(}\lambda\sum_{j=1}^{L}H_{j}\biggr{)}$ (2)

in the sense that $V(\lambda)\approx S_{1}({\lambda}/{r})^{r}$ for large $r$.
Suzuki systematically extended this formula to give a $(2k)$th-order
approximation $S_{2k}$, defined recursively by

$\displaystyle S_{2}(\lambda)$
$\displaystyle:=\prod_{j=1}^{L}\exp\biggl{(}\frac{\lambda}{2}H_{j}\biggr{)}\prod_{j=L}^{1}\exp\biggl{(}\frac{\lambda}{2}H_{j}\biggr{)}$
(3) $\displaystyle S_{2k}(\lambda)$
$\displaystyle:=S_{2k-2}(p_{k}\lambda)^{2}\,S_{2k-2}((1-4p_{k})\lambda)\,S_{2k-2}(p_{k}\lambda)^{2}$

with $p_{k}:=1/(4-4^{1/(2k-1)})$ [33]. Again we have $V(\lambda)\approx
S_{2k}({\lambda}/{r})^{r}$ for large $r$, and the approximation obtained with
a given value of $r$ improves as $k$ increases (albeit with a prefactor that
grows exponentially in $k$). We refer to all such formulas as _product
formulas_. When they are used for quantum simulation, $H$ is chosen to be the
Hamiltonian and $\lambda=-it$, where $t$ is the evolution time. Although other
approaches to quantum simulation have better proven asymptotic performance as
a function of various parameters [8, 4, 6, 7, 26, 27, 20], product formulas
perform well in practice [17] and are widely used in experimental
implementations [3, 12, 24] due to their simplicity and the fact that they do
not require any ancilla qubits.

The main challenge in applying product formulas to quantum simulation is to
choose the number of segments $r$ to ensure the simulation error is at most
some allowed threshold $\epsilon$. To simulate $H=\sum_{j=1}^{L}H_{j}$ for
time $t$, rigorous error analysis shows that

$r_{1,\text{det}}=O\biggl{(}\frac{(t\Lambda L)^{2}}{\epsilon}\biggr{)}$ (4)

suffices to ensure error at most $\epsilon$ for the first-order formula and

$r_{2k,\text{det}}=O\biggl{(}\frac{(t\Lambda
L)^{1+\frac{1}{2k}}}{\epsilon^{\frac{1}{2k}}}\biggr{)}$ (5)

suffices for $(2k)$th order [5], where $\Lambda:=\max_{j}\left\lVert
H_{j}\right\rVert$ is a spectral-norm upper bound on the summands and det
indicates that these formulas are constructed deterministically. However,
numerical simulations suggest that the product formula algorithm can perform
significantly better in practice than the best proven error bounds demonstrate
[2, 31, 32, 17]. Indeed, recent work suggests that it can even asymptotically
outperform more sophisticated simulation algorithms with better proven running
times [17]. This dramatic gap between the provable and the actual behavior of
product formula simulation suggests that it may be possible to significantly
improve their analysis, and thereby give more efficient algorithms for quantum
simulation.

It is sometimes possible to improve the analysis of product formulas using
further information about the form of the Hamiltonian. In particular, the cost
of simulation can be reduced when many pairs of summands commute [25, 2, 17].
However, this approach can only be applied for structured Hamiltonians that
contain many commuting summands. Furthermore, the best known bounds of this
type give only modest improvement, remaining orders of magnitude away from the
empirical performance even in cases where many summands commute [17].

Randomization can be a powerful tool for improving the performance of quantum
simulation algorithms. For example, Poulin et al. gave improved simulations of
time-dependent Hamiltonians by sampling the Hamiltonian at random times [29].
Closer in spirit to the present paper, Zhang studied the effect of randomizing
the ordering and/or duration of evolutions in a product formula, showing in
particular that randomly ordering the summands in the first-order formula in
either forward or reverse order can give an improved algorithm [37].

In this paper, we explore a closely related approach for higher-order product
formulas, which can achieve significantly better asymptotic performance.
Specifically, we analyze the effect of randomly permuting the summands. The
resulting algorithm is not much more complicated than a deterministic product
formula, but the savings in the simulation cost are substantial. For any
permutation $\sigma\in\operatorname{Sym}(L)$ of the $L$ summands, let

$\displaystyle S_{2}^{\sigma}(\lambda)$
$\displaystyle:=\prod_{j=1}^{L}\exp\biggl{(}\frac{\lambda}{2}H_{\sigma(j)}\biggr{)}\prod_{j=L}^{1}\exp\biggl{(}\frac{\lambda}{2}H_{\sigma(j)}\biggr{)}$
(6) $\displaystyle S_{2k}^{\sigma}(\lambda)$
$\displaystyle:=[S_{2k-2}^{\sigma}(p_{k}\lambda)]^{2}S_{2k-2}^{\sigma}((1-4p_{k})\lambda)[S_{2k-2}^{\sigma}(p_{k}\lambda)]^{2}.$

We show that the $(2k)$th-order randomized simulation has error

$\left\lVert\mathcal{V}(-it)-\biggl{(}\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\mathcal{S}_{2k}^{\sigma}\bigl{(}-it/r\bigr{)}\biggr{)}^{r}\right\rVert_{\diamond}=O\biggl{(}\frac{(\Lambda
tL)^{4k+2}}{r^{4k+1}}+\frac{(\Lambda t)^{2k+1}L^{2k}}{r^{2k}}\biggr{)}.$ (7)

where $\mathcal{V}(-it)$ and $\mathcal{S}_{2k}^{\sigma}(-it/r)$ are quantum
channels describing the unitary transformation $V(-it)$ and the random unitary
$S_{2k}^{\sigma}(-it/r)$, respectively, and
$\left\lVert\cdot\right\rVert_{\diamond}$ is the diamond norm (defined in
Section 2).

Our analysis uses a mixing lemma of Campbell and Hastings [22, 13] to bound
the diamond norm distance from the ideal evolution. (Even for the first-order
case, this improves over the analysis of Zhang, which uses similar methods but
only bounds the trace distance from the ideal final state [37], a metric that
does not account for entanglement with a reference system.) Informally, the
lemma of [22, 13] states that if we can approximate a desired operation as the
average over some set of operations, then the overall error depends linearly
on the error in the average operation but only quadratically on the error in
any individual operation. Standard error bounds for product formulas do not
depend on how the summands are ordered, but we show that randomizing the
ordering gives a more accurate average evolution. We motivate this approach in
Section 2, where we consider the effect of randomizing how the summands are
ordered in the simple case of the first-order formula. Assuming
$\Lambda:=\max_{j}\left\lVert H_{j}\right\rVert$ is constant, the randomized
first-order algorithm has gate complexity
$g_{1}^{\text{rand}}=O\big{(}t^{1.5}L^{2.5}/\epsilon^{0.5}\big{)}$, improving
over $g_{1}^{\det}=O\bigl{(}t^{2}L^{3}/\epsilon\bigr{)}$ in the deterministic
case.

Analyzing the effect of randomization on higher-order formulas is more
challenging. For terms of order at most $L$ in the Taylor expansion of a
product formula, the majority of the error comes from terms in which no
summands are repeated. We call such contributions _nondegenerate terms_. In
Section 3, we give a combinatorial argument to compute nondegenerate terms of
the average evolution
$\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}S_{2k}^{\sigma}(\lambda)$ in
closed form. (In fact, we prove a more general result that applies to the
average evolution as a special case.) As a corollary, we show that the
nondegenerate terms completely cancel in the randomized product formula.

Section 4 presents our main technical result, an upper bound on the error in a
randomized higher-order product formula simulation. This bound follows by
using the mixing lemma to combine an error bound for the average evolution
operator with standard product formula error bounds for the error of the
individual terms. Section 5 discusses the overall performance of the resulting
algorithm and compares it with deterministic approaches. For the $(2k)$th-
order product formula, assuming $\Lambda:=\max_{j}\left\lVert
H_{j}\right\rVert$ is constant, our randomized Hamiltonian simulation
algorithm has complexity

$g_{2k}^{\text{rand}}=\max\bigg{\\{}O\biggl{(}tL^{2}\biggl{(}\frac{tL}{\epsilon}\biggr{)}^{\frac{1}{4k+1}}\biggr{)},O\biggl{(}tL^{2}\biggl{(}\frac{t}{\epsilon}\biggr{)}^{\frac{1}{2k}}\biggr{)}\bigg{\\}},$
(8)

compared to
$g_{2k}^{\det}=O\bigl{(}tL^{2}(tL/\epsilon)^{\frac{1}{2k}}\bigr{)}$ in the
deterministic case. Thus our algorithm always improves the dependence on $L$
and sometimes achieves better dependence on $t$ and $\epsilon$ as well.

We also show in Section 5 that our bound can outperform a previous bound that
takes advantage of the structure of the Hamiltonian. Specifically, we compare
our randomized product formula algorithm with the deterministic algorithm
using the commutator bound of [17] for a one-dimensional Heisenberg model in a
random magnetic field. We find that over a significant range of parameters,
the randomized algorithm has better proven performance, despite using less
information about the form of the Hamiltonian.

In light of the large gap between proven and empirical performance of product
formulas, it is natural to ask whether randomized product formulas still offer
an improvement under the best possible error bounds. To address this question,
we present numerical comparisons of the deterministic and randomized product
formulas in Section 6. In particular, we show that the randomized approach can
sometimes outperform the deterministic approach even with respect to their
empirical performance.

Finally, we conclude in Section 7 with a brief discussion of the results and
some open questions.

##  2 The power of randomization

To see how randomness can improve a product formula simulation, consider a
simple Hamiltonian expressed as a sum of two operators, $H=H_{1}+H_{2}$. The
Taylor expansion of the first-order formula as a function of
$\lambda\in\mathbb{C}$ is

$S_{1}(\lambda)=\exp(\lambda H_{1})\exp(\lambda
H_{2})=I+\lambda(H_{1}+H_{2})+\frac{\lambda^{2}}{2}(H_{1}^{2}+2H_{1}H_{2}+H_{2}^{2})+O(\lambda^{3}),$
(9)

whereas the Taylor series of the ideal evolution is

$V(\lambda)=\exp((H_{1}+H_{2})\lambda)=I+\lambda(H_{1}+H_{2})+\frac{\lambda^{2}}{2}(H_{1}^{2}+H_{1}H_{2}+H_{2}H_{1}+H_{2}^{2})+O(\lambda^{3}).$
(10)

Using the triangle inequality, we can bound the spectral-norm error as

$\left\lVert
V(\lambda)-S_{1}(\lambda)\right\rVert\leq\left\lVert[H_{1},H_{2}]\right\rVert\frac{|\lambda|^{2}}{2}+O((\Lambda|\lambda|)^{3}),$
(11)

where $\Lambda:=\max\\{\left\lVert H_{1}\right\rVert,\left\lVert
H_{2}\right\rVert\\}$. Since $H_{1}$ and $H_{2}$ need not commute,
$S_{1}(\lambda)$ approximates $V(\lambda)$ to first order in $\lambda$, as
expected.

It is clearly impossible to approximate $V(\lambda)$ to second order using a
product of only two exponentials of $H_{1}$ and $H_{2}$: any such product can
have only one of the products $H_{1}H_{2}$ and $H_{2}H_{1}$ in its Taylor
expansion, whereas $V(\lambda)$ contains both of these products in its second-
order term. However, we can obtain both products by taking a uniform mixture
of $S_{1}(\lambda)$ and

$S_{1}^{\mathrm{rev}}(\lambda):=\exp(\lambda H_{2})\exp(\lambda H_{1}).$ (12)

Indeed, a simple calculation shows that

$\left\lVert
V(\lambda)-\frac{1}{2}\bigl{(}S_{1}(\lambda)+S_{1}^{\mathrm{rev}}(\lambda)\bigr{)}\right\rVert=O\bigl{(}(\Lambda|\lambda|)^{3}\bigr{)}.$
(13)

However, $\bigl{(}S_{1}(-it)+S_{1}^{\mathrm{rev}}(-it)\bigr{)}/2$ is not a
unitary operation in general. We could in principle implement a linear
combination of unitaries using the techniques of [6], but such an approach
would use ancillas and could have high cost, especially when the Hamiltonian
contains many summands. A simpler approach is to apply one of the two
operations $S_{1}(-it)$ and $S_{1}^{\mathrm{rev}}(-it)$ chosen uniformly at
random (as in Algorithm 2 of [37]), thereby implementing a quantum channel
that gives a good approximation to the desired evolution.

We now introduce some notation that is useful to analyze the performance of
randomized product formulas. Let $X$ be a matrix acting on a finite-
dimensional Hilbert space $\mathcal{H}$. We write $\left\lVert X\right\rVert$
for its spectral norm (the largest singular value) and $\left\lVert
X\right\rVert_{1}$ for its trace norm (the sum of its singular values, i.e.,
its Schatten 1-norm). Let $\mathcal{E}\colon X\mapsto\mathcal{E}(X)$ be a
linear map on the space of matrices on $\mathcal{H}$. The diamond norm of
$\mathcal{E}$ is

$\left\lVert\mathcal{E}\right\rVert_{\diamond}:=\max\\{\left\lVert(\mathcal{E}\otimes\mathds{1}_{\mathcal{H}})(Y)\right\rVert_{1}:\left\lVert
Y\right\rVert_{1}\leq 1\\},$ (14)

where the maximization is taken over all matrices $Y$ on
$\mathcal{H}\otimes\mathcal{H}$ satisfying $\left\lVert Y\right\rVert_{1}\leq
1$.

The following mixing lemma bounds how well we can approximate a unitary
operation using a random unitary channel. Specifically, the error is linear in
the distance between the target unitary and the average of the random
unitaries, and only quadratic in the distance between the target unitary and
each individual random unitary.

######  Lemma 1 (Mixing lemma [13, 22]).

Let $V$ and $\\{U_{j}\\}$ be unitary matrices, with associated quantum
channels $\mathcal{V}\colon\rho\mapsto V\rho V^{\dagger}$ and
$\mathcal{U}_{j}\colon\rho\mapsto U_{j}\rho U_{j}^{\dagger}$, and let
$\\{p_{j}\\}$ be a collection of positive numbers satisfying
$\sum_{j}p_{j}=1$. Suppose that

  1. (i)

$\left\lVert U_{j}-V\right\rVert\leq a$ for all $j$ and

  2. (ii)

$\bigl{\|}(\sum_{j}p_{j}U_{j})-V\bigr{\|}\leq b$.

Then the average evolution $\mathcal{E}:=\sum_{j}p_{j}\mathcal{U}_{j}$
satisfies $\left\lVert\mathcal{E}-\mathcal{V}\right\rVert_{\diamond}\leq
a^{2}+2b$.

To simulate the Hamiltonian $H=H_{1}+H_{2}$ for time $t$, we divide the
evolution into $r$ segments of duration $t/r$ and implement each segment via
the random unitary operation

$\frac{1}{2}\bigl{(}\mathcal{S}_{1}(-it/r)+\mathcal{S}_{1}^{\mathrm{rev}}(-it/r)\bigr{)}$
(15)

using one bit of randomness per segment, where $\mathcal{S}_{1}$ and
$\mathcal{S}_{1}^{\mathrm{rev}}$ are the quantum channels associated with
$S_{1}$ and $S_{1}^{\mathrm{rev}}$. Invoking the mixing lemma with
$a=O\bigl{(}{(\Lambda t)^{2}}/{r^{2}}\bigr{)}$ and $b=O\bigl{(}{(\Lambda
t)^{3}}/{r^{3}}\bigr{)}$, we find that

$\left\lVert\mathcal{V}(-it/r)-\frac{1}{2}\bigl{(}\mathcal{S}_{1}(-it/r)+\mathcal{S}_{1}^{\mathrm{rev}}(-it/r)\bigr{)}\right\rVert_{\diamond}=O\bigg{(}\frac{(\Lambda
t)^{3}}{r^{3}}\bigg{)}.$ (16)

Since the diamond norm distance between quantum channels is subadditive under
composition [35, p. 178], the error of the entire simulation is

$\left\lVert\mathcal{V}(-it)-\frac{1}{2^{r}}\bigl{(}\mathcal{S}_{1}(-it/r)+\mathcal{S}_{1}^{\mathrm{rev}}(-it/r)\bigr{)}^{r}\right\rVert_{\diamond}=O\bigg{(}\frac{(\Lambda
t)^{3}}{r^{2}}\bigg{)}.$ (17)

Thus the randomized first-order formula is effectively a second-order formula.

This approach easily extends to a sum of $L$ operators, again effectively
making the first-order formula accurate to second order (cf. [37], which shows
the same result with respect to trace distance of the output state). Keeping
track of all the prefactors, we find the following error bound for the
randomized first-order formula.

######  Theorem 1 (Randomized first-order error bound).

Let $\\{H_{j}\\}_{j=1}^{L}$ be Hermitian matrices. Let

$V(-it):=\exp\biggl{(}-it\sum_{j=1}^{L}H_{j}\biggr{)}$ (18)

be the evolution induced by the Hamiltonian $H=\sum_{j=1}^{L}H_{j}$ for time
$t\in\mathbb{R}$. Define

$\displaystyle S_{1}(\lambda)$ $\displaystyle:=\prod_{j=1}^{L}\exp(\lambda
H_{j})$ $\displaystyle and$ $\displaystyle S_{1}^{\mathrm{rev}}(\lambda)$
$\displaystyle:=\prod_{j=L}^{1}\exp(\lambda H_{j}).$ (19)

Let $r\in\mathbb{N}$ be a positive integer and $\Lambda:=\max\left\lVert
H_{j}\right\rVert$. Then

$\left\lVert\mathcal{V}(-it)-\frac{1}{2^{r}}\bigl{(}\mathcal{S}_{1}(-it/r)+\mathcal{S}_{1}^{\mathrm{rev}}(-it/r)\bigr{)}^{r}\right\rVert_{\diamond}\leq\frac{(\Lambda|t|L)^{4}}{r^{3}}\exp\biggl{(}2\frac{\Lambda|t|L}{r}\biggr{)}+\frac{2(\Lambda|t|L)^{3}}{3r^{2}}\exp\biggl{(}\frac{\Lambda|t|L}{r}\biggr{)}$
(20)

where, for $\lambda=-it$, we associate channels $\mathcal{V}(\lambda)$,
$\mathcal{S}_{1}(\lambda)$, and $\mathcal{S}_{1}^{\mathrm{rev}}(\lambda)$ with
the unitaries $V(\lambda)$, $S_{1}(\lambda)$, and
$S_{1}^{\mathrm{rev}}(\lambda)$, respectively.

To guarantee that the simulation error is at most $\epsilon$, we upper bound
the right-hand side of (20) by $\epsilon$ and solve for $r$. Assuming
$\Lambda:=\max_{j}\left\lVert H_{j}\right\rVert$ is constant, we find that it
suffices to choose
$r_{1}^{\mathrm{rand}}=O\bigl{(}(tL)^{1.5}/\epsilon^{0.5}\bigr{)}$, giving a
simulation algorithm with gate complexity
$g_{1}^{\mathrm{rand}}=O\bigl{(}t^{1.5}L^{2.5}/\epsilon^{0.5}\bigr{)}$. In
comparison, the gate complexity in the deterministic case is
$g_{1}^{\det}=O\bigl{(}t^{2}L^{3}/\epsilon\bigr{)}$. Therefore, the randomized
first-order product formula algorithm improves over the deterministic
algorithm with respect to all parameters of interest.

It is natural to ask whether a similar randomization strategy can improve
higher-order product formulas (as defined in (3)). While it turns out that
randomization does not improve the order of the formula, it does result in a
significant reduction of the error, and in particular, lowers the dependence
on the number of summands in the Hamiltonian. The more complicated structure
of higher-order formulas makes this analysis more involved than in the first-
order case (in particular, we randomly permute the $L$ summands instead of
simply choosing whether or not to reverse them, so we use $\Theta(L\log L)$
bits of randomness per segment instead of only a single bit). As discussed at
the end of Section 1, our proof is based on a randomization lemma (established
in the next section) that evaluates the dominant contribution to the Taylor
series of the randomized product formula in closed form.

##  3 Randomization lemma

In this section, we study the Taylor expansion of the average evolution
operator obtained by randomizing how the summands of a Hamiltonian are
ordered. We consider a formula of the form

$\displaystyle\exp(q_{1}\lambda H_{\pi_{1}(1)})\exp(q_{1}\lambda
H_{\pi_{1}(2)})\cdots\exp(q_{1}\lambda H_{\pi_{1}(L)})$ (21)
$\displaystyle\exp(q_{2}\lambda H_{\pi_{2}(1)})\exp(q_{2}\lambda
H_{\pi_{2}(2)})\cdots\exp(q_{2}\lambda H_{\pi_{2}(L)})$ $\displaystyle\cdots$
$\displaystyle\exp(q_{\kappa}\lambda
H_{\pi_{\kappa}(1)})\exp(q_{\kappa}\lambda
H_{\pi_{\kappa}(2)})\cdots\exp(q_{\kappa}\lambda H_{\pi_{\kappa}(L)})$

for real numbers $q_{1},\ldots,q_{\kappa}\in\mathbb{R}$, a complex number
$\lambda\in\mathbb{C}$, Hermitian matrices $H_{1},\ldots,H_{L}$, and
permutations $\pi_{1},\ldots,\pi_{\kappa}\in\operatorname{Sym}(L)$. By
choosing appropriate values of $q_{1},\ldots,q_{\kappa}\in\mathbb{R}$ and
ordering $H_{1},\ldots,H_{L}$ in both forward and backward directions, we can
write any product formula $S_{2k}(\lambda)$ in this form.

We now permute the summands to get the average evolution

$\displaystyle\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}$
$\displaystyle\exp(q_{1}\lambda H_{\sigma(\pi_{1}(1))})\exp(q_{1}\lambda
H_{\sigma(\pi_{1}(2))})\cdots\exp(q_{1}\lambda H_{\sigma(\pi_{1}(L))})$ (22)
$\displaystyle\exp(q_{2}\lambda H_{\sigma(\pi_{2}(1))})\exp(q_{2}\lambda
H_{\sigma(\pi_{2}(2))})\cdots\exp(q_{2}\lambda H_{\sigma(\pi_{2}(L))})$
$\displaystyle\cdots$ $\displaystyle\exp(q_{\kappa}\lambda
H_{\sigma(\pi_{\kappa}(1))})\exp(q_{\kappa}\lambda
H_{\sigma(\pi_{\kappa}(2))})\cdots\exp(q_{\kappa}\lambda
H_{\sigma(\pi_{\kappa}(L))}).$

In its Taylor expansion, we call the sum of the form

$\sum_{\begin{subarray}{c}m_{1},\ldots,m_{s}\\\ \text{pairwise
different}\end{subarray}}\alpha_{m_{1}\ldots m_{s}}\lambda^{s}H_{m_{1}}\cdots
H_{m_{s}},$ (23)

with coefficients $\alpha_{m_{1}\ldots m_{s}}\in\mathbb{C}$, the $s$th-order
_nondegenerate term_. This term contributes $\Theta(L^{s})$ to the $s$th-order
error, whereas the remaining (degenerate) terms only contribute $O(L^{s-1})$.

The following lemma shows how to compute the $s$th-order nondegenerate term
for an arbitrary average evolution.

######  Lemma 2 (Randomization lemma).

Define an average evolution operator as in (22) and let $s\leq L$ be a
positive integer. The $s$th-order nondegenerate term of this operator is

$\frac{[(q_{1}+\cdots+q_{\kappa})\lambda]^{s}}{s!}\sum_{\begin{subarray}{c}m_{1},\ldots,m_{s}\\\
\text{pairwise different}\end{subarray}}H_{m_{1}}\cdots H_{m_{s}}.$ (24)

###### Proof.

We take all possible products of $s$ terms from the Taylor expansion of (22).
Observe that the exponentials in (22) are organized in an array with $\kappa$
rows and $L$ columns. We use $\kappa_{1},\ldots,\kappa_{s}$ and
$l_{1},\ldots,l_{s}$ to label the row and column indices, respectively, of the
exponentials from which the terms are chosen. To avoid double counting, we
take terms with smaller row indices first (i.e.,
$\kappa_{1}\leq\cdots\leq\kappa_{s}$). Within each row, we take terms with
smaller column indices first. To get the $s$th-order nondegenerate term, we
require that $\pi_{\kappa_{1}}(l_{1}),\ldots,\pi_{\kappa_{s}}(l_{s})$ are
pairwise different. The $s$th-order nondegenerate term of (22) can then be
expressed as

$\displaystyle\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\sum_{\kappa_{1}\leq\cdots\leq\kappa_{s}}\sum_{\begin{subarray}{c}\pi_{\kappa_{1}}(l_{1}),\ldots,\pi_{\kappa_{s}}(l_{s})\\\
\text{pairwise different}\end{subarray}}(q_{\kappa_{1}}\lambda
H_{\sigma(\pi_{\kappa_{1}}(l_{1}))})\cdots(q_{\kappa_{s}}\lambda
H_{\sigma(\pi_{\kappa_{s}}(l_{s}))}).$ (25)

A direct calculation shows that

$\displaystyle\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\
\sum_{\kappa_{1}\leq\cdots\leq\kappa_{s}}\
\sum_{\begin{subarray}{c}\pi_{\kappa_{1}}(l_{1}),\ldots,\pi_{\kappa_{s}}(l_{s})\\\
\text{pairwise different}\end{subarray}}(q_{\kappa_{1}}\lambda
H_{\sigma(\pi_{\kappa_{1}}(l_{1}))})\cdots(q_{\kappa_{s}}\lambda
H_{\sigma(\pi_{\kappa_{s}}(l_{s}))})$ (26)
$\displaystyle\quad=\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\
\sum_{\kappa_{1}\leq\cdots\leq\kappa_{s}}\
\sum_{\begin{subarray}{c}\pi_{\kappa_{1}}(l_{1}),\ldots,\pi_{\kappa_{s}}(l_{s})\\\
\text{pairwise different}\end{subarray}}\
\sum_{\begin{subarray}{c}m_{1}=\sigma(\pi_{\kappa_{1}}(l_{1})),\ldots,\\\
m_{s}=\sigma(\pi_{\kappa_{s}}(l_{s}))\end{subarray}}(q_{\kappa_{1}}\lambda
H_{m_{1}})\cdots(q_{\kappa_{s}}\lambda H_{m_{s}})$
$\displaystyle\quad=\frac{1}{L!}\sum_{\begin{subarray}{c}m_{1},\ldots,m_{s}\\\
\text{pairwise different}\end{subarray}}\
\sum_{\kappa_{1}\leq\cdots\leq\kappa_{s}}\
\sum_{\begin{subarray}{c}\pi_{\kappa_{1}}(l_{1}),\ldots,\pi_{\kappa_{s}}(l_{s})\\\
\text{pairwise different}\end{subarray}}\
\sum_{\begin{subarray}{c}\sigma\in\operatorname{Sym}(L):\\\
\sigma(\pi_{\kappa_{1}}(l_{1}))=m_{1},\ldots,\\\
\sigma(\pi_{\kappa_{s}}(l_{s}))=m_{s}\end{subarray}}(q_{\kappa_{1}}\lambda
H_{m_{1}})\cdots(q_{\kappa_{s}}\lambda H_{m_{s}})$
$\displaystyle\quad=\frac{(L-s)!}{L!}\sum_{\begin{subarray}{c}m_{1},\ldots,m_{s}\\\
\text{pairwise
different}\end{subarray}}\bigg{[}\sum_{\kappa_{1}\leq\cdots\leq\kappa_{s}}\
\sum_{\begin{subarray}{c}\pi_{\kappa_{1}}(l_{1}),\ldots,\pi_{\kappa_{s}}(l_{s})\\\
\text{pairwise
different}\end{subarray}}(q_{\kappa_{1}}\lambda)\cdots(q_{\kappa_{s}}\lambda)\bigg{]}H_{m_{1}}\cdots
H_{m_{s}}.$

Now observe that the summand
$(q_{\kappa_{1}}\lambda)\cdots(q_{\kappa_{s}}\lambda)$ depends only on the row
indices. Letting $r_{1},\ldots,r_{\kappa}$ denote the number of terms picked
from row $1,\ldots,\kappa$, respectively, we can re-express this summand as
$(q_{1}\lambda)^{r_{1}}\cdots(q_{\kappa}\lambda)^{r_{\kappa}}$. We determine
the coefficient of this term as follows. The number of ways of choosing
$l_{1},\ldots,l_{s}$ pairwise different is $L(L-1)\cdots(L-s+1)$. However,
when we apply permutations $\pi_{\kappa_{1}},\ldots,\pi_{\kappa_{s}}$, we may
double count some terms. In particular, if $\kappa_{i}=\kappa_{i+1}$, we are
to pick terms from the same row $\kappa_{i}$ and we must have $l_{i}<l_{i+1}$.
This implies that the ordering of $\pi_{\kappa_{i}}(l_{i})$ and
$\pi_{\kappa_{i+1}}(l_{i+1})$ is uniquely determined. Altogether, we see that
we have overcounted by a factor of $(r_{1}!)\cdots(r_{\kappa}!)$. Therefore,
we have

$\displaystyle\sum_{\kappa_{1}\leq\cdots\leq\kappa_{s}}\sum_{\begin{subarray}{c}\pi_{\kappa_{1}}(l_{1}),\ldots,\pi_{\kappa_{s}}(l_{s})\\\
\text{pairwise
different}\end{subarray}}\\!(q_{\kappa_{1}}\lambda)\cdots(q_{\kappa_{s}}\lambda)$
$\displaystyle=\\!\sum_{\begin{subarray}{c}r_{1},\ldots,r_{\kappa}:\\\
r_{1}+\cdots+r_{\kappa}=s\end{subarray}}\\!\\!\\!\frac{L(L-1)\cdots(L-s+1)}{(r_{1}!)\cdots(r_{\kappa}!)}(q_{1}\lambda)^{r_{1}}\cdots(q_{\kappa}\lambda)^{r_{\kappa}}$
(27)
$\displaystyle=L(L-1)\cdots(L-s+1)\frac{[(q_{1}+\cdots+q_{\kappa})\lambda]^{s}}{s!},$

where the last equality follows by the multinomial theorem.

Substituting (27) into (26) completes the proof. ∎

As an immediate corollary, we compute the $s$th-order nondegenerate term of
the average evolution operator
$\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}S_{2k}^{\sigma}(\lambda)$.

######  Corollary 1.

Let $\\{H_{j}\\}_{j=1}^{L}$ be Hermitian operators; let
$\lambda\in\mathbb{C}$, $k,s\in\mathbb{N}$, and $s\leq L$. Then the $s$th-
order nondegenerate term of the average evolution
$\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}S_{2k}^{\sigma}(\lambda)$,
with $S_{2k}^{\sigma}(\lambda)$ defined in (6), is

$\frac{\lambda^{s}}{s!}\sum_{\begin{subarray}{c}m_{1},\ldots,m_{s}\\\
\text{pairwise different}\end{subarray}}H_{m_{1}}\cdots H_{m_{s}}.$ (28)

###### Proof.

The fact that $S_{2k}^{\sigma}(\lambda)$ is at least first-order accurate
implies that $q_{1}+\cdots+q_{\kappa}=1$ in (24). ∎

Observe that the $s$th-order nondegenerate term of
$V(\lambda)=\exp\bigl{(}\lambda\sum_{j=1}^{L}H_{j}\bigr{)}$ is also given by
(28). Therefore, the $s$th-order nondegenerate term completely cancels in

$V(\lambda)-\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}S_{2k}^{\sigma}(\lambda).$
(29)

##  4 Error bounds

In this section we establish our main result, an upper bound on the error of a
randomized product formula simulation. To apply the mixing lemma, we need to
bound the error of the average evolution. We now present an error bound for an
arbitrary fixed-order term in the Taylor expansion of the average evolution
operator.

######  Lemma 3.

Let $\\{H_{j}\\}_{j=1}^{L}$ be Hermitian operators; let $\lambda\in\mathbb{C}$
and $k,s\in\mathbb{N}$. Define the target evolution $V(\lambda)$ as in (2),
and define the permuted $(2k)$th-order formula $S_{2k}^{\sigma}(\lambda)$ as
in (6). Then the $s$th-order error of the approximation

$V(\lambda)-\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}S_{2k}^{\sigma}(\lambda)$
(30)

is at most

$\begin{cases}0&0\leq s\leq 2k,\\\ \frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{s}}{(s-2)!}L^{s-1}&s>2k,\end{cases}$ (31)

where $\Lambda:=\max\left\lVert H_{j}\right\rVert$.

The proof of this error bound uses the following estimate of a fixed-order
degenerate term in the average evolution operator.

######  Lemma 4.

Let $\\{H_{j}\\}_{j=1}^{L}$ be Hermitian operators with
$\Lambda:=\max_{j}\left\lVert H_{j}\right\rVert$; let
$q_{1},\ldots,q_{\kappa}\in\mathbb{R}$ with $\max_{k}|q_{k}|\leq 1$; and let
$s\leq L$ be a positive integer. Then the norm of the $s$th-order degenerate
term of the ideal evolution operator $V(\lambda)$ as in (2) is at most

$\frac{(\Lambda|\lambda|)^{s}}{s!}\big{[}L^{s}-L(L-1)\cdots(L-s+1)\big{]}$
(32)

and the norm of the $s$th-order degenerate term of the average evolution
operator as in (22) is at most

$\frac{(\kappa\Lambda|\lambda|)^{s}}{s!}\big{[}L^{s}-L(L-1)\cdots(L-s+1)\big{]}.$
(33)

###### Proof.

The $s$th-order term of $V(\lambda)$ is

$\frac{\big{(}\lambda\sum_{j=1}^{L}H_{j}\big{)}^{s}}{s!}=\frac{\lambda^{s}}{s!}\sum_{\begin{subarray}{c}m_{1},\ldots,m_{s}\end{subarray}}H_{m_{1}}\cdots
H_{m_{s}}$ (34)

and its nondegenerate term is

$\frac{\lambda^{s}}{s!}\sum_{\begin{subarray}{c}m_{1},\ldots,m_{s}\\\
\text{pairwise different}\end{subarray}}H_{m_{1}}\cdots H_{m_{s}}.$ (35)

We use the following strategy to bound the norms of these terms: (i) bound the
norm of a sum of terms by summing the norms of each term; (ii) bound the norm
of a product of terms by multiplying the norms of each term; (iii) bound the
norm of each summand by $\Lambda$; and (iv) replace $\lambda$ by $|\lambda|$.
Applying this strategy, we find that the norm of the $s$th-order term is at
most $(L\Lambda|\lambda|)^{s}/s!$, where the nondegenerate term contributes
precisely $L(L-1)\cdots(L-s+1)(\Lambda|\lambda|)^{s}/s!$. Taking the
difference gives the desired bound (32).

According to Lemma 2, the $s$th-order nondegenerate term of the average
evolution is

$\frac{[(q_{1}+\cdots+q_{\kappa})\lambda]^{s}}{s!}\sum_{\begin{subarray}{c}m_{1},\ldots,m_{s}\\\
\text{pairwise different}\end{subarray}}H_{m_{1}}\cdots H_{m_{s}}.$ (36)

Following the same strategy as for $V(\lambda)$ and also upper bounding the
norm of each $q_{k}$ by $1$ as part of step (iv), we find that the norm of
this term is at most

$\frac{(\kappa\Lambda|\lambda|)^{s}}{s!}L(L-1)\cdots(L-s+1).$ (37)

It remains to find an upper bound for the entire $s$th-order term of the
average evolution. To this end, we start with the average evolution (22) and
apply the following strategy: (i′) replace each summand of the Hamiltonian by
$\Lambda$; (ii′) replace each $q_{k}$ by $1$ and each $\lambda$ by
$|\lambda|$; and (iii′) expand all exponentials into their Taylor series and
extract the $s$th-order term. In other words, we extract the $s$th-order term
of $\sum_{\sigma\in\operatorname{Sym}(L)}\exp(\kappa L\Lambda|\lambda|)/L!$ to
get

$\frac{(\kappa L\Lambda|\lambda|)^{s}}{s!}.$ (38)

The equivalence of strategies (i)–(iv) and (i′)–(iii′) can be seen from [17,
Eq. (57)]. Finally, taking the difference between (38) and (37) gives the
desired bound (33). ∎

###### Proof of Lemma 3.

We first prove a stronger bound, namely that the $s$th-order error is at most

$\begin{cases}0&0\leq s\leq 2k,\\\ 2\frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{s}}{s!}[L^{s}-L(L-1)\cdots(L-s+1)]&2k<s\leq L,\\\
2\frac{(2\cdot 5^{k-1}\Lambda|\lambda|)^{s}}{s!}L^{s}&s>L.\end{cases}$ (39)

The first and third cases in this expression are straightforward. The formula
$S_{2k}^{\sigma}$ is exact for terms with order $0\leq s\leq 2k$ (this is what
it means for the formula to have order $2k$), so the error is zero in this
case. When $s>L$, the randomization lemma is not applicable and the error can
be bounded as in [17, Proof of Proposition F.3].

To handle the remaining case $2k<s\leq L$, we apply Lemma 4 with
$\kappa=2\cdot 5^{k-1}$. This choice of $\kappa$ follows from the definition
of the ($2k$)th-order formula (3). The norm of the $s$th-order degenerate
terms can be upper bounded by

$\frac{(\Lambda|\lambda|)^{s}}{s!}\big{[}L^{s}-L(L-1)\cdots(L-s+1)\big{]}+\frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{s}}{s!}\big{[}L^{s}-L(L-1)\cdots(L-s+1)\big{]}.$
(40)

According to Corollary 1, the $s$th-order nondegenerate term of (30) cancels,
which proves (39) for $2k<s\leq L$.

To finish the proof, we need a unified error expression for order $s>2k$. When
$2k<s\leq L$, we have

$\displaystyle L^{s}-L(L-1)\cdots(L-s+1)$ (41)
$\displaystyle\quad=\\#\big{\\{}(l_{1},\ldots,l_{s})\in[L]^{s}\big{\\}}-\\#\big{\\{}(l_{1},\ldots,l_{s})\in[L]^{s}:\forall
i,j,\ l_{i}\neq l_{j}\big{\\}}$
$\displaystyle\quad=\\#\big{\\{}(l_{1},\ldots,l_{s})\in[L]^{s}\big{\\}}-\\#\bigcap_{i<j}\big{\\{}(l_{1},\ldots,l_{s})\in[L]^{s}:\
l_{i}\neq l_{j}\big{\\}}$
$\displaystyle\quad=\\#\bigcup_{i<j}\big{\\{}(l_{1},\ldots,l_{s})\in[L]^{s}:\
l_{i}=l_{j}\big{\\}}$ $\displaystyle\quad\leq\binom{s}{2}L^{s-1},$

with $\\#\\{\cdot\\}$ denoting the size of a set and $[L]:=\\{1,\ldots,L\\}$,
where the inequality follows from the union bound. Therefore, we have

$\displaystyle 2\frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{s}}{s!}[L^{s}-L(L-1)\cdots(L-s+1)]$
$\displaystyle\leq\frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{s}}{s!}s(s-1)L^{s-1}$ (42)
$\displaystyle=\frac{(2\cdot 5^{k-1}\Lambda|\lambda|)^{s}}{(s-2)!}L^{s-1}.$

If $s>L\in\mathbb{N}$, we have $s(s-1)\geq(L+1)L\geq 2L$ and

$2\frac{(2\cdot 5^{k-1}\Lambda|\lambda|)^{s}}{s!}L^{s}\leq\frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{s}}{(s-2)!}L^{s-1}.$ (43)

This completes the proof. ∎

We also use the following standard tail bound on the exponential function [17,
Lemma F.2].

######  Lemma 5.

For any $x\in\mathbb{C}$ and $\kappa\in\mathbb{N}$, we have

$\bigg{|}\sum_{s=\kappa}^{\infty}\frac{x^{s}}{s!}\bigg{|}\leq\frac{|x|^{\kappa}}{\kappa!}\exp(|x|).$
(44)

We now establish the main theorem, which upper bounds the error of a higher-
order randomized product formula.

######  Theorem 2 (Randomized higher-order error bound).

Let $\\{H_{j}\\}_{j=1}^{L}$ be Hermitian matrices. Let

$V(-it):=\exp\biggl{(}-it\sum_{j=1}^{L}H_{j}\biggr{)}$ (45)

be the evolution induced by the Hamiltonian $H=\sum_{j=1}^{L}H_{j}$ for time
$t$. For any permutation $\sigma\in\operatorname{Sym}(L)$, define the permuted
$(2k)$th-order formula recursively by

$\displaystyle S_{2}^{\sigma}(\lambda)$
$\displaystyle:=\prod_{j=1}^{L}\exp\bigg{(}\frac{\lambda}{2}H_{\sigma(j)}\bigg{)}\prod_{j=L}^{1}\exp\bigg{(}\frac{\lambda}{2}H_{\sigma(j)}\bigg{)}$
(46) $\displaystyle S_{2k}^{\sigma}(\lambda)$
$\displaystyle:=[S_{2k-2}^{\sigma}(p_{k}\lambda)]^{2}S_{2k-2}^{\sigma}((1-4p_{k})\lambda)[S_{2k-2}^{\sigma}(p_{k}\lambda)]^{2},$

with $p_{k}:=1/(4-4^{1/(2k-1)})$ for $k>1$. Let $r\in\mathbb{N}$ and
$\Lambda:=\max\left\lVert H_{j}\right\rVert$. Then

$\displaystyle\Biggl{\|}\mathcal{V}(-it)-\biggl{(}\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\mathcal{S}_{2k}^{\sigma}(-it/r)\biggr{)}^{r}\Biggr{\|}_{\diamond}$
(47) $\displaystyle\quad\leq 4\frac{(2\cdot
5^{k-1}\Lambda|t|L)^{4k+2}}{\bigl{(}(2k+1)!\bigr{)}^{2}r^{4k+1}}\exp\biggl{(}4\cdot
5^{k-1}\frac{\Lambda|t|L}{r}\biggr{)}+2\frac{(2\cdot
5^{k-1}\Lambda|t|)^{2k+1}L^{2k}}{(2k-1)!r^{2k}}\exp\biggl{(}2\cdot
5^{k-1}\frac{\Lambda|t|L}{r}\biggr{)}$

where, for $\lambda=-it$, we associate quantum channels $\mathcal{V}(\lambda)$
and $\mathcal{S}_{2k}^{\sigma}(\lambda)$ with the unitaries $V(\lambda)$ and
$S_{2k}^{\sigma}(\lambda)$, respectively.

###### Proof.

We first prove that

$\displaystyle\left\lVert\mathcal{V}(\lambda)-\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\mathcal{S}_{2k}^{\sigma}(\lambda)\right\rVert_{\diamond}$
(48) $\displaystyle\quad\leq 4\frac{(2\cdot
5^{k-1}\Lambda|\lambda|L)^{4k+2}}{\bigl{(}(2k+1)!\bigr{)}^{2}}\exp\bigl{(}4\cdot
5^{k-1}\Lambda|\lambda|L\bigr{)}+2\frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{2k+1}L^{2k}}{(2k-1)!}\exp\bigl{(}2\cdot
5^{k-1}\Lambda|\lambda|L\bigr{)}.$

To this end, note that the $s$th-order error of
$V(\lambda)-S_{2k}^{\sigma}(\lambda)$ is at most

$\begin{cases}0&0\leq s\leq 2k,\\\ \frac{2(2\cdot
5^{k-1}\Lambda|\lambda|)^{s}}{s!}L^{s}&s>2k\end{cases}$ (49)

(as before, this follows as in [17, Proof of Proposition F.3]). Thus Lemma 5
gives

$\left\lVert V(\lambda)-S_{2k}^{\sigma}(\lambda)\right\rVert\leq
2\frac{(2\cdot 5^{k-1}\Lambda|\lambda|L)^{2k+1}}{(2k+1)!}\exp\bigl{(}2\cdot
5^{k-1}\Lambda|\lambda|L\bigr{)}.$ (50)

On the other hand, Lemma 3 implies that the $s$th-order error of
$V(\lambda)-\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}S_{2k}^{\sigma}(\lambda)$
is at most

$\begin{cases}0&0\leq s\leq 2k,\\\ \frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{s}}{(s-2)!}L^{s-1}&s>2k,\end{cases}$ (51)

so again Lemma 5 gives

$\left\lVert
V(\lambda)-\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}S_{2k}^{\sigma}(\lambda)\right\rVert\leq\frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{2k+1}L^{2k}}{(2k-1)!}\exp\bigl{(}2\cdot
5^{k-1}\Lambda|\lambda|L\bigr{)}.$ (52)

Equation (48) now follows from Lemma 1 by setting

$\displaystyle a$ $\displaystyle=2\frac{(2\cdot
5^{k-1}\Lambda|\lambda|L)^{2k+1}}{(2k+1)!}\exp\bigl{(}2\cdot
5^{k-1}\Lambda|\lambda|L\bigr{)},$ (53) $\displaystyle b$
$\displaystyle=\frac{(2\cdot
5^{k-1}\Lambda|\lambda|)^{2k+1}L^{2k}}{(2k-1)!}\exp\bigl{(}2\cdot
5^{k-1}\Lambda|\lambda|L\bigr{)}.$

To simulate the evolution for time $t$, we divide it into $r$ segments. The
error within each segment is obtained from (48) by setting $\lambda=-it/r$.
Then subadditivity of the diamond norm distance gives

$\left\lVert\mathcal{V}(-it)-\bigg{(}\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\mathcal{S}_{2k}^{\sigma}\big{(}-it/r\big{)}\bigg{)}^{r}\right\rVert_{\diamond}\leq
r\left\lVert\mathcal{V}(-it/r)-\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\mathcal{S}_{2k}^{\sigma}\big{(}-it/r\big{)}\right\rVert_{\diamond},$
(54)

which completes the proof. ∎

##  5 Algorithm performance and comparisons

We now analyze the complexity of our randomized product formula algorithm.
Assume that $k\in\mathbb{N}$ is fixed, $\Lambda=O(1)$ is constant, and $r>tL$.
By Theorem 2, the asymptotic error of the $(2k)$th-order randomized product
formula is

$\left\lVert\mathcal{V}(-it)-\bigg{(}\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\mathcal{S}_{2k}^{\sigma}\big{(}-it/r\big{)}\bigg{)}^{r}\right\rVert_{\diamond}\leq
O\bigg{(}\frac{(tL)^{4k+2}}{r^{4k+1}}+\frac{t^{2k+1}L^{2k}}{r^{2k}}\bigg{)}.$
(55)

To guarantee that the simulation error is at most $\epsilon$, we upper bound
the right-hand side of (55) by $\epsilon$ and solve for $r$. We find that it
suffices to use

$\displaystyle r_{2k}^{\mathrm{rand}}$
$\displaystyle=\max\bigg{\\{}O\bigg{(}\frac{(tL)^{\frac{4k+2}{4k+1}}}{\epsilon^{\frac{1}{4k+1}}}\bigg{)},O\bigg{(}\frac{t^{\frac{2k+1}{2k}}L}{\epsilon^{\frac{1}{2k}}}\bigg{)}\bigg{\\}}$
(56)
$\displaystyle=\max\bigg{\\{}O\bigg{(}tL\bigg{(}\frac{tL}{\epsilon}\bigg{)}^{\frac{1}{4k+1}}\bigg{)},O\bigg{(}tL\bigg{(}\frac{t}{\epsilon}\bigg{)}^{\frac{1}{2k}}\bigg{)}\bigg{\\}}$

segments, giving a simulation algorithm with

$g_{2k}^{\mathrm{rand}}=O(Lr_{2k}^{\mathrm{rand}})=\max\bigg{\\{}O\bigg{(}tL^{2}\bigg{(}\frac{tL}{\epsilon}\bigg{)}^{\frac{1}{4k+1}}\bigg{)},O\bigg{(}tL^{2}\bigg{(}\frac{t}{\epsilon}\bigg{)}^{\frac{1}{2k}}\bigg{)}\bigg{\\}}$
(57)

elementary gates.

For comparison, the error in the ($2k$)th-order deterministic formula
algorithm is at most [17, Proposition F.4]

$\left\lVert V(-it)-\bigl{[}S_{2k}(-it/r)\bigr{]}^{r}\right\rVert\leq
O\biggl{(}\frac{(tL)^{2k+1}}{r^{2k}}\biggr{)}.$ (58)

While this bound quantifies the simulation error in terms of the spectral-norm
distance, it can easily be adapted to the diamond-norm distance using either
Lemma 1 or [8, Lemma 7]. This translation introduces only constant-factor
overhead, so we have

$\left\lVert\mathcal{V}(-it)-\bigl{[}\mathcal{S}_{2k}(-it/r)\bigr{]}^{r}\right\rVert_{\diamond}\leq
O\biggl{(}\frac{(tL)^{2k+1}}{r^{2k}}\biggr{)}.$ (59)

Therefore, the number of segments that suffice to ensure error at most
$\epsilon$ satisfies

$r_{2k}^{\det}=O\biggl{(}tL\biggl{(}\frac{tL}{\epsilon}\biggr{)}^{\frac{1}{2k}}\biggr{)},$
(60)

giving an algorithm with

$g_{2k}^{\det}=O(Lr_{2k}^{\det})=O\biggl{(}tL^{2}\biggl{(}\frac{tL}{\epsilon}\biggr{)}^{\frac{1}{2k}}\biggr{)}$
(61)

elementary gates. Comparing to (57), we see that the randomized product
formula strictly improves the complexity as a function of $L$. Indeed, the
$(2k)$th-order randomized approach either provides an improvement with respect
to all parameters of interest over the $(2k)$th order deterministic approach
(if the first term of (57) obtains the maximum), or has better dependence on
the number of terms in the Hamiltonian than any deterministic formula (if the
second term dominates).

We can also compare our result to the commutator bound of [17], which depends
on the specific structure of the Hamiltonian. For concreteness, we consider a
one-dimensional nearest-neighbor Heisenberg model with a random magnetic
field, as studied in [17]. Specifically, let

$\displaystyle
H=\sum_{j=1}^{n}(\vec{\sigma}_{j}\cdot\vec{\sigma}_{j+1}+h_{j}\sigma_{j}^{z})$
(62)

with periodic boundary conditions (i.e.,
$\vec{\sigma}_{n+1}=\vec{\sigma}_{1}$), and $h_{j}\in[-h,h]$ chosen uniformly
at random, where
$\vec{\sigma}_{j}=(\sigma^{x}_{j},\sigma^{y}_{j},\sigma^{z}_{j})$ denotes a
vector of Pauli $x$, $y$, and $z$ matrices on qubit $j$. The $(2k)$th-order
deterministic formula with the commutator bound has error at most [17, Eq.
(146)]

$\left\lVert\mathcal{V}(-it)-\big{[}\mathcal{S}_{2k}\big{(}-it/r\big{)}\big{]}^{r}\right\rVert_{\diamond}\leq
O\biggl{(}\frac{(tL)^{2k+2}}{r^{2k+1}}+\frac{t^{2k+1}L^{2k}}{r^{2k}}\biggr{)},$
(63)

where we have again used Lemma 1 (or [8, Lemma 7]) to relate the spectral-norm
distance to the diamond-norm distance. To guarantee that the simulation error
is at most $\epsilon$, it suffices to choose

$\displaystyle r_{2k}^{\mathrm{comm}}$
$\displaystyle=\max\biggl{\\{}O\biggl{(}\frac{(tL)^{\frac{2k+2}{2k+1}}}{\epsilon^{\frac{1}{2k+1}}}\biggr{)},O\biggl{(}\frac{t^{\frac{2k+1}{2k}}L}{\epsilon^{\frac{1}{2k}}}\biggr{)}\biggr{\\}}$
(64)
$\displaystyle=\max\biggl{\\{}O\biggl{(}tL\biggl{(}\frac{tL}{\epsilon}\biggr{)}^{\frac{1}{2k+1}}\biggr{)},O\biggl{(}tL\biggl{(}\frac{t}{\epsilon}\biggr{)}^{\frac{1}{2k}}\biggr{)}\biggr{\\}}$

segments, giving an algorithm with

$g_{2k}^{\mathrm{comm}}=O(Lr_{2k}^{\mathrm{comm}})=\max\biggl{\\{}O\biggl{(}tL^{2}\biggl{(}\frac{tL}{\epsilon}\biggr{)}^{\frac{1}{2k+1}}\biggr{)},O\biggl{(}tL^{2}\biggl{(}\frac{t}{\epsilon}\biggr{)}^{\frac{1}{2k}}\biggr{)}\biggr{\\}}$
(65)

elementary gates. Comparing to the corresponding bound (57) for randomized
product formulas, we see that the only difference is that the exponent
$1/(2k+1)$ for the commutator bound becomes $1/(4k+1)$ in the randomized case.
Thus the randomized approach can provide a slightly faster algorithm despite
using less information about the structure of the Hamiltonian. More
specifically, the relationship between $t$ and $L$ determines whether the
randomized approach offers an improvement. If $t=\Omega(L^{2k})$, then the
second term of (65) achieves the maximum, and both approaches have asymptotic
complexity
$O\big{(}tL^{2}\bigl{(}\frac{t}{\epsilon}\big{)}^{\frac{1}{2k}}\bigr{)}$.
However, if $t=o(L^{2k})$, then the randomized formula is advantageous.

##  6 Empirical performance

While randomization provides a useful theoretical handle for establishing
better provable bounds, those bounds may still be far from tight. As described
in Section 1, our original motivation for considering randomization was the
observation that product formulas appear to perform dramatically better in
practice than the best available proven bounds would suggest. To investigate
the empirical behavior of product formulas, we numerically evaluate their
performance for simulations of the Heisenberg model (62) with $t=n$ and $h=1$,
targeting error $\epsilon=10^{-3}$, as previously considered in [17]. We
collect data for the first-, fourth-, and sixth-order formulas as the latter
two orders have the best performance in practice for small $n$ and the first-
order formula offers a qualitatively better theoretical improvement.

For the deterministic formula, we order the operators of the Hamiltonian in
the same way as [17], namely

$\displaystyle\sigma_{1}^{x}\sigma_{2}^{x},\ldots,\sigma_{n-1}^{x}\sigma_{n}^{x},\sigma_{n}^{x}\sigma_{1}^{x},\;\sigma_{1}^{y}\sigma_{2}^{y},\ldots,\sigma_{n-1}^{y}\sigma_{n}^{y},\sigma_{n}^{y}\sigma_{1}^{y},\;\sigma_{1}^{z}\sigma_{2}^{z},\ldots,\sigma_{n-1}^{z}\sigma_{n}^{z},\sigma_{n}^{z}\sigma_{1}^{z},\;\sigma_{1}^{z},\ldots,\sigma_{n}^{z}.$
(66)

We compute the error in terms of the spectral-norm distance and convert it to
the diamond-norm distance using Lemma 7 of [8] (i.e., we multiply by $2$). To
analyze the randomized formula, we would like to numerically evaluate the
diamond-norm distances

$\left\lVert\mathcal{V}(-it)-\frac{1}{2^{r}}\bigl{(}\mathcal{S}_{1}(-it/r)+\mathcal{S}_{1}^{\mathrm{rev}}(-it/r)\bigr{)}^{r}\right\rVert_{\diamond}$
(67)

and

$\left\lVert\mathcal{V}(-it)-\biggl{(}\frac{1}{L!}\sum_{\sigma\in\operatorname{Sym}(L)}\mathcal{S}_{2k}^{\sigma}\bigl{(}-it/r\bigr{)}\biggr{)}^{r}\right\rVert_{\diamond}.$
(68)

While the diamond norm can be computed using a semidefinite program [34],
direct computation is prohibitive as the channel contains $(L!)^{r}$ Kraus
operators. Instead, we use Lemma 1 to estimate the error. We randomly choose
the ordering of the summands in each of the $r$ segments, exponentiate each
individual operator, and construct a unitary operator by concatenating the
exponentials according to the given product formula. We follow this procedure
to obtain a Monte Carlo estimate of the average error

$\left\lVert
V(-it)-\frac{1}{M}\sum_{m=1}^{M}S_{2k}^{\sigma_{m,r}}\bigl{(}-it/r\bigr{)}\cdots
S_{2k}^{\sigma_{m,1}}\bigl{(}-it/r\bigr{)}\right\rVert$ (69)

for the ($2k$)th-order formula and similarly for the first-order case. Here,
$M$ is the number of samples in the Monte Carlo estimation, which can be
increased to get more accurate estimate. In practice, we find that it suffices
to take only three samples, as the standard deviations are already negligibly
small (about $10^{-5}$). We then invoke Lemma 1 to bound the diamond-norm
error in (68). To the extent that the bound of Lemma 1 is loose, we expect the
empirical performance to be better in practice.

$6$ $7$ $8$ $9$ $10$ $10^{3}$$10^{4}$$10^{5}$$10^{6}$$n$$r$First
orderDeterministicRandomized

$6$ $7$ $8$ $9$ $10$ $10^{1}$$10^{2}$$10^{3}$$n$$r$Fourth
orderDeterministicRandomized

$6$ $7$ $8$ $9$ $10$ $10^{1}$$10^{2}$$n$$r$Sixth orderDeterministicRandomized

Figure 1:  Comparison of the values of $r$ between deterministic and
randomized product formulas. Error bars are omitted when they are negligibly
small on the plot. Straight lines show power-law fits to the data.

Using five randomly generated instances for each value of $n$, we apply binary
search to determine the smallest number of segments $r$ that suffices to give
error at most $10^{-3}$. Figure 1 shows the resulting data for the first-,
fourth-, and sixth-order formulas, which are well-approximated by power laws.
Fitting the data, we estimate that

$\displaystyle r_{1}^{\mathrm{remp}}$ $\displaystyle=300.0n^{1.806}$
$\displaystyle r_{4}^{\mathrm{remp}}$ $\displaystyle=5.458n^{1.439}$
$\displaystyle r_{6}^{\mathrm{remp}}$ $\displaystyle=2.804n^{1.152}$ (70)

segments should suffice to give error at most $10^{-3}$. We thus observe that
the empirical complexity of the randomized algorithm is still significantly
better than the provable performance

$\displaystyle r_{1}^{\mathrm{rand}}$ $\displaystyle=O(n^{3})$ $\displaystyle
r_{4}^{\mathrm{rand}}$ $\displaystyle=O(n^{2.25})$ $\displaystyle
r_{6}^{\mathrm{rand}}$ $\displaystyle=O(n^{2.17}).$ (71)

For comparison, analogous empirical fits for deterministic formulas give the
comparable values

$\displaystyle r_{1}^{\mathrm{demp}}$ $\displaystyle=4143n^{2.066}$
$\displaystyle r_{4}^{\mathrm{demp}}$ $\displaystyle=5.821n^{1.471}$
$\displaystyle r_{6}^{\mathrm{demp}}$ $\displaystyle=2.719n^{1.160},$ (72)

(cf. [17, Eq. (147)], but note that we have generated new data using [8, Lemma
7] to bound the diamond-norm distance in terms of the spectral-norm distance),
whereas the rigorous commutator bound gives the larger exponents [17]

$\displaystyle r_{1}^{\mathrm{comm}}$ $\displaystyle=O(n^{3})$ $\displaystyle
r_{4}^{\mathrm{comm}}$ $\displaystyle=O(n^{2.4})$ $\displaystyle
r_{6}^{\mathrm{comm}}$ $\displaystyle=O(n^{2.28}).$ (73)

We see that the randomized bound offers significantly better empirical
performance at first order, consistent with the observation that randomization
improves the order of approximation in this case. The fourth-order formula
slightly improves both the exponent and the constant factor. While this
improvement is small, it is nevertheless notable since it involves only a
minor change to the algorithm. At sixth order we see negligible improvement.
Since the proven bounds give less improvement with each successive order, it
is perhaps not surprising to see that the empirical performance shows similar
behavior.

$10$ $20$ $30$ $40$ $50$ $60$ $70$ $80$ $90$ $100$
$10^{5}$$10^{6}$$10^{7}$$10^{8}$$10^{9}$$10^{10}$$10^{11}$$10^{12}$$10^{13}$$n$number
of exponentialsDeterministic PF4 (min)Deterministic PF4 (com)Randomized PF4
(min)Deterministic PF4 (emp)Randomized PF4 (emp)Deterministic PF6
(min)Randomized PF6 (min)Deterministic PF6 (emp)Randomized PF6 (emp)

Figure 2: Comparison of the total number of elementary exponentials for
product formula simulations of the Heisenberg model using deterministic and
randomized product formulas of fourth and sixth order with both rigorous and
empirical error bounds. Note that since the empirical performance of
deterministic and randomized sixth-order product formulas is almost the same,
the latter data points are obscured by the former.

To illustrate the effect of using different formulas and different error
bounds to simulate larger systems, Figure 2 compares the cost of simulating
our model system for sizes up to $n=100$ with deterministic and randomized
formulas of orders $4$ and $6$, using both proven error bounds and the above
empirical estimates. (We omit the first-order formula since it is not
competitive even at such small sizes.) We give rigorous bounds for
deterministic formulas using the minimized bound of [17], and for fourth order
we also show the result of using the commutator bound. We see that
randomization gives a significant improvement over the deterministic formula
using the minimized bound, although the commutator bound outperforms the
randomized bound at the system sizes shown here. For sufficiently large $n$,
the randomized bound gives lower complexity, but this requires a fairly large
$n$ since the difference in exponents is small and the commutator bound
achieves a favorable constant prefactor. Empirical estimates of the error
improve the performance by several orders of magnitude, with randomization
giving a small advantage for the fourth-order formula as indicated above.
However, for systems of size larger than about $n=25$, the sixth-order bound
prevails, and in this case randomization no longer offers a significant
advantage.

##  7 Discussion

We have shown that randomization can be used to establish better performance
for quantum simulation algorithms based on product formulas. By simply
randomizing how the summands in the Hamiltonian are ordered, we introduce
terms in the average evolution that could not appear in any deterministic
product formula approximation of the same order, and thereby give a more
efficient algorithm. Indeed, this approach can outperform the commutator bound
even though that method uses more information about the structure of the
Hamiltonian. A randomized product formula simulation algorithm is not much
more complicated than the corresponding deterministic formula, using only
$O(L\log L)$ bits of randomness per segment and no ancilla qubits.
Furthermore, we showed that randomization can even offer improved empirical
performance in some cases.

While randomization has allowed us to make some progress on the challenge of
proving better bounds on the performance of product formulas, our strengthened
bounds remain far from the apparent empirical performance. We expect that
other ideas will be required to improve the product-formula approach [15, 28].
Although our bounds have better asymptotic $n$-dependence than the previous
commutator bound, they only offer an improvement if the system is sufficiently
large. It could be fruitful to establish bounds for randomized product
formulas that take advantage of the structure of the Hamiltonian, perhaps
offering better performance both asymptotically and for small system sizes.
More generally, it may be of interest to investigate other scenarios in which
random choices can be used to improve the analysis of quantum simulation [14,
10] and other quantum algorithms.

## Acknowledgments

We thank Guoming Wang for helpful discussions during the initial stages of
this work and anonymous referees for their helpful comments on our manuscript.

This work was supported in part by the Army Research Office (MURI award
W911NF-16-1-0349), the Canadian Institute for Advanced Research, the
Department of Energy (grant 17-020469), and the National Science Foundation
(grant 1526380).

## References

  * Aharonov and Ta-Shma [2003] Dorit Aharonov and Amnon Ta-Shma.  Adiabatic quantum state generation and statistical zero knowledge.  In _Proceedings of the 35th ACM Symposium on Theory of Computing_ , pages 20–29, 2003.  doi: 10.1145/780542.780546.  arXiv:quant-ph/0301023.
  * Babbush et al. [2015] Ryan Babbush, Jarrod McClean, Dave Wecker, Alán Aspuru-Guzik, and Nathan Wiebe.  Chemical basis of Trotter-Suzuki errors in quantum chemistry simulation.  _Physical Review A_ , 91:022311, 2015.  doi: 10.1103/PhysRevA.91.022311.  arXiv:1410.8159.
  * Barends et al. [2015] R. Barends, L. Lamata, J. Kelly, L. García-Álvarez, A. G. Fowler, A Megrant, E Jeffrey, T. C. White, D. Sank, J. Y. Mutus, B. Campbell, Yu Chen, Z. Chen, B. Chiaro, A. Dunsworth, I.-C. Hoi, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, A. Vainsencher, J. Wenner, E. Solano, and John M. Martinis.  Digital quantum simulation of fermionic models with a superconducting circuit.  _Nature Communications_ , 6:7654, 2015.  doi: 10.1038/ncomms8654.  arXiv:1501.07703.
  * Berry and Childs [2012] Dominic W. Berry and Andrew M. Childs.  Black-box Hamiltonian simulation and unitary implementation.  _Quantum Information and Computation_ , 12(1-2):29–62, 2012.  arXiv:0910.4157.
  * Berry et al. [2007] Dominic W. Berry, Graeme Ahokas, Richard Cleve, and Barry C. Sanders.  Efficient quantum algorithms for simulating sparse Hamiltonians.  _Communications in Mathematical Physics_ , 270(2):359–371, 2007.  doi: 10.1007/s00220-006-0150-x.  arXiv:quant-ph/0508139.
  * Berry et al. [2014] Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma.  Exponential improvement in precision for simulating sparse Hamiltonians.  In _Proceedings of the 46th ACM Symposium on Theory of Computing_ , pages 283–292, 2014.  doi: 10.1145/2591796.2591854.  arXiv:1312.1414.
  * Berry et al. [2015a] 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(9):090502, 2015a.  doi: 10.1103/PhysRevLett.114.090502.  arXiv:1412.4687.
  * Berry et al. [2015b] Dominic W. Berry, Andrew M. Childs, and Robin Kothari.  Hamiltonian simulation with nearly optimal dependence on all parameters.  In _Proceedings of the 56th IEEE Symposium on Foundations of Computer Science_ , pages 792–809, 2015b.  doi: 10.1109/FOCS.2015.54.  arXiv:1501.01715.
  * Berry et al. [2017] Dominic W. Berry, Andrew M. Childs, Aaron Ostrander, and Guoming Wang.  Quantum algorithm for linear differential equations with exponentially improved dependence on precision.  _Communications in Mathematical Physics_ , 356:1057–1081, 2017.  doi: 10.1007/s00220-017-3002-y.  arXiv:1701.03684.
  * Berry et al. [2019] Dominic W. Berry, Andrew M. Childs, Yuan Su, Xin Wang, and Nathan Wiebe.  Time-dependent Hamiltonian simulation with $L^{1}$-norm scaling, 2019\.  arXiv:1906.07115.
  * Brandao and Svore [2017] Fernando G. S. L. Brandao and Krysta M. Svore.  Quantum speed-ups for solving semidefinite programs.  In _Proceedings of the 58th IEEE Symposium on Foundations of Computer Science_ , pages 415–426, 2017.  doi: 10.1109/FOCS.2017.45.  arXiv:1609.05537.
  * Brown et al. [2006] Kenneth R. Brown, Robert J. Clark, and Isaac L. Chuang.  Limitations of quantum simulation examined by simulating a pairing Hamiltonian using nuclear magnetic resonance.  _Physical Review Letters_ , 97:050504, 2006.  doi: 10.1103/PhysRevLett.97.050504.  arXiv:quant-ph/0601021.
  * Campbell [2017] Earl Campbell.  Shorter gate sequences for quantum computing by mixing unitaries.  _Physical Review A_ , 95:042306, Apr 2017.  doi: 10.1103/PhysRevA.95.042306.  arXiv:1612.02689.
  * Campbell [2019] Earl Campbell.  Random compiler for fast Hamiltonian simulation.  _Physical Review Letters_ , 123:070503, Aug 2019.  doi: 10.1103/PhysRevLett.123.070503.  arXiv:1811.08017.
  * Childs and Su [2019] Andrew M. Childs and Yuan Su.  Nearly optimal lattice simulation by product formulas.  _Physical Review Letters_ , 123:050503, Aug 2019.  doi: 10.1103/PhysRevLett.123.050503.  arXiv:1901.00564.
  * Childs et al. [2003] Andrew M. Childs, Richard Cleve, Enrico Deotto, Edward Farhi, Sam Gutmann, and Daniel A. Spielman.  Exponential algorithmic speedup by quantum walk.  In _Proceedings of the 35th ACM Symposium on Theory of Computing_ , pages 59–68, 2003.  doi: 10.1145/780542.780552.  arXiv:quant-ph/0209131.
  * Childs et al. [2018] 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(38):9456–9461, 2018.  doi: 10.1073/pnas.1801723115.  arXiv:1711.10980.
  * Farhi et al. [2008] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum algorithm for the Hamiltonian NAND tree.  _Theory of Computing_ , 4(1):169–190, 2008.  doi: 10.4086/toc.2008.v004a008.
  * Feynman [1982] Richard P. Feynman.  Simulating physics with computers.  _International Journal of Theoretical Physics_ , 21(6-7):467–488, 1982.  doi: 10.1007/BF02650179.
  * Haah et al. [2018] Jeongwan Haah, Matthew B. Hastings, Robin Kothari, and Guang Hao Low.  Quantum algorithm for simulating real time evolution of lattice Hamiltonians.  In _2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 350–360, Oct 2018.  doi: 10.1109/FOCS.2018.00041.  arXiv:1801.03922.
  * Harrow et al. [2009] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  _Physical Review Letters_ , 103(15):150502, 2009\.  doi: 10.1103/PhysRevLett.103.150502.  arXiv:0811.3171.
  * Hastings [2017] Matthew B. Hastings.  Turning gate synthesis errors into incoherent errors.  _Quantum Information and Computation_ , 17(5-6):488–494, 2017.  arXiv:1612.01011.
  * Jordan et al. [2012] Stephen P. Jordan, Keith S. M. Lee, and John Preskill.  Quantum algorithms for quantum field theories.  _Science_ , 336(6085):1130–1133, 2012.  doi: 10.1126/science.1217069.  arXiv:1111.3633.
  * Lanyon et al. [2011] B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Gerritsma, F. Zähringer, P. Schindler, J. T. Barreiro, M. Rambach, G. Kirchmair, M. Hennrich, P. Zoller, R. Blatt, and C. F. Roos.  Universal digital quantum simulation with trapped ions.  _Science_ , 334(6052):57–61, 2011.  doi: 10.1126/science.1208001.  arXiv:1109.1512.
  * Lloyd [1996] Seth Lloyd.  Universal quantum simulators.  _Science_ , 273(5278):1073–1078, 1996.  doi: 10.1126/science.273.5278.1073.
  * Low and Chuang [2017] Guang Hao Low and Isaac L. Chuang.  Optimal Hamiltonian simulation by quantum signal processing.  _Physical Review Letters_ , 118:010501, 2017.  doi: 10.1103/PhysRevLett.118.010501.  arXiv:1606.02685.
  * Low and Chuang [2019] Guang Hao Low and Isaac L. Chuang.  Hamiltonian Simulation by Qubitization.  _Quantum_ , 3:163, July 2019.  doi: 10.22331/q-2019-07-12-163.  arXiv:1610.06546.
  * Low et al. [2019] Guang Hao Low, Vadym Kliuchnikov, and Nathan Wiebe.  Well-conditioned multiproduct Hamiltonian simulation, 2019.  arXiv:1907.11679.
  * Poulin et al. [2011] David Poulin, Angie Qarry, Rolando D. Somma, and Frank Verstraete.  Quantum simulation of time-dependent Hamiltonians and the convenient illusion of Hilbert space.  _Physical Review Letters_ , 106(17):170501, 2011\.  doi: 10.1103/PhysRevLett.106.170501.  arXiv:1102.1360.
  * Poulin et al. [2015] David Poulin, Matthew B. Hastings, Dave Wecker, Nathan Wiebe, Andrew C. Doherty, and Matthias Troyer.  The Trotter step size required for accurate quantum simulation of quantum chemistry.  _Quantum Information and Computation_ , 15(5-6):361–384, 2015.  arXiv:1406.4920.
  * Raeisi et al. [2012] Sadegh Raeisi, Nathan Wiebe, and Barry C. Sanders.  Quantum-circuit design for efficient simulations of many-body quantum dynamics.  _New Journal of Physics_ , 14:103017, 2012.  doi: 10.1088/1367-2630/14/10/103017.  arXiv:1108.4318.
  * Reiher et al. [2017] Markus Reiher, Nathan Wiebe, Krysta M. Svore, Dave Wecker, and Matthias Troyer.  Elucidating reaction mechanisms on quantum computers.  _Proceedings of the National Academy of Sciences_ , 114(29):7555–7560, 2017.  doi: 10.1073/pnas.1619152114.  arXiv:1605.03590.
  * Suzuki [1991] Masuo Suzuki.  General theory of fractal path integrals with applications to many-body theories and statistical physics.  _Journal of Mathematical Physics_ , 32(2):400–407, 1991.  doi: 10.1063/1.529425.
  * Watrous [2013] John Watrous.  Simpler semidefinite programs for completely bounded norms.  _Chicago Journal of Theoretical Computer Science_ , 2013(8), 2013.  doi: 10.4086/cjtcs.2013.008.
  * Watrous [2018] John Watrous.  _The Theory of Quantum Information_.  Cambridge University Press, 2018.  doi: 10.1017/9781316848142.
  * Wecker et al. [2014] Dave Wecker, Bela Bauer, Bryan K. Clark, Matthew B. Hastings, and Matthias Troyer.  Gate count estimates for performing quantum chemistry on small quantum computers.  _Physical Review A_ , 90:022305, 2014.  doi: 10.1103/PhysRevA.90.022305.  arXiv:1312.1695.
  * Zhang [2012] Chi Zhang.  Randomized algorithms for Hamiltonian simulation.  In Leszek Plaskota and Henryk Woźniakowski, editors, _Monte Carlo and Quasi-Monte Carlo Methods 2010_ , pages 709–719, Berlin, Heidelberg, 2012. Springer Berlin Heidelberg.  ISBN 978-3-642-27440-4.  doi: 10.1007/978-3-642-27440-4_42.