††thanks<EMAIL_ADDRESS>P. K. F. and M. S. have
contributed equally.

# Randomizing multi-product formulas for improved Hamiltonian simulation

P. K. Faehrmann Dahlem Center for Complex Quantum Systems, Freie Universität
Berlin, 14195 Berlin, Germany  M. Steudtner Dahlem Center for Complex Quantum
Systems, Freie Universität Berlin, 14195 Berlin, Germany  R. Kueng  Institute
for Integrated Circuits, Johannes Kepler University Linz, Austria  M.
Kieferová Centre for Quantum Computation and Communication Technology, Centre
for Quantum Software and Information, University of Technology Sydney, NSW
2007, Australia  J. Eisert Dahlem Center for Complex Quantum Systems, Freie
Universität Berlin, 14195 Berlin, Germany  Helmholtz-Zentrum Berlin für
Materialien und Energie, Hahn-Meitner-Platz 1, 14109 Berlin, Germany

(August 27, 2024)

###### Abstract

Quantum simulation, the simulation of quantum processes on quantum computers,
suggests a path forward for the efficient simulation of problems in condensed-
matter physics, quantum chemistry and materials science. While the majority of
quantum simulation algorithms is deterministic, a recent surge of ideas has
shown that randomization can greatly benefit algorithmic performance. In this
work, we introduce a scheme for quantum simulation that unites the advantages
of randomized compiling on the one hand and higher-order multi-product
formulas as they are used for example in linear-combination-of-unitaries (LCU)
algorithms or quantum error mitigation on the other hand. In doing so, we
propose a framework of randomized sampling that is expected to be useful for
programmable quantum simulators and present two new multi-product formula
algorithms tailored to it. Our framework greatly reduces the circuit depth by
circumventing the need for oblivious amplitude amplification required by the
implementation of multi-product formulas using standard LCU methods, rendering
it especially useful for near-term quantum computing. Our algorithms achieve a
simulation error that shrinks exponentially with the circuit depth. To
corroborate their functioning, we prove rigorous performance bounds as well as
the concentration of the randomized sampling procedure. We demonstrate the
functioning of the approach for several physically meaningful examples of
Hamiltonians, including fermionic systems and the Sachdev–Ye–Kitaev model, for
which the method provides a favorable scaling in effort.

##  I Introduction

The simulation of quantum processes on quantum computers is one of the most
eagerly anticipated use-cases for quantum computing. The ability to simulate a
system’s time evolution promises to provide insights into the dynamics of
interacting quantum systems in situations where approximate classical
simulation methods fail and constitutes one of the cornerstones of quantum
technologies [1].

This work aims at substantially improving algorithms that simulate the
dynamics of expectation values of observables. The need for developing such a
machinery stems from the observation that state-of-the-art quantum devices are
still rather limited in their realizable circuit depths and control. We
therefore assume only access to a quantum-oracle machine which implements
single-qubit state preparation, controlled time evolution and quantum
measurements and strive for minimizing the required depth of suitable quantum
algorithms. Such a setting explicitly allows for the use of product formulas,
which are the earliest algorithms proposed for the simulation of time-
independent local Hamiltonians [2].

Such Trotter-Suzuki methods, as they are called, have evolved from comparably
simple prescriptions for local Hamiltonians to sophisticated schemes able to
capture more general sparse time-independent Hamiltonians [3, 4, 5] as well as
time-dependent Hamiltonians [6, 7] and open quantum systems [8, 9]. In spite
of their relative simplicity, product formulas are still at the forefront of
Hamiltonian simulation techniques. A numerical study has shown that product
formulas can in practice outperform more complex techniques [10] and their
complexity is better than initial estimates [11] suggest. In fact, product
formulas are nearly optimal for lattice model simulations [12].

Multi-product formulas introduced in the _linear-combinations-of-unitaries
(LCU)_ [13] approach have been built upon previous results of Trotter-Suzuki
methods and have recently been improved [14]. The idea of linearly combining
unitaries has led to quantum algorithms with an exponential speedup in
precision [15, 16, 17, 18] which have been recently gaining popularity.
Besides quantum simulation in the LCU framework, the mathematical construct of
multi-product formulas (see Fig. 1) can also be used in quantum error
mitigation [19]. Although these methods have optimal asymptotic error scaling,
they inherently require deep quantum circuits for their implementation and are
thus not suitable for medium-term applications. It is therefore crucial to
find algorithms that require shorter circuits and fewer digital gates.

Multi-product formulasQuantum error mitigation [19]Linear-combination-of-
unitaries [13, 14]Randomized sampling [here] Figure 1: Multi-product formulas,
although introduced with use in the LCU framework in mind, can also be used in
different frameworks, such as quantum error mitigation or randomized sampling,
which is the main contribution of this work.

Recently, yet a new element has been introduced to aid this search: the
element of _classical randomness_. The idea of randomization in Hamiltonian
simulation has heralded a renaissance of product formula methods [20, 21, 22,
23, 24]. For single steps of such techniques, a rigorous understanding has
recently been reached [24]. Randomized algorithms can also be considered when
one is only interested in estimating expectation values. For such a task we do
not need to prepare the time-evolved state (from which the observable will be
measured) perfectly. Instead, a lower-effort randomized algorithm can be used
such that the correct expectation value is obtained only after averaging the
measurement outcomes.

The aims of this work are twofold: First, we strive for combining the
advantages of higher-order multi-product formulas with those of schemes of
randomized compiling, to create a novel framework in which multi-product
formulas can be put to good use, which we dub “randomized sampling”.
Specifically, we adopt the above scenario of computing expectation values of
observables and propose to sample product formulas from multi-product
formulas, by which we circumvent the need to use the LCU framework and its
corresponding methods such as block encodings. Instead, we implement multi-
product formulas on average through random sampling. This results in a notable
reduction in algorithmic depth at the expense of additional circuit
evaluations. In this way, we see how notions of randomized compiling and
higher-order multi-product formulas – when suitably brought together – allow
for more resource-efficient notions of quantum simulation amenable to
programmable devices. Second, we develop alternative multi-product formulas
tailored to this new framework, which promise to outperform the accuracy of
the multi-product formulas introduced by Childs and Wiebe in Ref. [13] in the
regime of short simulation times. We show that using a quantum device limited
to the aforementioned operations can already yield substantial improvements to
fully analog approaches. It is therefore especially useful when algorithms
with a pure focus on noisy-intermediate-scale-quantum devices [25], such as
variational quantum algorithms, reach their limits, but fully digital
algorithms are not yet feasible.

The remainder of this work is organized as follows: After a short review of
multi-product formulas in Section II, we present the framework of randomized
sampling and a short summary of our main results in Section III. Section IV
then contains a detailed analysis of the proposed randomized sampling
framework and multi-product formulas and gives the proofs to the results
discussed in the previous section. To conclude, we compare the performance of
our formulas to that of Trotter-Suzuki and Childs and Wiebe type multi-product
formulas [13] in Section V, discussing both their error bounds and actual
performance on a number of physically plausible and interesting Hamiltonian
models of strongly correlated quantum systems, for which we find a favorable
performance over known schemes of quantum simulation.

##  II Multi-product formulas

Multi-product formulas constitute the main ingredient of our work, and hence
we will briefly review the underlying ideas, starting with product formulas.
The goal of quantum simulations is to approximate the quantum dynamics of a
complex quantum system described by a many-body Hamiltonian decomposed as

$H=\sum_{k=1}^{L}h_{k}$ (1)

composed of Hermitian Hamiltonian terms $\\{h_{k}\\}$ of neither necessarily
small nor geometrically local support defined on a quantum lattice equipped
with a Hilbert space ${\cal H}=(\mathbb{C}^{d})^{\otimes n}$.

Here, $n$ is the number of degrees of freedom of finite local dimension $d$.
To access the Hamiltonian, we assume to have oracles $\mathbb{O}_{k}(t)$
implementing time evolutions under each term in the Hamiltonian

$\mathbb{O}_{k}(t)\ket{\psi}=e^{-\mathrm{i}h_{k}t}\ket{\psi}$ (2)

for any $k\in[1,L]$ and $t\in\mathbb{R}$. In digital quantum simulation, these
oracles are built from Clifford gates and phase rotations, but we do not
assume that such a decomposition is available to us, as the oracles could be
implemented by the time evolution of a programmable device. We only require to
have control over the implementation, i.e., that we can apply the oracle
$\mathbb{O}_{k}(t)$ depending on the state of a subsystem encoded as a qubit
as

$\displaystyle\ket{0}\\!\\!\bra{0}\otimes\mathbb{I}+\ket{1}\\!\\!\bra{1}\otimes\mathbb{O}_{k}(t)\,.$
(3)

Now, quantum simulation aims at making reliable predictions of the expectation
values of observables $O$ (which in the ideal case are local with a small
support on the lattice, but again, this is not a necessity) at later times
$T>0$

$\langle O(T)\rangle\coloneqq{\rm tr}(U(T)\rho U^{\dagger}(T)O)$ (4)

for initial quantum states $\rho$, where $U$ is the time evolution operator.
In practice, this commonly means to establish ways to approximate the
corresponding time evolution operator

$U(T)=\mathrm{e}^{-\mathrm{i}HT}\,.$ (5)

Although it is worth stressing that a-priori information about the time $T$
and properties of the initial state $\rho$, as well as features of locality of
the underlying Hamiltonian $H$ can be exploited to come up with highly
specialized approximations, we do not assume any underlying structure or
limitations in this work.

A simple, yet effective approach to approximating the time evolution operator
is that of using product formulas, which make use of $N$ consecutive
evolutions of individual Hamiltonian terms $h_{k_{j}}$ by associated short
time intervals $\alpha_{j}t$, defined with real coefficients $\alpha_{j}$ such
that $\sum_{j}\alpha_{j}=1$. Partial backwards evolutions are explicitly
allowed, i.e. $\alpha_{j}<0$ for some $j$, even though $t>0$. In general, such
a formula is defined as a product of time evolutions

$\mathrm{e}^{-\mathrm{i}Ht}\approx\prod_{j=1}^{N}\mathrm{e}^{-\mathrm{i}\alpha_{j}\,h_{k_{j}}t}\,$
(6)

where each single term evolution can be implemented via $\mathbb{O}_{k}(t)$.
Here, we distinguish between a direct evolution by time $t$ and a repeated
evolution of short time slices $t/r<1$ such that

$\mathrm{e}^{-\mathrm{i}Ht}\approx\left(\prod_{j=1}^{N}\mathrm{e}^{-\mathrm{i}\alpha_{j}\,h_{k_{j}}\frac{t}{r}}\right)^{r}\,.$
(7)

The most accurate known formulas of this type are _Trotter-Suzuki formulas_
[26]. They are recursively defined for any positive integer $\chi$ and any
time $t$ by

$\displaystyle S_{1}(t)=\prod_{k=1}^{L}\mathrm{e}^{-\mathrm{i}h_{k}t}\,,$ (8)
$\displaystyle
S_{2}(t)=\prod_{k=1}^{L}\mathrm{e}^{-\mathrm{i}h_{k}t/2}\prod_{k=L}^{1}\mathrm{e}^{-\mathrm{i}h_{k}t/2}\,,$
(9) $\displaystyle
S_{2\chi}(t)=\Big{(}S_{2\chi-2}(s_{2\chi-2}t)\Big{)}^{2}S_{2\chi-2}\Big{(}[1-4s_{2\chi-2}]t\Big{)}\Big{(}S_{2\chi-2}(s_{2\chi-2}t)\Big{)}^{2}$
(10)

with $s_{2\chi}\coloneqq(4-4^{1/(2\chi+1)})^{-1}$ for any positive integer
$\chi$. This specific choice of $s_{2\chi}$ ensures that the Taylor series of
$S_{2\chi}(t)$ matches that of the actual time evolution $U(t)$ up to
$\mathcal{O}(t^{2\chi+1})$, which makes it a good approximation for $t\ll 1$.
It is important to stress that constructing an $\mathcal{O}(t^{2\chi+1})$
approximation requires the application of $2\cdot 5^{\chi-1}(L-1)+1$
individual oracles $\mathbb{O}_{k}(t)$, where $L$ is again the number of
Hamiltonian terms. More generally, we would represent them by oracle calls and
we will refer to the number of exponentials $N_{\rm exp}$. Using repetitions
as in Eq. (7), the number of expoentials of the algorithm is of
$\mathcal{O}(rL5^{\chi-1})$ to approximate the actual time evolution $U(T)$ up
to an accuracy of $\varepsilon\in\mathcal{O}((T/r)^{2\chi+1})$ in operator
distance. The exponential dependence of the circuit depth on $\chi$ can pose a
challenge for real-world implementations, and especially so within the NISQ
regime. Altogether, the number of oracle calls needed to simulate the
evolution up to an error $\varepsilon\leq 1$ is upper bounded by

$N_{\rm
oracle}\leq\frac{2L5^{2\chi}\left(L\tau\right)^{1+1/2\chi}}{\varepsilon^{1/2\chi}},$
(11)

where $\tau:=\sum_{k=1}^{L}\lVert h_{k}\rVert T$ and $\lVert\cdot\rVert$
denotes the operator norm [4]. An alternative complexity bound using
commutator relations of the individual Hamiltonian terms can be found in Ref.
[11].

A known technique to decrease the number of required short time evolutions,
i.e., oracle calls, is the use of multi-product formulas [27, 28]. While the
Trotter-Suzuki approximation cancels erroneous contributions of higher-order
terms by adding backwards evolutions, multi-product formulas achieve the same
cancellations by superposing different product formulas. Conventionally, one
employs multi-product formulas that describe the same time evolution (up to a
fixed order of $t$), but whose erroneous higher-order contributions are of
different strength and can thus be made to cancel. This is achieved by
approximating the time evolution to the same Trotter-Suzuki order, but
considering product formulas different in the number of time slices [27, 29].
The Childs and Wiebe multi-product formula discussed in Ref. [13] is of the
form

$M_{K,2\chi}(t)\coloneqq\sum_{q=1}^{K+1}C_{q}S_{2\chi}(t/\ell_{q})^{\ell_{q}},$
(12)

where $K$ is an integer defining a cutoff and $\\{\ell_{q}\\}$ a set of
pairwise different integers. The coefficients $\\{C_{q}\\}$ are determined via

$\left(\begin{matrix}1&1&\,&\cdots&\,&1\\\
\ell_{1}^{-2\chi}&\ell_{2}^{-2\chi}&&\cdots&&\ell_{K+1}^{-2\chi}\\\
\ell_{1}^{-2\chi-2}&\ell_{2}^{-2\chi-2}&&\cdots&&\ell_{K+1}^{-2\chi-2}\\\
\vdots&\vdots&&\ddots&&\vdots\\\
\ell_{1}^{-2(K+\chi-1)}&\ell_{2}^{-2(K+\chi-1)}&&\cdots&&\ell_{K+1}^{-2(K+\chi-1)}\end{matrix}\right)\left(\begin{matrix}C_{1}\\\
C_{2}\\\ C_{3}\\\ C_{4}\\\ \vdots\\\ C_{K+1}\\\
\end{matrix}\right)=\left(\begin{matrix}1\\\ 0\\\ 0\\\ 0\\\ \vdots\\\
0\end{matrix}\right),$ (13)

ensuring the error terms in the multi-product formula vanish up to
$\mathcal{O}(t^{2(K+\chi)})$, resulting in

$\left\lVert
M_{K,\chi}-U(t)\right\rVert=\mathcal{O}\left(t^{2(K+\chi)+1}\right).$ (14)

Unlike Childs and Wiebe, we will henceforth use the simplest version achieved
by setting $\ell_{q}=q$ for all $q$.

In Definitions 2 and 3 and Section IV, we employ a different approach for
multi-product formulas and develop two techniques whose errors scale with
$\mathcal{O}(t^{2\chi R+1})$, and where $R$ is comparable to $K+1$.

Multi-product formulas were firstly used for quantum simulation by Childs and
Wiebe [13], who developed the _linear-combinations-of-unitaries (LCU)_
approach to directly implement multi-product formulas on a quantum system.
Note that sums of unitaries are not inherently physical operations since the
unitary group is not closed under addition. Childs and Wiebe use a non-
deterministic approach to implement multi-product formulas that can lead to
large algorithmic depth. Berry et al. [15] have implemented LCU for a
truncated Taylor series nearly perfectly but their use of oblivious amplitude
amplification requires an additional register and a complex state preparation
procedure. Recently, Ref. [14] has improved the condition number of multi-
product formulas and thus extended the use of LCU by amplitude amplification.

Using _randomized sampling_ such as proposed in Algorithm 1, we can circumvent
the need for these potentially deep circuits. Randomized compiling of multi-
product formulas entails sampling among the individual terms of Eq. (12)
requiring the implementation of just twice the number of oracle calls, rather
than the entire linear combination and the overhead for its implementation.

##  III Randomized sampling using multi-product formulas

LCU-type algorithms, despite having a finite failure probability, are
deterministic when successfully applied. Apart from such deterministic
implementations, there is a recent interest in randomized algorithms [21, 22,
23, 24]. While these novel results improve the performance of product formulas
by randomizing the order of the short time evolution, we propose the
alternative setting of randomized compiling. With the goal of reducing the
circuit depth, we introduce the setting of randomized sampling: Given an
observable $O$ and a quantum system in the initial state $\rho$, our goal is
to find the expectation value of the observable after a time evolution
$U=\mathrm{e}^{-\mathrm{i}Ht}$ governed by the Hamiltonian
$H=\sum_{k=1}^{L}h_{k}$. That is, we wish to compute

$\langle O\rangle={\rm tr}(OU\rho U^{\dagger})\,.$ (15)

In the following, we describe a randomized algorithm for such a task, give
convergence guarantees and present novel multi-product formulas suited for
this framework. Note that we will require one part of our system to be encoded
as a qubit, while the rest acts as the simulator.

Figure 2: Quantum circuit for randomized sampling of some observable $O$
according to Algorithm 1. When $V_{\bullet}$ and $V_{\circ}$ are chosen from
$\set{V_{1},\ldots,V_{M}}$ with probabilities $\set{p_{1},\ldots,p_{M}}$ such
that $\sum_{k=1}^{M}p_{k}V_{k}=V$, the expectation value of the measurement
outcome is equal to $\operatorname{tr}({OV\rho V^{\dagger}})$.

######  Algorithm 1 (Randomized sampling).

Given an observable $O$ and an ensemble
$\mathcal{V}=\set{(V_{k},p_{k})}_{k=1}^{M}$ of $M$ unitaries $\\{V_{k}\\}$ and
corresponding probabilities $\\{p_{k}\\}$, we consider $N$ independent
repetitions of the following protocol:

  1. 1.

Prepare $\ket{+}\\!\\!\bra{+}\otimes\rho\,$.

  2. 2.

Sample $V_{\circ}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{V}$ and
apply the anti-controlled unitary $\ket{0}\\!\\!\bra{0}\otimes
V_{\circ}+\ket{1}\\!\\!\bra{1}\otimes\mathbb{I}\,$.

  3. 3.

Sample $V_{\bullet}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{V}$ and
apply the controlled unitary
$\ket{0}\\!\\!\bra{0}\otimes\mathbb{I}+\ket{1}\\!\\!\bra{1}\otimes
V_{\bullet}\,$.

  4. 4.

Perform a single shot of the POVM measurement associated with the observable
$X\otimes O\,$.

  5. 5.

Store the measurement outcomes $\\{o_{j}\\}$.

The corresponding quantum circuit is shown in Fig. 2. As we show in more
detail in Section IV, we find the following result for infinitely many runs of
the algorithm:

######  Corollary 1 (Sample mean convergence).

Let $O$ be an observable with $\lVert O\rVert\leq 1$, let $\rho$ be an initial
state and let $\mathcal{V}=\set{(V_{k},p_{k})}_{k=1}^{M}$ be an ensemble of
unitaries $\set{V_{i}}$ and corresponding probabilities $\set{p_{i}}$, such
that

$\mathbb{E}_{\mathcal{V}}[V_{k}]=\sum_{k=1}^{M}p_{k}V_{k}=V\,.$ (16)

Then, the sample mean of Algorithm 1 converges to the expectation value of the
random measurement outcomes $\\{o_{j}\\}_{j}$, given by

$\lim_{N\to\infty}\frac{1}{N}\sum_{j=1}^{N}o_{j}=\mathrm{tr}\left(OV\rho
V^{\dagger}\right)\,.$ (17)

Furthermore, since an infinite number of measurements is infeasible, we give
an expression for the number $N$ of repetitions of Algorithm 1 sufficient for
achieving a desired accuracy and confidence for the goal of randomized
sampling. Using Hoeffding’s inequality [30], we can formulate the following
result:

######  Theorem 1 (Randomized implementation of sums of unitaries).

Let $O$ be an observable with $\lVert O\rVert\leq 1$, let $\rho$ an initial
state and let $\mathcal{V}=\set{(V_{k},p_{k})}_{k=1}^{M}$ be an ensemble of
unitaries $\set{V_{k}}$ and corresponding probabilities $\set{p_{k}}$, such
that

$\mathbb{E}_{\mathcal{V}}[V_{k}]=\sum_{k=1}^{M}p_{k}V_{k}=V\,.$ (18)

Then, for a fixed accuracy $\varepsilon\in(0,1)$ and confidence
$\delta\in(0,1)$, a total of

$N\geq\frac{2\log\left(2/\delta\right)}{\varepsilon^{2}}$ (19)

repetitions of Algorithm 1 suffice to accurately approximate the expectation
value $\mathrm{tr}\left(OV\rho V^{\dagger}\right)$ using the sampling mean
estimator, i.e.,

$\left\lvert\frac{1}{N}\sum_{j=1}^{N}o_{j}-\mathrm{tr}\left(OV\rho
V^{\dagger}\right)\right\rvert\leq\varepsilon$ (20)

with probability at least $1-\delta$.

So far, this framework is fairly general: Given a set $\mathcal{V}$, the
algorithm samples from the linear transform with the unitary $V$. For quantum
simulation, $V$ (a decent approximation to $U$) and $\mathcal{V}$ still have
to be found. This is where multi-product formulas come into play. Multi-
product formulas approximate $U$ by a weighted sum of product formulas, and so
the featured product formulas and their weights could make up the set
$\mathcal{V}$. The problem is that we require $\\{p_{k}\\}$ to form a
probability distribution, i. e., $p_{k}>0$ and $\sum_{k=1}^{M}p_{k}=1$. This
disqualifies us from identifying the sets $\\{p_{k}\\}$ with $\\{C_{q}\\}$
straightforwardly. While $\sum_{q}C_{q}=1$ holds for any multi-product formula
per construction (see the first component of Eq. (13)), the sign of some
$C_{q}$ will be negative [31]. We thus have to absorb the signs of these
negative $C_{q}$ into their unitaries $V_{q}$, but then we find that the
operation is not normalized as $\sum_{q}\lvert C_{q}\rvert>1$. At this point
we introduce the quantity $\Xi\coloneqq\sum_{q}|C_{q}|$, that we call
resolution factor. The resolution is practically used to define a probability
distribution via $p_{k}=\left|C_{k}\right|/\Xi$ such that Algorithm 1 samples
from

$\frac{1}{\Xi^{2}}\mathrm{tr}(OU\rho U^{\dagger})\,.$ (21)

As the time evolution is now approximated by $\Xi V$ rather than $V$, $\Xi$
factors into the number of required circuit evaluations. With $\Xi>1$
requiring to run the algorithm more often to achieve an error comparable with
the situation where $\Xi=1$, this resolution factor can be regarded as a
penalty. However, especially in near-term and medium-term applications, a
deeper circuit may pose a bottleneck much tighter than a slightly higher
number of circuit evaluations does.

We can make the following guarantees for approximating time evolutions in
randomized sampling with a nontrivial resolution and a finite number of
measurements:

######  Theorem 2 (Approximating unitaries with a finite number of
measurements,).

Let $U$ be a unitary, $O$ an observable with $\lVert O\rVert\leq 1$, $\rho$
some initial state and let $\mathcal{V}=\set{(V_{k},p_{k})}_{k=1}^{M}$ be an
ensemble of unitaries $\set{V_{k}}$ with a corresponding probability
distribution $\set{p_{k}}$ such that
$\mathbb{E}_{\mathcal{V}}[V_{k}]=\sum_{k=1}^{M}p_{k}V_{k}=V$. Let there be a
constant factor $\Xi\in\mathbb{R}_{+}$ such that

$\lVert\Xi V-U\rVert\leq\varepsilon/3\,.$ (22)

Then, it is sufficient to run

$N\geq 8\ln\left(2/\delta\right)\left(\frac{\Xi}{\varepsilon}\right)^{2}$ (23)

repetitions of Algorithm 1 to achieve

$\left\lvert\frac{1}{N}\sum_{j=1}^{N}o_{j}-\mathrm{tr}(OU\rho
U^{\dagger})\right\rvert\leq(1+\Xi)\,\varepsilon\,,$ (24)

with probability at least $1-\delta$.

The multi-product formulas of Childs and Wiebe can easily be adapted to this
randomized sampling scheme. Since we can ignore constraints encountered in LCU
techniques, our sole interest lies in a low algorithmic depth and a moderate
resolution and so one can set $\ell_{q}=q$ for all $q=1,\dots,K+1$ in Eq. (12)
and (13). However, those multi-product formulas make relatively little use of
the Trotter-Suzuki order $2\chi$ of its components. In fact the only reason
not to minimize the depth by using $S_{2}(\cdot)$ blocks is the high
resolution factor. To ease the impact of the resolution at higher orders of
the time evolution, we present a family of novel multi-product formulas with
an improved scaling in the Trotter-Suzuki order. While Childs and Wiebe’s
approach manipulates higher-order error terms, we employ a construction to
modulate the entire Taylor expansion.

######  Definition 1 (Linear combination of time evolution operators).

Let $\chi\geq 1$ and $R\geq 1$ be integers and $S_{2\chi}(t)$ a product
formula approximation to $U(t)$ with $\left\lVert
S_{2\chi}(t)-U(t)\right\rVert\in\mathcal{O}(t^{2\chi+1})$. Then, for any
$t\in\mathbb{R}$ and arbitrary
$\boldsymbol{b}=(b_{0},\,b_{1},\,\dots\,,\,b_{2\chi R})^{\top}$,
$\boldsymbol{\nu}=(\nu_{0},\,\nu_{2},\,\dots\,,\,\nu_{2\chi
R})^{\top}\in\mathbb{R}^{2\chi R+1}$ , we define the multi-product formula
$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu},\boldsymbol{b},t)$ as

$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu},\boldsymbol{b},t)\coloneqq\sum_{q=0}^{2\chi
R}C_{q}(\boldsymbol{\nu},\boldsymbol{b})\,S_{2\chi}\\!\left(b_{q}t\right),$
(25)

where

$\left(\begin{matrix}1&1&\,&\cdots&\,&1\\\ b_{0}&b_{1}&&\cdots&&b_{2\chi R}\\\
b_{0}^{2}&b_{1}^{2}&&\cdots&&b_{2\chi R}^{2}\\\
b_{0}^{3}&b_{1}^{3}&&\cdots&&b_{2\chi R}^{3}\\\
\vdots&\vdots&&\ddots&&\vdots\\\ \,b_{0}^{2\chi R}&\,b_{1}^{2\chi
R}&&\cdots&&\,b_{2\chi R}^{2\chi
R}\end{matrix}\right)\left(\begin{matrix}C_{0}\\\ C_{1}\\\ C_{2}\\\ C_{3}\\\
\vdots\\\ C_{2\chi R}\end{matrix}\right)=\left(\begin{matrix}\nu_{0}\\\
\nu_{1}\\\ \nu_{2}\\\ \nu_{3}\\\ \vdots\\\ \nu_{2\chi R}\end{matrix}\right)$
(26)

for some $(C_{0},\,C_{1},\,\dots\,,\,C_{2\chi R})^{\top}\in\mathbb{R}^{2\chi
R+1}$.

Here, the coefficients $\\{C_{q}\\}$ are related to the parameters
$\boldsymbol{b}$ and $\boldsymbol{\nu}$ via the inverse of the corresponding
Vandermonde matrix in Eq. (26) as is discussed in detail in Section IV. The
parameters $\boldsymbol{\nu}$ will have an important role in the following
construction, whereas $\boldsymbol{b}$ are tuned to strike a balance between
$\Xi$ and $\varepsilon$. Note that a similar construction has recently been
used to approximate energy measurements in the fully analog setting [32]. We
can now turn towards the definition of the first new multi-product formula.

######  Definition 2 (Matching multi-product formula).

Let $\chi\geq 1$ and $R\geq 1$ be integers and
$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu},\boldsymbol{b},t)$ as defined in
Definition 1. Then, for any $t\in\mathbb{R}$, we define the multi-product
formula $\widetilde{M}_{2\chi,R}^{(\mathrm{m})}(t)$ as

$\widetilde{M}_{2\chi,R}^{(\mathrm{m})}(t)\coloneqq\prod_{r=1}^{R}\mathcal{L}_{2\chi,R}\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t\right)$
(27)

with $\nu_{0}^{(r)}=1$ and $\nu_{k}^{(r)}=0$ for $k>2\chi$ and all $r=1\dots
R$. The remaining $\left\\{\nu_{q}^{(r)}\right\\}$ are fixed by

$\sum_{k_{1}+k_{2}+\ldots+k_{R}=k}\frac{\nu_{k_{1}}^{(1)}\nu_{k_{2}}^{(2)}\cdots\nu_{k_{R}}^{(R)}}{k_{1}!\,k_{2}!\,\cdots\,k_{R}!}=\frac{1}{k!}$
(28)

for all $0\leq k\leq 2\chi R$, whereas the parameters
$\\{\boldsymbol{b}^{(r)}\\}$ can be chosen arbitrarily.

Multiplying the coefficients $C_{q}$ of the individual building blocks
$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t)$ of the
matching multi-product formula results in a resolution factor of

$\Xi^{(\mathrm{m})}\coloneqq\prod_{r}\left(\sum_{q}\left\lvert
C_{q}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)}\\!\right)\right\rvert\right).$
(29)

For the time evolution to work, the set of vectors
$\left\\{\boldsymbol{\nu}^{(r)}\right\\}$ has to be found satisfying the
constraint in Eq. (28) for different $R$ and $2\chi$. When such a
decomposition is unavailable, we can employ a second version of the multi-
product formula in which the $\left\\{\boldsymbol{\nu}^{(r)}\right\\}$ are
already determined at the cost of a higher $\Xi$.

######  Definition 3 (Closed-form multi-product formula).

Let $\chi\geq 1$ and $R\geq 1$ be integers and
$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu},\boldsymbol{b},t)$ as defined in
Definition 1. Then for any $t\in\mathbb{R}$, we define the multi-product
formula $\widetilde{M}_{2\chi,R}^{(\mathrm{cf})}(t)$ as

$\widetilde{M}_{2\chi,R}^{(\mathrm{cf})}(t)\coloneqq\sum_{r=1}^{R}\left(\mathcal{L}_{2\chi,R}\\!\left(\boldsymbol{\nu}^{(0)},\boldsymbol{b}^{(0)},t\right)\right)^{r-1}\mathcal{L}_{2\chi,R}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t\right)$
(30)

with

$\displaystyle\nu^{(0)}_{k}$
$\displaystyle=\begin{cases}1,&\text{for}\;k=2\chi\\\
0,&\text{else},\end{cases}$ (31) $\displaystyle\nu^{(1)}_{k}$
$\displaystyle=\begin{cases}1,&\text{for}\;k\leq 2\chi\\\
0,&\text{else},\end{cases}$ (32) $\displaystyle\nu^{(n)}_{k}$
$\displaystyle=\begin{cases}\frac{k!\left((2\chi)!\right)^{n-1}}{(2\chi(n-1)+k)!},&\text{for}\;0<k\leq
2\chi\\\ 0,&\text{else},\end{cases}$ (33)

for all $\;0\leq k\leq 2\chi R\;$ and $\;1<n\leq R$. The parameters
$\\{\boldsymbol{b}^{(r)}\\}$ can be chosen arbitrarily.

Multiplying the coefficients of the individual
$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t)$,
results in a resolution factor for the closed-form multi-product formula of

$\Xi^{(\mathrm{cf})}:=\sum_{r}\left(\sum_{q}\left\lvert
C_{q}\\!\left(\boldsymbol{\nu}^{(0)},\boldsymbol{b}^{(0)}\\!\right)\right\rvert\right)^{r-1}\left(\sum_{q}\left\lvert
C_{q}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)}\\!\right)\right\rvert\right).$
(34)

For these multi-product formulas, we can prove the following two results:

######  Corollary 2 (Low depth unitary approximation).

Let $\widetilde{M}_{2\chi,R}(t)$ be a multi-product formula constructed
according to Definition 2 or 3. We then find that

$\left\lVert
U(t)-\widetilde{M}_{2\chi,R}(t)\right\rVert=\mathcal{O}\left(t^{2\chi
R+1}\right).$ (35)

######  Theorem 3 (Error bound for custom multi-product formulas).

For $\widetilde{M}_{2\chi,R}(t)$ being a multi-product formula constructed
according to Definitions 2 or 3, we have

$\left\lVert U(t)-\widetilde{M}_{2\chi,R}(t)\right\rVert\leq\left(1+\zeta
g_{\chi}^{2\chi R+1}\right)\frac{\left(\Lambda t\right)^{2\chi
R+1}}{\left(2\chi R+1\right)!}\,,$ (36)

with $g_{\chi}\ \coloneqq\frac{4\chi}{5}\left(\frac{5}{3}\right)^{\chi-1}$,
$\Lambda\coloneqq\sum\limits_{k}\left\lVert h_{k}\right\rVert$ and
corresponding

$\displaystyle\zeta^{(\mathrm{m})}$
$\displaystyle\coloneqq\sum_{q_{1}=0}^{2\chi R}\sum_{q_{2}=0}^{2\chi
R}\cdots\sum_{q_{r}=0}^{2\chi R}\left\lvert
C_{q_{1}}^{(1)}C_{q_{2}}^{(2)}\cdots C_{q_{r}}^{(r)}\right\rvert\left(\lvert
b_{q_{1}}^{(1)}\rvert+\ldots+\lvert b_{q_{r}}^{(r)}\rvert\right)$ (37)
$\displaystyle\zeta^{(\mathrm{cf})}$
$\displaystyle\coloneqq\sum_{r=1}^{R}\sum_{q_{1}=0}^{2\chi
R}\sum_{q_{2}=0}\cdots\sum_{q_{r}=0}^{2\chi R}\left\lvert
C_{q_{1}}^{(0)}C_{q_{2}}^{(0)}\cdots
C_{q_{r-1}}^{(0)}C_{q_{r}}^{(r)}\right\rvert\left(\lvert
b_{q_{1}}^{(0)}\rvert+\ldots+\lvert b_{q_{r-1}}^{(0)}\rvert+\lvert
b_{q_{r}}^{(r)}\rvert\right),$ (38)

for matching and closed-form multi-product formulas respectively.

Consequently, these formulas can be used for randomized sampling via Theorem
2. A numerical comparison of the error bounds of all presented formulas can be
found in Fig. 4, while an overview of circuit depths, error scaling and error
bounds is given in Table 1.

The matching and closed-form multi-product formulas both rely on products of
$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t)$. For
Algorithm 1, these products can either be expanded and the resulting operators
and coefficients be identified with the sets $\\{V_{k}\\}$, $\\{p_{k}\\}$, or
sets of
$(V_{\circ}^{(1)},V_{\bullet}^{(1)}),\dots,(V_{\circ}^{(R)},V_{\bullet}^{(R)})$
be drawn and applied as in Fig. 3. For the matching multi-product formula, all
$V_{\circ}^{(r)}$ and $V_{\bullet}^{(r)}$ are drawn independently from each
other while for the closed-form multi-product formula, the set from which
$V_{\circ/\bullet}^{(r+1)}$ is drawn depends on the set from which
$V_{\circ/\bullet}^{(r)}$ has been sampled.

Figure 3: Quantum circuit for randomized sampling with a multi-product formula
$\prod_{r=1}^{R}\left(\sum_{k}p^{(r)}_{k}V^{(r)}_{k}\right)$. The unitaries
$V_{\circ}^{(r)},V_{\bullet}^{(r)}$ are drawn independently at random from
$\\{V^{(r)}_{k}\\}_{k}$ according to the probabilities $\\{p^{(r)}_{k}\\}_{k}$
for all $r=1\dots R$, such that the circuit samples from $\mathbb{E}(\langle
O\rangle)$.

Furthermore, since the above framework randomizes only over product formulas,
we can also think of a doubly-stochastic version in which the short-term
evolutions comprising the product formulas are sampled randomly as well. In
that sense, it is allowed to construct the multi-product formulas with
building blocks of

$\displaystyle\widehat{S}_{2}(t)=\frac{1}{2}\left(S_{1}^{\phantom{\dagger}}(t)+S_{1}^{\dagger}(-t)\right)\,,$
(39)

rather than (9). This is not possible for Childs and Wiebe’s multi-product
formulas, which require symmetric building blocks.

##  IV Proofs

In this section, we provide proofs of the statements presented above. In
Section IV.1, we are proving our statements regarding the errors and
uncertainty of the randomized sampling procedure. We then turn our attention
to the matching and closed-form multi-product formula, providing the intuition
behind their construction and verifying their error scaling in Section IV.2.
We finally analyze the upper bounds for their error (which is the error of the
average time evolution they describe) in Section IV.3.

###  IV.1 Randomized sampling

Let us start by following Algorithm 1 step by step. The controlled and anti-
controlled applications of $V_{\circ}$ and $V_{\bullet}$ are defined as

$\displaystyle\mathrm{\overline{C}}V_{\circ}$
$\displaystyle\coloneqq\ket{0}\\!\\!\bra{0}\otimes
V_{\circ}+\ket{1}\\!\\!\bra{1}\otimes\mathbb{I},$ (40)
$\displaystyle\mathrm{C}V_{\bullet}$
$\displaystyle\coloneqq\ket{0}\\!\\!\bra{0}\otimes\mathbb{I}+\ket{1}\\!\\!\bra{1}\otimes
V_{\bullet}\,,$ (41)

and the initial state $\rho$ will be transformed to

$\displaystyle\widetilde{\rho}$
$\displaystyle=\mathrm{C}V_{\bullet}\mathrm{\overline{C}}V_{\circ}\,\left(\left|+\right\rangle\\!\\!\left\langle+\right|\otimes\rho\right)\,(\mathrm{\overline{C}}V_{\circ})^{\dagger}(\mathrm{C}V_{\bullet})^{\dagger}$
(42) $\displaystyle=\frac{1}{2}\left(\ket{0}\\!\\!\bra{0}\otimes
V^{\phantom{\dagger}}_{\circ}\rho
V_{\circ}^{\dagger}+\ket{0}\\!\\!\bra{1}\otimes
V^{\phantom{\dagger}}_{\circ}\rho
V_{\bullet}^{\dagger}+\ket{1}\\!\\!\bra{0}\otimes
V^{\phantom{\dagger}}_{\bullet}\rho
V_{\circ}^{\dagger}+\ket{1}\\!\\!\bra{1}\otimes
V^{\phantom{\dagger}}_{\bullet}\rho V_{\bullet}^{\dagger}\right)$ (43)

after the third step of the protocol. The expectation value for the succeeding
measurement of $X\otimes O$ is then given by

$\mathrm{tr}\left(X\otimes
O\widetilde{\rho}\right)=\frac{1}{2}\mathrm{tr}\left(O(V^{\phantom{\dagger}}_{\circ}\rho
V_{\bullet}^{\dagger}+V^{\phantom{\dagger}}_{\bullet}\rho
V_{\circ}^{\dagger})\right).$ (44)

We now use that $V_{\circ}$ and $V_{\bullet}$ are sampled independently from
the same ensemble

$V_{\circ},V_{\bullet}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathcal{V}=\set{(V_{k},p_{k})}_{k=1}^{M}$ (45)

such that

$\mathbb{E}_{\mathcal{V}}[V_{\circ}]=\mathbb{E}_{\mathcal{V}}[V_{\bullet}]=\sum_{k=1}^{M}p_{k}V_{k}=V\,.$
(46)

Then, we can combine the quantum average form Eq. (44) with the classical
average from Eq. (46) to conclude that the expectation value of the random
single-shot measurement outcome $o$ is given by

$\displaystyle\mathbb{E}[o]$
$\displaystyle=\mathbb{E}_{\mathcal{V}}\left[\mathrm{tr}\left(X\otimes
O\widetilde{\rho}\right)\right]=\mathrm{tr}\left(OV\rho V^{\dagger}\right),$
(47)

which proves Corollary 1.

For real world applications we will never achieve a perfect expectation value.
It is therefore vital to inspect the single-shot behavior of Algorithm 1 and
give an estimate for the number of iterations required to achieve the desired
precision. Note that

$\displaystyle O=\sum_{j=1}^{l}\omega_{j}P_{j}$ (48)

is an observable with eigenvalues $\lvert\omega_{j}\rvert\leq 1$ for $\lVert
O\rVert\leq 1$ and $\\{P_{j}\\}$ is a set of projectors. We thus find that
according to Born’s rule, the POVM measurement defined by $X\otimes O$ has
$2l$ possible measurement outcomes $o_{j}$ given by $\omega_{j}\in[-1,1]$.
Consequently, by applying Hoeffding’s inequality, we can prove the validity of
Theorem 1.

###### Proof of Theorem 1.

Since the single-shot measurement outcomes $\\{o_{j}\\}$ obtained from $N$
iterations of Algorithm 1 are individually bounded by the interval $[-1,1]$,
Hoeffding’s inequality states

$\mathrm{Prob}\left(\left\lvert\frac{1}{N}\sum_{j=1}^{N}o_{j}-\mathrm{tr}\left(OV\rho
V^{\dagger}\right)\right\rvert\geq\varepsilon\right)\leq
2\exp{\left(-\frac{N\varepsilon^{2}}{2}\right)}.$ (49)

For a fixed accuracy $\varepsilon\in(0,1)$ and confidence $\delta\in(0,1)$, it
is therefore sufficient to perform

$N\geq\frac{2\ln\left(2/\delta\right)}{\varepsilon^{2}}$ (50)

repetitions of Algorithm 1 to achieve the desired accuracy and confidence,
concluding the proof. ∎

Given recent developments, for example in Ref. [33], one might wonder whether
substantial improvements are possible by using a more refined estimator, most
notably median of means. Unfortunately, this is not very realistic. For
observables that obey $O^{2}=\mathbb{I}$, such as local and global Pauli
observables, the variance of a single-shot outcome $o_{j}\in\left[-1,1\right]$
becomes
$\mathrm{Var}\left[o_{j}\right]=1-\mathbb{E}\left[o_{j}\right]=\mathcal{O}(1)$.
In this case, the variance is of the same order as the magnitude and it is
impossible to (asymptotically) improve over Hoeffding’s inequality (asymptotic
normality) [34]. Now, we additionally assume that $V$ times the resolution
$\Xi$ approximates the time evolution operator $U$, i.e.,

$\lVert\Xi V-U\rVert\leq\hat{\varepsilon}\,.$ (51)

Furthermore, we will use the following result.

######  Lemma 1 (Closeness of expectation values).

Let $U$ be unitary and $V$ and $\Xi$ as given above such that Eq. (51) holds.
Furthermore, fix a state $\rho$ and observable $O$ with $\lVert O\rVert\leq
1$. Then,

$\left\lvert\mathrm{tr}\left(OU\rho
U^{\dagger}\right)-\Xi^{2}\mathrm{tr}\left(OV\rho
V^{\dagger}\right)\right\rvert\leq 3\hat{\varepsilon}\,.$ (52)

###### Proof.

We begin by defining $\widetilde{U}$ as the difference between the exact and
approximated time evolution

$\Xi V=\Xi\sum_{k=1}^{M}p_{k}V_{k}=U+\widetilde{U}.$ (53)

According to Eq. (51) we find that
$\lVert\widetilde{U}\rVert\leq\hat{\varepsilon}$. Consequently,

$\displaystyle\left\lvert\mathrm{tr}\left(OU\rho
U^{\dagger}\right)-\Xi^{2}\mathrm{tr}\left(OV\rho
V^{\dagger}\right)\right\rvert$
$\displaystyle=\left\lvert\mathrm{tr}\left(\widetilde{U}^{\dagger}OU\rho\right)+\mathrm{tr}\left(U^{\dagger}O\widetilde{U}\rho\right)+\mathrm{tr}\left(\widetilde{U}^{\dagger}O\widetilde{U}\rho\right)\right\rvert$
$\displaystyle\leq 2\lVert\widetilde{U}\rVert\lVert U\rVert\lVert
O\rVert+\lVert\widetilde{U}\rVert^{2}\lVert O\rVert$ (54) $\displaystyle\leq
2\hat{\varepsilon}+\hat{\varepsilon}^{2}$ (55) $\displaystyle\leq
3\hat{\varepsilon}.$ (56)

∎

By combining the insights from the previous discussion, we can now tackle
Theorem 2:

###### Proof of Theorem 2.

Add and subtract $\mathrm{tr}\left(OV\rho V^{\dagger}\right)$ to find

$\left\lvert\frac{\Xi^{2}}{N}\sum_{j=1}^{N}o_{j}-\mathrm{tr}\left(OU\rho
U^{\dagger}\right)\right\rvert\leq\Xi^{2}\left\lvert\frac{1}{N}\sum_{j=1}^{N}o_{j}-\mathrm{tr}\left(OV\rho
V^{\dagger}\right)\right\rvert+\left\lvert\vphantom{\sum_{j=1}^{N}}\Xi^{2}\mathrm{tr}\left(OV\rho
V^{\dagger}\right)-\mathrm{tr}\left(OU\rho U^{\dagger}\right)\right\rvert.$
(57)

By applying Theorem 1 with accuracy $\varepsilon/\Xi$ to the first term and
using Lemma 1 with $\varepsilon=3\hat{\varepsilon}$, we find that this upper
bound is given by $(1+\Xi)\varepsilon$. ∎

###  IV.2 Construction of the new multi-product formulas

In Definitions 2 and 3, we propose two alternatives to Childs and Wiebe’s
multi-product formulas, which reduce the number of circuit evaluations
required for a randomized implementation and have improved error scaling. In
the following, we will motivate their construction and stress the advantages
they provide. First of all, note that every approximation of the exact time
evolution $\widetilde{U}(t)\approx\mathrm{e}^{-\mathrm{i}tH}$ can be written
as a sum of operators $\hat{A}_{k}$ such that

$\displaystyle\widetilde{U}(t)=\sum_{k=0}^{\infty}\hat{A}_{k}\,t^{k}\,,$ (58)

with $\hat{A}_{0}=\mathbb{I}$. For any approximated time evolution with an
error of $\mathcal{O}(t^{m+1})$, we find that $\hat{A}_{k}=(iH)^{k}/{k!}$ for
all $k\leq m$. In other words, the Taylor expansion of the exact time
evolution and its approximation has the same Taylor expansion for the first
$m$ non-trivial terms.

For the standard Trotter-Suzuki formula with $\widetilde{U}(t)=S_{2\chi}(t)$,
we have

$\hat{A}_{k}=(iH)^{k}/{k!}$ (59)

for all $k\leq 2\chi$. For $k>2\chi$, these operators $\hat{A}_{k}$ resemble
uncontrolled, erroneous operators. In Definition 1, we propose a superposition
of approximations with differently scaled times $S_{2\chi}(b_{n}t)$,
introducing controllable modulation parameters $\nu_{k}$ up to a Taylor order
of $2\chi R$. Its Taylor expansion is given by

$\displaystyle\mathcal{L}_{2\chi,R}(\boldsymbol{\nu},\boldsymbol{b},t)$
$\displaystyle=\sum_{q=0}^{2\chi
R}C_{q}(\boldsymbol{\nu},\boldsymbol{b})\,S_{2\chi}(b_{q}t)$ (60)
$\displaystyle=\sum_{k=0}^{\infty}\left(\sum_{q=0}^{2\chi
R}C^{\phantom{k}}_{q}\\!(\boldsymbol{\nu},\boldsymbol{b})\,b_{q}^{k}\right)\hat{A}_{k}t^{k}$
(61) $\displaystyle=\sum_{k=0}^{2\chi
R}\nu_{k}\,\hat{A}_{k}t^{k}\;+\;\mathcal{O}\left(t^{2\chi R+1}\right),$ (62)

where $C_{q},b_{q}\in\mathbb{R}$ and
$\nu_{k}=(\sum_{q}C^{\phantom{k}}_{q}b_{q}^{k})$ is fixed via the linear
transformation of the coefficients

$\boldsymbol{C}=(C_{0},C_{1},\dots,C_{2\chi R})^{\top}$

with the Vandermonde matrix $B_{j,k}=b_{k}^{j-1}$,
$B\boldsymbol{C}=\boldsymbol{\nu}$ as in Eq. (26). The condition number of the
Vandermonde matrix as the product of its Hilbert-Schmidt norm and the norm of
the pseudo-inverse can be bounded from above and below by explicit expressions
involving the vector defining the Vandermonde matrix [35]. Choosing the vector
$\boldsymbol{b}$, we can calculate the coefficients $\boldsymbol{C}$ for a
fixed solution vector $\boldsymbol{\nu}$ using the inverse Vandermonde matrix
$B^{-1}$, which is found exactly to be [36, 37]

$\displaystyle\left(B^{-1}\right)_{j,k}=\frac{(-1)^{k-1}}{\prod\limits_{m\,\in\,\boldsymbol{\mu}(j)}\left(b_{m}-b_{j}\right)}\sum_{\scriptsize\quad\begin{matrix}\boldsymbol{a}\in\mathbb{F}_{2}^{{}^{2\chi
R}}\\\ |\boldsymbol{a}|=2\chi R-k\end{matrix}}\prod_{i=0}^{2\chi
R-1}\left(b_{\mu_{i}(j)}\right)^{a_{i}}$ (63)
$\displaystyle\text{with}\quad\boldsymbol{\mu}(j)=(0,1,\ldots,j-1,j+1,\ldots,2\chi
R)\,,$

where the sum runs over all binary strings
$\boldsymbol{a}=(a_{0}\,a_{1}\,\dots\,a_{2\chi R-1})$ of length $2\chi R$ and
Hamming weight $2\chi R-k$. However, we have found that for numerical
purposes, a matrix inversion of $B$ clearly outmatches the analytical
computation of $B^{-1}$ in terms of runtime.

Using Eq. (62), we can now manipulate the Taylor expansion of a Trotter-Suzuki
block $S_{2\chi}(t)$ by replacing it with some
$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu},\boldsymbol{b},t)$. The key insight
here is that while the operators $\hat{A}_{k}$ are only correct for Taylor
orders $k\leq 2\chi$, we gain control of the prefactors up to order $2\chi R$.
Consequently, we can eliminate all erroneous $\hat{A}_{k}$ for $2\chi<k\leq
2\chi R$ by setting the corresponding $\nu_{k}$ to zero. This allows us to
construct the matching and closed-form multi-product formulas using blocks of
$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu},\boldsymbol{b},t)$ with different
$\boldsymbol{\nu}$ (and possibly $\boldsymbol{b}$).

####  IV.2.1 Matching multi-product formula

The first version of our proposed multi-product formula builds upon the
multiplication of $R$ multi-product building blocks
$\mathcal{L}_{2\chi,R}(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t)$ that
according to Eq. (62) can be written as

$\mathcal{L}_{2\chi,R}\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t\right)=\sum_{k=0}^{2\chi}\frac{\nu^{(r)}_{k}}{k!}\left(-\mathrm{i}Ht\right)^{k}+\sum_{k=2\chi+1}^{2\chi
R}\nu_{k}\hat{A}_{k}+\mathcal{O}\left(t^{2\chi R+1}\right)\,,$ (64)

where we can eliminate the second sum by setting $\nu_{k}=0$ for $2\chi<k\leq
2\chi R$. Their product now yields

$\prod_{r=1}^{R}\mathcal{L}_{2\chi,R}\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t\right)=\sum_{k=0}^{2\chi
R}\mu_{k}\left(-\mathrm{i}Ht\right)^{k}+\mathcal{O}(t^{2\chi R+1})$ (65)

with

$\mu_{k}\,=\sum_{i_{1}+i_{2}+\ldots+i_{R}=k}\frac{\nu_{i_{1}}^{(1)}\nu_{i_{2}}^{(2)}\cdots\nu_{i_{R}}^{(R)}}{i_{1}!i_{2}!\cdots
i_{R}!}\,.$ (66)

To mimic the exact time evolution up to $\mathcal{O}(t^{2\chi R+1})$, we
require $\mu_{k}=1/k!$, leading to Definition 2.

####  IV.2.2 Closed-form multi-product formula

The closed-form version of the presented multi-product formula sums products
of building blocks $\mathcal{L}$ . To motivate the specific construction, it
is useful to take a look at the specific case of $2\chi=4$ and $R=3$ for some
choice of $\boldsymbol{b}$ and $t$ using the shorthand

$\mathcal{L}\left(\begin{smallmatrix}\nu_{0}\\\ :\\\ \nu_{2\chi
R}\end{smallmatrix}\right):=\mathcal{L}(\boldsymbol{\nu},\boldsymbol{b},t),$
(67)

to get

$\widetilde{M}^{(\mathrm{cf})}_{4,3}=\mathcal{L}\left(\begin{matrix}1\\\
\boldsymbol{1}\\\ \boldsymbol{1}\\\ \boldsymbol{1}\\\ \boldsymbol{1}\\\ 0\\\
0\\\ \vdots\end{matrix}\right)+4!\,\mathcal{L}\left(\begin{matrix}0\\\
\boldsymbol{0}\\\ \boldsymbol{0}\\\ \boldsymbol{0}\\\ \boldsymbol{1}\\\ 0\\\
0\\\ \vdots\end{matrix}\right)\mathcal{L}\left(\begin{matrix}0\\\
\boldsymbol{1!/5!}\\\ \boldsymbol{2!/6!}\\\ \boldsymbol{3!/7!}\\\
\boldsymbol{4!/8!}\\\ 0\\\ 0\\\
\vdots\end{matrix}\right)+(4!)^{2}\mathcal{L}\left(\begin{matrix}0\\\
\boldsymbol{0}\\\ \boldsymbol{0}\\\ \boldsymbol{0}\\\ \boldsymbol{1}\\\ 0\\\
0\\\ \vdots\end{matrix}\right)\mathcal{L}\left(\begin{matrix}0\\\
\boldsymbol{0}\\\ \boldsymbol{0}\\\ \boldsymbol{0}\\\ \boldsymbol{1}\\\ 0\\\
0\\\ \vdots\end{matrix}\right)\mathcal{L}\left(\begin{matrix}0\\\
\boldsymbol{1!/9!}\\\ \boldsymbol{2!/10!}\\\ \boldsymbol{3!/11!}\\\
\boldsymbol{4!/12!}\\\ 0\\\ 0\\\ \vdots\end{matrix}\right)\,,$ (68)

where the orders $1,\dots,2\chi$ have been visually highlighted for clarity.
In this example the eleventh order in $t$ is revealed by multiplying the
corresponding terms from the Taylor expansion of $\mathcal{L}$ in Eq. (60),
with the corresponding coefficients in (68)

$(4!)^{2}\times\frac{(-\mathrm{i}Ht)^{4}}{4!}\times\frac{(-\mathrm{i}Ht)^{4}}{4!}\times
3!/11!\times\frac{(-\mathrm{i}Ht)^{3}}{3!}=\frac{(-\mathrm{i}Ht)^{11}}{11!}\,.$
(69)

The first term of Eq. (68) takes care of the zeroth and first $2\chi$ order of
$U$. The second term multiplies all terms of orders $t^{1}\dots t^{2\chi}$
with $t^{2\chi}$ and so takes care of the next $2\chi$ terms of the expansion.
The third term multiplies with a $\mathcal{O}(t^{2\chi})$ term twice, taking
care of the subsequent $2\chi$ terms – a pattern emerges. In a general setting
(with arbitrary $2\chi$ and $R$) we can write this sum as

$\displaystyle\widetilde{M}_{2\chi,R}^{(\mathrm{cf})}(t)$
$\displaystyle\coloneqq\sum_{r=1}^{R}\left(\mathcal{L}_{2\chi,R}\left(\boldsymbol{\nu}^{(0)},\boldsymbol{b}^{(0)},t\right)\right)^{r-1}\mathcal{L}_{2\chi,R}\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t\right)$
$\displaystyle=\sum_{r=1}^{R}\left(\sum_{j=0}^{2\chi
R}\nu_{j}^{(0)}\frac{\left(-\mathrm{i}Ht\right)^{j}}{j!}+\mathcal{O}\left(t^{2\chi
R+1}\right)\right)^{r-1}\left(\sum_{k=0}^{2\chi
R}\nu_{k}^{(r)}\frac{\left(-\mathrm{i}Ht\right)^{k}}{k!}+\mathcal{O}\left(t^{2\chi
R+1}\right)\right)\,$ (70)
$\displaystyle=\sum_{r=1}^{R}\left(\frac{(-\mathrm{i}Ht)^{2\chi}}{(2\chi)!}\right)^{r-1}\sum_{k=0}^{2\chi
R}\nu_{k}^{(r)}\frac{\left(-\mathrm{i}Ht\right)^{k}}{k!}\;+\;\mathcal{O}\left(t^{2\chi
R+1}\right)$ (71)

where we have used Eq. (62) in Eq. (70) and $\nu_{j}^{(0)}=\delta_{2\chi,j}$
from Definition 3 to collapse the first factor into
$(-\mathrm{i}Ht)^{2\chi(r-1)}/(2\chi)!$ in Eq. (71). Considering also that
$\nu_{j}^{(0)}=1$ for $0\leq j\leq 2\chi$, the $r=1$ term (that we recognize
are the first $2\chi$ orders of the time evolution) can be separated from the
sum. Discarding all terms that vanish due to $\nu_{k}^{(r)}=0$ we rewrite Eq.
(71) to

$\displaystyle\widetilde{M}_{2\chi,R}^{(\mathrm{cf})}(t)\;$
$\displaystyle=\;\sum_{j=0}^{2\chi}\frac{(-\mathrm{i}Ht)^{j}}{j!}\;+\;\sum_{r=2}^{R}\sum_{k=1}^{2\chi}\nu_{k}^{(r)}\frac{\left(-\mathrm{i}Ht\right)^{2\chi(r-1)+k}}{\left((2\chi)!\right)^{r-1}k!}\;+\;\mathcal{O}\left(t^{2\chi
R+1}\right)\,,$ (72)

for which we consult Definition 3, a last time resolving the remaining
$\nu_{k}^{(r)}$. This leaves us with the correct time evolution up to order
$2\chi R+1$, thus proving the definition.

###  IV.3 Error of the averaged operators

Following upon the insights from Section IV.2 we find Corollary 2 already
proven by Eq. (64) and Eq. (72) as long as $\widetilde{M}_{2\chi,R}$ is
constructed according to Definition 2 or 3. The proof of Theorem 3 requires a
more involved error analysis.

###### Proof of Theorem 3.

We first need to bound the remainder terms of the Taylor series expansions of
$U(t)$ and $\widetilde{M}^{(\mathrm{m})}$. In the following,
$\mathrm{\textbf{R}}_{\ell}(f)$ denotes the remainder term of the Taylor
series of an operator-valued function $f$ truncated at $\ell$-th order in $t$.
We thus find that

$\displaystyle\left\lVert\vphantom{\prod_{r}^{R}}U(t)-\widetilde{M}^{(\mathrm{m})}_{2\chi,R}\right\rVert$
$\displaystyle=\left\lVert\mathrm{e}^{-\mathrm{i}tH}-\prod_{r=1}^{R}\left[\sum_{q=0}^{2\chi
R}C_{q}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)}\\!\right)\,S_{2\chi}\left(b^{(r)}_{q}t\right)\right]\right\rVert$
$\displaystyle\leq\left\lVert\vphantom{\mathrm{\textbf{R}}_{2\chi
R}\left(\prod_{r=1}^{R}\left[\sum_{q=0}^{2\chi
R}C_{q}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)}\\!\right)\,S_{2\chi}\left(b^{(r)}_{q}t\right)\right]\right)}\mathrm{\textbf{R}}_{2\chi
R}\left(\mathrm{e}^{-\mathrm{i}tH}\right)\right\rVert+\left\lVert\mathrm{\textbf{R}}_{2\chi
R}\left(\prod_{r=1}^{R}\left[\sum_{q=0}^{2\chi
R}C_{q}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)}\\!\right)\,S_{2\chi}\left(b^{(r)}_{q}t\right)\right]\right)\right\rVert.$
(73)

Moving forward, we employ some recently established results on the theory of
Trotter errors. Specifically, we make use of the ‘Trotter error with 1-norm
scaling’ lemma of Ref. [11]. Building upon these insights in this fresh
context, we find that the exponential remainder of a product formula such as
in Eq. (6) can be bounded by

$\displaystyle\left\lVert\mathrm{\textbf{R}}_{\ell}\left(\prod_{j=1}^{N}\mathrm{e}^{-\mathrm{i}\alpha_{j}h_{k_{j}}t}\right)\right\rVert$
$\displaystyle=\frac{t^{\ell+1}}{(\ell+1)!}\left\lVert\left(\frac{\partial}{\partial
t}\right)^{\ell}\prod_{j=1}^{N}\mathrm{e}^{-\mathrm{i}\alpha_{j}h_{k_{j}}t}\right\rVert$
(74)
$\displaystyle=\frac{t^{\ell+1}}{(\ell+1)!}\left\lVert\sum_{x_{1}+\dots+x_{N}=\ell+1}\frac{(\ell+1)!}{x_{1}!\cdots
x_{N}!}\prod_{j=1}^{N}\left(\frac{\partial}{\partial
t}\right)^{x_{n}}\mathrm{e}^{-\mathrm{i}\alpha_{j}h_{k_{j}}t}\right\rVert$
$\displaystyle=\frac{t^{\ell+1}}{(\ell+1)!}\left\lVert\sum_{x_{1}+\dots+x_{N}=\ell+1}\frac{(\ell+1)!}{x_{1}!\cdots
x_{N}!}\prod_{j=1}^{N}\left(-\mathrm{i}\alpha_{j}\,h_{k_{j}}\,t\right)^{x_{j}}\mathrm{e}^{-\mathrm{i}\alpha_{j}h_{k_{j}}t}\right\rVert$
$\displaystyle\leq\frac{t^{\ell+1}}{(\ell+1)!}\sum_{x_{1}+\dots+x_{N}=\ell+1}\frac{(\ell+1)!}{x_{1}!\cdots
x_{N}!}\prod_{j=1}^{N}\left\lVert\alpha_{j}\,h_{k_{j}}\,t\right\rVert^{x_{j}}\cdot\underbrace{\left\lVert\mathrm{e}^{-\mathrm{i}\alpha_{j}h_{k_{j}}t}\right\rVert}_{=1}$
$\displaystyle=\frac{\left(\sum_{j=1}^{N}\left\lVert\alpha_{j}\,h_{k_{j}}\right\rVert
t\right)^{\ell+1}}{(\ell+1)!}\,.$

This bound can now be applied to (73). For the first term we obtain

$\left\lVert\mathrm{\textbf{R}}_{2\chi
R}\left(\mathrm{e}^{-\mathrm{i}tH}\right)\right\rVert\leq\frac{\left(\Lambda
t\right)^{2\chi R+1}}{(2\chi R+1)!},$ (75)

while the second term can be bounded by

$\displaystyle\left\lVert\mathrm{\textbf{R}}_{2\chi
R}\left(\prod_{r=1}^{R}\left[\sum_{q=0}^{2\chi
R}C_{q}^{(r)}S_{2\chi}\left(b^{(r)}_{q}t\right)\right]\right)\right\rVert$
$\displaystyle\leq\sum_{q_{1}=0}^{2\chi R}\cdots\sum_{q_{r}=0}^{2\chi
R}\left\lvert C_{q_{1}}^{(1)}\cdots
C_{q_{r}}^{(r)}\right\rvert\left\lVert\mathrm{\textbf{R}}_{2\chi
R}\left(S_{2\chi}\left(b^{(1)}_{q_{1}}t\right)\cdots
S_{2\chi}\left(b^{(r)}_{q_{r}}t\right)\right)\right\rVert$
$\displaystyle\leq\frac{\left(g_{\chi}\Lambda t\right)^{2\chi R+1}}{(2\chi
R+1)!}\sum_{q_{1}=0}^{2\chi R}\cdots\sum_{q_{r}=0}^{2\chi R}\left\lvert
C_{q_{1}}^{(1)}\cdots C_{q_{r}}^{(r)}\right\rvert\left(\lvert
b_{q_{1}}^{(1)}\rvert+\ldots+\lvert b_{q_{r}}^{(r)}\rvert\right),$ (76)

where we have introduced the factor $g_{\chi}$ from the Trotter-Suzuki
decomposition [13], which is defined as

$g_{\chi}\coloneqq\frac{4\chi}{5}\left(\frac{5}{3}\right)^{\chi-1}.$ (77)

At the same time, we have defined $\Lambda\coloneqq\sum_{k}\lVert
h_{k}\rVert$. Consequently, we can bound the error of the matching multi-
product formula via

$\left\lVert
U(t)-\widetilde{M}^{(\mathrm{m})}_{2\chi,R}(t)\right\rVert\leq\left(1+\zeta^{(\mathrm{m})}g_{\chi}^{2\chi
R+1}\right)\frac{\left(\Lambda t\right)^{2\chi R+1}}{\left(2\chi
R+1\right)!}\,,$ (78)

with

$\zeta^{(\mathrm{m})}\coloneqq\sum_{q_{1}=0}^{2\chi R}\sum_{q_{2}=0}^{2\chi
R}\cdots\sum_{q_{r}=0}^{2\chi R}\left\lvert
C_{q_{1}}^{(1)}C_{q_{2}}^{(2)}\cdots C_{q_{r}}^{(r)}\right\rvert\left(\lvert
b_{q_{1}}^{(1)}\rvert+\ldots+\lvert b_{q_{r}}^{(r)}\rvert\right),$ (79)

concluding the proof.

The error analysis of the closed-form multi-product formula follows along
similar lines: Again, we bound the remainder terms of the Taylor series
expansions of $U(t)$ and $\widetilde{M}^{(\mathrm{cf})}$, finding

$\displaystyle\left\lVert
U(t)-\widetilde{M}^{(\mathrm{cf})}_{2\chi,R}\right\rVert=\left\lVert\mathrm{e}^{-\mathrm{i}tH}-\sum_{r=1}^{R}\left(\mathcal{L}_{2\chi,R}\left(\boldsymbol{\nu}^{(0)},\boldsymbol{b}^{(0)},t\right)\right)^{r-1}\mathcal{L}_{2\chi,R}\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)},t\right)\right\rVert$
$\displaystyle\leq\left\lVert\vphantom{\mathrm{\textbf{R}}_{2\chi
R}\left(\sum_{r=1}^{R}\left[\sum_{q=0}^{2\chi
R}C_{q}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)}\\!\right)\,S_{2\chi}\left(b^{(r)}_{q}t\right)\right]^{r-1}\left[\sum_{q=0}^{2\chi
R}C_{q}^{(r)}\,S_{2\chi}\left(b_{q}t\right)\right]\right)}\mathrm{\textbf{R}}_{2\chi
R}\left(\mathrm{e}^{-\mathrm{i}tH}\right)\right\rVert$
$\displaystyle\quad+\;\left\lVert\mathrm{\textbf{R}}_{2\chi
R}\left(\sum_{r=1}^{R}\left[\sum_{q=0}^{2\chi
R}C_{q}\\!\left(\boldsymbol{\nu}^{(0)},\boldsymbol{b}^{(0)}\\!\right)\,S_{2\chi}\\!\left(b_{q}^{(0)}t\right)\right]^{r-1}\left[\sum_{q=0}^{2\chi
R}C_{q}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)}\\!\right)\,\,S_{2\chi}\\!\left(b_{q}^{(r)}t\right)\right]\right)\right\rVert\,,$
(80)

where the second term can be bounded by

$\displaystyle\left\lVert\mathrm{\textbf{R}}_{2\chi
R}\left(\sum_{r=1}^{R}\left[\sum_{q=0}^{2\chi
R}C_{q}\\!\left(\boldsymbol{\nu}^{(0)},\boldsymbol{b}^{(0)}\\!\right)\,S_{2\chi}\\!\left(b_{q}^{(0)}t\right)\right]^{r-1}\left[\sum_{q=0}^{2\chi
R}C_{q}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)}\\!\right)\,\,S_{2\chi}\\!\left(b_{q}^{(r)}t\right)\right]\right)\right\rVert$
$\displaystyle\;\leq\;\sum_{r=1}^{R}\sum_{q_{1}=0}^{2\chi
R}\cdots\sum_{q_{r}=0}^{2\chi R}\left\lvert
C_{q_{1}}(\boldsymbol{\nu}^{(0)},\boldsymbol{b}^{(0)})\cdots
C_{q_{r-1}}(\boldsymbol{\nu}^{(0)},\boldsymbol{b}^{(0)})\,C_{q_{r}}(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)})\right\rvert$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times\left\lVert\mathrm{\textbf{R}}_{2\chi
R}\left(S_{2\chi}\left(b_{q_{1}}^{(0)}\right)\cdots
S_{2\chi}\left(b_{q_{r-1}}^{(0)}\right)\,S_{2\chi}\left(b_{q_{r}}^{(r)}\right)\right)\right\rVert$
(81)

Defining

$\displaystyle\zeta^{(\mathrm{cf})}$ $\displaystyle\coloneqq$
$\displaystyle\sum_{r=1}^{R}\sum_{q_{1}=0}^{2\chi
R}\sum_{q_{2}=0}\cdots\sum_{q_{r}=0}^{2\chi R}\left\lvert
C_{q_{1}}^{(0)}C_{q_{2}}^{(0)}\cdots
C_{q_{r-1}}^{(0)}C_{q_{r}}^{(r)}\right\rvert\left(\lvert
b_{q_{1}}^{(0)}\rvert+\ldots+\lvert b_{q_{r-1}}^{(0)}\rvert+\lvert
b_{q_{r}}^{(r)}\rvert\right),$ (82) $\displaystyle\Lambda$
$\displaystyle\coloneqq$ $\displaystyle\sum_{k}\lVert h_{k}\rVert\,,$ (83)

we find

$\left\lVert
U(t)-\widetilde{M}^{(\mathrm{cf})}_{2\chi,R}(t)\right\rVert\leq\left(1+\zeta^{(\mathrm{cf})}g_{\chi}^{2\chi
R+1}\right)\frac{\left(\Lambda t\right)^{2\chi R+1}}{\left(2\chi
R+1\right)!}\,,$ (84)

concluding the proof. ∎

##  V Comparison and numerical validation

Algorithm |  Trotter-Suzuki [26] |  Childs & Wiebe [13] |  Novel multi-product formulas [here]
---|---|---|---
Formula | $S_{2\chi}(t/r)^{r}$ | $\sum_{q=1}^{K+1}C_{q}S_{2\chi}(t/q)^{q}$ | $\prod_{r=1}^{R}\sum_{q=0}^{2\chi R}C_{q}\\!\left(\boldsymbol{\nu}^{(r)},\boldsymbol{b}^{(r)}\right)\;S_{2\chi}\\!\left(b_{q}^{(r)}t\right)$
Max. depth | $r\cdot 2\cdot 5^{\chi-1}$ | $(K+1)\cdot 2\cdot 5^{\chi-1}$ | $R\cdot 2\cdot 5^{\chi-1}$
Error scaling | $\mathcal{O}((t/r)^{2\chi+1})$ | $\mathcal{O}(t^{2(\chi+K)+1})$ | $\mathcal{O}(t^{2\chi R+1})$
Error bound | $2\frac{\left(g_{\chi}\Lambda t/r\right)^{2\chi+1}}{(2\chi+1)!}$ |  | $\left(1+g_{\chi}^{2(\chi+K)+1}\sum\limits_{q}\left|C_{q}\right|\right)$
---
$\times\;\frac{\left(\Lambda t\right)^{2(\chi+K)+1}}{(2(\chi+K)+1)!}$
$\left(1+g_{\chi}^{2\chi R+1}\zeta\right)\frac{\left(\Lambda t\right)^{2\chi
R+1}}{\left(2\chi R+1\right)!}$
Table 1: Comparison of the standard Trotter-Suzuki formula, the multi-product
formula used by Childs and Wiebe and the closed-form / matching multi-product
formula introduced in this work. Here, $g_{\chi}=(5/3)^{\chi-1}4\chi/5$
originates in the exponential tail of the Trotter-Suzuki terms and
$\Lambda=\sum_{k}\left\lVert h_{k}\right\rVert$. The formula dependent
$\zeta^{(\mathrm{m})}$ and $\zeta^{(\mathrm{cf})}$ are defined and derived in
equations (79) and (82), respectively, and depend on the used $C_{q}$ and
$b_{q}$. The maximum circuit depth is given for their implementation using the
randomized sampling framework introduced in Algorithm 1. Figure 4: Comparing
error bounds for Trotter-Suzuki product formulas as well as several multi-
product formulas for different $\tau=\Lambda t$. The formulas for these bounds
are given in Table 1. We have chosen $2\chi=4$ and $r=R=K+1=3$ to fix the
depth of all methods to be 30 oracle calls. The parameters
$\\{\boldsymbol{b}^{(r)}\\}$ have been numerically optimized with an initial
guess of $\boldsymbol{b}=\set{1,-1,2,-2,\ldots,7}$. We find the resolution
factors for the three multi-product formulas to be
$\Xi^{(\mathrm{CW})}\\!\approx 3.13$, $\Xi^{(\mathrm{m})}\\!\approx 1.22$ and
$\Xi^{(\mathrm{cf})}\\!\approx 1.36$. In the black box optimization, we have
used objective functions $\propto(\Xi^{\mathrm{(m)}})^{20}$ and
$\propto(\Xi^{\mathrm{(cf)}})^{10}$ to ensure reasonable resolution factors.

We will now compare the Trotter-Suzuki product formula algorithm and multi-
product formulas in the randomized sampling framework. To make comparisons as
fair as possible, we will assume that each algorithm uses at most $R$-fold
sequences of $S_{2\chi}(\cdot)$ blocks. Product formulas may use repetitions
as in Eq. (7) achieving an error of $\mathcal{O}((t/R)^{2\chi+1})$, which is a
reliable way to approximate time evolutions for longer times $R\gg\tau>1$.
However, repetitions improve the accuracy of shorter time evolutions only
minimally. The situation is different if $R$ is an integer power of five,
which allows us to build the next higher Trotter-Suzuki order and approximate
the time evolution up to a leading order of $2(\chi+\log_{5}R)+1$ in $t$. This
means improving the leading power of $t$ comes at an exponential cost for the
circuit depth. The situation can be remedied by sampling from a multi-product
formula. Remarkably, the statement for sampling observable eigenvalues with
product formulas is very similar to Theorem 2.

This allows us to disregard the sampling error in the comparison, and only
think about asymptotic limits. The main difference between product and multi-
product formulas is that Trotter-Suzuki algorithms do not have a resolution
factor. This is equivalent to $\Xi=1$ with consequences for sampling
complexity and error bounds. This advantage is, however, quickly outweighed as
Childs and Wiebe’s product formula clearly delivers an improved approximation
that is exact up to a leading order of $2(\chi+R)-1$ in $t$. Modifying the
leading power of $t$ is now possible by adding exponentially fewer terms. Note
that they require the same circuit depth for $R=K+1$ in the randomized
sampling framework due to the final term in Eq. (12). With matching and
closed-form multi-product formulas, the approximation can be further improved
to leading order $2\chi R+1$ in $t$. This means that only one additional
$S_{2\chi}(\cdot)$ block improves the order by $2\chi$ rather than $2$ as for
multi-product formulas of the prior art. While the scaling in $t$ of the novel
multi-product formulas is superior to those of Childs and Wiebe, their
theoretical error bound is not as tight, leading to a crossover of error
bounds at $\sum|h_{k}|t<1$.

We have plotted error bounds of all formulas for fixed circuit depth between
all methods with respect to one particular optimization of the multi-product
formula parameters in Fig. 4. Optimizing the set of parameters
$\\{\boldsymbol{b}^{(r)}\\}$, we tend to achieve noticeably lower resolution
factors than obtained using Childs and Wiebe’s multi-product formula. We
generally find the matching multi-product formula to have a smaller resolution
factor than the closed-form multi-product formula. While the matching multi-
product formula has a better resolution, it requires an additional layer of
classical optimization.

It is important to note that since we did not find useful bounds relating
$\boldsymbol{b}$ to any of the resolution factors, further improvements of the
bounds in Fig. 4 seem possible. In the absence of relations $\Xi(b)$, it is
necessary to numerically optimize all $\boldsymbol{b}^{(r)}$ with respect to a
chosen loss function. We found that global optimization using basin hopping
combined with Nelder-Mead optimization yields the best results, since the
optimization landscape exhibits a large number of local minima. The
optimization is also sensitive to the initial guess for those parameters, with
an equal spread of positive and negative integers, i.e.,
$\boldsymbol{b}_{\text{init}}=(1,-1,2,-2,3,-3,\ldots,\chi R+1)$, leading to
better results. It is also fruitful to vary the loss function used for
optimization. While using the error bound for a fixed $\tau$ yields the lowest
error, the corresponding resolution factors are too large to be practical.
Amplifying the impact of the resolution factor by using $\Xi^{p}$ for some
power $p$ in the loss function yielded better results, though other useful
loss functions could be explored further. We have also found that bounding
each parameter $b$ by the total maximum to lead to more robust results.

To corroborate the functioning of the newly proposed multi-product formulas
beyond theoretical bounds, we also compare the actual performance with that of
the conventional multi-product formula and Trotter-Suzuki. Here, we compare
the actual operator distance to the ideal time evolution for the following
five physically plausible and meaningful Hamiltonians, with the results shown
in Fig. 5:

Figure 5: Numerical comparison of Trotter-Suzuki product formulas, Childs and
Wiebe’s multi-product formula and the multi-product formulas proposed in this
work for different $\tau=\sum_{k}||h_{k}||t$. Here, we approximate the time
evolution operator and plot the operator distance between the approximation
and the ideal evolution operator for a number of physically plausible and
interesting local, but not necessarily geometrically local, Hamiltonians,
defined in equations (87), (86), (88) and (89). The simulated distances
(thick, dashed lines) are much smaller than the theoretical bounds (thin,
solid lines), but do have remarkably similar features overall.

First, we consider a standard _Heisenberg Hamiltonian_ with periodic boundary
conditions, described by

$H_{\mathrm{Heisenberg}}=-\sum_{\langle
i,j\rangle}\left(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}\right)+2\sum_{i}X_{i}.$ (85)

This is a Hamiltonian that plays an important role in condensed matter
physics, as a prototypical model capturing ferromagnetism. Trotter-Suzuki and
Childs-Wiebe formulas perform extremely well, as the Hamiltonian comprises
many commuting terms. Since their performance guarantees rely on nested
commutators [11], this behavior is to be expected and demonstrates the
superior performance of these well studies formulas for lattice Hamiltonians,
which is not reached by the newly proposed multi-product formulas. However,
they are less optimal for Hamiltonians with fewer commuting terms. Motivated
by these findings, we now turn to investigating a Hamiltonian comprising
mutually anti-commuting terms, defined as

$H_{\mathrm{anti}}=\sum_{i=0}^{6}Z^{\otimes i}\otimes(X+Y)+Z^{\otimes 7}.$
(86)

Anti-commuting terms play an important role when using quantum simulation
algorithms to investigate _fermionic quantum models_. Here, we find a notable
and in fact significant advantage of the newly proposed multi-product formulas
over Trotter-Suzuki and Childs-Wiebe already for quite large parameters of
$\tau=\sum|h_{k}|t$. Consequently, we expect our formulas to work well and
exceed previous methods for fermionic systems, which once transformed into
qubit Hamiltonians by virtue of the Jordan-Wigner transformation, will have
fewer commuting terms than standard lattice spin Hamiltonians.

This behavior is confirmed and further corroborated by our results of the
_Sachdev–Ye–Kitaev (SYK)_ model as defined in Ref. [38], whose Hamiltonian is
given by

$H_{\mathrm{SYK}}=\frac{1}{4\cdot
4!}\sum_{p,q,r,s=0}^{N-1}J_{p,q,r,s}\gamma_{p}\gamma_{q}\gamma_{r}\gamma_{s},$
(87)

where $N$ is the number of Majorana fermion mode operators $\gamma_{p}$ and
the $J_{p,q,r,s}$ are real-valued scalars drawn randomly from a normal
distribution with variance $\sigma^{2}=3!/N^{3}$. This is an intricate local,
but not geometrically local, model that is believed to provide insights into
instances of strongly correlated quantum materials. It is used in the study of
scrambling dynamics and has a close relation with discrete models that capture
aspects of holography in the black hole context.

Again invoking the Jordan-Wigner transformation, $N$ Majorana fermion mode
operators can be mapped onto $N/2$ qubits. For our simulations, we thus chose
$N=10$ and $N=14$, leading to a five and seven qubit model, respectively.
Here, we again find a notable and in instances substantial advantage for large
$\tau$ and even for large times $t$. This is a physically highly plausible and
interesting model for which our new simulation methods fare well.

As a further fermionic system, we look at the spinful Hubbard model on two by
two sites, defined as [39]

$H_{\mathrm{Hub}}=-t\sum_{\langle
i,j\rangle,\sigma}(a^{\dagger}_{i,\sigma}a_{j,\sigma}+a^{\dagger}_{j,\sigma}a_{i,\sigma})+U\sum_{i}a^{\dagger}_{i,\uparrow}a_{i,\uparrow}a^{\dagger}_{i,\downarrow}a_{i,\downarrow}-\mu\sum_{i}\sum_{\sigma}a^{\dagger}_{i,\sigma}a_{i,\sigma}-h\sum_{i}(a^{\dagger}_{i,\uparrow}a_{i,\uparrow}-a^{\dagger}_{i,\downarrow}a_{i,\downarrow}),$
(88)

with spin $\sigma$, tunneling amplitude $t=2$, Coulomb potential $U=2$,
magnetic field $h=0.5$ and chemical potential $\mu=0.25$. Also, $a$ and
$a^{\dagger}$ represent fermionic annihilation and creation operators. These
are comparably small system sizes, but already show the substantial potential
of the proposed simulation method.

Finally, as a last family of examples, in order to gauge the performance of
our methods for larger system sizes, we investigate a system of 200 non-
interacting (“free”) fermions with nearest neighbor interactions and periodic
boundary conditions

$H_{\mathrm{ff}}(h)=\sum_{i,j}h_{i,j}a^{\dagger}_{i}a^{\dagger}_{j},$ (89)

with $h_{i,j}=1$ if the respective fermions are nearest neighbors and
$h_{i,j}=0$ otherwise. Note that for gauging the performance for free
fermions, we do not compare the time evolution operator
$U(t)=\mathrm{e}^{-\mathrm{i}Ht}$ to its approximation, but the Greens
function propagator $G(t)=\mathrm{e}^{-\mathrm{i}ht}$ to its approximation.

All comparisons are done in the randomized sampling framework for a Trotter-
Suzuki order of $2\chi=4$, the number of repetitions $R=3$ and corresponding
parameters for all other algorithms, ensuring an equal depth measured in the
number of the required oracle calls. Although the straightforward comparison
of their bounds in Fig. 4 suggests an advantage of our multi-product formulas
at about $\tau=0.1$, we find that this is the case for much larger $\tau$
already on actual systems as shown in Fig. 5. Note also that since the
Hamiltonians we consider are not necessarily geometrically local, known
classical simulation techniques will be heavily challenged even for comparably
short simulation times

While the performance of the simple multi-product formula and Trotter-Suzuki
is also much better than their bounds indicate, we find that the presented
bounds for our proposed multi-product formulas are comparably looser. It is
also worth to note that performance of the matching version is slightly better
than that of the closed-form version, although it comes with the penalty of an
additional, classical optimization loop. We also find that advantages of our
multi-product formulas are visible even for $\tau>1$; in the case of the SYK
model, this holds true even for actual simulation times $t>1$. Additionally,
we find that this advantage can also be maintained for larger system sizes, as
indicated by their performance on the model of free fermions. The presented
numerical studies therefore provide strong arguments for the functioning of
our proposed multi-product formulas and their advantage in the presented
regimes. It is also important to note that the corresponding resolution
factors required for the randomized sampling scheme of
$\Xi^{\mathrm{(cf)}}\approx 1.36$ and $\Xi^{\mathrm{(m)}}\approx 1.22$ are
significantly better than the $\Xi^{\mathrm{(CW)}}\approx 3.1$ of the multi-
product formula proposed by Childs and Wiebe in Ref. [13].

##  VI Discussion and conclusion

In this work, we have brought together two main ingredients of methods of
quantum simulation. The results are substantially more resource efficient ways
of performing short-time Hamiltonian simulation. These are on the one hand
higher-order multi-product formulas [13, 14], on the other an element which
has long been underappreciated but recently been of high interest: This is the
element of _randomness_ [20, 21, 22, 23, 24, 24]. Overcoming the prejudice
that the time evolution has to be completed in each run of a compilation, we
have introduced a novel framework for implementing multi-product formulas,
with the goal of estimating expectation values of time evolved observables.

Concretely, we have proposed a randomized sampling approach that focuses on
the time evolution of the observable, not the state. When implementing multi-
product formulas in a randomized fashion rather than via block encodings in
the LCU framework, we can circumvent the need for additional amplitude
amplification or post-selection. The results presented here have been obtained
by only requiring access to a quantum-oracle machine which implements single-
qubit state preparation, controlled time evolution and quantum measurements.
They are thus especially relevant for pre-digital settings, in which NISQ
algorithms reach their limits but where full-fledged, digital, long-time
evolution algorithms are not yet available. Consequently, this work may be
seen as targeting a regime in between the digital and analog setting, where we
have some form of parametric control over a simulator system allowing us to
compile the target time evolution with sequences of the simulator’s time
evolutions [40, 41]. This programmable regime then constitutes a departure
from the analog setting with relatively little control over the simulator, and
is not as strict as digital simulation, where the control over the quantum
system is strong enough to fashion its interactions into quantum gates.

Within this randomized sampling framework, we have proposed two new multi-
product formulas. These schemes have been equipped with full rigorous
performance guarantees. Furthermore, we have included a detailed estimation of
the number of circuit evaluations that are required, a vital metric for
randomized approaches. Comparing error bounds of these newly introduced
algorithms with Trotter-Suzuki product formulas and Childs and Wiebe’s multi-
product formula, we find that they outperform the latter for a fixed circuit
depth in a practically relevant regime. While their performance for long
simulation times would scale exponential in time, the base of this
exponential, the resolution factor in our case, while always strictly larger
than one, could be optimized further and even for our toy example was in the
range of $1.2-1.4$. Also, since we do not solely consider geometrically local
Hamiltonians, even comparably short simulation times are a highly difficult
task for known classical simulation techniques.

Benchmarking them on five different Hamiltonians, all of which physically well
motivated and each interesting in its own right, we have found that this
advantage can be expected already at comparably large simulation times. While
lattice Hamiltonians with many commuting terms are most likely best
approximated using Trotter-Suzuki or Childs-Wiebe formulas, the newly proposed
multi-product formulas show a clear advantage for fermionic Hamiltonians and
those with a small number of commuting term. This insight points to the
direction that there might not be a universally optimal quantum simulation
algorithm for digital quantum simulation. Instead, some algorithms could be
better suited to capture the specifics of a given local Hamiltonian model. The
downside of the proposed methods is that more measurements are required to
reach the desired precision through the resolution factor. This resolution
factor can be optimized using a classical black-box optimization, which is de-
facto a requirement for the functioning of the proposed multi-product
formulas.

The present work is essentially bridging the gap between analog and (perhaps
error-corrected), fully digital quantum technology: Not only do we expect
there to be other randomized sampling schemes in digital quantum simulation,
but once one is able to replace the element of randomness with block
encodings, one can switch from these expectation value based algorithm to
algorithms based on quantum phase estimation. Overall, the method introduced
gives rise to a less resource demanding way of performing Hamiltonian
simulation, while also remaining conceptually and technologically simpler than
for instance qubitization [18], bringing such ideas to an extent closer to the
realm of near-term quantum computing.

Looking ahead, it remains an open problem to relate the parameters
$\boldsymbol{b}$ in Definition 1 to the resolution factor in a way that would
eliminate the need for black-box optimization. So far, we can only connect the
two quantities analytically, and that involves the complicated product with
the Vandermonde matrix (26). Alternatively, one could improve the optimization
rather than replacing it. We have used only simple optimizers and loss
functions, and expect possible improvements for more involved loss functions
and optimization algorithms. Furthermore, the presented constructions are just
two of the plethora of new multi-product formulas that could be constructed
with Definition 1 and might exhibit better error bounds and resolution.

##  VII Acknowledgments

This work has been supported by the DFG (CRC 183 project B01 and A03, EI
519/21-1). This work has also received funding from the European Unions
Horizon 2020 research and innovation program under grant agreement No. 817482
(PASQuanS), specifically dedicated to programmable quantum simulators. It has
also been supported by the BMWi (PlanQK) and the BMBF (DAQC on notions of
digital-analog quantum simulation and FermiQP on fermionic quantum
processors).

M. K. acknowledges funding from ARC Centre of Excellence for Quantum
Computation and Communication Technology (CQC2T), project number CE170100012.
The authors endorse Scientific CO2nduct [42] and provide a CO2 emission table
in the appendix.

## References

  * Acin  _et al._ [2018] A. Acin, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, J. Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm, The European quantum technologies roadmap, New J. Phys. 20, 080201 (2018).
  * Lloyd [1996] S. Lloyd, Universal quantum simulators, Science 273, 1073 (1996).
  * Aharonov and Ta-Shma [2003] D. Aharonov and A. Ta-Shma, Adiabatic quantum state generation and statistical zero knowledge, arXiv:quant-ph/0301023 (2003).
  * Berry  _et al._ [2007] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Efficient quantum algorithms for simulating sparse Hamiltonians, Commun. Math. Phys. 270, 359 (2007).
  * Wiebe  _et al._ [2010] N. Wiebe, D. Berry, P. Høyer, and B. C. Sanders, Higher order decompositions of ordered operator exponentials, J. Phys. A 43, 065203 (2010).
  * Wiebe  _et al._ [2011] N. Wiebe, D. W. Berry, P. Høyer, and B. C. Sanders, Simulating quantum dynamics on a quantum computer, J. Phys. A 44, 445308 (2011).
  * Poulin  _et al._ [2011] D. Poulin, A. Qarry, R. Somma, and F. Verstraete, Quantum simulation of time-dependent Hamiltonians and the convenient illusion of Hilbert space, Phys. Rev. Lett. 106, 170501 (2011).
  * Kliesch  _et al._ [2011] M. Kliesch, T. Barthel, C. Gogolin, M. Kastoryano, and J. Eisert, Dissipative quantum Church-Turing theorem, Phys. Rev. Lett. 107, 120501 (2011).
  * Sweke  _et al._ [2016] R. Sweke, M. Sanz, I. Sinayskiy, F. Petruccione, and E. Solano, Digital quantum simulation of many-body non-markovian dynamics, Phys. Rev. A 94, 022317 (2016).
  * Childs  _et al._ [2018] A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Toward the first quantum simulation with quantum speedup, Proc. Natl. Ac. Sc. 115, 9456 (2018).
  * Childs  _et al._ [2020] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, A theory of Trotter error, arXiv:1912.08854 (2020).
  * Childs and Su [2019] A. M. Childs and Y. Su, Nearly optimal lattice simulation by product formulas, Phys. Rev. Lett. 123, 050503 (2019).
  * Childs and Wiebe [2012] A. M. Childs and N. Wiebe, Hamiltonian simulation using linear combinations of unitary operations, Quant. Inf. Comp. 12, 901 (2012).
  * Low  _et al._ [2019] G. H. Low, V. Kliuchnikov, and N. Wiebe, Well-conditioned multiproduct Hamiltonian simulation, arXiv:1907.11679 (2019).
  * Berry  _et al._ [2015a] D. W. Berry, A. M. Childs, and R. Kothari, Hamiltonian simulation with nearly optimal dependence on all parameters, in  _2015 IEEE 56th Annual Symposium on Foundations of Computer Science_ (IEEE, 2015) pp. 792–809.
  * Berry  _et al._ [2014] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Exponential improvement in precision for simulating sparse Hamiltonians, in  _Proceedings of the forty-sixth annual ACM symposium on Theory of computing_ (2014) pp. 283–292.
  * Berry  _et al._ [2015b] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Simulating Hamiltonian dynamics with a truncated Taylor series, Phys. Rev. Lett. 114, 090502 (2015b).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).
  * Endo  _et al._ [2021] S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, Hybrid quantum-classical algorithms and quantum error mitigation, J. Phys. Soc. Jap. 90 (2021).
  * Campbell [2017] E. T. Campbell, Shorter gate sequences for quantum computing by mixing unitaries, Phys. Rev. A 95, 042306 (2017).
  * Campbell [2019] E. T. Campbell, A random compiler for fast Hamiltonian simulation, Phys. Rev. Lett. 123, 070503 (2019).
  * Childs  _et al._ [2019] A. M. Childs, A. Ostrander, and Y. Su, Faster quantum simulation by randomization, Quantum 3, 182 (2019).
  * Ouyang  _et al._ [2020] Y. Ouyang, D. R. White, and E. T. Campbell, Compilation by stochastic Hamiltonian sparsification, Quantum 4, 235 (2020).
  * Chen  _et al._ [2020] C.-F. Chen, H.-Y. Huang, R. Kueng, and J. A. Tropp, Quantum simulation via randomized product formulas: Low gate complexity with accuracy guarantees, arXiv:2008.11751 [quant-ph] (2020).
  * Preskill [2018] J. Preskill, Quantum computing in the NISQ era and beyond, arXiv:1801.00862 [cond-mat, physics:quant-ph] 10.22331/q-2018-08-06-79 (2018).
  * Suzuki [1991] M. Suzuki, General theory of fractal path integrals with applications to many‐body theories and statistical physics, J. Math. Phys. 32, 400 (1991).
  * Blanes  _et al._ [1999] S. Blanes, F. Casas, and J. Ros, Extrapolation of symplectic integrators, Cel. Mech. Dyn. Astr. 75, 149 (1999).
  * Chin [2010] S. A. Chin, Multi-product splitting and Runge-Kutta-Nyström integrators, Celestial Mechanics and Dynamical Astronomy 106, 10.1007/s10569-010-9255-9 (2010).
  * Yoshida [1990] H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A 150, 262 (1990).
  * Hoeffding [1963] W. Hoeffding, Probability inequalities for sums of bounded random variables, J. Am. Stat. Ass. 58, 13 (1963).
  * Sheng [1989] Q. Sheng, Solving linear partial differential equations by exponential splitting, IMA J. Num. An. 9, 99 (1989).
  * Bespalova and Kyriienko [2020] T. A. Bespalova and O. Kyriienko, Hamiltonian operator approximation for energy measurement and ground state preparation, arXiv:2009.03351 (2020).
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nat. Phys. 16, 1050––1057 (2020).
  * Le Cam [1960] L. Le Cam, Locally asymptotically normal families of distributions. Certain approximations to families of distributions and their use in the theory of estimation and testing hypotheses, Univ. California Publ. Statist. 3, 37 (1960).
  * Bazan [2006] F. S. V. Bazan, Conditioning of rectangular Vandermonde matrices with nodes in the unit disk, SIAM J. Mat. An. App. 21, 679 (2006).
  * El-Mikkawy [2003] M. E. A. El-Mikkawy, Explicit inverse of a generalized Vandermonde matrix, Appl. Math. Comp. 146, 643 (2003).
  * Knuth [1997] D. E. Knuth,  _The art of computer programming: Volume 1: fundamental algorithms_ (Addison-Wesley, Boston, 1997).
  * Babbush  _et al._ [2019] R. Babbush, D. W. Berry, and H. Neven, Quantum simulation of the Sachdev-Ye-Kitaev model by asymmetric qubitization, Physical Review A 99, 040301 (2019).
  * McClean  _et al._ [2020] J. R. McClean, K. J. Sung, I. D. Kivlichan, Y. Cao, C. Dai, E. S. Fried, C. Gidney, B. Gimby, P. Gokhale, T. Häner, T. Hardikar, V. Havlíček, O. Higgott, C. Huang, J. Izaac, Z. Jiang, X. Liu, S. McArdle, M. Neeley, T. O’Brien, B. O’Gorman, I. Ozfidan, M. D. Radin, J. Romero, N. Rubin, N. P. D. Sawaya, K. Setia, S. Sim, D. S. Steiger, M. Steudtner, Q. Sun, W. Sun, D. Wang, F. Zhang, and R. Babbush, Openfermion: The electronic structure package for quantum computers, Quantum Science and Technology 5, 034014 (2020).
  * Trotzky  _et al._ [2012] S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert, and I. Bloch, Probing the relaxation towards equilibrium in an isolated strongly correlated one-dimensional Bose gas, Nature Physics 8, 325 (2012).
  * Parra-Rodriguez  _et al._ [2020] A. Parra-Rodriguez, P. Lougovski, L. Lamata, E. Solano, and M. Sanz, Digital-analog quantum computation, Phys. Rev. A 101, 022305 (2020).
  * con [2021] Scientific co2nduct, online (2021).

##  Appendix A CO2 emission table

Numerical simulations |
---|---
Total Kernel Hours [$\mathrm{h}$] | $\approx 800$
Thermal Design Power Per Kernel [$\mathrm{W}$] | 5.75
Total Energy Consumption Simulations [$\mathrm{kWh}$] | 4.6
Average Emission Of CO2 In Germany [$\mathrm{kg/kWh}$] | 0.56
Total CO2-Emission For Numerical Simulations [$\mathrm{kg}$] | 2.6
Were The Emissions Offset? | Yes
Transport |
Total CO2-Emission For Transport [$\mathrm{kg}$] | 0
Total CO2-Emission [$\mathrm{kg}$] | 2.6