# Well-conditioned multi-product formulas for hardware-friendly Hamiltonian
simulation

Almudena Carrera Vazquez  IBM Quantum, IBM Research Europe – Zurich  ETH
Zurich  Daniel J. Egger  IBM Quantum, IBM Research Europe – Zurich  David
Ochsner  IBM Quantum, IBM Research Europe – Zurich  ETH Zurich  Stefan Woerner
<EMAIL_ADDRESS>IBM Quantum, IBM Research Europe – Zurich

###### Abstract

Simulating the time-evolution of a Hamiltonian is one of the most promising
applications of quantum computers. Multi-Product Formulas (MPFs) are well
suited to replace standard product formulas since they scale better with
respect to time and approximation errors. Hamiltonian simulation with MPFs was
first proposed in a fully quantum setting using a linear combination of
unitaries. Here, we analyze and demonstrate a hybrid quantum-classical
approach to MPFs that classically combines expectation values evaluated with a
quantum computer. This has the same approximation bounds as the fully quantum
MPFs, but, in contrast, requires no additional qubits, no controlled
operations, and is not probabilistic. We show how to design MPFs that do not
amplify the hardware and sampling errors, and demonstrate their performance.
In particular, we illustrate the potential of our work by theoretically
analyzing the benefits when applied to a classically intractable spin-boson
model, and by computing the dynamics of the transverse field Ising model using
a classical simulator as well as quantum hardware. We observe an error
reduction of up to an order of magnitude when compared to a product formula
approach by suppressing hardware noise with Pauli Twirling, pulse efficient
transpilation, and a novel zero-noise extrapolation based on scaled cross-
resonance pulses. The MPF methodology reduces the circuit depth and may
therefore represent an important step towards quantum advantage for
Hamiltonian simulation on noisy hardware.

##  1 Introduction

A major goal of quantum mechanics is to accurately describe how a given
quantum system evolves with time. This was the original incentive argued by
Feynman in 1981 to develop quantum computers [Feynman(1982)]. More
specifically, Hamiltonian simulation seeks to implement the operator
$e^{-iHt}$ as accurately as possible with an efficient quantum algorithm,
which allows to study the time dynamics of a corresponding quantum system.
Here, $H$ is a Hermitian operator and efficient means that the running time
should be polynomial in the size of the corresponding system, i.e., ${\rm
polylog}(N)$ for a $N\times N$ Hamiltonian. For example, in the context of the
spin-boson model this would enable the study of decoherence models at scale
[Leggett et al.(1987)Leggett, Chakravarty, Dorsey, Fisher, Garg, and Zwerger,
Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli].

Hamiltonian simulation can be done with system specific machines in which the
hardware natively implements $H$, as exemplified by cold atoms [Gross and
Bloch(2017)], or with general purpose gate-based quantum computers. System
specific machines allow large-scale simulation but only for the subset of
Hamiltonians they natively support. Accurate Hamiltonian simulation on gate-
based quantum computers is however challenging since $e^{-iHt}$ has to be
implemented on hardware with a native gate set and limited connectivity
[Weidenfeller et al.(2022)Weidenfeller, Valor, Gacon, Tornow, Bello, Woerner,
and Egger].

Many quantum algorithms have been proposed to simulate Hamiltonians since
Feynman’s talk [Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli]. The
choice of the algorithm depends on the structure of the Hamiltonian [Childs
and Kothari(2010), Childs(2010)], e.g., whether it is sparse or dense, time-
dependent [Berry et al.(2020)Berry, Childs, Su, Wang, and Wiebe, An et
al.(2022)An, Fang, and Lin] or time-independent, and whether the time
evolution is real or imaginary [Endo et al.(2020)Endo, Sun, Li, Benjamin, and
Yuan, Kamakari et al.(2022)Kamakari, Sun, Motta, and Minnich].

Here, we focus on real-time simulations of a time-independent sparse
Hamiltonian. Then, the two main approaches are variational quantum algorithms
(VQA) [Li and Benjamin(2017), Cîrstoiu et al.(2020)Cîrstoiu, Holmes, Iosue,
Cincio, Coles, and Sornborger, Yuan et al.(2019)Yuan, Endo, Zhao, Li, and
Benjamin, Zoufal et al.(2021)Zoufal, Sutter, and Woerner] and product formulas
(PF) [Lloyd(1996), Aharonov and Ta-Shma(2003), Berry et al.(2007)Berry,
Ahokas, Cleve, and Sanders, Childs and Kothari(2011), Childs et
al.(2003)Childs, Cleve, Deotto, Farhi, Gutmann, and Spielman, Ahokas(2004)].
VQAs do not implement the operator $e^{-iHt}$ but instead create a
parameterized variational ansatz $\ket{\psi(\boldsymbol{\theta})}$ which
approximates the state $\ket{\Phi(t)}$ governed by the Schrödinger’s equation
[Li and Benjamin(2017)]. The quantum computer is then used to determine the
coefficients in an equation of motion (EOM) that is classically integrated.
The depth and structure of the circuit determine how well the variational
approach can approximate the true state evolution. A priori, this is not
clear, and we have to rely on error bounds that can be evaluated a posteriori
to improve the circuit design [Zoufal et al.(2021)Zoufal, Sutter, and
Woerner]. Further, VQAs must compute the coefficients in the EOM to high
accuracy. This requires a large number of circuit executions rendering the VQA
approach in its current implementation unlikely to lead to a quantum advantage
for time-dependent problems [Miessen et al.(2021)Miessen, Ollitrault, and
Tavernelli]. First steps towards overcoming this computational cost could
include approximating the Quantum Fisher Information matrix using simultaneous
perturbation stochastic approximation techniques [Gacon et al.(2021)Gacon,
Zoufal, Carleo, and Woerner]. In addition, the cost of estimating complex
observables can be reduced with positive operator valued measures [Fischer et
al.(2022)Fischer, Miller, Tacchino, Barkoutsos, Egger, and Tavernelli, García-
Pérez et al.(2021)García-Pérez, Rossi, Sokolov, Tacchino, Barkoutsos, Mazzola,
Tavernelli, and Maniscalco] or classical shadows [Huang et al.(2020)Huang,
Kueng, and Preskill, Hadfield et al.(2022)Hadfield, Bravyi, Raymond, and
Mezzacapo]. In contrast, PFs are analytical techniques which perform well in
practice and implement $e^{-iHt}$ as an operator by breaking the Hamiltonian
into sub-Hamiltonians, each simulated as a circuit for a small time-step. The
circuit depth grows polynomially with the inverse approximation error and the
simulation time. Interestingly, the existing error bounds for PFs can be
several orders of magnitude larger than measured errors, even for small
systems [Childs et al.(2018)Childs, Maslov, Nam, Ross, and Su, Hao Low et
al.(2019)Hao Low, Kliuchnikov, and Wiebe, Childs et al.(2021)Childs, Su, Tran,
Wiebe, and Zhu, Layden(2022)], which suggests that further investigation is
needed.

A promising candidate to replace PFs are multi-product formulas (MPF) [Blanes
et al.(1999)Blanes, Casas, and Ros, Childs and Wiebe(2012)], which can be seen
as a form of Richardson extrapolation [Richardson(1911), Butcher(2016),
Sidi(2003)]. Here, we show how to use MPFs to mitigate algorithmic errors in
PFs by _classically_ combining expectation values instead of implementing a
quantum Linear Combination of Unitaries (LCU) [Childs and Wiebe(2012)], which
was the original quantum algorithm proposed in 2012. We run a PF several times
on a gate-based quantum computer with a different number of Trotter steps to
estimate the same expectation value. Then, we linearly combine these estimates
on a classical computer into a better estimate of the expectation value, which
would otherwise only be obtainable with a large number of Trotter steps, i.e.,
significantly deeper circuits.

Our work makes both an algorithmic and a hardware contribution to Hamiltonian
simulation with MPFs. We now describe how we build on existing references. The
idea to classically combine expectation values to mitigate algorithmic errors
was first introduced in 2018 in Ref. [Endo et al.(2019)Endo, Zhao, Li,
Benjamin, and Yuan] applied to Hamiltonian simulation. Next, Ref. [Vazquez
Carrera et al.(2022)Vazquez Carrera, Hiptmair, and Woerner] showed that such a
classical combination reduces the resources needed to implement the HHL
algorithm [Harrow et al.(2009)Harrow, Hassidim, and Lloyd]. In particular,
Ref. [Endo et al.(2019)Endo, Zhao, Li, Benjamin, and Yuan] noticed that
Richardson extrapolation can mitigate both physical [Temme et al.(2017)Temme,
Bravyi, and Gambetta, Li and Benjamin(2017)] and algorithmic errors, and
showed promise in numerical experiments performed with two and three
extrapolation points. In 2019 Ref. [Hao Low et al.(2019)Hao Low, Kliuchnikov,
and Wiebe] noticed that a non-careful choice of sequences of Trotter steps
leads to exponentially small probabilities of success for MPFs implemented as
an LCU. Our work generalizes the classical combination which Ref. [Vazquez
Carrera et al.(2022)Vazquez Carrera, Hiptmair, and Woerner] applied to
tridiagonal Toeplitz symmetric matrices to arbitrary Hamiltonians. We further
improve upon Ref. [Vazquez Carrera et al.(2022)Vazquez Carrera, Hiptmair, and
Woerner], which focused on the HHL algorithm and thus a Hamiltonian simulation
embedded in the quantum phase estimation algorithm, by showing that the
classical combination is applicable to observables arising in a general
Hamiltonian simulation. We also solve the vanishing success probability issue
raised by Ref. [Hao Low et al.(2019)Hao Low, Kliuchnikov, and Wiebe]. We also
formalize the classical combination of expectation values for Hamiltonian
simulation of Ref. [Endo et al.(2019)Endo, Zhao, Li, Benjamin, and Yuan], our
contribution is an algorithm to create well-conditioned MPFs implemented as a
classical combination of expectation values. Here, we explain how ill-
conditioning, presented in Ref. [Hao Low et al.(2019)Hao Low, Kliuchnikov, and
Wiebe], increases non-algorithmic errors when classically combining
expectation values. The result of our work is a hardware-friendly MPF. Its
circuit depth is quadratically lower than MPFs implemented as an LCU [Hao Low
et al.(2019)Hao Low, Kliuchnikov, and Wiebe] since the individual product
formulas in the MPF are independently run on the hardware. The ill-
conditionning is avoided by a careful choice of Trotter exponents. We compare
the different algorithms in Table 1.

Crucially, we implement for the first time Hamiltonian simulation by MPF on
quantum hardware. We observe an algorithmic error reduction of up to an order
of magnitude while suppressing hardware noise [Kim et al.(2023)Kim, Wood,
Yoder, Merkel, Gambetta, Temme, and Kandala] with Pauli Twirling, pulse
efficient transpilation [Earnest et al.(2021)Earnest, Tornow, and Egger], and
a novel zero-noise extrapolation (ZNE) [Temme et al.(2017)Temme, Bravyi, and
Gambetta, Li and Benjamin(2017)] based on scaled cross-resonance pulses
[Stenger et al.(2021)Stenger, Bronn, Egger, and Pekker]. This last method is
similar to digital ZNE [Dumitrescu et al.(2018)Dumitrescu, McCaskey, Hagen,
Jansen, Morris, Papenbrock, Pooser, Dean, and Lougovski, He et al.(2020)He,
Nachman, de Jong, and Bauer, Giurgica-Tiron et al.(2020)Giurgica-Tiron, Hindy,
LaRose, Mari, and Zeng] but allows for a continuous range of stretch factors
without the need for costly pulse-calibration. We illustrate our work by
computing the dynamics in simulation and on quantum hardware of the transverse
field Ising model; a model ubiquitous in condensed matter and statistical
mechanics [Sachdev(2011)]. The Ising Hamiltonian is

$H=-J\sum_{\langle i,j\rangle}Z_{i}Z_{j}-h\sum_{i}X_{i},$ (1)

where $J$ denotes the interaction strength between nearest neighboring spins
$\langle i,j\rangle$, $h$ the transverse magnetic field, and $X$ and $Z$ are
Pauli operators. This model is often studied on both gate-based [Kim et
al.(2023)Kim, Wood, Yoder, Merkel, Gambetta, Temme, and Kandala, Chen et
al.(2022)Chen, Burdick, Yao, Orth, and Iadecola] and analogue quantum
simulations [Browaeys and Lahaye(2020), Bernien et al.(2017)Bernien, Schwartz,
Keesling, Levine, Omran, Pichler, Choi, Zibrov, Endres, Greiner, Vuletić, and
Lukin, Zhang et al.(2017)Zhang, Pagano, Hess, Kyprianidis, Becker, Kaplan,
Gorshkov, Gong, and Monroe]. Additionally, we extend the scaling predictions
made by Ref. [Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli] of the
resources needed to simulate a classically intractable spin-boson model by a
first-order PF and VQA. We show that in this setting our method requires
exponentially less resources with respect to the system size than a first-
order PF.

The main limitation of PF based approaches, including this work, is that the
number of gates depends linearly on the number of terms $L$ in the
Hamiltonian, which can grow very large for chemistry applications [Wecker et
al.(2014)Wecker, Bauer, Clark, Hastings, and Troyer]. Using randomization, the
qDRIFT protocol [Campbell(2019)] addresses this problem by trading accuracy
for an $L$-independent algorithm. Hybrid approaches that interpolate between
qDRIFT and PFs improve the asymptotic complexity with respect to both methods
[Ouyang et al.(2020)Ouyang, White, and Campbell, Hagan and Wiebe(2022)].
However, they do not improve the scaling with respect to the simulation time
nor the precision. Randomizing the implementation of the LCU has also been
explored. For example, Ref. [Faehrmann et al.(2022)Faehrmann, Steudtner,
Kueng, Kieferova, and Eisert] shows a drastic reduction in circuit depth by
implementing the LCU on average and increasing the number of circuits.
However, this method still requires two controlled PFs per circuit.
Alternatively, Ref. [Wan et al.(2022)Wan, Berta, and Campbell] presents a
randomized phase estimation algorithm based partially on the LCU approach and
has an $L$-independent complexity. Although only one ancilla is required and
most of the complexity is handled by the sampling, the circuits still require
$\mathcal{O}(\epsilon^{-2})$ controlled Pauli rotations, where $\epsilon$ is
the desired precision. The MPF approach we present here is suitable for
simulating general dynamics, in contrast to Ref. [Wan et al.(2022)Wan, Berta,
and Campbell], and for applications that require higher precision, since we
achieve more accurate results by increasing the number of shallow circuits we
run.

This paper is structured as follows. In Sec. 2 we review the theory of PFs and
then show how to implement well-conditioned MPFs by _classically_ combining
the outcome of multiple PFs. Further, in Sec. 3 we update the resource
estimate for spin-boson models in Ref. [Miessen et al.(2021)Miessen,
Ollitrault, and Tavernelli] with our MPF approach to illustrate that we can
achieve the same approximation accuracy as a first-order PF using shallower
circuits. In Sec. 4, we discuss error mitigation strategies for quantum
hardware and then illustrate our MPF method on fixed-frequency superconducing
transmon qubits by simulating an Ising model. We conclude our work in Sec. 5
and provide suggestions for further research.

Table 1:  Second-order PF compared to different MPFs. The last row corresponds to the algorithm presented in this paper. $l$ denotes the number of extrapolation points, $k_{1}<\ldots<k_{l}$ to the Trotter exponents, and $\varepsilon^{\prime}$ to the non-algorithmic errors such as physical or sampling noise. We assume that the MPFs are built from $S_{2}$, that $l\leq 15$, and that the simulation time satisfies $t\leq 1$. The harmonic Trotter sequence corresponds to $[k_{1},k_{1}+1,\ldots,k_{1}+l-1]$ and exemplifies a generic ill-conditioned sequence that may arise from a non-careful choice of Trotter exponents, such as proposed in Refs. [Childs and Wiebe(2012), Endo et al.(2019)Endo, Zhao, Li, Benjamin, and Yuan], where no particular sequences were specified. The constant factor of the circuit depth scaling inside the $\order{\cdot}$-notation for LCU implementations can be very large once the circuit is transpiled to the hardware due to the limited connectivity and the many controlled gates.  Algo. | Combination | Trotter seq. | #circuits | Max. circ. depth | Success prob. | Approx. error | $\|\vec{a}\|_{1}$
---|---|---|---|---|---|---|---
$S_{2}$ | N/A | $[k]$ | 1 | $\order{k}$ | 1 | $\order{\frac{1}{k^{2}}+\varepsilon^{\prime}}$ | N/A
MPF |  quantum (LCU) |  harmonic [Childs and Wiebe(2012)] | $1$ | $\order{l^{2}}$ | $\order{\|\vec{a}\|_{1}^{-2}}$ | $\order{\frac{1}{k_{1}^{2l+1}}+\varepsilon^{\prime}}$ | $e^{\Omega(l)}$
well-conditioned [Hao Low et al.(2019)Hao Low, Kliuchnikov, and Wiebe] | $\order{l^{2}\log(l)}$ | $\leq 2$
classical |  harmonic [Endo et al.(2019)Endo, Zhao, Li, Benjamin, and Yuan] | $l$ | $\order{l}$ | 1 | $\order{\frac{1}{k_{1}^{2l+1}}+\|\vec{a}\|_{1}\varepsilon^{\prime}}$ | $e^{\Omega(l)}$
well-conditioned (our work) | $\order{l\log(l)}$ | $\leq 2$

##  2  Multi-product formulas

Product formulas are used to build a unitary matrix approximating $e^{-iHt}$
by (i) writing $H$ as $\sum_{j}H_{j}$ (e.g. where each $H_{j}$ is a $1$-sparse
matrix [Childs et al.(2003)Childs, Cleve, Deotto, Farhi, Gutmann, and
Spielman, Ahokas(2004)]) such that (ii) each $e^{-iH_{j}t}$ has an efficient
quantum circuit representation, and finally (iii), if the matrices $H_{j}$ do
not commute with each other, use an approximation. The first-order PF—also
known as the Lie-Trotter formula—is

$S_{1}(t)\coloneqq\prod_{j=1}^{J}e^{-iH_{j}t}$ (2)

and has a quadratic error term in $t$, i.e.
$S_{1}(t)=e^{-iHt}+\mathcal{O}(t^{2})$. Higher order PFs converge faster and
are defined recursively as follows.

######  Definition 2.1 (Lie-Trotter-Suzuki, [Suzuki(1991)]).

Let $H\in\mathbb{C}^{N\times N}$ be a Hermitian matrix with
$H=\sum_{j=1}^{J}H_{j}$, and $\chi\in\mathbb{N}$. Then, the $2\chi$-order
symmetric formula approximating $e^{-iHt}$ is recursively defined by

$\displaystyle S_{2}(t)$ $\displaystyle\coloneqq$
$\displaystyle\prod_{j=1}^{J}e^{-iH_{j}t/2}\prod_{j=J}^{1}e^{-iH_{j}t/2},$
$\displaystyle S_{2\chi}(t)$ $\displaystyle\coloneqq$ $\displaystyle
S_{2\chi-2}\left(s_{\chi}t\right)^{2}S_{2\chi-2}\left[\left(1-4s_{\chi}\right)t\right]S_{2\chi-2}\left(s_{\chi}t\right)^{2},$

where $s_{p}=\left(4-4^{1/(2p-1)}\right)^{-1}$ for $p\in\mathbb{N}$.

When the duration $t$ of the time evolution is long we approximate $e^{-iHt}$
by splitting $t$ into $k$ equal _Trotter steps_ of length $t/k$ and apply a
product formula $S_{\chi}$ in each interval such that

$S_{\chi}^{k}\left({\textstyle\frac{t}{k}}\right)=\left[S_{\chi}\left({\textstyle\frac{t}{k}}\right)\right]^{k}=e^{-iHt}+\order{t\left({\textstyle\frac{t}{k}}\right)^{\chi}}.$
(4)

The last equality follows from the triangle inequality. As the order $\chi$ of
the PF increases the number of Trotter steps $k$ required to reach a target
error decreases. However, the depth of the quantum circuit scales
exponentially with $\chi$ [Suzuki(1991), Childs and Wiebe(2012)]. Therefore,
choosing the optimal formula for a given problem is a non-trivial task further
complicated by the lack of tight analytical bounds [Childs et al.(2018)Childs,
Maslov, Nam, Ross, and Su]. Moreover, even if one could find the PF that
minimizes the resources needed for a given Hamiltonian, it would still require
deep circuits, rendering the PF approach very challenging to implement on
noisy quantum hardware.

Multi-product formulas, which can be seen as a form of Richardson
extrapolation, are a promising candidate to replace product formulas [Blanes
et al.(1999)Blanes, Casas, and Ros, Childs and Wiebe(2012)]. MPFs only require
a sum of polynomially many unitary operations to achieve the same accuracy as
PFs with exponentially many unitaries would. They achieve this by linearly
combining PFs of the same order $\chi$ but varying Trotter exponents $k$ to
cancel lower order error terms.

######  Definition 2.2 (Multi-product formula).

For $l\in\mathbb{N}$, we define the Trotter exponents as
$\vec{k}=(k_{1},\dots,k_{l})$ where the $k_{i}$’s are distinct natural
numbers, and the extrapolation weights as $\vec{a}=(a_{1},\dots,a_{l})$ with
$a_{i}\in\mathbb{R}$. Then, the multi-product formula $M_{l,\chi}(t)$ is

$M_{l,\chi}(t)=\sum_{j=1}^{l}a_{j}S^{k_{j}}_{\chi}\left({\textstyle\frac{t}{k_{j}}}\right),$
(5)

where $\chi,S_{\chi}$ are as in Def. 2.1. To cancel the first $l-1$ error
terms, the extrapolation weights must satisfy

$\displaystyle\sum_{j=1}^{l}a_{j}=1,\quad\text{and}\quad\sum_{j=1}^{l-1}a_{j}/k_{j}^{\eta}=0,$
(6)

where $\eta=\chi+2n$ for symmetric product formulas and $\eta=\chi+n$
otherwise, for $n=0,...,l-2$. For more details, see Appendix A.

###  2.1 Circuit implementation

Ref. [Childs and Wiebe(2012)] introduced MPFs as a quantum algorithm to
approximate Hamiltonian simulation. The linear combination in Eq. (5) is
translated to a quantum circuit by the Linear Combination of Unitaries (LCU)
method, illustrated in Fig. 1. Asymptotically, this method with $k_{j}=j$,
i.e. the minimal sequence with $l$ different natural numbers, implements
$\order{l^{2}}$ Trotter steps and approximates the Hamiltonian better than PFs
with as many Trotter steps. More precisely, MPFs require exponentially fewer
Trotter steps [Chin(2010)] to achieve the same accuracy. However, the LCU
implementation is probabilistic with a success probability that decreases
exponentially as $e^{-\Omega(l)}$. Furthermore, even the shallowest circuits
are too deep to implement on noisy hardware since the LCU method sums Trotter
steps of varying order as controlled unitaries, requiring additional qubits
and controlled operations. For example, simulating a five linearly connected
spin Ising model with an MPF with $\vec{k}=[k_{1},k_{2},k_{3}]$ implemented by
an LCU requires $108(k_{1}+k_{2})+40k_{3}+4$ CNOT gates, see Appendix B.

While Eq. (5) is a linear combination of operators, in practice we measure
operator expectation values after a quantum circuit. Crucially, as we show
here and illustrate in Fig. 1, the linear combination in MPFs can thus be done
classically. This allows to parallelize the circuits implementing the PF for
different Trotter steps, resulting in shallow circuits with $\order{l}$
Trotter steps, which is better for noisy hardware. We therefore implement each
$S_{\chi}^{k_{j}}$ in Eq. (5) as an independent quantum circuit. As shown in
Appendix C we can substitute the $S_{\chi}$’s in Eq. (5) by their expectation
values. More precisely, computing an expectation value $E$ after approximately
simulating $e^{-iHt}$ with the PF $S^{k_{j}}_{2\chi}(t/k_{j})$ yields an
approximation $E_{j}$ of the expectation value $E$ for each Trotter exponent
$k_{j}$ with an approximation error related to $\chi$ and $k_{j}$ by Eq. (4).
Finally, combining the $E_{j}$’s as

$\sum_{j=1}^{l}a_{j}E_{j}=E+\mathcal{O}\left(\frac{t^{2(\chi+l)+1}}{k_{1}^{2(\chi+l)}}\right),$
(7)

gives an estimation of $E$ where the remaining asymptotic error depends on the
base PF, see Appendix A.

Figure 1: Implementation of MPFs. The LCU circuit probabilistically implements
the linear combination $\sum_{j}S_{1}^{k_{j}}(t/k_{j})$ and is successful if
the measurement outcome of the top two qubits is $0$. Our method implements
each of the PFs as a separate quantum circuit and calculates the linear
combination classically. Running the PFs separately instead of in sequence
circumvents the need of auxiliary qubits and significantly reduces the circuit
depth, especially as the number of terms in the linear combination increases.

###  2.2  Choosing Trotter exponents

MPFs can achieve exponentially lower algorithmic errors than PFs while
requiring similar circuit depths. However, a poor choice of the Trotter
exponents and extrapolation weights amplifies other error sources such as
physical or sampling noise so that the error in the final quantity is larger
than the largest algorithmic error in the MPF. This choice-of-sequences issue
was raised in Ref. [Hao Low et al.(2019)Hao Low, Kliuchnikov, and Wiebe] where
the authors noted that choosing the arithmetic progression for the Trotter
exponents $k_{j}=j$, which gives the shortest circuits, leads to an ill-
conditioned MPF, i.e. $\norm{\vec{a}}_{1}=e^{\Omega(l)}$ where
$\norm{\vec{a}}_{1}\coloneqq\sum_{j}^{l}\absolutevalue{a_{j}}$. When
implemented as a quantum circuit as an LCU this ill-conditioning translates
into an exponentially small success probability scaling as
$\mathcal{O}\big{(}\norm{\vec{a}}_{1}^{-2}\big{)}$. Since our approach is not
probabilistic, i.e., we classically compute the linear combination in Eq. (7)
rather than implementing $M_{l,\chi}$ as an LCU, the condition number
$\norm{\vec{a}}_{1}$ amplifies any error that is not the Trotter error by a
factor of $\norm{\vec{a}}_{1}$. More precisely, suppose that in Eq. (7) we
calculated the expectation values $E_{j}$ with a quantum computer and we have
several error sources. If $\varepsilon(k_{j})$ denotes the algorithmic error
of the Trotter exponent $k_{j}$ and other sources arising from, e.g., sampling
or physical noise are labelled $\varepsilon^{\prime}$, then

$E_{j}=E+\varepsilon(k_{j})+\varepsilon^{\prime}.$ (8)

Now, the result of the MPF is

$\displaystyle\sum_{j=1}^{l}a_{j}\left(E+\varepsilon(k_{j})+\varepsilon^{\prime}\right)$
(9)
$\displaystyle=E+\mathcal{O}\left(\frac{t^{2(l+\chi)+1}}{k_{1}^{2(l+\chi)}}\right)+\mathcal{O}\left(\norm{\vec{a}}_{1}\varepsilon^{\prime}\right).$
(10)

Eq. (10) shows that if $\norm{\vec{a}}_{1}$ is large, then
$\varepsilon^{\prime}$ is amplified and as a result, the sum in Eq. (9) can be
a worse approximation to $E$ than the $E_{j}$’s. Appendix D shows a toy
example illustrating the asymptotic behaviour, where $\varepsilon^{\prime}$ is
the sampling noise of a Bernoulli variable. Here, the increase in non-
algorithmic errors results in an increase in the sampling overhead.

Ref. [Hao Low et al.(2019)Hao Low, Kliuchnikov, and Wiebe] provides choices of
$\vec{k}$ for which $\norm{\vec{a}}_{1}$ remains “small” to avoid
exponentially low success probabilities, i.e., ill-conditioned MPFs. The
authors minimize either $\norm{\vec{a}}_{1}\lVert\vec{k}\rVert_{1}$ or
$\norm{\vec{a}}_{1}$ and provide the corresponding values for $\chi=2$ and
$4$. Their method requires $\order{l^{2}\log l}$ Trotter steps to implement
the LCU instead of $\order{l^{2}}$ for $k_{j}=j$. By contrast, our classical
combination of expectation values needs only $\order{l\log l}$ Trotter steps
for the longest circuit. Importantly, the condition number decreases
exponentially from $\norm{\vec{a}}_{1}=e^{\Omega(l)}$ to
$\norm{\vec{a}}_{1}=\mathcal{O}(\log l)$. Furthermore, Tab. I in Ref. [Hao Low
et al.(2019)Hao Low, Kliuchnikov, and Wiebe] shows that for an MPF based on
$S_{2}$ and a practical value of $l\leq 15$, the condition number
$\norm{\vec{a}}_{1}<1.7$ and is therefore negligible. To avoid increasing non-
algorithmic errors, see Eq. (10), we chose the Trotter exponents by minimizing
$\norm{\vec{a}}_{1}$ under the constraints in Eq. (6) for a given possible
range of the $k_{j}$’s. The algorithm is described in Appendix E.

We illustrate how to choose the Trotter exponents with the Ising model for
five linearly connected spins, see Eq. (1), with
$H_{1}=-J\sum_{i=0}^{3}Z_{i}Z_{i+1}$ and $H_{2}=-h\sum_{i=0}^{4}X_{i}$.
Therefore, $e^{-iH_{1}t}$ and $e^{-iH_{2}t}$ are implemented with layers of
two-qubit $R_{ZZ}$ and single-qubit $R_{X}$ rotations, respectively. As
benchmark we compute the exact local magnetization $\langle Z_{0}\rangle^{*}$
for $h=1$, $J=0.5$ and $t=0.5$ by numerically exponentiating $H_{1}+H_{2}$.
Next, we compute the local magnetization $\langle Z_{0}(k_{j})\rangle$ with a
first-order PF and $k_{j}$ Trotter steps. We calculate the relative error
$\absolutevalue{\langle Z_{0}(k_{j})\rangle-\langle Z_{0}\rangle^{*}}/|\langle
Z_{0}\rangle^{*}|$ which we fit to $k_{j}^{-1}$, see the blue dots and dashed
line in Fig. 2. The condition number of an MPF based on $S_{1}$ grows faster
than for $S_{2}$, see Appendix D. Therefore, we empirically set the threshold
defining a well-conditioned MPF to $\norm{\vec{a}}_{1}\leq 3$ instead of 1.7.
Thus, we compute well- and ill-conditioned MPFs, with $\norm{\vec{a}}_{1}\leq
3$ and $\norm{\vec{a}}_{1}>3$, respectively, with different Trotter exponents
by classically combining the expectation values $\langle Z_{0}(k_{j})\rangle$.
The values of $\vec{a}$ are given in Appendix E. In the absence of non-
algorithmic errors, i.e., $\varepsilon^{\prime}=0$, all MPFs reduce the error
compared to the benchmark, see triangles in Fig. 2. Next, we artificially add
a fixed perturbation $\varepsilon^{\prime}=10^{-3}$ following $\langle
Z_{0}(k_{j})\rangle\to\langle Z_{0}(k_{j})\rangle+{\rm
sign}(a_{j})\varepsilon^{\prime}$ before performing the linear combination.
The well-conditioned MPFs still reduce the error compared to ill-conditioned
sequences which are worse than the benchmark, see stars in Fig. 2. For
example, the noisy $S_{1}$-based MPF $[2,4]$ whose deepest circuit requires
$k_{2}=4$ Trotter steps has a lower error than the PF with $k=8$ Totter steps.
This is a factor of two reduction in circuit depth.

Figure 2: Local magnetization of the five-spin Ising model. The Ising
simulation is done with a first-order Trotter formula. The y-axis shows the
relative error $\absolutevalue{\langle Z_{0}\rangle-\langle
Z_{0}\rangle^{*}}/\absolutevalue{\langle Z_{0}\rangle^{*}}$ between the
simulated and exact values. The extrapolated values are computed using well-
conditioned (green markers) and ill-conditioned (red markers) MPFs with
different Trotter exponents (shown in the brackets $[k_{1},k_{2},...]$) by
combining the expectation values classically. The calculation is done without
non-algorithmic errors ($\varepsilon^{\prime}=0$, triangular markers) and with
non-algorithmic errors ($\varepsilon^{\prime}=10^{-3}$, star markers).

##  3  Spin-boson model

The spin-boson model [Leggett et al.(1987)Leggett, Chakravarty, Dorsey,
Fisher, Garg, and Zwerger] is a model of quantum dissipation. It gives
insights into the transition between coherent and incoherent behaviour as well
as a phase transition suppressing quantum tunneling. Its dynamics can be
simulated with analogue simulators [Egger and Wilhelm(2013), Indrajeet et
al.(2020)Indrajeet, Wang, Hutchings, Taketani, Wilhelm, LaHaye, and Plourde]
or gate-based quantum computers [Miessen et al.(2021)Miessen, Ollitrault, and
Tavernelli]. In Ref. [Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli]
the authors compare the performance of VQAs based on McLachlan’s variational
principle [McLachlan(1964)] against a first-order PF to simulate the dynamics
of the spin-boson model. The authors made scaling predictions for a
classically intractable system for both, VQA and first-order PFs. They
conclude that VQA in its current form is unlikely to lead to a quantum
advantage due to the large number of circuit executions, although it is
assumed to have shallower circuits than first-order PFs. Further, they show
that first-order PFs that achieve the desired accuracy lead to very deep
circuits, which might render them infeasible as well. We now extend their
results by adding second-order PFs and MPFs to the comparison, showing a
significant circuit depth reduction.

More precisely, the spin-boson Hamiltonian

$H=\sum^{M}_{k=1}\omega_{k}a_{k}^{\dagger}a_{k}+\frac{\omega_{s}}{2}Z+\Delta
X+\sum_{k=1}^{M}g_{k}X\left(a^{\dagger}_{k}+a_{k}\right)$ (11)

describes a two-level system, such as an atom, coupled to a bath of $M$ bosons
[Frisk Kockum et al.(2019)Frisk Kockum, Miranowicz, De Liberato, Savasta, and
Nori, Di Paolo et al.(2020)Di Paolo, Barkoutsos, Tavernelli, and Blais]. Here,
the bosonic operators $a^{\dagger}_{k}$ $\left(a_{k}\right)$ create
(annihilate) harmonic basis states with eigenfrequencies $\omega_{k}$. The
Pauli matrices $X$, $Z$ act on the state of the spin with eigenfrequency
$\omega_{s}$ and tunneling rate $\Delta$. The coupling strength between the
spin and the $k^{\text{th}}$ bosonic mode is $g_{k}$. As in Ref. [Miessen et
al.(2021)Miessen, Ollitrault, and Tavernelli] we consider the resonant case
$\omega_{k}\equiv\omega$, $g_{k}\equiv g$ and ultrastrong coupling, i.e.
$g/\omega=0.5$ with $\omega_{s}=-1$ and $\Delta=0$. We consider different
spin-boson models in which the $M$ bosonic modes are truncated at different
occupation numbers $n^{\max}$.

As in Ref. [Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli], we are
interested in $t=10$. Since we suspected that the algorithm would also work
for sequences with $t/k_{1}>1$, we first computed several sequences of Trotter
exponents with $l=2,3$ and $4$. These sequences and their corresponding
extrapolation weights satisfied $\norm{\vec{a}}_{1}\leq 2$ to ensure that
errors arising from noise or sampling are amplified by at most a factor of
$2$. Next, we computed the approximation error $\varepsilon$ for each MPF and
selected the one which required less resources to achieve the target accuracy
$\varepsilon_{t}$. We observed that this method works well in practice even
when $t/k_{1}>1$, confirming that the analytical error bounds on PFs are not
tight. Furthermore, we found that the second-order PF always requires
shallower circuits than the PFs $S_{4}$ and $S_{6}$ when used in an MPF. We
also found that the PF $S_{2}$ alone always requires shallower circuits than
$S_{1}$. For small spin-boson systems we see a reduction of up to a factor of
eight in circuit depth, see Appendix F.

We derive the scaling of the MPF approach taking into account the dependence
of the error on the number of qubits $N_{q}$ and the number of different
Trotter exponents $l$. First, the circuit depth of the spin-boson model grows
linearly with $N_{q}$ [Miessen et al.(2021)Miessen, Ollitrault, and
Tavernelli]. Next, as a worst case estimate we set $\alpha=t$ such that
$t/k_{j}\leq 1$ for $j=1,\cdots,l$ so that we can bound the error in the MPF
by $\mathcal{O}\left(t^{2l+1}/(2l+1)!\right)$ [Hao Low et al.(2019)Hao Low,
Kliuchnikov, and Wiebe, Vazquez Carrera et al.(2022)Vazquez Carrera, Hiptmair,
and Woerner]. A short calculation detailed in Appendix F shows that our hybrid
scheme requires the deepest circuit to have

$\displaystyle\mathcal{O}\left(\frac{t\log^{2}(N_{q}/\varepsilon_{t})}{\log^{2}\log(N_{q}/\varepsilon_{t})}\right)$
(12)

Trotter steps instead of the
$\mathcal{O}\big{(}\sqrt{\smash[b]{t^{3}N_{q}^{2}/\varepsilon_{t}}}\big{)}$
required by $S_{2}$. A more careful choice of Trotter exponents and $\alpha$
may allow better error bounds. Although the sampling overhead of
$\varepsilon^{-2}_{t}$, not shown in Eq. (12), cannot be avoided in either
approach, the depth of the MPF circuits grows logarithmically instead of
polynomially in $\varepsilon^{-1}_{t}$. Our method is advantageous for large
systems as it reduces the linear dependency of errors on $N_{q}$ to a
logarithmic one. Ref. [Miessen et al.(2021)Miessen, Ollitrault, and
Tavernelli] concluded that first-order Trotter formulas are a better approach
than VQA to simulate spin-boson models of intractable system sizes. However,
$S_{1}$ exhibits a scaling linear in $N_{q}$, making such a simulation
challenging to implement on noisy hardware. MPFs implemented as a linear
combination of expectation values are therefore promising to simulate
classically intractable quantum systems.

##  4  Ising model computation on hardware

Quantum computing hardware has increased in scale and quality [Jurcevic et
al.(2021)Jurcevic, Javadi-Abhari, Bishop, Lauer, Bogorin, Brink, Capelluto,
Günlük, Itoko, Kanazawa, and _et al._]. However, gate errors on cross-
resonance based hardware are non-negligible. For example, a typical CNOT gate
on IBM Quantum hardware has an average error ranging from 0.5% to 1%. This
makes it challenging to implement MPFs on noisy quantum hardware as hardware-
related errors $\varepsilon^{\prime}$ must be kept comparable to or smaller
than the total Trotter error which is close to 10%, see Eq. (10) and Fig. 2.
We can however suppress hardware errors in the measured observables with noise
mitigation methods which we describe in Sec. 4.1 – 4.3 in the same order as we
apply them as Qiskit transpiler passes [Qiskit contributors(2023)]. Hardware
results are in Sec. 4.4. They show a factor of six reduction in circuit depth
of an MPF implementation over a PF.

Figure 3:  Transpiler steps that implement error mitigation. (a) A single
Trotter step on four qubits with $H_{1}$ and $H_{2}$ implemented as single-
and two-qubit gates, respectively. (b) Pauli Twirling applied to the $R_{ZZ}$
gates. The Pauli gates $U_{j}$ applied to the qubits before the
$R_{ZZ}(\theta)$ gates are chosen uniformly from the set $\mathbb{G}_{ZZ}$,
which commutes with the twirled gate. The same Pauli gates are applied after
the $R_{ZZ}(\theta)$ gate so that the net operation is unchanged. (c) A pulse-
efficient transpiler pass scales the cross-resonance pulses of the backend-
provided calibrated CNOT gates to the intended rotation angle $\theta$. The
$R_{ZZ}$ gates are implemented as echoed cross-resonance gates shown in the
inset. (d) We identify all $R_{ZX}(-\theta/2)XR_{ZX}(\theta/2)$ gate sequences
in the first Trotter layer and append the sequence
$R_{ZX}(-\theta_{c})R_{ZX}(\theta_{c})$ to perform ZNE. (e) Pulse sequence
implementing $R_{ZX}(\theta_{c})$ with $\theta_{c}=\pi$. The width $w$ of the
flat-top of the CR pulses and the corresponding rotation angle $\theta_{c}$
are chosen to produce a target stretch factor $c$.

###  4.1 Pauli Twirling

Pauli twirling (PT) converts an arbitrary noise channel into a stochastic
Pauli error channel, suppressing off-diagonal coherent error contributions
[Wallman and Emerson(2016), Kim et al.(2023)Kim, Wood, Yoder, Merkel,
Gambetta, Temme, and Kandala]. We first apply PT to the $R_{ZZ}$ gates, the
dominant error source, to convert their noise channels into a stochastic Pauli
error channel. PT is implemented by sandwiching each $R_{ZZ}$ gate between
randomly sampled Pauli gates chosen so that the net operation of the circuit
is unaffected in the absence of noise. The gate noise is reduced to a
stochastic form by averaging the results of $N_{\text{PT}}\in\mathbb{N}$
logically equivalent random circuits [Wallman and Emerson(2016), Cai and
Benjamin(2019), Hashim et al.(2021)Hashim, Naik, Morvan, Ville, Mitchell,
Kreikebaum, Davis, Smith, Iancu, O’Brien, Hincks, Wallman, Emerson, and
Siddiqi, Chen et al.(2022)Chen, Burdick, Yao, Orth, and Iadecola, Kim et
al.(2023)Kim, Wood, Yoder, Merkel, Gambetta, Temme, and Kandala]. As in Ref.
[Kim et al.(2023)Kim, Wood, Yoder, Merkel, Gambetta, Temme, and Kandala], we
randomly sample twirling gates from the set $\mathbb{G}_{ZZ}=\\{I\otimes
I,X\otimes X,Y\otimes Y,Z\otimes Z,X\otimes Y,Y\otimes X,Z\otimes I,I\otimes
Z\\}$ that commutes with $R_{ZZ}(\theta)$. We apply the same Pauli gates
before and after the $R_{ZZ}(\theta)$ gate, see Fig. 3(a) and (b). Each
quantum circuit is executed with $N/N_{\text{PT}}$ shots on that hardware
where $N$ is the total shot budget.

###  4.2 Pulse efficient transpilation and calibration

Second, we use a pulse-efficient circuit transpilation which replaces the
$R_{ZZ}(\theta)$ gates by $R_{ZX}$ gates with echos exposed in the circuit as
presented in Ref. [Stenger et al.(2021)Stenger, Bronn, Egger, and Pekker,
Earnest et al.(2021)Earnest, Tornow, and Egger], see Fig. 3(c). This avoids
the onerous double CNOT implementation of $R_{ZZ}(\theta)$. The $R_{ZX}$ gates
are implemented by scaling the hardware native cross-resonance (CR) pulses
[Sheldon et al.(2016)Sheldon, Magesan, Chow, and Gambetta, Sundaresan et
al.(2020)Sundaresan, Lauer, Pritchett, Magesan, Jurcevic, and Gambetta] found
in the calibrated schedules of the CNOT gates. Since we run a Trotter
simulation where the rotation angles $\theta$ are small the flat-top of the
square pulses with Gaussian flanks vanish. We must therefore implement
$R_{ZX}(\theta)$ with Gaussian CR pulses of varying amplitude. To overcome
non-linear effects we calibrate the amplitude of the CR pulses, see Appendix
H.

###  4.3 Zero-noise extrapolation with scaled cross-resonance gates

The noise in quantum computers can be mitigated with zero-noise extrapolation
(ZNE) which runs logically equivalent circuits with different noise levels
[Temme et al.(2017)Temme, Bravyi, and Gambetta, Li and Benjamin(2017)].
Ideally, such circuits are implemented by reducing the amplitude and
increasing the duration of all the underlying pulses to increase the noise
[Kandala et al.(2019)Kandala, Temme, Córcoles, Mezzacapo, Chow, and Gambetta].
This requires extensive gate calibration, which is impractical and reduces the
availability of a cloud-based quantum computing system. To overcome this
limitation digital ZNE replaces each CNOT gate with $2n+1$ CNOT gates
[Dumitrescu et al.(2018)Dumitrescu, McCaskey, Hagen, Jansen, Morris,
Papenbrock, Pooser, Dean, and Lougovski, Stamatopoulos et
al.(2020)Stamatopoulos, Egger, Sun, Zoufal, Iten, Shen, and Woerner, Carrera
Vazquez and Woerner(2021)]. This, however, limits the possible stretch factors
to $c=2n+1$. Similarly, the Mitiq package [LaRose et al.(2022)LaRose, M.,
Kaiser, Karalekas, Alves, Czarnik, Mandouh, Gordon, Hindy, Robertson, T.,
Shammah, and Zeng] inserts pairs of CNOT gates but allows for fractional
stretch factors by repeating selected gate layers instead of repeating all
CNOT gates.

Here, we implement ZNE by adding the gate sequence
$R_{ZX}(-\theta_{c})R_{ZX}(\theta_{c})$ after the first layer of
$R_{ZZ}(\theta)$ gates, see Fig. 3(d). Crucially, the rotation angle
$\theta_{c}\neq\theta$ which allows us to create a _continuous_ range of
stretch factors by linearly scaling the width $w$ of the echoed CR pulses in
the $R_{ZX}(\pm\theta_{c})$ gates as discussed in Ref. [Stenger et
al.(2021)Stenger, Bronn, Egger, and Pekker, Earnest et al.(2021)Earnest,
Tornow, and Egger], see Fig. 3(e). Finally, we calculate the effective stretch
factor $c$ of each circuit as the fraction between the duration of the
stretched circuit and the duration of the original circuit, both without
measurements. The smallest possible stretch factor $c_{\min}$ we can
implement, indicated in Fig. 4(b), occurs when the flat-top $w$ of the CR
pulses vanishes, see Fig. 3(e). By scaling only the $R_{ZX}$ gates rather than
each gate in the circuit, it is possible to use a wider range of stretch
factors for the same $T_{1}$ times of the hardware. Each observable we measure
is then extrapolated to the zero-noise limit with an exponential fit
$ae^{-bc}+d$ where $c$ is the stretch factor and $a$, $b$, and $d$ are fit
parameters, see Fig. 4(b). We compare our ZNE method to the Mitiq package in
Appendix I.

###  4.4 Results

Figure 4:  Global magnetization of a five-spin Ising model. (a) The global
magnetization is the average local magnetization of qubits 13, 14, 16, 19, and
20. (b) Zero-noise extrapolation of $\langle Z(k_{j})\rangle$ for $k_{j}=2$
Trotter steps. The black dots represent $\langle Z(2)\rangle$ at different
circuit stretch factors $c$. (c) Error on the simulated global magnetization
relative to the exact value $\langle Z\rangle^{*}$ obtained numerically by
exponentiating $H$, i.e., $\absolutevalue{\langle Z\rangle-\langle
Z\rangle^{*}}/\absolutevalue{\langle Z\rangle^{*}}$. The blue dots show a
noiseless state-vector simulation. The corresponding hardware measured ZNE
expectation values are depicted as green stars. The error on the global
magnetization obtained from MPFs based on the state-vector simulation and
hardware measurements are shown as blue pluses and green crosses,
respectively. The Trotter exponents underlying the MPF are shown in square
brackets. The results were obtained between May 19${}^{\text{th}}$ and June
2${}^{\text{nd}}$ 2022.

We compute the dynamics of the Ising model with five linearly connected spins
and $J=0.5$, $h=1$, and an evolution time of $t=0.5$. To select five good
qubits on _ibmq_montreal_ we perform quantum process tomography [Mohseni et
al.(2008)Mohseni, Rezakhani, and Lidar, Bialczak et al.(2010)Bialczak,
Ansmann, Hofheinz, Lucero, Neeley, O’Connell, Sank, Wang, Wenner, Steffen,
Cleland, and Martinis] of scaled-CR $R_{ZX}(\theta_{c})R_{ZX}(-\theta_{c})$
gates and retain the five linearly connected qubits 13, 14, 16, 19, and 20 as
their gate fidelity decay visually resembles their $T_{1}$ decay, see Appendix
J. We calculate the global magnetization $\langle Z\rangle$ by measuring the
local magnetization $\langle Z_{i}\rangle$ of each qubit and take the average
$\langle Z\rangle=\frac{1}{5}\sum_{i=1}^{5}\langle Z_{i}\rangle$, see Fig.
4(a). As benchmark we compute the first-order Trotter formula $S_{1}$ with
exponents $k_{j}\in\\{1,...,10\\}$ on a noiseless state-vector simulator and
compare it to the exact analytical value of $e^{-iHt}$, see Eq. (2) and the
blue dots in Fig. 4(c). Next, we compute $S_{1}$ on _ibmq_montreal_ with the
error mitigation strategies described in Sec. 4.1-4.3. The properties of
_ibmq_montreal_ are shown in Appendix J. We also ran the corresponding
circuits with and without dynamical decoupling but did not see a meaningful
effect. We elaborate on this in Appendix G. For each Trotter step we create
eight random PT circuits by dressing the $R_{ZZ}$ gates with gates uniformly
drawn from $\mathbb{G}_{ZZ}$. Each of these eight circuits is executed 20
times with different stretch-factors $\theta_{c}$ linearly spaced in the
interval $[\theta_{c}^{\min},\theta_{c}^{\max}]$. Here, $\theta_{c}^{\max}$ is
chosen such that the pulse-schedule of ten Trotter steps does not exceed the
shortest $T_{1}$ time. For every $\theta_{c}$ each of the eight circuits is
evaluated with $10^{5}/8$ shots to keep the physical and sampling noise below
$\varepsilon^{\prime}=10^{-2}$. This process is repeated ten times to gather
statistics. Therefore, a zero-noise extrapolated global magnetization $\langle
Z(k_{j})\rangle$, indicated as green stars in Fig. 4(b) and (c), is computed
with a total of $10\times 20\times 10^{5}=2\cdot 10^{7}$ shots. Before
computing each $\langle Z(k_{j})\rangle$ we calibrate the $R_{ZX}$ gates using
the procedure in Appendix H. All the measured expectation values are corrected
for readout errors with a tensor product of individual qubit calibration
matrices [Bravyi et al.(2020)Bravyi, Gosset, König, and Tomamichel, Kim et
al.(2023)Kim, Wood, Yoder, Merkel, Gambetta, Temme, and Kandala]. The noise
mitigation significantly increases the accuracy in the measured $\langle
Z(k_{j})\rangle$, see Fig. 4(b). For up to four Trotter steps the measured
$\langle Z(k_{j})\rangle$ agree well with the noiseless benchmark. At Trotter
steps $k_{j}>4$ the $\langle Z(k_{j})\rangle$ have a lower error than
predicted by the noiseless benchmark which is a coincidence.

We build the MPF estimation of the zero-noise extrapolated values of $\langle
Z\rangle$ using Eq. (7). We select the Trotter exponents $k_{j}$ to combine
and the corresponding weights $a_{j}$ following the optimization algorithm
described in Appendix E, which also contains the values of $\vec{k}$ and
$\vec{a}$, with the threshold $\norm{\vec{a}}_{1}\leq 3$. Crucially, we
observe that the MPF significantly reduces the error. For example, the
combinations [1, 2], [1, 3], and [2, 4] all have the same accuracy as a first-
order formula with 24 Trotter steps, see green crosses in Fig. 4(c). This is a
reduction in circuit depth by a factor of twelve to six over the product
formula.

##  5  Conclusion

We present a hybrid quantum-classical approach to MPFs that requires neither
additional qubits nor controlled operations as do previous approaches based on
an LCU [Childs and Wiebe(2012), Hao Low et al.(2019)Hao Low, Kliuchnikov, and
Wiebe]. Our method approximates a Hamiltonian evolution with the same accuracy
as PFs or MPFs implemented with an LCU but with significantly shallower
circuits by implementing each unitary in the LCU as a standalone circuit.
Therefore, even though the total number of matrix exponentials in both
approaches is the same, the reduction in gate count for the deepest circuit
can be substantial as illustrated by the $l=3$ example for which our approach
has 22 times less two qubit gates than an LCU implementation when
$\vec{k}=[1,2,7]$. Reducing the gate and qubit count makes it possible to
execute MPFs on noisy hardware that achieve accurate results, such as those
shown in Fig. 4, that would not be possible if implemented as an LCU.

The fast convergence of MPFs to the exact operator $e^{-iHt}$ also means that
the approximation error soon becomes smaller than the noise levels on the
hardware. Error mitigation is thus crucial to suppress hardware noise. As the
hardware and error mitigation techniques improve, the advantages of our MPF
method will become increasingly noticeable. To mitigate errors we used PT and
pulse-efficient transpilation and introduced an implementation of ZNE that
provides a continuous range of stretch factors while simultaneously avoiding
onerous pulse calibration. Crucially, this method overcomes the
$c\in\\{1,3,5,...\\}$ limitation of the digital ZNE which rapidly becomes
infeasible for deep circuits. With these error mitigation methods we
demonstrated on hardware that MPFs can achieve the same accuracy as a PF but
with up to a factor of 12 reduction in circuit depth for a five spin Ising
Hamiltonian.

Our approach has shallower circuits than MPFs based on an LCU for $t\leq 1$.
Future work may investigate whether there is a regime for large values of $t$
in which the LCU circuits are shallower than the classical combination of
expectation values we propose. For example, since an MPF implemented as an LCU
is one quantum circuit, it may be possible to repeat $t/\Delta t$ times the
circuit for a time $\Delta t$. This may use less Trotter steps than a
classical combination of expectation values. In addition, tighter bounds could
be derived, since as mentioned in Sec. 3, our experiments showed that the
scheme still works even when $t/k_{1}>1$. Regarding ZNE as presented in Sec.
4.3: a deeper study of the optimal number of stretch factors is warranted.
Furthermore, we empirically found that inserting the gate sequence
$R_{ZX}(-\theta_{c})R_{ZX}(\theta_{c})$ after the gates in the first layer of
$R_{ZZ}$ gates allowed us to mitigate errors. However, the choices of the
sequence and where to insert them are not unique. While we tried several
possibilities, a more rigorous study is left as further work. We obtained our
ZNE expectation values with an exponential fit using the effective stretch
factors as the independent variable. These choices could be either studied
empirically in more detail or backed up with theory.

We conjecture that our scheme is optimal when used with $S_{2}$, where optimal
means requiring the shallowest circuits to achieve a given approximation
error. In summary, the linear combinations of expectation values makes it
possible to use the MPF approach on noisy hardware. This allows us to reduce
circuit depth by orders of magnitude which we believe is crucial to reach a
quantum advantage for Hamiltonian simulation on noisy quantum hardware.

##  6  Acknowledgements

The authors thank Nate Earnest, Andreas Fuhrer, Youngseok Kim, David Layden
and Alexander Miessen for useful discussions. This work was supported as a
part of NCCR QSIT, a National Centre of Competence (or Excellence) in
Research, funded by the Swiss National Science Foundation (grant number
51NF40-185902).

## References

  * [ibm()] IBM Quantum.  https://quantum-computing.ibm.com/.
  * [Aharonov and Ta-Shma(2003)] D. Aharonov and A. Ta-Shma.  Adiabatic quantum state generation and statistical zero knowledge.  In _Proc. 35th Annu. ACM Symp. Theory Comput._ , page 20–29, 2003\.  doi: 10.1145/780542.780546.
  * [Ahokas(2004)] G. Ahokas.  Improved algorithms for approximate quantum Fourier transforms and sparse Hamiltonian simulations.  2004\.  doi: 10.11575/PRISM/22839.
  * [Alexander et al.(2020)Alexander, Kanazawa, Egger, Capelluto, Wood, Javadi-Abhari, and McKay] T. Alexander, N. Kanazawa, D. J. Egger, L. Capelluto, C. J. Wood, A. Javadi-Abhari, and D. C. McKay.  Qiskit pulse: programming quantum computers through the cloud with pulses.  _Quantum Sci. Technol._ , 5(4):044006, aug 2020\.  doi: 10.1088/2058-9565/aba404.
  * [An et al.(2022)An, Fang, and Lin] D. An, D. Fang, and L. Lin.  Time-dependent Hamiltonian simulation of highly oscillatory dynamics and superconvergence for Schrödinger equation.  _Quantum_ , 6:690, 2022.  doi: 10.22331/q-2022-04-15-690.
  * [Bernien et al.(2017)Bernien, Schwartz, Keesling, Levine, Omran, Pichler, Choi, Zibrov, Endres, Greiner, Vuletić, and Lukin] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin.  Probing many-body dynamics on a 51-atom quantum simulator.  _Nature_ , 551(7682):579–584, Nov 2017.  ISSN 1476-4687.  doi: 10.1038/nature24622.
  * [Berry et al.(2007)Berry, Ahokas, Cleve, and Sanders] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders.  Efficient quantum algorithms for simulating sparse Hamiltonians.  _Commun. Math. Phys._ , 270(2):359–371, Mar 2007\.  ISSN 1432-0916.  doi: 10.1007/s00220-006-0150-x.
  * [Berry et al.(2020)Berry, Childs, Su, Wang, and Wiebe] D. W. Berry, A. M. Childs, Y. Su, X. Wang, and N. Wiebe.  Time-dependent Hamiltonian simulation with $L^{1}$-norm scaling.  _Quantum_ , 4:254, 2020.  doi: 10.22331/q-2020-04-20-254.
  * [Bialczak et al.(2010)Bialczak, Ansmann, Hofheinz, Lucero, Neeley, O’Connell, Sank, Wang, Wenner, Steffen, Cleland, and Martinis] R. C. Bialczak, M. Ansmann, M. Hofheinz, E. Lucero, M. Neeley, A. D. O’Connell, D. Sank, H. Wang, J. Wenner, M. Steffen, A. N. Cleland, and J. M. Martinis.  Quantum process tomography of a universal entangling gate implemented with Josephson phase qubits.  _Nat. Phys._ , 6(6):409–413, Jun 2010.  ISSN 1745-2481.  doi: 10.1038/nphys1639.
  * [Blanes et al.(1999)Blanes, Casas, and Ros] S. Blanes, F. Casas, and J. Ros.  Extrapolation of symplectic Integrators.  _Celestial Mech. Dyn. Astron._ , 75:149–161, 1999.  doi: 10.1023/A:1008364504014.
  * [Bravyi et al.(2020)Bravyi, Gosset, König, and Tomamichel] S. Bravyi, D. Gosset, R. König, and M. Tomamichel.  Quantum advantage with noisy shallow circuits.  _Nat. Phys._ , 16(10):1040–1045, Oct 2020.  ISSN 1745-2481.  doi: 10.1038/s41567-020-0948-z.
  * [Browaeys and Lahaye(2020)] A. Browaeys and T. Lahaye.  Many-body physics with individually controlled Rydberg atoms.  _Nat. Phys._ , 16(2):132–142, Feb 2020.  ISSN 1745-2481.  doi: 10.1038/s41567-019-0733-z.
  * [Butcher(2016)] J. C. Butcher.  _Numerical Methods for Ordinary Differential Equations_.  Wiley, 3 edition, 7 2016.  ISBN 9781119121534.  doi: 10.1002/9781119121534.  Chapter 33.
  * [Cai and Benjamin(2019)] Z. Cai and S. C. Benjamin.  Constructing smaller Pauli twirling sets for arbitrary error channels.  _Sci. Rep._ , 9(1):11281, Aug 2019.  ISSN 2045-2322.  doi: 10.1038/s41598-019-46722-7.
  * [Campbell(2019)] E. Campbell.  Random compiler for fast hamiltonian simulation.  _Phys. Rev. Lett._ , 123:070503, Aug 2019.  doi: 10.1103/PhysRevLett.123.070503.
  * [Carrera Vazquez and Woerner(2021)] A. Carrera Vazquez and S. Woerner.  Efficient state preparation for quantum amplitude estimation.  _Phys. Rev. Applied_ , 15:034027, Mar 2021.  doi: 10.1103/PhysRevApplied.15.034027.
  * [Chen et al.(2022)Chen, Burdick, Yao, Orth, and Iadecola] I. Chen, B. Burdick, Y. Yao, P. P. Orth, and T. Iadecola.  Error-mitigated simulation of quantum many-body scars on quantum computers with pulse-level control.  _Phys. Rev. Res._ , 4:043027, Oct 2022.  doi: 10.1103/PhysRevResearch.4.043027.
  * [Chiesa et al.(2019)Chiesa, Tacchino, Grossi, Santini, Tavernelli, Gerace, and Carretta] A. Chiesa, F. Tacchino, M. Grossi, P. Santini, I. Tavernelli, D. Gerace, and S. Carretta.  Quantum hardware simulating four-dimensional inelastic neutron scattering.  _Nat. Phys._ , 15(5):455–459, May 2019.  ISSN 1745-2481.  doi: 10.1038/s41567-019-0437-4.
  * [Childs(2010)] A. M. Childs.  On the relationship between continuous- and discrete-time quantum walk.  _Commun. Math. Phys._ , 294(2):581–603, Mar 2010\.  ISSN 1432-0916.  doi: 10.1007/s00220-009-0930-1.
  * [Childs and Kothari(2010)] A. M. Childs and R. Kothari.  Limitations on the simulation of non-sparse Hamiltonians.  _Quantum Info. Comput._ , 10(7):669–684, July 2010.  ISSN 1533-7146.  doi: 10.5555/2011373.2011380.
  * [Childs and Kothari(2011)] A. M. Childs and R. Kothari.  Simulating sparse Hamiltonians with star decompositions.  In _Proc. 5th Conf. Theory Quantum Comput. Commun. Cryptogr._ , pages 94–103, 2011.  ISBN 3-642-18072-8, 978-3-642-18072-9.  doi: 10.5555/1946127.19461355.
  * [Childs and Su(2019)] A. M. Childs and Y. Su.  Nearly optimal lattice simulation by product formulas.  _Phys. Rev. Lett._ , 123:050503, Aug 2019.  doi: 10.1103/PhysRevLett.123.050503.
  * [Childs and Wiebe(2012)] A. M. Childs and N. Wiebe.  Hamiltonian simulation using linear combinations of unitary operations.  _Quantum Info. Comput._ , 12(11–12):901–924, November 2012.  ISSN 1533-7146.  doi: 10.5555/2481569.2481570.
  * [Childs et al.(2003)Childs, Cleve, Deotto, Farhi, Gutmann, and Spielman] A. M. Childs, R. Cleve, E. Deotto, E. Farhi, S. Gutmann, and D. A. Spielman.  Exponential algorithmic speedup by a quantum walk.  In _Proc. 35th Annu. ACM Symp. Theory Comput._ , pages 59–68, 2003\.  ISBN 1-58113-674-9.  doi: 10.1145/780542.780552.
  * [Childs et al.(2018)Childs, Maslov, Nam, Ross, and Su] A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su.  Toward the first quantum simulation with quantum speedup.  _PNAS_ , 115(38):9456–9461, 2018.  ISSN 0027-8424.  doi: 10.1073/pnas.1801723115.
  * [Childs et al.(2021)Childs, Su, Tran, Wiebe, and Zhu] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu.  Theory of Trotter error with commutator scaling.  _Phys. Rev. X_ , 11(1):011020, February 2021.  doi: 10.1103/PhysRevX.11.011020.
  * [Chin(2010)] S. A. Chin.  Multi-product splitting and Runge-Kutta-Nyström integrators.  _Celestial Mech. Dyn. Astron._ , 106(4):391–406, Apr 2010.  ISSN 1572-9478.  doi: 10.1007/s10569-010-9255-9.
  * [Cîrstoiu et al.(2020)Cîrstoiu, Holmes, Iosue, Cincio, Coles, and Sornborger] C. Cîrstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger.  Variational fast forwarding for quantum simulation beyond the coherence time.  _npj Quantum Inf._ , 6(1):82, Sep 2020.  ISSN 2056-6387.  doi: 10.1038/s41534-020-00302-0.
  * [Corless et al.(1996)Corless, Gonnet, Hare, Jeffrey, and Knuth] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth.  On the Lambert W function.  _Adv. Comput. Math._ , 5(4):329–359, 1996.  doi: 10.1007/BF02124750.
  * [Di Paolo et al.(2020)Di Paolo, Barkoutsos, Tavernelli, and Blais] A. Di Paolo, P. Kl. Barkoutsos, I. Tavernelli, and A. Blais.  Variational quantum simulation of ultrastrong light-matter coupling.  _Phys. Rev. Research_ , 2:033364, Sep 2020.  doi: 10.1103/PhysRevResearch.2.033364.
  * [Dumitrescu et al.(2018)Dumitrescu, McCaskey, Hagen, Jansen, Morris, Papenbrock, Pooser, Dean, and Lougovski] E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean, and P. Lougovski.  Cloud quantum computing of an atomic nucleus.  _Phys. Rev. Lett._ , 120:210501, May 2018.  doi: 10.1103/PhysRevLett.120.210501.
  * [Earnest et al.(2021)Earnest, Tornow, and Egger] N. Earnest, C. Tornow, and D. J. Egger.  Pulse-efficient circuit transpilation for quantum applications on cross-resonance-based hardware.  _Phys. Rev. Research_ , 3:043088, Oct 2021.  doi: 10.1103/PhysRevResearch.3.043088.
  * [Egger and Wilhelm(2013)] D. J. Egger and F. K. Wilhelm.  Multimode circuit quantum electrodynamics with hybrid metamaterial transmission lines.  _Phys. Rev. Lett._ , 111:163601, Oct 2013.  doi: 10.1103/PhysRevLett.111.163601.
  * [Endo et al.(2019)Endo, Zhao, Li, Benjamin, and Yuan] S. Endo, Q. Zhao, Y. Li, S. Benjamin, and X. Yuan.  Mitigating algorithmic errors in a Hamiltonian simulation.  _Phys. Rev. A_ , 99:012334, Jan 2019.  doi: 10.1103/PhysRevA.99.012334.
  * [Endo et al.(2020)Endo, Sun, Li, Benjamin, and Yuan] S. Endo, J. Sun, Y. Li, S. C. Benjamin, and X. Yuan.  Variational quantum simulation of general processes.  _Phys. Rev. Lett._ , 125:010501, Jun 2020.  doi: 10.1103/PhysRevLett.125.010501.
  * [Faehrmann et al.(2022)Faehrmann, Steudtner, Kueng, Kieferova, and Eisert] P. K. Faehrmann, M. Steudtner, R. Kueng, M. Kieferova, and J. Eisert.  Randomizing multi-product formulas for Hamiltonian simulation.  _Quantum_ , 6:806, September 2022.  ISSN 2521-327X.  doi: 10.22331/q-2022-09-19-806.
  * [Feynman(1982)] R. P. Feynman.  Simulating physics with computers.  _Int. J. Theor. Phys._ , 21(6):467–488, Jun 1982\.  ISSN 1572-9575.  doi: 10.1007/BF02650179.
  * [Fischer et al.(2022)Fischer, Miller, Tacchino, Barkoutsos, Egger, and Tavernelli] L. E. Fischer, D. Miller, F. Tacchino, P. Kl. Barkoutsos, D. J. Egger, and I. Tavernelli.  Ancilla-free implementation of generalized measurements for qubits embedded in a qudit space.  _Phys. Rev. Research_ , 4:033027, Jul 2022.  doi: 10.1103/PhysRevResearch.4.033027.
  * [Frisk Kockum et al.(2019)Frisk Kockum, Miranowicz, De Liberato, Savasta, and Nori] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori.  Ultrastrong coupling between light and matter.  _Nat. Rev. Phys._ , 1(1):19–40, Jan 2019.  ISSN 2522-5820.  doi: 10.1038/s42254-018-0006-2.
  * [Gacon et al.(2021)Gacon, Zoufal, Carleo, and Woerner] J. Gacon, C. Zoufal, G. Carleo, and S. Woerner.  Simultaneous Perturbation Stochastic Approximation of the Quantum Fisher Information.  _Quantum_ , 5:567, October 2021.  ISSN 2521-327X.  doi: 10.22331/q-2021-10-20-567.
  * [García-Pérez et al.(2021)García-Pérez, Rossi, Sokolov, Tacchino, Barkoutsos, Mazzola, Tavernelli, and Maniscalco] G. García-Pérez, M. A. C. Rossi, B. Sokolov, F. Tacchino, P. Kl. Barkoutsos, G. Mazzola, I. Tavernelli, and S. Maniscalco.  Learning to measure: Adaptive informationally complete generalized measurements for quantum algorithms.  _PRX Quantum_ , 2:040342, Nov 2021.  doi: 10.1103/PRXQuantum.2.040342.
  * [Giurgica-Tiron et al.(2020)Giurgica-Tiron, Hindy, LaRose, Mari, and Zeng] T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng.  Digital zero noise extrapolation for quantum error mitigation.  In _IEEE Int. Conf. on QCE_ , pages 306–316, 2020.  doi: 10.1109/QCE49297.2020.00045.
  * [Gross and Bloch(2017)] C. Gross and I. Bloch.  Quantum simulations with ultracold atoms in optical lattices.  _Science_ , 357(6355):995–1001, 2017.  doi: 10.1126/science.aal3837.
  * [Hadfield et al.(2022)Hadfield, Bravyi, Raymond, and Mezzacapo] C. Hadfield, S. Bravyi, R. Raymond, and A. Mezzacapo.  Measurements of quantum Hamiltonians with locally-biased classical shadows.  _Commun. Math. Phys._ , 391(3):951–967, May 2022\.  ISSN 1432-0916.  doi: 10.1007/s00220-022-04343-8.
  * [Hagan and Wiebe(2022)] M. Hagan and N. Wiebe.  Composite quantum simulations.  2022\.  URL https://arxiv.org/abs/2206.06409.
  * [Hao Low et al.(2019)Hao Low, Kliuchnikov, and Wiebe] G. Hao Low, V. Kliuchnikov, and N. Wiebe.  Well-conditioned multiproduct Hamiltonian simulation.  2019\.  URL https://arxiv.org/abs/1907.11679.
  * [Harrow et al.(2009)Harrow, Hassidim, and Lloyd] A. W. Harrow, A. Hassidim, and S. Lloyd.  Quantum algorithm for linear systems of equations.  _Phys. Rev. Lett._ , 103:150502, Oct 2009.  doi: 10.1103/PhysRevLett.103.150502.
  * [Hashim et al.(2021)Hashim, Naik, Morvan, Ville, Mitchell, Kreikebaum, Davis, Smith, Iancu, O’Brien, Hincks, Wallman, Emerson, and Siddiqi] A. Hashim, R. K. Naik, A. Morvan, J. L. Ville, B. Mitchell, J. M. Kreikebaum, M. Davis, E. Smith, C. Iancu, K. P. O’Brien, I. Hincks, J. J. Wallman, J. Emerson, and I. Siddiqi.  Randomized compiling for scalable quantum computing on a noisy superconducting quantum processor.  _Phys. Rev. X_ , 11:041039, Nov 2021.  doi: 10.1103/PhysRevX.11.041039.
  * [He et al.(2020)He, Nachman, de Jong, and Bauer] A. He, B. Nachman, W. A. de Jong, and C. W. Bauer.  Zero-noise extrapolation for quantum-gate error mitigation with identity insertions.  _Phys. Rev. A_ , 102:012426, Jul 2020.  doi: 10.1103/PhysRevA.102.012426.
  * [Hoorfar and Hassani(2008)] A. Hoorfar and M. Hassani.  Inequalities on the Lambert W function and hyperpower function.  _JIPAM_ , 9, 01 2008.  URL http://eudml.org/doc/130024.
  * [Huang et al.(2020)Huang, Kueng, and Preskill] H.-Y. Huang, R. Kueng, and J. Preskill.  Predicting many properties of a quantum system from very few measurements.  _Nat. Phys._ , 16(10):1050–1057, Oct 2020.  ISSN 1745-2481.  doi: 10.1038/s41567-020-0932-7.
  * [Indrajeet et al.(2020)Indrajeet, Wang, Hutchings, Taketani, Wilhelm, LaHaye, and Plourde] S. Indrajeet, H. Wang, M.D. Hutchings, B.G. Taketani, F. K. Wilhelm, M.D. LaHaye, and B.L.T. Plourde.  Coupling a superconducting qubit to a left-handed metamaterial resonator.  _Phys. Rev. Applied_ , 14:064033, Dec 2020.  doi: 10.1103/PhysRevApplied.14.064033.
  * [Jurcevic et al.(2021)Jurcevic, Javadi-Abhari, Bishop, Lauer, Bogorin, Brink, Capelluto, Günlük, Itoko, Kanazawa, and _et al._] P. Jurcevic, A. Javadi-Abhari, L. S. Bishop, I. Lauer, D. F. Bogorin, M. Brink, L. Capelluto, O. Günlük, T. Itoko, N. Kanazawa, and _et al._ Demonstration of quantum volume 64 on a superconducting quantum computing system.  _Quantum Sci. Technol._ , 6(2):025020, mar 2021\.  doi: 10.1088/2058-9565/abe519.
  * [Kamakari et al.(2022)Kamakari, Sun, Motta, and Minnich] H. Kamakari, S. N. Sun, M. Motta, and Austin J. Minnich.  Digital quantum simulation of open quantum systems using quantum imaginary–time evolution.  _PRX Quantum_ , 3:010320, Feb 2022.  doi: 10.1103/PRXQuantum.3.010320.
  * [Kanazawa et al.(2023)Kanazawa, Egger, Ben-Haim, Zhang, Shanks, Aleksandrowicz, and Wood] Naoki Kanazawa, Daniel J. Egger, Yael Ben-Haim, Helena Zhang, William E. Shanks, Gadi Aleksandrowicz, and Christopher J. Wood.  Qiskit experiments: A python package to characterize and calibrate quantum computers.  _J. Open Source Softw._ , 8(84):5329, 2023.  doi: 10.21105/joss.05329.
  * [Kandala et al.(2019)Kandala, Temme, Córcoles, Mezzacapo, Chow, and Gambetta] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta.  Error mitigation extends the computational reach of a noisy quantum processor.  _Nature_ , 567(7749):491–495, Mar 2019.  ISSN 1476-4687.  doi: 10.1038/s41586-019-1040-7.
  * [Kim et al.(2023)Kim, Wood, Yoder, Merkel, Gambetta, Temme, and Kandala] Y. Kim, C. J. Wood, T. J. Yoder, S. T. Merkel, J. M. Gambetta, K. Temme, and A. Kandala.  Scalable error mitigation for noisy quantum circuits produces competitive expectation values.  _Nat. Phys._ , 19(5):752–759, May 2023.  doi: 10.1038/s41567-022-01914-3.
  * [LaRose et al.(2022)LaRose, M., Kaiser, Karalekas, Alves, Czarnik, Mandouh, Gordon, Hindy, Robertson, T., Shammah, and Zeng] R. LaRose, A. M., S. Kaiser, P. J. Karalekas, A. A. Alves, P. Czarnik, M. El Mandouh, M. H. Gordon, Y. Hindy, A. Robertson, P. T., N. Shammah, and W. J. Zeng.  Mitiq: A software package for error mitigation on noisy quantum computers.  _Quantum_ , 6:774, August 2022.  ISSN 2521-327X.  doi: 10.22331/q-2022-08-11-774.
  * [Layden(2022)] D. Layden.  First-order Trotter error from a second-order perspective.  _Phys. Rev. Lett._ , 128:210501, May 2022.  doi: 10.1103/PhysRevLett.128.210501.
  * [Leggett et al.(1987)Leggett, Chakravarty, Dorsey, Fisher, Garg, and Zwerger] A. J. Leggett, S. Chakravarty, A. T. Dorsey, Matthew P. A. Fisher, A. Garg, and W. Zwerger.  Dynamics of the dissipative two-state system.  _Rev. Mod. Phys._ , 59:1–85, Jan 1987.  doi: 10.1103/RevModPhys.59.1.
  * [Li and Benjamin(2017)] Y. Li and S. C. Benjamin.  Efficient variational quantum simulator incorporating active error minimization.  _Phys. Rev. X_ , 7:021050, Jun 2017.  doi: 10.1103/PhysRevX.7.021050.
  * [Lidar(2014)] D. A. Lidar.  _Review of Decoherence-Free Subspaces, Noiseless Subsystems, and Dynamical Decoupling_ , pages 295–354.  John Wiley & Sons, Ltd, 2014.  ISBN 9781118742631.  doi: https://doi.org/10.1002/9781118742631.ch11.
  * [Lloyd(1996)] S. Lloyd.  Universal quantum simulators.  _Science_ , 273(5278):1073–1078, 1996.  doi: 10.1126/science.273.5278.1073.
  * [Magesan and Gambetta(2020)] E. Magesan and J. M. Gambetta.  Effective Hamiltonian models of the cross-resonance gate.  _Phys. Rev. A_ , 101:052308, May 2020.  doi: 10.1103/PhysRevA.101.052308.
  * [McLachlan(1964)] A. D. McLachlan.  A variational solution of the time-dependent Schrodinger equation.  _Mol. Phys._ , 8(1):39–44, 1964.  doi: 10.1080/00268976400100041.
  * [Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli] A. Miessen, P. J. Ollitrault, and I. Tavernelli.  Quantum algorithms for quantum dynamics: A performance study on the spin-boson model.  _Phys. Rev. Research_ , 3:043212, Dec 2021.  doi: 10.1103/PhysRevResearch.3.043212.
  * [Mohseni et al.(2008)Mohseni, Rezakhani, and Lidar] M. Mohseni, A. T. Rezakhani, and D. A. Lidar.  Quantum-process tomography: Resource analysis of different strategies.  _Phys. Rev. A_ , 77:032322, Mar 2008.  doi: 10.1103/PhysRevA.77.032322.
  * [Ouyang et al.(2020)Ouyang, White, and Campbell] Y. Ouyang, D. R. White, and E. T. Campbell.  Compilation by stochastic Hamiltonian sparsification.  _Quantum_ , 4:235, February 2020.  ISSN 2521-327X.  doi: 10.22331/q-2020-02-27-235.
  * [Pokharel et al.(2018)Pokharel, Anand, Fortman, and Lidar] B. Pokharel, N. Anand, B. Fortman, and D. A. Lidar.  Demonstration of fidelity improvement using dynamical decoupling with superconducting qubits.  _Phys. Rev. Lett._ , 121:220502, Nov 2018.  doi: 10.1103/PhysRevLett.121.220502.
  * [Qiskit contributors(2023)] Qiskit contributors.  Qiskit: An open-source framework for quantum computing, 2023.  URL https://doi.org/10.5281/zenodo.2562111.
  * [Richardson(1911)] L. F. Richardson.  The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam.  _Philos. Trans. R. Soc. London, Ser. A_ , 210:307–357, 1911\.  ISSN 02643952.  URL http://www.jstor.org/stable/90994.
  * [Sachdev(2011)] S. Sachdev.  _Quantum Phase Transitions_.  Cambridge University Press, 2 edition, 2011.  doi: 10.1017/CBO9780511973765.
  * [Sheldon et al.(2016)Sheldon, Magesan, Chow, and Gambetta] S. Sheldon, E. Magesan, J. M. Chow, and J. M. Gambetta.  Procedure for systematically tuning up cross-talk in the cross-resonance gate.  _Phys. Rev. A_ , 93:060302, Jun 2016.  doi: 10.1103/PhysRevA.93.060302.
  * [Sidi(2003)] A. Sidi.  _Practical Extrapolation Methods: Theory and Applications_.  Cambridge University Press, 2003.  ISBN 9780511546815.  doi: 10.1017/CBO9780511546815.
  * [Stamatopoulos et al.(2020)Stamatopoulos, Egger, Sun, Zoufal, Iten, Shen, and Woerner] N. Stamatopoulos, D. J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen, and S. Woerner.  Option Pricing using Quantum Computers.  _Quantum_ , 4:291, July 2020.  ISSN 2521-327X.  doi: 10.22331/q-2020-07-06-291.
  * [Stenger et al.(2021)Stenger, Bronn, Egger, and Pekker] J. P. T. Stenger, N. T. Bronn, D. J. Egger, and D. Pekker.  Simulating the dynamics of braiding of Majorana zero modes using an IBM quantum computer.  _Phys. Rev. Research_ , 3:033171, Aug 2021.  doi: 10.1103/PhysRevResearch.3.033171.
  * [Sundaresan et al.(2020)Sundaresan, Lauer, Pritchett, Magesan, Jurcevic, and Gambetta] N. Sundaresan, I. Lauer, E. Pritchett, E. Magesan, P. Jurcevic, and J. M. Gambetta.  Reducing unitary and spectator errors in cross resonance with optimized rotary echoes.  _PRX Quantum_ , 1:020318, Dec 2020.  doi: 10.1103/PRXQuantum.1.020318.
  * [Suzuki(1991)] M. Suzuki.  General theory of fractal path integrals with applications to many-body theories and statistical physics.  _J. Math. Phys._ , 32(2):400–407, 1991.  doi: 10.1063/1.529425.
  * [Temme et al.(2017)Temme, Bravyi, and Gambetta] K. Temme, S. Bravyi, and J. M. Gambetta.  Error mitigation for short-depth quantum circuits.  _Phys. Rev. Lett._ , 119:180509, Nov 2017.  doi: 10.1103/PhysRevLett.119.180509.
  * [The IBM Decision Optimization on Cloud team(2022)] The IBM Decision Optimization on Cloud team.  DOCPLEX, 2022.  URL https://pypi.org/project/docplex/.
  * [Tornow et al.(2022)Tornow, Kanazawa, Shanks, and Egger] C. Tornow, N. Kanazawa, W. E. Shanks, and D. J. Egger.  Minimum quantum run-time characterization and calibration via restless measurements with dynamic repetition rates.  _Phys. Rev. Applied_ , 17:064061, Jun 2022.  doi: 10.1103/PhysRevApplied.17.064061.
  * [Tripathi et al.(2022)Tripathi, Chen, Khezri, Yip, Levenson-Falk, and Lidar] V. Tripathi, H. Chen, M. Khezri, K. W. Yip, E. M. Levenson-Falk, and D. A. Lidar.  Suppression of crosstalk in superconducting qubits using dynamical decoupling.  _Phys. Rev. Appl._ , 18:024068, Aug 2022.  doi: 10.1103/PhysRevApplied.18.024068.
  * [Vazquez Carrera et al.(2022)Vazquez Carrera, Hiptmair, and Woerner] A. Vazquez Carrera, R. Hiptmair, and S. Woerner.  Enhancing the quantum linear systems algorithm using richardson extrapolation.  _ACM Trans. Quantum Comput._ , 3(1), jan 2022.  ISSN 2643-6809.  doi: 10.1145/3490631.
  * [Wallman and Emerson(2016)] J. J. Wallman and J. Emerson.  Noise tailoring for scalable quantum computation via randomized compiling.  _Phys. Rev. A_ , 94:052325, Nov 2016.  doi: 10.1103/PhysRevA.94.052325.
  * [Wan et al.(2022)Wan, Berta, and Campbell] K. Wan, M. Berta, and E. T. Campbell.  Randomized quantum algorithm for statistical phase estimation.  _Phys. Rev. Lett._ , 129:030503, Jul 2022.  doi: 10.1103/PhysRevLett.129.030503.
  * [Wecker et al.(2014)Wecker, Bauer, Clark, Hastings, and Troyer] D. Wecker, B. Bauer, B. K. Clark, M. B. Hastings, and M. Troyer.  Gate-count estimates for performing quantum chemistry on small quantum computers.  _Phys. Rev. A_ , 90:022305, Aug 2014.  doi: 10.1103/PhysRevA.90.022305.
  * [Weidenfeller et al.(2022)Weidenfeller, Valor, Gacon, Tornow, Bello, Woerner, and Egger] J. Weidenfeller, L. C. Valor, J. Gacon, C. Tornow, L. Bello, S. Woerner, and D. J. Egger.  Scaling of the quantum approximate optimization algorithm on superconducting qubit based hardware.  _Quantum_ , 6:870, December 2022.  ISSN 2521-327X.  doi: 10.22331/q-2022-12-07-870.
  * [Yongxin(2022)] Y. Yongxin.  Quantum many-body scars on QC, 2022.  URL https://gitlab.com/QANED/rydberg_mfim.
  * [Yuan et al.(2019)Yuan, Endo, Zhao, Li, and Benjamin] X. Yuan, S. Endo, Q. Zhao, Ying Li, and S. Benjamin.  Theory of variational quantum simulation.  _Quantum_ , 3:191, 2019.  doi: 10.22331/q-2019-10-07-191.
  * [Zhang et al.(2017)Zhang, Pagano, Hess, Kyprianidis, Becker, Kaplan, Gorshkov, Gong, and Monroe] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe.  Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator.  _Nature_ , 551(7682):601–604, Nov 2017.  ISSN 1476-4687.  doi: 10.1038/nature24654.
  * [Zoufal et al.(2021)Zoufal, Sutter, and Woerner] C. Zoufal, D. Sutter, and S. Woerner.  Error bounds for variational quantum time evolution.  2021\.  URL https://arxiv.org/abs/2108.00022.

##  Appendix A Error cancellations in multi-product formulas

We illustrate how errors in multi-product formulas cancel with the first-order
formula. The example is easily adapted to the $2\chi$-order. We write
$S^{k_{j}}_{1}(t/k_{j})$ as a power series in terms of the error [Chin(2010)]

$S^{k_{j}}_{1}\left({\textstyle\frac{t}{k_{j}}}\right)=e^{-iHt}+\sum_{n=1}^{\infty}A_{n}t^{n+1}/k_{j}^{n}.$
(13)

Here, the $A_{n}$ are matrices that depend only on commutators of the $H_{j}$.
In particular, they are independent of the time $t$ and the Trotter exponents
$k_{j}$. Therefore, the lower-order error terms can be cancelled out by a
careful choice of extrapolation weights $a_{j}$ that are independent of the
Hamiltonian. The $A_{n}$ are thus irrelevant in subsequent calculations. For
the MPF $M_{l,1}(t)$ to approximate $e^{-iHt}$ without bias, the extrapolation
weights must satisfy $\sum_{j=1}^{l}a_{j}=1$. Furthermore, to cancel the first
$l-1$ error terms in $S_{1}$ we impose the constraint
$\sum_{j=1}^{l-1}a_{j}/k_{j}^{1+n}=0$ for $n=0,...,l-2$. Therefore, the
resulting error is given by

$\displaystyle M_{l,1}(t)$
$\displaystyle=\sum_{j=1}^{l}a_{j}S^{k_{j}}_{1}\left({\textstyle\frac{t}{k_{j}}}\right)$
$\displaystyle=\sum_{j=1}^{l}a_{j}\left(e^{-iHt}+\sum_{n=1}^{l-1}A_{n}\frac{t^{n+1}}{k_{j}^{n}}+\sum_{n=l}^{\infty}A_{n}\frac{t^{n+1}}{k_{j}^{n}}\right)$
$\displaystyle=e^{-iHt}+\mathcal{O}\left(\frac{t^{l+1}}{k_{1}^{l}}\right).$
(14)

##  Appendix B Comparison with LCU

Here, we compare the circuit depth of MPFs implemented by classically
combining expectation values to MPFs implemented as LCUs with $l=3$
extrapolation points as example. The LCU method implements an operation
proportional to the linear combination $\kappa U+V$ for some unitary gates
$U,V$ and $\kappa\geq 0$. A linear combination of three unitaries is
implemented recursively as $\alpha_{1}(\alpha_{2}U_{1}+U_{2})+U_{3}$.
Similarly, we can implement an operation proportional to
$\sum_{j=1}^{3}a_{j}S_{1}^{k_{j}}(t/k_{j})$ by choosing the appropriate
coefficients $\alpha_{1},\alpha_{2}\geq 0$, shown in Fig. 5. We count the
number of CNOT gates for the Ising model described in Sec. 4.4 required by
each approach. The first-order Trotter formula implementing a five qubit Ising
model requires four $R_{ZZ}$ gates and five $R_{X}$ gates. Thus, the deepest
circuit of a classical combination will contain $4k_{3}$ $R_{ZZ}$ gates. The
number of different gates required by an LCU implementation is given in Table
2. The table also shows the number of CNOT gates required to implement each
gate, as transpiled by Qiskit into the native gate set.

Table 2: Number of times each gate is required by an MPF with $l=3$ implemented as an LCU or as a classical linear combination. The second column shows the number of CNOT gates required to implement each gate. Here we do not include pulse-efficient transpilation. Here, we do not explicitly include pulse-efficient transpilation but treat the cost of the $R_{Z}Z$ gate as a single CNOT gate. Gate | CNOTs | LCU | Class.
---|---|---|---
$R_{ZZ}$ | 1 | 0 | $4k_{3}$
1-controlled $R_{ZZ}$ | 8 | $4k_{3}$ | 0
1-controlled single-qubit $U$ | 2 |  $4k_{3}$ \+ 2 | 0
2-controlled $R_{ZZ}$ | 20 | $4(k_{1}+k_{2})$ | 0
2-controlled $R_{X}$ | 7 | $4(k_{1}+k_{2})$ | 0

The LCU approach would therefore require
$4(k_{1}+k_{2})\cdot(20+7)+4k_{3}\cdot(8+2)+2\cdot
2=108(k_{1}+k_{2})+40k_{3}+4$ CNOTs and two additional qubits. By contrast,
our approach only requires $4k_{3}$ CNOTs. For $l=3$, this is thus a reduction
of the number of CNOT gates by a factor of
$[108(k_{1}+k_{2})+40k_{3}+4]/4k_{3}$. With, e.g., $\vec{k}=[1,2,7]$ this
corresponds to $22\times$ less CNOT gates than an LCU implementation.

Figure 5: Circuit implementing an MPF with $l=3$ as an LCU. For a five spin
Ising model the LCU method therefore requires seven qubits for $l=3$. Here,
each first-order PF $S_{1}^{k_{j}}$ requires $4k_{j}$ $R_{ZZ}$ gates.

##  Appendix C  Extrapolation properties

A crucial step in our scheme is to compute expectation values with the
$S_{\chi}^{k}$’s in Eq. (5) instead of $M_{l,\chi}$, since it allows us to
directly compute the extrapolated observable from the measured expectation
values without measuring cross-terms in the MPF. Here, cross-terms correspond
to $\bra{\psi}S_{\chi}^{k_{l}}OS_{\chi}^{k_{j}}\ket{\psi}$ for $k_{j}\neq
k_{l}$ when computing the expectation of an observable $O$ with Eq. (5). The
proof included below for completeness follows similar arguments as those from
Refs. [Endo et al.(2019)Endo, Zhao, Li, Benjamin, and Yuan, Vazquez Carrera et
al.(2022)Vazquez Carrera, Hiptmair, and Woerner].

Let $S_{\chi}^{k_{j}}=U+E_{j}$, where $U=e^{-iHt}$ is the matrix exponential
we want to approximate via a PF and
$E_{j}=\sum_{i}\alpha_{\eta_{i}}\frac{t^{\eta_{i}+1}}{k_{j}^{\eta_{i}}}$ is
the error incurred in the approximation, where $\alpha_{k}$ denotes matrix
commutators. Let $O$ denote the observable of interest and $a_{1},...,a_{l}$
the extrapolation weights that satisfy Eq. (6). We will prove that if

$\displaystyle\sum_{j=1}^{l}a_{j}S_{\chi}^{k_{j}}$
$\displaystyle=\sum_{j=1}^{l}a_{j}U+\sum_{j=1}^{l}a_{j}\sum_{i}\alpha_{\eta_{i}}\frac{t^{\eta_{i}+1}}{k_{j}^{\eta_{i}}}$
$\displaystyle=U+\mathcal{O}\left(k_{1}^{-n}\right),$ (15)

where $\eta_{l-1}<n\in\mathbb{Z}$, then

$\displaystyle\sum_{j=1}^{l}a_{j}\bra{\psi}\left(S_{\chi}^{k_{j}}\right)^{\dagger}OS_{\chi}^{k_{j}}\ket{\psi}$
$\displaystyle=\bra{\psi}U^{\dagger}OU\ket{\psi}$
$\displaystyle\phantom{{}=}+\mathcal{O}\left(k_{1}^{-n}\right).$ (16)

Eq. (C) shows that by ignoring cross-terms we still obtain an approximation
with an $\mathcal{O}(k_{1}^{-n})$ error.

With $S_{\chi}^{k_{j}}=U+E_{j}$ the expectation value $\langle
O\rangle=\bra{\psi}\left(S_{\chi}^{k_{j}}\right)^{\dagger}OS_{\chi}^{k_{j}}\ket{\psi}$
is

$\displaystyle\langle O\rangle$
$\displaystyle=\bra{\psi}U^{\dagger}OU\ket{\psi}+\bra{\psi}U^{\dagger}OE_{j}\ket{\psi}$
$\displaystyle+\bra{\psi}E_{j}^{\dagger}OU\ket{\psi}+\bra{\psi}E_{j}^{\dagger}OE_{j}\ket{\psi}.$
(17)

Next, the first order in $E_{j}$ error terms is

$\displaystyle\epsilon$
$\displaystyle=\sum_{j=1}^{l}a_{j}\bra{\psi}U^{\dagger}OE_{j}\ket{\psi}$
$\displaystyle=\sum_{j=1}^{l}a_{j}\sum_{i}\bra{\psi}U^{\dagger}O\alpha_{\eta_{i}}\frac{t^{\eta_{i}+1}}{k_{j}^{\eta_{i}}}\ket{\psi}$
$\displaystyle=\sum_{i}\bra{\psi}U^{\dagger}O\alpha_{\eta_{i}}t^{\eta_{i}+1}\ket{\psi}\sum_{j=1}^{l}\frac{a_{j}}{k_{j}^{\eta_{i}}}$
$\displaystyle=\mathcal{O}\left(k_{1}^{-n}\right),$ (18)

where the last equality comes from combining Eq. (6), i.e. the first $l-1$
error terms cancel out, and Eq. (C), i.e.
$\alpha_{n}\frac{t^{n+1}}{k_{j}^{n}}=\mathcal{O}(k_{1}^{-n})$. A similar
result can be shown for $\bra{\psi}E_{j}^{\dagger}OU\ket{\psi}$ and
$\bra{\psi}E_{j}^{\dagger}OE_{j}\ket{\psi}$. Thus, Eq. (C) holds and
extrapolation can be applied directly to the measured observables.

##  Appendix D Sampling noise of a Bernoulli variable

We illustrate the asymptotic behaviour of the error term
$\mathcal{O}\left(\norm{\vec{a}}_{1}\varepsilon^{\prime}\right)$ from Eq. (9)
restated below

$\displaystyle\sum_{j=1}^{l}a_{j}\left(E_{j}+\varepsilon(k_{j})+\varepsilon^{\prime}\right)$
$\displaystyle=E+\mathcal{O}\left(\frac{t^{2(l+\chi)+1}}{k_{1}^{2(l+\chi)}}\right)+\mathcal{O}\left(\norm{\vec{a}}_{1}\varepsilon^{\prime}\right).$

We consider as example the task of sampling a Bernoulli discrete random
variable $\mathcal{X}$ and computing the probability $p$ that $\mathcal{X}=1$
by taking $M$ samples. We thus have $E_{j}=E=\mathbb{E}[\mathcal{X}]$ and
$\varepsilon(k_{j})=0$ for all $j$. The number of samples allows us to control
the sampling noise as $\varepsilon^{\prime}(M)=1/\sqrt{M}$ on average. We use
the extrapolation weights $\vec{a}$ of an ill-conditioned MPF built from the
PF $S_{2}$ and by taking $k_{j}=j$ as the sequence of Trotter exponents.
Indeed, this setting uniquely defines the sequences $\vec{a}$ by Eq. (6) whose
1-norm scales exponentially with the number of extrapolation points $l$, i.e.,
$\norm{\vec{a}}_{1}\in e^{\Omega(l)}$ [Hao Low et al.(2019)Hao Low,
Kliuchnikov, and Wiebe], see Fig. 6(a). We expect the sampling error to behave
as

$\sum_{j=1}^{l}a_{j}\left(E+\frac{1}{\sqrt{M}}\right)=E+\mathcal{O}\left(\frac{\norm{\vec{a}}_{1}}{\sqrt{M}}\right),$
(19)

i.e. ill-conditioned sequences have an error that grows exponentially with
$l$. Adding extrapolation points therefore increases the error, see Fig. 6(b).
Here, each line corresponds to a different $l$ and we estimate
$E_{j}(M)\coloneqq E_{j}+\varepsilon^{\prime}(M)$ for $j=1,\cdots,l$ by
sampling $M$ times from $\mathcal{X}$. Thus, each point in the plot
corresponds to the value

$E^{l}(M)\coloneqq\sum_{j=1}^{l}a_{j}E_{j}(M).$ (20)

This example illustrates how a non-careful choice of Trotter exponents
exponentially increases the sampling error of a quantum algorithm. Note that
the example above can be seen as a weighted average of $l$ times $M$ samples
with weights $\vec{a}$. A simple average is recovered with $a_{j}=l^{-1}$.

Figure 6:  Illustration of an ill-conditioned MPF. (a) Condition number of
ill-conditioned MPFs with base PFs $S_{1}$ and $S_{2}$ when taking Trotter
exponents $k_{j}=j$. (b) Error $\absolutevalue{E-E^{l}(M)}$ obtained from
sampling a Bernoulli variable with $p=0.3$ using Eq. (20). Different values of
$l$ are plotted as different lines.

##  Appendix E Optimal Trotter exponents

We obtain sequences of Trotter exponents by solving an optimization problem in
Python that minimizes $\norm{\vec{a}}_{1}$. We set Eq. (6) as the constraints
for the optimization problem. Then, from all the possible sequences of Trotter
exponents we accept those for which the condition number returned by the
optimization algorithm is below a set threshold. All the results presented
here satisfy $\norm{\vec{a}}_{1}\leq 3$. The Trotter exponents and
extrapolation weights used in the Ising simulation experiments are given in
Table 3. They were obtained with DOCPLEX [The IBM Decision Optimization on
Cloud team(2022)] using the code shown below.

⬇

1from docplex.mp.model import Model

2

3# npts := number of extrapolation points

4num_rows = npts - 1

5h = np.ones((num_rows, npts))

6# 2nd order

7h[0, :] = [1/f/f for f in factors]

8# 1st order

9h[0, :] = [1/f for f in factors]

10for i in range(1, num_rows):

11 h[i, :] = h[0, :]**(i+1)

12

13# create empty model

14mdl = Model()

15

16# create decision variables

17a = mdl.continuous_var_list(

18 [f’a{i}’ for i in range(npts)],

19 lb=-np.inf,

20 )

21x = mdl.continuous_var_list(

22 [f’x{i}’ for i in range(npts)],

23 )

24

25# minimize sum_i x[i]

26mdl.minimize(mdl.sum(x))

27

28# x[i] >= a[i] and x[i] >= -a[i]

29for i in range(num_points):

30 mdl.add_constraint(x[i] >= a[i])

31 mdl.add_constraint(x[i] >= -a[i])

32

33# sum(a) == 1

34mdl.add_constraint(mdl.sum(a) == 1)

35

36# for each i: sum_j a_i*h_ij = 0

37for i in range(num_rows):

38 row = [a[j] * h[i, j]

39 for j in range(npts)]

40 row_sum = mdl.sum(row)

41 mdl.add_constraint(row_sum == 0)

Table 3: Trotter exponents and extrapolation weights used to build the $S_{1}$-based MPFs to approximate the Ising Hamiltonian. We define well-conditioned MPFs those that satisfy $\norm{\vec{a}}_{1}\leq 3$. Trotter exponents $\vec{k}$ |  Extrapolation weights $\vec{a}$ | $\norm{\vec{a}}_{1}$
---|---|---
Well-conditioned MPF
$[1,2]$ | $\quad[-1.0,2.0]$ | $\quad 3$
$[1,3]$ | $\quad[-0.5,1.5]$ | $\quad 2$
$[2,4]$ | $\quad[-1.0,2.0]$ | $\quad 3$
$[2,5]$ | $\quad[2/3,5/3]$ | $\quad 2.33$
$[1,2,6]$ | $\quad[0.2,-1,1.8]$ | $\quad 3$
$[1,2,7]$ | $\quad[1/6,-4/5,49/30]$ | $\quad 2.6$
Ill-conditioned MPF
$[6,7]$ | $\quad[-6,7]$ | $\quad 13$
$[3,4,5,6,7]$ | $\quad[27/8,-128/3,375/4,-216,2401/24]$ | $\quad 518.33$
$[1,2,3,4,5,6,7]$ | $\quad[1/720,-8/15,243/16,-1024/9,15625/48,-3888/10,117649/720]$ | $\quad 1007.22$

##  Appendix F Quantum circuit scaling of the spin-boson model

Table 4:  Total number of repetitions of the different circuits in the deepest Trotter circuit needed to simulate spin-boson systems of different sizes $(M,n^{\max})$ to achieve the target accuracy $\varepsilon_{t}\in\\{10^{-2},10^{-3},10^{-4}\\}$. Here, $M$ is the number of bosons and $n^{\max}$ is the maximum occupation number in each bosonic mode. We refer to Ref. [Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli] for details on the quantum implementation. However, compared to Ref. [Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli] we count the number of different circuits $e^{-iH_{j}t}$ instead of the number of first-order Trotter layers. For example, if $H=H_{1}+H_{2}$ we count $e^{-iH_{1}t}e^{-iH_{2}t}$ as two repetitions while Miessen _et al._ count it as one. The number in parenthesis in the $S_{2}$ and $S_{4}$ columns represent $l$ in the MPF, i.e. in column 4, row 7, $13(2)$ means that an MPF with $S_{2}$ requires two circuits with the deepest circuit repeating $e^{-iH_{j}t}$ $13$ times to achieve an error below $10^{-2}$. The values for $S_{1}$ and VQA are taken from Ref. [Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli]. |  | $\quad\qquad\varepsilon_{t}=10^{-2}$ | $\qquad\qquad\varepsilon_{t}=10^{-3}$ | $\qquad\qquad\varepsilon_{t}=10^{-4}$ |
---|---|---|---|---|---
$N_{q}$ |  ($M$, $n^{\max}$) | $\quad\quad S_{1}$ | $\quad\quad S_{2}$ | $\quad\quad S_{4}$ | $\quad\quad\quad S_{1}$ | $\quad\quad S_{2}$ | $\quad\quad S_{4}$ | $\quad\quad\quad S_{1}$ | $\quad\quad S_{2}$ | $\quad\quad S_{4}$ |  VQA
3 | (1,1) | 48 | 11(2) | 31(2) | 170 | 23(2) | 41(2) | 552 | 43(2) | 81(2) | 2
4 | (1,2) | 69 | 21(2) | 61(3) | 240 | 29(3) | 101(2) | 777 | 45(3) | 121(2) | 3
5 | (1,3) | 69 | 25(3) | 61(3) | 219 | 41(2) | 101(2) | 690 | 113(1) | 101(2) | 6
(2,1) | 54 | 13(2) | 51(2) | 196 | 23(2) | 71(3) | 646 | 57(2) | 81(3) | 4
7 | (2,2) | 102 | 13(2) | 121(2) | 333 | 77(2) | 181(2) | 1068 | 89(2) | 201(2) | 9
(3,1) | 64 | 21(2) | 51(2) | 228 | 21(2) | 51(2) | 750 | 21(2) | 181(3) | 6
11 | (2,4) | 117 | 41(3) | 101(4) | 372 | 81(1) | 161(2) | 1167 | 145(1) | 401(1) | 12
(5,1) | 96 | 13(2) | 71(2) | 314 | 49(2) | 91(2) | 1008 | 135(1) | 211(3) | 10
Table 5:  Total number of repetitions of the different circuits in the deepest Trotter circuit needed to simulate spin-boson systems of different sizes $(M,n^{\max})$ to achieve the target accuracy $\varepsilon_{t}\in\\{10^{-2},10^{-3},10^{-4}\\}$. The setting is the same as in Table 4 but here we compare only different PFs. |  | $\quad\qquad\varepsilon_{t}=10^{-2}$ | $\qquad\qquad\varepsilon_{t}=10^{-3}$ | $\qquad\qquad\varepsilon_{t}=10^{-4}$
---|---|---|---|---
$N_{q}$ |  ($M$, $n^{\max}$) | $\quad\quad S_{1}$ | $\quad\quad S_{2}$ | $\quad\quad S_{4}$ | $\quad\quad\quad S_{1}$ | $\quad\quad S_{2}$ | $\quad\quad S_{4}$ | $\quad\quad\quad S_{1}$ | $\quad\quad S_{2}$ | $\quad\quad S_{4}$
3 | (1,1) | 48 | 27 | 111 | 170 | 45 | 121 | 552 | 79 | 171
4 | (1,2) | 69 | 37 | 221 | 240 | 93 | 241 | 777 | 165 | 261
5 | (1,3) | 69 | 37 | 221 | 219 | 65 | 121 | 690 | 113 | 261
(2,1) | 54 | 37 | 111 | 196 | 63 | 121 | 646 | 111 | 171
7 | (2,2) | 102 | 61 | 161 | 333 | 109 | 221 | 1068 | 189 | 361
(3,1) | 64 | 37 | 81 | 228 | 63 | 121 | 750 | 111 | 201
11 | (2,4) | 117 | 57 | 221 | 372 | 81 | 241 | 1167 | 145 | 401
(5,1) | 96 | 45 | 101 | 314 | 77 | 141 | 1008 | 135 | 241

Mapping the spin-boson Hamiltonian to a qubit-lattice shows that the
approximation error of PFs grows linearly with the number of qubits $N_{q}$ as
$\mathcal{O}\left(N_{q}t^{2}/k_{j}\right)$ [Chiesa et al.(2019)Chiesa,
Tacchino, Grossi, Santini, Tavernelli, Gerace, and Carretta, Childs and
Su(2019), Miessen et al.(2021)Miessen, Ollitrault, and Tavernelli]. Now, we
calculate for the spin-boson model the scaling of the deepest circuit needed
by an MPF that classically combines expectation values. We start from the
Taylor expansion of a symmetric PF

$S_{2\chi}^{k_{j}}\left({\textstyle\frac{t}{k_{j}}}\right)=e^{-iHt}+\sum_{n=\chi}^{\infty}A_{2n}t^{2n+1}/k_{j}^{2n},$
(21)

where the $A_{n}$ are matrices that depend on the commutators of the $H_{j}$
[Chin(2010), Hao Low et al.(2019)Hao Low, Kliuchnikov, and Wiebe]. Using this
notation we bound the MPF error by
$\mathcal{O}\left(\frac{\norm{A_{2l+1}}t^{2l+1}}{(2l+1)!}\right)$ [Hao Low et
al.(2019)Hao Low, Kliuchnikov, and Wiebe, Vazquez Carrera et al.(2022)Vazquez
Carrera, Hiptmair, and Woerner]. Ref. [Childs and Su(2019)] shows that
$\norm{A_{n}}\leq N_{q}$ by the triangle inequality and by grouping the
Hamiltonian terms $H_{j}$ in an even-odd pattern. We therefore update the MPF
error bound to $\mathcal{O}\left(N_{q}/(2l+1)!\right)$ for which we also
assumed $t\leq 1$. Thus, for a target error tolerance $\varepsilon_{t}$ the
asymptotic scaling of $l$ satisfies the inequality

$\frac{N_{q}}{(2l+1)!}<\varepsilon_{t}.$ (22)

For simplicity, we write $y\coloneqq N_{q}/\varepsilon_{t}$ and
$x\coloneqq(2l+1)$ such that Eq. (22) reads $x!>y$. By Stirling’s
approximation, i.e. $x!\sim\sqrt{2\pi x}(x/e)^{x}$, the inequality $x!>y$
holds if we set $y=(x/e)^{x}$. Rearranging this equation yields

$\frac{1}{x}y^{\frac{1}{x}}=e^{-1},$ (23)

which has the solution

$2l+1=x=\frac{\ln y}{W\left(\frac{\ln y}{e}\right)},$ (24)

where $W(\cdot)$ denotes the Lambert function [Corless et al.(1996)Corless,
Gonnet, Hare, Jeffrey, and Knuth]. For $z>e$ it holds that $\ln z-\ln\ln
z<W(z)<\ln z$ [Hoorfar and Hassani(2008)]. We can thus substitute $W(e^{-1}\ln
y)$ by $\ln(e^{-1}\ln y)-\ln\ln(e^{-1}\ln y)=\mathcal{O}(\ln\ln y)$ to find
that the number of extrapolation points $l$ scales as

$\displaystyle\mathcal{O}\left(\frac{\log(N_{q}/\varepsilon_{t})}{\log\log(N_{q}/\varepsilon_{t})}\right).$
(25)

Equation (25) holds for $t\leq 1$. For an arbitrary time, similarly to the
main text, we rescale $\vec{k}$ by $\alpha=t$ to ensure that $t/k_{j}\leq 1$.
Furthermore, for $t<1$ the number of Trotter steps of the deepest circuit
$k_{l}$ is bounded by $l^{2}$ for the sequences in Ref. [Hao Low et
al.(2019)Hao Low, Kliuchnikov, and Wiebe] which becomes $k_{l}\leq\alpha
l^{2}$ when we rescale $\vec{k}$ to accommodate for $t>1$. We thus conclude
that $k_{l}$ scales as

$\mathcal{O}\left(\frac{t\log^{2}(N_{q}/\varepsilon_{t})}{\log^{2}\log(N_{q}/\varepsilon_{t})}\right)$
(26)

as presented in the main text. We further emphasize that Eq. (26) is a worst
case upper bound which may yet be improved, for instance, by better choices of
$\alpha$.

##  Appendix G Dynamical Decoupling

Dynamical decoupling (DD) applies carefully chosen gate sequences on idle
qubits to suppress, e.g., decoherence [Pokharel et al.(2018)Pokharel, Anand,
Fortman, and Lidar] and cross-talk [Tripathi et al.(2022)Tripathi, Chen,
Khezri, Yip, Levenson-Falk, and Lidar]. These sequences are chosen to cancel
system-environment interactions to a given order in time-dependent
perturbation theory [Pokharel et al.(2018)Pokharel, Anand, Fortman, and Lidar,
Lidar(2014)]. Since, in practice, the duration of the two-qubit cross-
resonance gates on IBM hardware varies across qubit pairs and because the
circuit structure has layers of one- and two-qubit gates there are often
periods where some qubits idle and therefore suffer decoherence and energy
relaxation. The simplest DD sequence replaces sufficiently long idle periods
of length $\tau$ with the gate sequence $\frac{\tau^{\prime}}{4}$ – $X$ –
$\frac{\tau^{\prime}}{2}$ – $X$ – $\frac{\tau^{\prime}}{4}$, as done in Ref.
[Jurcevic et al.(2021)Jurcevic, Javadi-Abhari, Bishop, Lauer, Bogorin, Brink,
Capelluto, Günlük, Itoko, Kanazawa, and _et al._, Chen et al.(2022)Chen,
Burdick, Yao, Orth, and Iadecola]. Here, the idle time $\tau^{\prime}$ is the
idle time without DD less the duration of the $X$ gate. We did not see a
significant impact of DD on the hardware results and therefore we did not use
it. We conjecture this is because single-qubit gates have a $4\cdot 10^{-4}$
error and that the idle times are small. Indeed, the duration of the average
delay which can accommodate the DD sequence in the circuit without the added
ZNE gates is $252~{}{\rm ns}$. Under a $T_{1}$ time of $125~{}\mu{\rm s}$, the
average $T_{1}$ across _ibmq_montreal_ , this delay can be seen as an error of
$1-e^{-0.292/125}\approx 2\cdot 10^{-3}$, i.e. comparable to twice the single-
qubit gate error. If the idle times in the circuit were longer we believe that
it would have been advantageous to use DD.

##  Appendix H Calibration of $R_{ZX}$ pulses

Computing the dynamics of a quantum system with a Trotter expansion requires
many high accuracy small rotations. To avoid decoherence the $R_{ZX}$ gates
are built from scaled cross-resonance gates following the linear methodology
in Ref. [Stenger et al.(2021)Stenger, Bronn, Egger, and Pekker, Earnest et
al.(2021)Earnest, Tornow, and Egger]. Scaling calibrated CNOT gates to
implement $R_{ZX}(\theta)$ with small values of $\theta$, i.e. up to $\sim
50~{}{\rm mrad}$, results in Gaussian cross-resonance pulses without a flat-
top and a $\theta$-dependent amplitude. Since the rotation angle $\theta$
implemented by the cross-resonance drive depends non-linearly on the pulse
amplitude [Magesan and Gambetta(2020), Alexander et al.(2020)Alexander,
Kanazawa, Egger, Capelluto, Wood, Javadi-Abhari, and McKay, Earnest et
al.(2021)Earnest, Tornow, and Egger], we mitigate rotation errors with a fine
amplitude calibration experiment [Tornow et al.(2022)Tornow, Kanazawa, Shanks,
and Egger]. Here, the gate sequence $[R_{ZX}(\theta)]^{n}\ket{00}$ is repeated
for a variable $n$ and a desired $\theta$ after which the target qubit is
measured and both qubits are reset, see Fig. 7(a). The underlying pulse
sequences, implemented by linearly scaling the cross-resonance amplitude,
creates an effective rotation $[R_{ZX}({\rm d}\theta+\theta)]^{n}$ in the
measured population, see Fig. 7(b) which we fit to the function

$\displaystyle y(n)=\frac{A}{2}e^{-n/\tau}\cos\left(2\pi n[{\rm
d}\theta+\theta]+\pi\right)+B.$ (27)

We chose the maximum $n$ as $3/\theta$ to have approximately three full
oscillations. The small rotation angles we are dealing with thus require pulse
sequences with a duration which can exceed the $T_{1}$ times of the qubits. To
produce good fits, as measured by $\chi^{2}$, we add the damping factor
$e^{-n/\tau}$ in Eq. (27). We run the fine amplitude calibration for all qubit
pairs and 20 values of $\theta$ linearly spaced between 0 and the largest
$\theta$ for which the scaled CR gates have no flat-top. The measured ${\rm
d}\theta$ yields a correction factor $\theta/({\rm d}\theta+\theta)$ that we
multiply with the amplitude of the scaled CR gates, see Fig. 8. We observe
that this amplitude calibration improves the value of the measured expectation
values, see Fig. 9.

Figure 7: Fine amplitude calibration of a $R_{ZX}$ gate with a small rotation
angle. (a) Pulse schedule of a single $R_{ZX}$ gate. (b) Gate sequence
designed to amplify rotation errors. (c) Population of the target qubit as a
function of the number of gates fit to a damped cosinusoidal oscillation.
Figure 8: Small $R_{ZX}$ rotation angle amplitude correction factors
$\theta/({\rm d}\theta+\theta)$ obtained with fine amplitude calibrations on
_ibmq_montreal_ on May 19${}^{\text{th}}$ 2022. Figure 9: Comparison of the
local ZNE magnetization $\langle Z_{i}\rangle$ using amplitude calibrated
(blue markers) versus uncalibrated (red markers) $R_{ZX}$ gates. That data is
for the five-spin Ising model described in Sec. 4.4 simulated on
_ibmq_montreal_. Pannels (a) and (b) correspond to qubit 16 and 19,
respectively. Each point is the average of ten different runs and the error
bars show twice the standard deviation.

##  Appendix I Comparison with Mitiq

Ref. [Chen et al.(2022)Chen, Burdick, Yao, Orth, and Iadecola] simulates an
Ising model using a similar approach to error mitigation, that is, Pauli
Twirling, a pulse efficient transpilation of the $R_{ZZ}$ gates and dynamical
decoupling. By contrast, we do not apply dynamical decoupling and implement
ZNE by linearly scaling cross-resonance pulses. Ref. [Chen et al.(2022)Chen,
Burdick, Yao, Orth, and Iadecola] uses the Mitiq package [LaRose et
al.(2022)LaRose, M., Kaiser, Karalekas, Alves, Czarnik, Mandouh, Gordon,
Hindy, Robertson, T., Shammah, and Zeng] to implement random gate folding,
i.e. repeating random layers of $R_{ZZ}$ gates. Here, we use the example
described in Sec. 4.4 as a benchmark to compare both approaches in Fig. 10. We
observe that ZNE implemented by folding gates improves the data with respect
to the unmitigated results in some cases. Instead, implementing ZNE by scaling
CR pulses improves the results in almost all cases.

Figure 10: Comparison of two ZNE methodologies on _ibmq_montreal_. The top row
compares simulation, hardware with no error mitigation (hardware: no EM) and
hardware with readout error mitigation, PT, DD and ZNE as implemented by Mitiq
using the code publicly available from Ref. [Yongxin(2022)]. In the bottom row
the hardware error mitigated results were obtained using readout error
mitigation, PT and the ZNE approach described in Sec. 4.3. In all cases the
$R_{ZZ}$ gates are obtained by pulse-efficient transpilation as described in
Sec. 4.2.

##  Appendix J  Hardware properties

Figure 11: Layout of _ibmq_montreal_. We choose five qubits in a line such
that the fidelity decay resembles the $T_{1}$ decay. We chose to work with
qubits $\\{13,14,16,19,20\\}$. Figure 12: Fidelity of
$R_{ZX}(\theta)R_{ZX}(-\theta)$ on _ibmq_montreal_ on May 29${}^{\text{th}}$
2022 for different rotation angles $\theta\in(0.5,21)$. Each panel corresponds
to a different qubit pair $(q_{1},q_{2})$. The dots represent the fidelity of
the gates, while the solid lines the $T_{1}$ decay of the qubits. Here we
chose to work with qubits $\\{13,14,16,19,20\\}$ and the corresponding panels
are highlighted in blue.

Error mitigation techniques can significantly decrease the physical noise in
the measured observables. However, their efficacy is limited and error rates
are mostly hardware related. In particular, choosing _good_ qubits is a simple
but worthwhile strategy towards achieving accurate results. Here _good_ qubits
refers to low error rates in the readout, single- and two-qubit gates. Most of
these values are provided by the IBM Quantum backends [ibm()]. However, our
circuits comprise mostly $R_{ZX}(\theta)$ gates, which are the main source of
noise. The fidelity of these gates is not provided by the backends as they are
built from scaled-CR gates [Earnest et al.(2021)Earnest, Tornow, and Egger].
Therefore, we perform quantum process tomography of scaled-CR
$R_{ZX}(\theta)R_{ZX}(-\theta)$ gates using Qiskit Experiments [Kanazawa et
al.(2023)Kanazawa, Egger, Ben-Haim, Zhang, Shanks, Aleksandrowicz, and Wood]
on _ibmq_montreal_ , see Fig. 11.

For each qubit pair in the device we calculate the fidelity of the gate
sequence for 20 values of $\theta\in(0.5,20)$, shown as dots in Fig. 12, and
plot it together with the $T_{1}$ decay of the two qubits, represented as
solid lines. We then inspect the plots and define as _good_ qubit pairs those
for which the gate fidelity decay resembles the $T_{1}$ decay. Finally, we
select 4 good pairs that form five linearly connected qubits, depicted in Fig.
11. The properties of _ibmq_montreal_ are shown in Table 6.

Table 6: Properties of the qubits and gates used to simulate the Ising Hamiltonian as reported by _ibmq_montreal_ on May 29${}^{\text{th}}$ 2022\. The $T_{1}$ and $T_{2}$ times are reported in $\mu s$. The CNOT error in each row is the error between the qubit of that row and the qubit of the next row. |  |  | Error rate $(\%)$ of
---|---|---|---
Qubit | $T_{1}$ | $T_{2}$ | $\sqrt{X}$ | CNOT | Measure
13 | 118 | 122 | 0.018 | 0.57 | 0.74
14 | 133 | 256 | 0.019 | 0.79 | 1.09
16 | 82 | 22 | 0.052 | 1.68 | 1.83
19 | 160 | 167 | 0.023 | 0.86 | 1.21
20 | 85 | 133 | 0.035 | – | 1.50
average | 115.6 | 140 | 0.029 | 0.97 | 1.27