# Composite QDrift-Product Formulas for Quantum and Classical Simulations in
Real and Imaginary Time

Matthew Pocrnic  Department of Physics, 60 Saint George St., University of
Toronto, Toronto, Ontario, M5S 1A7, Canada  Matthew Hagan  Department of
Physics, 60 Saint George St., University of Toronto, Toronto, Ontario, M5S
1A7, Canada  Juan Carrasquilla  Department of Physics, 60 Saint George St.,
University of Toronto, Toronto, Ontario, M5S 1A7, Canada  Vector Institute,
MaRS Centre, Toronto, Ontario, M5G 1M1, Canada  Department of Physics and
Astronomy, University of Waterloo, Ontario, N2L 3G1, Canada  Dvira Segal
Department of Chemistry and Centre for Quantum Information and Quantum
Control, University of Toronto, 80 Saint George St., Toronto, Ontario, M5S
3H6, Canada  Department of Physics, 60 Saint George St., University of
Toronto, Toronto, Ontario, M5S 1A7, Canada  Nathan Wiebe  Department of
Computer Science, University of Toronto, Toronto ON, Canada  Pacific Northwest
National Laboratory, Richland Wa, USA  Canadian Institute for Advanced
Research, Toronto ON, Canada

(June 27, 2023)

###### Abstract

Recent work has shown that it can be advantageous to implement a composite
channel that partitions the Hamiltonian $H$ for a given simulation problem
into subsets $A$ and $B$ such that $H=A+B$, where the terms in $A$ are
simulated with a Trotter-Suzuki channel and the $B$ terms are randomly sampled
via the QDrift algorithm. Here we show that this approach holds in imaginary
time, making it a candidate classical algorithm for quantum Monte-Carlo
calculations. We upper-bound the induced Schatten-$1\to 1$ norm on both
imaginary-time QDrift and Composite channels. Another recent result
demonstrated that simulations of Hamiltonians containing geometrically-local
interactions for systems defined on a finite lattices can be improved by
decomposing $H$ into subsets that contain only terms supported on that subset
of the lattice using a Lieb-Robinson argument. Here, we provide a quantum
algorithm by unifying this result with the composite approach into “local
composite channels” and we upper bound the diamond distance. We provide exact
numerical simulations of algorithmic cost by counting the number of gates of
the form $e^{-iH_{j}t}$ and $e^{-H_{j}\beta}$ to meet a certain error
tolerance $\epsilon$. We show constant factor advantages for a variety of
interesting Hamiltonians, the maximum of which is a $\approx 20$ fold speedup
that occurs for a simulation of Jellium.

##  I Introduction

Quantum simulations of quantum systems, suggested by Feynman in 1982 [1],
provides what is perhaps the most natural application for quantum computers
with expected exponential advantages over the best classical algorithms [2].
The quantum simulation problem can be summarized as follows: Given a
Hamiltonian $H$, a time $t$, and an error tolerance $\epsilon$, we wish to
implement, on a quantum computer, an approximate function $f$ to the time
evolution operator $e^{-iHt}$ such that
$\norm{e^{-iHt}-f(iHt)}_{\kappa}\leq\epsilon$ for some distance norm $\kappa$.
The significance of this operator comes from quantum mechanics where
$e^{-iHt}$ is the solution to the Schrödinger Equation with a time-independent
Hamiltonian: $i\partial_{t}\ket{\Psi}=H\ket{\Psi}$ in units where $\hbar=1$.
Motivation for quantum simulation is drawn from the fact that it is a
necessary subroutine for phase estimation, which allows one to learn the
eigenvalues of $H$ [3, 4]. This of particular interest in quantum chemistry
and materials [5, 6, 7] where Hamiltonians are large and analytic solutions to
the eigenvalue problem are not known. In addition, quantum simulation also
allows for the study of dynamical quantum observables [8]. Classically,
implementing the time evolution operator is a difficult problem due to how
rapidly the matrices grow with respect to the size of the system of interest.
As well, quantum effects such as interference are hard to simulate
classically. This problem has received extensive attention over the last few
decades, and many different solutions have been proposed. These solutions can
in some sense be broken down into two families; those built on product
formulas including Trotter-Suzuki [2, 9, 10] formulas and QDrift [11, 12, 13,
14] and those implementing more intricate Hamiltonian transformations such as
Linear Combinations of Unitaries [15, 16], and Qubitization [17, 18, 19]. Put
simply, if we have a Hamiltonian that can be expressed as a sum of terms
$H=\sum_{j}H_{j}$ with a time evolution operator $e^{-it\sum_{j}H_{j}}$, a
product formula is an approximation to $e^{-iHt}\approx\prod_{j}e^{-iH_{j}t}$
which is not exact in general due to non-zero commutators in the Taylor
series. In the absence of fault-tolerant quantum computers, powerful methods
have also been discovered to study the properties of $H$ on classical
computers, particularly in the field of Quantum Monte-Carlo [20, 21]. Here,
the propagator of interest is often the time evolution operator which has been
Wick transformed to imaginary time $e^{-\beta H}:it\rightarrow\beta$. This is
implemented both as a projection method to study the ground state of $H$, and
to study the partition function $Z=\Tr e^{-\beta H}$ from statistical physics
where $\beta$ is the inverse temperature of the system. Trotter-Suzuki
formulas found some of their first applications in Monte Carlo calculations
prior to becoming a popular ancilla free quantum simulation algorithm. In this
paper, we take the reverse approach and consider how Composite Trotter-QDrift
simulation strategies [22] from quantum computing carry over to imaginary
time. We also extend these composite strategies to account for locality of $H$
based on Lieb-Robinson bounds [23] in terms of their application for quantum
simulation (in real time).

Within this manuscript, we provide three contributions. First, we show that
QDrift and Composite QDrift-Trotter algorithms [22] can be extended to
imaginary time with the same asymptotic scaling as in real time. The basis of
these algorithms is a partitioning of the Hamiltonian such that $H=A+B$, where
we simulate the subset $A$ with Trotter-Suzuki formulas and $B$ with QDrift,
to take advantage of the strengths of each. This promotes them to a candidate
classical algorithm for quantum Monte-Carlo calculations. In doing so, we
provide new bounds on QDrift and Trotter-Suzuki formulas in imaginary time.
Next, we show that Composite algorithms can be applied to geometrically local
lattices, by first decomposing the lattice into local subsets [23], and then
partitioning the Hamiltonian terms supported on said subsets into localized
QDrift and Trotter simulations. In comparison to standard Composite
simulations, the Local Composite simulations recover the form of the original
error bounds with an additional error term that is exponentially small.
Finally, we provide a library to numerically evaluate the performance of
Composite channels in both real and imaginary time, as well as the performance
of the Local Composite quantum algorithm in real time. Here, we also
numerically introduce new optimization schemes for conducting this
partitioning based on heuristics and intuition that comes from the error
bounds of each respective algorithm. Rather than computing bounds, we evaluate
the exact exponential gate cost (the number of exponentials
$e^{-iH_{j}t^{\prime}}$ with time slice $t^{\prime}$ that appear in the
product formula) at simulation time $t$ for implementing an approximation to
the time evolution operator. Equivalent counts are performed in imaginary time
for $e^{-\beta H}$ at inverse temperature $\beta$. This is done for a series
of interesting Hamiltonians which are drawn from quantum chemistry and spin
systems, the former of which are generated using OpenFermion [24]. We find
significant constant factor advantages using our composite approaches, with
the maximum being a factor of 18.8 for a 6-site Jellium Hamiltonian. We make
the observation that composite methods seem to work better for Hamiltonians
containing a large number of terms, and for those with sharply peaked spectral
norm distributions.

The paper is organized as follows: We begin with detailed preliminary material
in Section II that introduces the algorithms in question as well as some
details regarding how they are analyzed analytically and numerically. Here we
also provide some background on the Lieb-Robinson bound and define
geometrically-local Hamiltonians for which they apply. In Section III, we show
how to extend composite channels to imaginary time and we provide bounds on
the error for said channels. In doing so, we show that the QDrift algorithm
holds in imaginary time and that the same $\epsilon$ scaling can be achieved,
suggesting it can also be used effectively for classical computation. We also
provide a quantum algorithm that allows us to include locality as means for
partitioning in the “local composite” channel, and show that doing so
introduces exponentially-small errors. In Section V, we outline the details of
our numerical simulations and describe the general features of our library. We
then provide plots of the gate cost, or the number of exponential operators
needed to meet a chosen error tolerance. We conclude in Sec. VI with a
discussion and consider open questions.

##  II Preliminaries

In this section, we introduce the general physical setup, notation and
simulation approaches that will be used throughout the paper. We then review
and summarize the algorithmic cost of the Trotter-Suzuki and the QDrift
product formulas of real time time evolution, as well as their hybrid
approach, which was introduced in Ref. 22.

###  II.1 General Notation and Assumptions

Before introducing the simulation algorithms, it is important to outline the
notation that will be used throughout the paper, as well as the assumptions
that remain relevant for calculations and proofs. The following assumptions
will remain valid unless otherwise stated. We assume that Hamiltonian
operators $H$ are time-independent, bounded, and that they live in a finite
$2^{n}$ dimensional Hilbert space $\mathscr{H}$. We further assume that these
operators can be written as a finite sum of terms, $H=\sum_{i=1}^{L}H_{i}$
with each $H_{i}\in\mathbb{C}^{{2^{n}\times 2^{n}}}$. Later, in Section IV, we
will discuss algorithms for geometrically-local Hamiltonians in which case
each term in this sum will only be supported on a finite subset of
$\mathscr{H}$; we save those details for later. Note that we can normalize the
terms in the above sum by factoring out their respective spectral norms
(Schatten infinity norm on matrices), $h_{i}=||H_{i}||$ such that
$H=\sum_{i}^{L}h_{i}H_{i}$. We also impose $h_{i}\geq 0$ by pulling any phase
factors into $H_{i}$; this can always be done without loss of generality. The
Schatten p-norms will be a useful concept throughout, and they are defined as
$\norm{H}_{p}=\left(\sum_{i}s_{i}^{p}\right)^{\frac{1}{p}}$, where $s_{i}$ are
the singular values of the matrix $H$. From this definition, we observe that
the Schatten infinity norm coincides with the spectral norm. Other norms will
be discussed throughout the paper. Wherever a norm $\norm{\cdot}$ is lacking a
subscript, it will be assumed to be the Schatten infinity norm. All other
norms will be subscripted appropriately. In what follows, Hamiltonian
operators are taken to be time independent, and they generate dynamics through
the time evolution operator, $U=e^{-iHt}$.

Throughout, we make use of the density matrix formalism to describe mixed
quantum states, with quantum channels being the object that evolve states,
rather than the evolution operator $e^{-iHt}$. By quantum channels we mean a
mapping $\mathcal{A}(\rho)=\sum_{i}E_{i}\rho E_{i}^{\dagger}$ that is
completely positive and preserves the trace, $\Tr(\rho)$. For this to be true,
the operators $E_{i}$ must satisfy the condition
$\sum_{i}E_{i}E_{i}^{\dagger}=\openone$. These conditions are required so that
information regarding the channel operation is preserved and not lost to an
external environment. For example, consider performing a measurement on the
output density matrix $\mathcal{A}(\rho)$. If we find that
$\Tr(\mathcal{A}(\rho))\leq 1$, we conclude that the mapping does not provide
a complete description of the processes involved, since there must be some
other process(es) and therefore measurement outcomes that occur with non-zero
probability [4]. For easy distinction, we reserve calligraphic fonts for
channels $\mathcal{A}$ while operators $A$ are written in the standard way.
Given a density matrix $\rho(t_{0})$ that describes the state at some initial
time $t_{0}$, the time evolution channel is then defined in the following way

$\mathcal{U}{(\rho(t_{0}),t)}:=e^{-iH(t-t_{0})}\rho(t_{0})e^{iH(t-t_{0})}\rightarrow\rho(t).$
(1)

When not specifically defining evolutionary channels, we will often use
$t_{0}=0$ and write $\rho(t=0)=\rho$ to be the initial state. If time
evolution is broken into $m$ time intervals $\\{\Delta t_{1},\Delta
t_{2},...,\Delta t_{m}\\}$,

the state vector evolution looks like $e^{-iH\Delta t_{m}}...e^{-iH\Delta
t_{2}}e^{-iH\Delta t_{1}}\ket{\psi}$. However, in the channel description we
do not take products of channels, rather we apply a channel multiple times in
a channel composition which we denote with $\circ$. For example,
$\mathcal{U}{(t_{i+1})}\circ\mathcal{U}{(t_{i})}(\rho)=e^{-iHt_{i+1}}e^{-iHt_{i}}\rho
e^{iHt_{i}}e^{iHt_{i+1}}$.

Note that in this expression, $t_{i}$ corresponds to the time interval during
which the channel acts.

A significant portion of this paper is dedicated to the investigation of
imaginary-time channels
$\mathcal{U}{(\rho,\beta)}=\frac{1}{\mathscr{N}}e^{-\beta H}\rho e^{-\beta H}$
and the algorithms that closely approximate them. Upon first glance, the
notion of imaginary time seems to violate the trace preserving property of a
quantum channel, so we introduce $\mathscr{N}$ as a normalization factor. To
avoid confusion, we use $\beta$ throughout these discussions instead of $t$,
motivated by statistical mechanics where $\beta$ is considered an inverse
temperature. In these cases, we are not considering a dynamical evolution of
the system in the physical sense, but are rather interested in calculating
properties of the Hamiltonian $H$, such as its partition function $Z=\Tr
e^{-\beta H}$ or its ground state. This is because time evolution of a quantum
system is isomorphic to the cooling of a statistical ensemble via a Wick
rotation. We wish to make clear that all discussions pertaining to algorithms
in imaginary time are not being proposed for digital quantum simulation, but
as classical algorithms for use in calculations like Quantum Monte-Carlo. For
a study of imaginary time evolution on a quantum computer, see [25]

###  II.2 Algorithmic Cost Model

Given that we are dealing with channels containing only product formulas, a
natural cost model is to count the number of exponential gates of the form
$e^{-iH_{j}t^{\prime}}$, required to achieve some error $\epsilon$ in a
simulation. Here, $H_{i}$ is a specific term chosen from a sum of Hamiltonian
terms and $t^{\prime}$ is an appropriately-chosen time slice of the total
simulation time $t$; the following sections on the algorithms of interest
explain how to choose these time slices to meet a desired simulation accuracy
$\epsilon$. Throughout, the aforementioned cost will be denoted $C$ and often
subscripted with the algorithm of interest. This cost model is convenient
analytically as the algorithms we explore contain an exact amount of
iterations or samples of terms that can be easily equated to a count of
exponential gates. This situation is different from common forms of cost in
the literature such as a query complexity, which requires an assumed access to
an oracle that necessarily has an implementation overhead. In this case, the
overhead would come within the implementation of each $e^{-iH_{j}t^{\prime}}$,
which depends on the form of the Hamiltonian terms. In the case of common
Pauli and Fermionic Hamiltonians, circuits are well known and relatively
simple to analyze [26]. For example, it is common in Heisenberg-like models to
have interactions that are a products of Pauli operators such as
$\sigma^{\nu}_{i}\sigma^{\nu}_{i+1}\sigma^{\nu}_{i+2}$. The time evolution of
these operators
$e^{-it(\sigma^{\nu}_{i}\sigma^{\nu}_{i+1}\sigma^{\nu}_{i+2})}$ can be simply
implemented using rotation gates
$R_{\nu}(\theta)=exp\left(-i\frac{\theta}{2}\sigma_{\nu}\right)$. In fact,
using the following circuit we can implement the time evolution for any
product of Pauli operators using just $R_{z}(\theta)$:

$\ket{j_{\idx}}$$R_{z}(\theta)$ Figure 1: A simple quantum circuit showing an
implementation of an exponentiated produce of Pauli spin operators.

This circuit effectively implements
$\exp(-i\frac{\theta}{2}\sigma^{z}_{i}\sigma^{z}_{i+1}\sigma^{z}_{i+2})$, and
we can build a product of any length by continuing this pattern. Switching
this to a product of $\sigma_{x}$ or $\sigma_{y}$ is straightforward given
that $R_{x}=\hat{H}R_{z}\hat{H}$ and $R_{y}=T^{2}R_{z}(T^{2})^{\dagger}$. We
can then adjust the simulation time $t$ through our choice of angle $\theta$.
Fermionic Hamiltonians follow in a similar way, given that we can map them to
a Pauli Hamiltonian using the Jordan-Wigner transformation, which will be
outlined in Section V.1.1. The exampled circuit also reinforces the need for
product formula approximations since while we can simply implement the above
mentioned products of operators, implementing $e^{-it(\sigma^{x}+\sigma^{z})}$
is highly non-trivial, and we approximate it as
$e^{-it\sigma^{x}}e^{-it\sigma^{z}}$ as we detail in the following sections.
In our cost model, this approximation would be counted as 2 gates.
Numerically, counting exponential gates is also a convenient cost model since
it allows for straightforward counting of gates as they are applied.

###  II.3 Trotter-Suzuki Formulas

Given an operator that can be written as a finite sum of terms
$A=\sum_{i}^{L}A_{i}$, the simplest approximation for the corresponding
operator exponential is then
$e^{\sum_{i}^{L}A_{i}}\approx\prod_{i}^{L}e^{A_{i}}$. While this holds as an
equality for exponential functions for variables, it does not in general hold
for operator exponentials since individual terms of $A$ do not necessarily
commute. This can be easily seen by considering the Taylor series of a simple
exponential approximation $e^{(A+B)x}\approx e^{Ax}e^{Bx}$ to second order in
$x$:

$\displaystyle e^{(A+B)x}-e^{Ax}e^{Bx}$
$\displaystyle=\left[\openone+Ax+Bx+\frac{(A+B)^{2}x^{2}}{2!}+...\right]-\left[(\openone+Ax+\frac{A^{2}x^{2}}{2!}+...)(\openone+Bx+\frac{B^{2}x^{2}}{2^{!}}+...)\right]$
(2) $\displaystyle=\frac{(AB+BA)x^{2}}{2}-ABx^{2}+...$ (3)
$\displaystyle=\frac{[B,A]x^{2}}{2}+...$ (4)

An error clearly arises due to a non-vanishing commutator, meaning the
approximation is exact for commuting algebras. However, the nature of this
error opens up the possibility for improving this approximation. If one can
symmetrize this formula, the second order commutator will vanish. In reference
to the previous example, a perhaps obvious symmetrization is
$e^{\frac{Ax}{2}}e^{B}e^{\frac{Ax}{2}}$, which annuls the error term
$\frac{[B,A]x^{2}}{2}$, leaving only error terms in
$\mathcal{O}\left(x^{3}\right)$. More generally, for an operator
$A=\sum_{i}A_{i}$, the 2nd order Trotter-Suzuki decomposition is given by

${S_{2}{\left(Ax\right)}}=\prod_{i}^{L}e^{\frac{A_{i}x}{2}}\prod^{i}_{L}e^{\frac{A_{i}x}{2}}$
(5)

where the flipping of indices in the two products is indicative of reversing
the ordering of products. One can readily imagine that it is then possible to
continue this strategy and build higher order formulas and cancel out
arbitrarily high commutator errors in the product formula approximation.
Indeed, we can build these highly non-trivial decompositions in the following
way [27, 28]:

######  Definition 1 (Trotter-Suzuki Product Formula).

Given a linear operator $A$ acting on a finite dimensional Hilbert space that
can be represented as a finite sum of linear operators $A=\sum_{i}^{L}A_{i}$,
a Trotter-Suzuki product formula approximation ${S_{2k}{\left(Ax\right)}}$ for
$A$ of order $2k$ ($k>1$) can be constructed in the following recursive
fashion:

${S_{2k}{\left(Ax\right)}}:=({S_{2k-2}{\left(Axs_{2k}\right)}})^{2}{S_{2k-2}{\left((1-4s_{2k})Ax\right)}}({S_{2k-2}{\left(Axs_{2k}\right)}})^{2},$
(6)

where $s_{2k}=\frac{1}{4-4^{\frac{2k-1}{2}}}$ .

Hence, in order to obtain a desired product formula of order $2k$, we start
with the second order formula ${S_{2}{\left(Ax\right)}}$ given in Equation (5)
and recursively build the higher order formulas. To reiterate, the advantage
of doing so is that for a product formula of order $2k$, the error of the
approximation is $\in\mathcal{O}\left(t^{2k+1}\right)$ with $t$ here as the
time interval [28].

When the general operator $e^{Ax}$ is replaced by the Hamiltonian time
evolution operator $e^{-iHt}$, where $H=\sum_{i}^{L}H_{i}$, it becomes obvious
how product formulas can be used to evolve quantum states $\ket{\psi}$.
However, we also wish to generalize this to evolve probabilistic mixtures of
states, which are represented by density matrices
$\rho=\sum_{i}p_{i}\outerproduct{\psi_{i}}{\psi_{i}}$. A Trotter-Suzuki
evolution formula for a channel can be written as follows:

######  Definition 2 (Trotter-Suzuki Channel).

Given a Hamiltonian $H$, a density matrix $\rho$, times $t$ and $t_{0}$
($t>t_{0}\geq 0$), and order 2k, then a Trotter-Suzuki channel
$\mathcal{T}^{2k}{(\rho(t_{0}),\Delta t)}$ performs the operation
$\mathcal{T}^{2k}{(\rho(t_{0}),\Delta t)}\rightarrow\rho(t)$ and can be
defined as:

$\mathcal{T}^{2k}{(\rho(t_{0}),t)}:={S_{2k}{\left(-iH\Delta
t\right)}}\rho{S_{2k}{\left(-iH\Delta t\right)}}^{\dagger}$ (7)

Where ${S_{2k}{\left(-iHt\right)}}$ represent the product formulae from
Equation (5) and $\Delta t=t-t_{0}$.

Since the error of this approximation accumulates considerably for long
iteration times, we introduce the iteration parameter $r$, and apply the
channel $r$ times with time intervals $\nicefrac{{t}}{{r}}$. As a result,
longer time simulations become more expensive in the number of iterations $r$
required to manage the error that accumulates to do the simulation time $t$.
The composition of $r$ Trotter-Suzuki channels is straightforward, and
applying $r$ channels each for time $\nicefrac{{t}}{{r}}$ can be written as

$\mathcal{T}^{2k}{(\rho,t)}^{\circ
r}=\prod^{r}_{j=1}{S_{2k}{\left(-iH\frac{t}{r}\right)}}\rho\prod^{r}_{j=1}{S_{2k}{\left(-iH\frac{t}{r}\right)}}^{\dagger}.$
(8)

The cost of this algorithm for first order decomposition is simply $rL$ where
$L$ is the number of Hamiltonian terms. For order $2k$, the cost is $\Upsilon
Lr$ where $\Upsilon=2\times 5^{k-1}$.

Looking at cost formulas for different orders $2k$, we desire information
about $r$, specifically, the minimum value of $r$ required to meet some error
tolerance $\epsilon$ where
$\norm{\mathcal{T}^{2k}{(\rho,\nicefrac{{t}}{{r}})}^{\circ
r}-\mathcal{U}{(\rho,t)}}_{\kappa}\leq\epsilon$ for some norm $\kappa$. In
Ref. 22, the diamond norm is selected, which is the completely bounded trace
norm between channels tensored with the environment. The diamond norm is
defined as follows:
$\left|\left|\mathcal{U}(t)-\mathcal{V}(t)\right|\right|_{\diamond}:=\sup_{\rho:\norm{\rho}\leq
1}\norm{\left(\mathcal{U}(\rho,t)-\mathcal{V}(\rho,t)\right)\otimes\openone}_{1}$.
This is a good analytical norm to choose, given that it maintains the
interpretation of distinguishability between states from the trace distance
(see Section B.1.2). However, it is a bound on channels (maximized over
possible input states $\rho$) and thus yields the interpretation of
distinguishability between quantum channels. This norm is also sub-additive
and sub-multiplicative, which makes it convenient to work with. The procedure
in Ref. 22 upper bounds this norm, sets the value equal to $\epsilon$, and
solves the expression for $r$ to arrive at the following first and $2k$th
order costs

$\displaystyle C_{TS}^{1}(H,t,\epsilon,1)$ $\displaystyle=Lr\leq
L\left\lceil\frac{t^{2}}{2\epsilon}\sum_{i,j}h_{i}h_{j}\norm{[H_{i},H_{j}]}_{\infty}\right\rceil$
(9) $\displaystyle C_{TS}^{2k}(H,t,\epsilon,2k)$ $\displaystyle=\Upsilon
Lr\leq\Upsilon L\left\lceil\frac{(\Upsilon
t)^{1+1/2k}}{\epsilon^{1/2k}}\left(\frac{4\alpha_{comm}(H,2k)}{2k+1}^{1/2k}\right)\right\rceil,$
(10)

where $C$, to reiterate, is the number of operator exponentials that need be
applied rather than an exact gate count where gates are drawn from some
universal gate set. Here,

$\alpha_{comm}(H,2k):=\sum_{\gamma_{i}\in\left\\{1,\ldots
L\right\\}}\left(\prod
h_{\gamma_{i}}\right)\norm{[H_{\gamma_{2k+1}},[H_{\gamma_{2k}},\ldots[H_{\gamma_{2}},H_{\gamma_{1}}]\ldots]}_{\infty}.$

This is a nested commutator error that arises from $2k$ symmetrizations in the
product formula. The main takeaway from this formula is the structure of the
Hamiltonian that the simulation cost depends on. Trotter-Suzuki channels
clearly depend on the number of terms in the Hamiltonian $L$ and on the
commutator structure of the Hamiltonian. The simulation cost does not directly
depend on the magnitude of the terms $\norm{H_{i}}$, only indirectly through
their commutators. This is important to keep in mind when examining the next
algorithm and its cost dependence.

###  II.4 QDrift Random Compiler

The QDrift algorithm was discovered by Campbell [11]. It generates random
product formulas allowing for the evolution of a system to stochastically
drift (through the Hilbert space) towards the true evolution with high
probability. Using the introduced quantum channel notation, the algorithm
works as follows:

######  Definition 3 (QDrift Channel).

Given a Hamiltonian $H=\sum_{i=1}^{L}h_{i}H_{i}$, density matrix $\rho$, times
$t$ and $t_{0}$ ($t>t_{0}\geq 0$), and let $p_{i}=\frac{h_{i}}{\lambda}$ with
$\lambda=\sum_{i=1}^{L}h_{i}$ be a probability distribution from which
Hamiltonian terms are drawn, then a QDrift channel
$\mathcal{Q}{(\rho(t_{0}),\Delta t)}\rightarrow\rho(t)$ can be defined as:

$\mathcal{Q}{(\rho(t_{0}),t)}:=\sum_{i=1}^{L}p_{i}e^{-iH_{i}\lambda\Delta
t}\rho e^{iH_{i}\lambda\Delta t},$ (11)

where $\Delta t=t-t_{0}$.

This channel mixes unitaries with a single sample from the distribution
$p_{i}$. Similarly to the Trotter-Suzuki approach, the accuracy of this
approach improves with the number of samples $N$ that go into our “random
product formula”. To write the channel for multiple samples, it is useful to
think of sampling a vector $\mathbf{j}=\\{j_{1},j_{2},...j_{N}\\}$ of length
$N$, i.i.d from the distribution
$P_{\mathbf{j}}=\lambda^{-N}\prod_{k=1}^{N}h_{j_{k}}$. This corresponds to
implementing the following unitary
$\mathbf{V_{j}}=\prod^{N}_{k=1}e^{-iH_{j_{k}}\tau}$ with
$\tau=\nicefrac{{t\lambda}}{{N}}$. With this in mind, we can neatly write down
the expression for the channel with arbitrary samples $N$ as

$\mathcal{Q}{(\rho,t)}^{\circ
N}=\sum_{\mathbf{j}}P_{\mathbf{j}}\mathbf{V_{j}}\rho\mathbf{V_{j}}^{\dagger},$
(12)

where the sum is performed over all possible vectors $\mathbf{j}$ of length
$N$. In terms of preparing these gates with a quantum circuit, this algorithm
can be thought of as a linear combination of unitaries under classical
control. Similarly to the Trotter-Suzuki formula, we are also interested in
the algorithmic cost of QDrift. A diamond distance upper bound is provided in
Ref. 11, which is then used via the same procedure to provide a cost function
in Ref. 22, by using $N=4\lambda^{2}t^{2}/\epsilon$ and restricting
$\epsilon\in(0,\lambda t\ln(2)/2)$,

$C_{QD}(H,t,\epsilon)=\frac{4\lambda^{2}t^{2}}{\epsilon}.$ (13)

The important take home message here is the difference in the structure of the
algorithmic cost between the Trotter-Suzuki formula [Eq. (9)] and QDrift [Eq.
(13)]. In QDrift, the cost does not depend on the number of Hamiltonian terms
$L$, or the commutator structure. In contrast, it depends on the size of the
terms in the sum of the spectral norms $\lambda$. It is important to keep this
in mind when constructing composite channels.

###  II.5 Composite Simulation Formulas

Equipped with the Trotter-Suzuki and QDrift channel formulas, we wish to
hybridize these approaches to form a composite time evolution channel. The
eventual goal is to do so in such a way that we can minimize the “side
effects” of each algorithm while taking advantage of the strengths of each
approach. Following Ref. [22], we define the composite channel as follows. We
consider the Hamiltonian $H=A+B$, where $A=\sum_{j=1}^{L_{A}}A_{j}$ and
$B=\sum_{j=1}^{L_{B}}B_{j}$. Other terms of interest, such as the number of
summands in a set $L$ are written with a subscript $A$ to indicate that they
belong to a set of Hamiltonian terms that will be simulated by a Trotter-
Suzuki channel, and the same goes for the $B$ terms, for a QDrift channel. The
composite channel takes the following form

$\mathcal{X}^{2k}(\rho,t)=\mathcal{Q}_{B}{(t)}\circ\mathcal{T}_{A}^{2k}{(t)}(\rho),$
(14)

where the channel subscripts indicate the subset of the Hamiltonian that it is
simulating. Note that there is substantial freedom in the construction of this
channel. The first obvious freedom is the partitioning into the sets $A$ and
$B$. Section V.2 is devoted to this task. The next choice lies in the Trotter-
Suzuki order $2k$. Further, we can construct an outer loop as in Ref. 22 which
uses a Trotter-Suzuki style symmetrization of the channels themselves. For
example, a composite channel with outer order $2k=2$ looks like:

$\mathcal{X}^{2k}(\rho,t)=\mathcal{Q}_{B}{(t/2)}\circ\mathcal{T}_{A}^{2k}{(t)}\circ\mathcal{Q}_{B}{(t/2)}(\rho).$
(15)

Such an “outer-loop” with a respective “outer-order” was introduced as a
strategy to make more detailed analysis of simulation costs. The outer loop
can be thought of as a product formula for composite channels, symmetrizing
the channels and allocating the respective time slice, while the inner Trotter
order does so on the operator exponentials as before. However, when
introducing the outer-order in Ref. 22 it was enforced to match the same order
as the inner factorization. In this paper, one of our goals is to numerically
evaluate the exact costs of these composite algorithms, and this gives us the
freedom to independently set the inner and outer orders. However, an actual
application of the outer loop is situational as it requires specific knowledge
about the commutator structure for it to be of use. We further discuss this in
later sections where we will resort to the notation $\mathcal{X}^{2k,2l}$ for
$2k$ inner, $2l$ outer orders respectively. If only a single superscript is
given, it is understood to be the inner order with the outer order fixed at 1.

We now restate the first order cost result for composite channels from Ref.
22, where the $A$ terms are placed in a Trotter channel and the $B$ terms are
placed in QDrift:

$\displaystyle
C_{Comp}(A,B,t,\epsilon)=(L_{A}+N_{B})r=(L_{A}+N_{B})\left\lceil\frac{t^{2}}{\epsilon}\left(\frac{1}{2}\sum_{i,j}a_{i}a_{j}\norm{[A_{i},A_{j}]}_{\infty}+\frac{1}{2}\sum_{i,j}a_{i}b_{j}\norm{[A_{i},B_{j}]}_{\infty}+\frac{4\lambda_{B}^{2}}{N_{B}}\right)\right\rceil.$
(16)

We only restate the first order Trotter-Suzuki result due to the fact that
this is the expression from which we will draw most of intuition to build our
heuristics for partitioning. As well, the second order cost formula is of a
similar structure with extra $\alpha_{comm}$ and $\Upsilon$ terms, but it is
based on channels with matching inner-order and outer-order, which our
numerical work does not follow in general. Inspecting Equation (16), it
clearly inherits the structure from Trotter in its dependence on the
commutators and number of terms $L_{A}$ in the set $A$ (first term), and on
QDrift in the sum of the spectral norms $\lambda_{B}$ and number of samples
$N_{B}$ of the terms in the set $B$ (third term). The second term adds to the
cost due to the hybridization of the Trotter and QDrift channels. Clearly, in
order to gain as much advantage as possible using this algorithm over the non-
composite approaches, we desire that the Hamiltonian of interest has a
structure that we can exploit via a good choice of partitioning. An
interesting structure might be of an $H$ that contains terms of largely
varying magnitude, such that there might be a sharp contrast between the small
and large terms, rather than their spectral norms following say a Gaussian
distribution. Having small magnitude terms that largely outnumber the large-
magnitude terms is also likely desirable. This understanding naturally leads
into the next topic of how to choose the partitioning. This task can be
motivated by the structure of $H$ and its commutators, if known a-priori. For
example, inspecting Eq. (16), if we have a small number of large commuting
terms and numerous small non-commuting terms, it seems natural to place the
large terms into the Trotter channel, but sample the smaller magnitude terms
with QDrift. The question of choosing effective partitioning will be further
addressed in Section V.2.

###  II.6 Local Lattice Hamiltonians

Given the observation of commutator-dependent error in previous sections, an
interesting system to study are geometrically-local Hamiltonians. These
systems are defined on a lattice, and they have only local interactions, such
as nearest-neighbour terms. It turns out, via a Lieb-Robinson bound (Section
II.7), that commutators of interaction terms fall off exponentially with the
separation distance in the lattice. This motivates the search for a local
product decomposition to simulate the time evolution of local systems. In
fact, computational advantages have been demonstrated for local Hamiltonians
in Ref. 23. We follow this work in defining geometrically-local Hamiltonians
as

######  Definition 4.

Given a D-dimensional lattice $\Gamma\in\mathbb{Z}^{D}$, a local Hamiltonian
can be written in the following way:

$H=\sum_{X\subseteq\Gamma}H_{X},$ (17)

Where $H_{X}$ is only supported on the subset $X$ that contains only
immediately adjacent lattice points, and it acts as $\openone$ on
$\mathbb{Z}^{D}\in\Gamma\backslash X$. Equivalently, $H_{X}=0$ if $diam(X)>1$,
where $diam(X)=\max_{x,x^{\prime}\in X}dist(x,x^{\prime})$ and $dist$ is the
graph distance between the lattice indices. Each term in the Hamiltonian is
also normalized such that $\norm{H_{X}}\leq 1$.

In our analysis below, we will explicitly state when we are referring to a
local Hamiltonian. These Hamiltonians can be used to describe a wide variety
of physical systems, perhaps most famously is the quantum Heisenberg model,
which will be investigated numerically in Section V. If locality is not
explicitly mentioned then the Hamiltonian is of the aforementioned more
general structure.

###  II.7 Lieb-Robinson Bound

A Lieb-Robinson bound is essentially a bound on the speed at which information
propagates through a quantum system that has interactions governed by a local
Hamiltonian. This result is particularly interesting due to the fact that it
holds for non-relativistic quantum systems, meaning that in no part of the
system is the finite speed of light $c$ enforced. Instead, the locality of the
system’s interactions, as well as the geometry of the lattice leads to the
emergence of a Lieb-Robinson velocity $v_{LR}$, which bounds the speed of
causality. The central idea behind this bound is that the commutator between
the time evolution of an operator $A_{X}$ and another operator $B_{Y}$, such
that the operators are supported on disjoint sets $X\cap Y=\varnothing$, have
an exponentially small commutator. More formally, a general Lieb-Robinson
bound is often given in the following form [29, 30]:

$||[e^{iH_{\Gamma}t}Ae^{-iH_{\Gamma}t},B]||\leq
2||A||||B||Fe^{-\mu(dist(X,Y)-v_{LR}|t|)},$ (18)

where $\mu$ is a constant that depends on the lattice and $F$ depends on the
cardinality of the sets $X$ and $Y$. It is also important to note how the
distance between sets is defined, and it appears in the exponent above. The
distance between the sets $X$ and $Y$ is $dist(X,Y)=\min_{x\in X,y\in
Y}dist(x,y)$, where $dist(x,y)$ is the distance between elements in the set;
in a lattice this measure is just the absolute value of the difference of
their indices. From a simulation viewpoint, and specifically regarding our
Trotter-Suzuki product formulas, we know that the error in these
approximations depend on commutators of terms. Therefore, carefully-chosen
product decompositions may have exponentially-small error. This understanding
motivated the algorithm in Ref. 23, and in Sec. IV we combine this algorithm
with the composite framework.

##  III Imaginary Time Channels Algorithmic Analysis

Within this section, we show that the composite algorithms proposed by Hagan
and Wiebe in Ref. [22], for the purposes of digital real-time quantum
simulation, can be applied onto the imaginary time case with similar error
bounds and asymptotics. As mentioned previously, the algorithms in this
section are formulated with the intention of calculating properties of the
Hamiltonian on classical computers. These algorithms can therefore be
considered as quantum-inspired classical algorithms, a topic unrelated to the
problem of quantum imaginary time evolution [25]. Within this section, we
introduce the imaginary-time QDrift, Trotter-Suzuki, and the composite
channels, and we bound their distance norm with respect to the ideal unitary
evolution.

We first introduce the notion of an imaginary time evolution channel, given in
the definition below, and demonstrate how it can be viewed as a ground state
preparation, or “cooling” procedure.

######  Definition 5 (Imaginary Time Evolution Channel).

Given a Hamiltonian $H$ and input state $\rho$, we define the imaginary time
evolved state as the action of the following map

$\frac{\mathcal{U}{(\rho,\beta)}}{\Tr\mathcal{U}{(\rho,\beta)}}=\frac{e^{-\beta
H}\rho e^{-\beta H}}{\Tr e^{-\beta H}\rho e^{-\beta H}}.$ (19)

This can be seen as an imaginary time evolution channel most straightforwardly
when considering quantum state vectors as opposed to density matrix. The time-
independent Schödinger equation tells us
$\ket{\psi(t)}=e^{-iHt}\ket{\psi(0)}$. If we take this expression and perform
a Wick rotation that sends $t\to-i\beta$ we get that
$\ket{\psi(\beta)}=e^{-\beta H}\ket{\psi(0)}$. If we consider the action of
this matrix $e^{-\beta H}$ on the density matrix of the state $\ket{\psi(0}$
we get $e^{-\beta H}\outerproduct{\psi(0)}{\psi(0)}e^{-\beta H}$ due to the
fact that $e^{-\beta H}$ is a Hermitian operator. Relaxing the input from a
pure state to a density matrix and normalizing by the trace yields the channel
provided in Definition 5.

Often when imaginary time is discussed in quantum mechanics the state vector
representation is used, in which case the imaginary time exponential operator
has the following property

$\lim_{\beta\to\infty}\frac{e^{-\beta H}\ket{\psi_{i}}}{\norm{e^{-\beta
H}\ket{\psi_{i}}}_{2}}=\ket{\psi}_{GS}.$ (20)

This clearly only holds whenever our initial state $\ket{\psi}_{i}$ has
nonzero overlap with the ground state $\ket{\psi}_{GS}$, which is expressed as
$\innerproduct{\psi_{i}}{\psi_{GS}}\neq 0$. In the density matrix picture, we
have the following equivalent property:

$\lim_{\beta\to\infty}\frac{e^{-\beta H}\rho_{i}e^{-\beta H}}{\Tr e^{-\beta
H}\rho_{i}e^{-\beta H}}=\rho_{GS},$ (21)

such that
$\operatorname{supp}(\rho_{i})\cap\operatorname{supp}(\rho_{GS})\neq\varnothing$.
This property can be seen by expanding $\rho$ into an arbitrary mixed state in
the energy eigenbasis. Now that the maps of interest have been made clear, we
can consider algorithms that closely approximate them and analyze their
performance analytically by bounding their error. An important detail is the
choice of distance norm to quantify this error. In comparison to real time
composite algorithms, the aforementioned diamond norm is not a good choice for
imaginary time, given that this norm is about distinguishing quantum
operations on a system and ancillary space, whereas here we are simply
performing classical calculations. For this reason, we bound the following
induced Schatten $p\rightarrow q$ norm on a map $\Phi$:
$\norm{\Phi}_{p\rightarrow q}=\max_{\norm{\rho}_{q}=1}\norm{\Phi(\rho)}_{p}$.
Here, $\Phi$ will be the difference between the ideal map and that which is
generated by the composite algorithm, and we will investigate the case of the
induced trace norm where $p=q=1$. Despite the fact that the maps being
investigated in this section are non-linear, due to the trace operation in the
denominator, the outputs of the maps are Schatten-class linear operators, and
thus the induced Schatten norm is well defined. Further, by analyzing $1\to 1$
norms for two density matrices, such as
$\max_{\norm{\rho}_{1}=1}\norm{\rho-\sigma}_{1}\leq\epsilon$, we can guarantee
from properties of the trace distance that the measurement statistics for
$\rho$ and $\sigma$ will deviate by at most $\epsilon$.

###  III.1 Imaginary Time QDrift

Introduced in Section II.4, QDrift was introduced as a quantum simulation
algorithm that, in real-time simulations, has a unique property in that its
error does not depend on the number Hamiltonian summands, nor on the
commutators between them, only on the total size of these operators. To the
best of our knowledge, this algorithm has not been applied to imaginary time
calculations. We recover all the characteristic properties of the QDrift
channel in imaginary time with renormalizing trace operations in the following
Theorem.

######  Theorem III.1 (Imaginary Time QDrift).

Given a Hamiltonian $H=\sum_{j}h_{j}H_{j}$ with $\lambda=\sum_{j}h_{j}$,
imaginary time $\beta$, number of samples $N$, and a density matrix initial
state $\rho$, then the induced Schatten $1\rightarrow 1$ norm of the
difference between the imaginary time QDrift channel and exact imaginary time
evolutionary channel has the following bound:

$\norm{\frac{\mathcal{U}(\rho,\beta)}{\Tr\mathcal{U}(\rho,\beta)}-\left(\frac{\mathcal{Q}{(\rho,\nicefrac{{\beta}}{{N}})}}{\Tr\mathcal{Q}{(\rho,\nicefrac{{\beta}}{{N}})}}\right)^{\circ
N}}_{1\rightarrow
1}\leq\frac{4\beta^{2}\lambda^{2}}{N}\frac{e^{\nicefrac{{2\lambda\beta}}{{N}}}}{2-e^{\nicefrac{{2\beta\lambda}}{{N}}}}+2\frac{\frac{4\beta^{4}\lambda^{4}}{N^{3}}e^{\nicefrac{{6\lambda\beta}}{{N}}}}{1-\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}}+2\frac{2\beta^{2}\lambda^{2}}{N}e^{\nicefrac{{4\lambda\beta}}{{N}}},$
(22)

given that $|e^{\nicefrac{{2\beta\lambda}}{{N}}}-1|<1$ which is satisfied if
$N>\frac{2\beta\lambda}{\ln 2}$. Further, If the constraint
$\frac{\lambda}{N}\leq 0.01$ is satisfied, then the bound simplifies to the
following:

$\norm{\frac{\mathcal{U}(\rho,\beta)}{\Tr\mathcal{U}(\rho,\beta)}-\left(\frac{\mathcal{Q}{(\rho,\nicefrac{{\beta}}{{N}})}}{\Tr\mathcal{Q}{(\rho,\nicefrac{{\beta}}{{N}})}}\right)^{\circ
N}}_{1\rightarrow 1}\leq\frac{C\beta^{2}\lambda^{2}}{N},$ (23)

where the constant $C\approx 29.71747$.

For a proof of this theorem see Appendix A. The proof consists of repeated
applications of the triangle inequality and sub-multiplicative property of the
diamond norm, as well as manipulations with the geometric series and a tail
bound on a Taylor series. The geometric series arises when one uses a Taylor
expansion on the terms in the denominator and takes the trace of the zeroth
order term, which is just $\Tr\rho=1$. A constant theme throughout the
analytics in this paper are multiple exponential factors that arise due to the
techniques used to bound trace terms. Multiple techniques are used to
accomplish this throughout, and the exponential factors are unavoidable in
imaginary time, given that the exponential operator is no longer unitary and
has an exponentially decaying norm. For clearer interpretation of the bounds
we use linearization techniques to either show that they disappear under some
constraints when proving an inequality, or showing they become
$\mathcal{O}\left(1\right)$ in some limit for an asymptotic bound. The
important takeaway from this theorem is that the QDrift error scales
essentially the same was in imaginary time as it does in real time, thus
making it a candidate algorithm to be applied in Quantum-Monte Carlo
calculations. As well, this promise on the ratio $\frac{\lambda}{N}$ is
reasonable to make as this algorithm excels in the regime of small terms.

###  III.2 Imaginary Time Trotter-Suzuki

Different from QDrift, work has already been done to provide bounds on non-
unitary product formulas [10]. However, these bounds are given in terms of
operator norms rather than channel distances, as well as they are given
without renormalizing trace operations. However, these existing results are
integral in proving the theorem below.

######  Theorem III.2 (Imaginary-Time Trotter-Suzuki Channels).

Given a Hamiltonian $H=\sum_{j}h_{j}H_{j}$ with $\lambda=\sum_{j}h_{j}$,
imaginary time $\beta$, and a density matrix initial state $\rho$, then the
induced Schatten $1\rightarrow 1$ norm of the difference between the imaginary
time Trotter-Suzuki channel $\mathcal{T}^{2k}{(\rho,\beta)}$ and the exact
imaginary time evolutionary $\mathcal{U}{(\rho,\beta)}$ channel has the
following bound:

$\norm{\frac{\mathcal{U}{(\rho,\beta)}}{\Tr(\mathcal{U}{(\rho,\beta)})}-\left(\frac{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}\right)^{\circ
r}}_{1\rightarrow 1}\leq
2\Upsilon^{2k+1}\frac{\alpha_{comm}(H,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}e^{\nicefrac{{4\Upsilon\beta\lambda}}{{r}}}\left(e^{\nicefrac{{4\beta\lambda}}{{r}}}e^{\nicefrac{{2\beta\norm{H}}}{{r}}}+e^{\nicefrac{{2\beta\norm{H}}}{{r}}}\right),$
(24)

which yields the following asymptotic bound:

$\norm{\frac{\mathcal{U}{(\rho,\beta)}}{\Tr(\mathcal{U}{(\rho,\beta)})}-\left(\frac{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}\right)^{\circ
r}}_{1\rightarrow
1}\in\mathcal{O}\left(\Upsilon^{2k+1}\frac{\alpha_{comm}(H,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}\right).$
(25)

Here $\Upsilon$ is the number of stages of the product formula, $2k$ is the
order of the product formula, and $\lambda$ is the sum of the spectral norms
of the Hamiltonian summands. Here $\mathcal{O}\left(\cdot\right)$ is
understood in the infinite limit of its arguments.

For a proof of this theorem, see Appendix A. The proof follows similar
strategies to that of imaginary-time QDrift, however, one cannot simply expand
a geometric series in the denominator here given that we have a product
instead of a sum of terms. So we instead bound the trace terms using a
combination of von Neumann’s trace inequality and Weyl’s inequality regarding
the singular values of matrices. Similar to the QDrift case, we obtain a
similar looking bound to that obtained in [22] for the real time case, with
some additional exponential factors as expected. Considering asymptotics, the
bounds behave the same in both cases where we again see the dependence of the
error on a nested commutator sum.

###  III.3 Imaginary Time Composite Channels

Now equipped with error bounds on both QDrift and Trotter-Suzuki channels in
imaginary time, we can proceed with the analysis of channels composed of the
two. Given the difficulty of analyzing the renormalizing operations in the
trace terms, we will again make use of asymptotic notation as was done in the
Trotter-Suzuki analysis in the previous section. We will once again write a
partitioning of the Hamiltonian as $H=A+B$, where
$A=\sum_{i=1}^{L_{A}}a_{i}A_{i}$ are the terms placed in the Trotter-Suzuki
channel and $B=\sum_{i=1}^{L_{B}}b_{i}B_{i}$ are placed in the QDrift channel.

######  Theorem III.3 (Imaginary Time Composite Channels).

Given a partitioned Hamiltonian $H=\sum_{j}a_{j}A_{j}+\sum_{i}b_{i}B_{i}$ with
$\lambda_{A}=\sum_{j}a_{j}$, $\lambda_{B}=\sum_{j}b_{j}$, imaginary time
$\beta$, QDrift samples $N_{B}$ and a density matrix initial state $\rho$,
then the induced Schatten $1\rightarrow 1$ norm of the difference between the
imaginary time Composite channel $\mathcal{X}^{2k}(\rho,\beta)$ of order $2k$
and the exact imaginary time evolutionary channel $\mathcal{U}{(\rho,\beta)}$
has the following bound:

$\norm{\frac{\mathcal{U}{(\rho,\beta)}}{\Tr\mathcal{U}{(\rho,\beta)}}-\frac{\mathcal{X}^{2k}(\rho,\beta)}{\Tr\mathcal{X}^{2k}(\rho,\beta)}}_{1\rightarrow
1}\in\mathcal{O}\left(\frac{\beta^{2}\lambda_{B}^{2}}{N_{B}r}+\Upsilon^{2k+1}\frac{\alpha_{comm}(A,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}+\norm{[A,B]}\frac{\beta^{2}}{r}\right)$
(26)

given that $r\geq\Upsilon\beta\lambda$. Here, all parameters in this bound
correspond to the same quantities in Theorems III.1 and III.2, where the
subscripts $A$ and $B$ indicate their belonging corresponding set.

For a proof of this bound see Appendix A. The condition on $r$ is not overly
restrictive, nor arbitrary. In general, we do not expect to simulate the
imaginary time evolution with $r$ sub-linear in either the imaginary time
$\beta$ or the norms of the Hamiltonian $\lambda$, and $\Upsilon$ is an
artefact of using higher order Trotter-Suzuki formulas. This result yields
three error terms that we might expect. It is a sum of the QDrift and Trotter
error, plus the error that arises from the initial partitioning of the
Hamiltonian into two subsets. Therefore, it becomes apparent that if we have
information about the terms in the Hamiltonian, we can potentially exploit the
form of this error bound by choosing a good partition. By good partition, we
mean one that is close to the optimal one that minimizes the error. We later
show in Section V.2 how to possibly go about choosing a partition using
intuition from physics, as well as provide a machine learning regressor to
find the optimal partition. The later is more for a proof of principle rather
than something you would run in practice due to the large overhead.

###  III.4 Examples and Implementation

Within this section we illustrate an implementation of the imaginary time
composite channel with an example. Approximating the imaginary time propagator
$e^{-\beta H}$, which is a function of a Hamiltonian $H$ that is hard to
exponentiate in general, can be achieved through the algorithm written in
pseudo-code below. This algorithm returns a list of propagators (of which
their product approximates $e^{-\beta H}$) that are instead easily
diagonalizable.

Input: Simulation parameters: Trotter-Suzuki inner-order $2k$, QDrift samples
$N_{B}$, iterations $r$, and imaginary time $\beta$. A list of Hamiltonian
terms $H=\sum_{j}h_{j}H_{j}$, two deterministic classical functions:
PARTITION(_$H$_) which returns two lists of Hamiltonian terms
$A=\sum_{q}a_{q}A_{q}$ and $B=\sum_{p}b_{p}B_{p}$, and
TROTTER(_$A,2k,\beta,r$_) that returns a Trotter-Suzuki product formula of
order $2k$ with time-slice $\frac{\beta}{r}$. We also require a classical
oracle function SAMPLE(__) which returns a value $j$ from the probability
distribution $p_{j}=\frac{b_{j}}{\sum_{l}b_{l}}$.

Output: A composite imaginary time product formula in the form of a list of
propagators, with inner-order $2k$ and outer-order $1$, which closely
approximates the imaginary time propagator $e^{-\beta H}$.

Function  _CompositeListCompilation( $H,\beta,r,2k,N_{B}$)_:

$A,B\leftarrow$ PARTITION(_$H$_);

$m\leftarrow 1$;

$\mathcal{X}\leftarrow$ { }(empty ordered list);

$\mathcal{T}\leftarrow$ TROTTER(_$A,2k,\beta,r$_);

while  _$m\leq r$_ do

$\mathcal{Q}\leftarrow$ { }(empty ordered list);

$\lambda_{B}\leftarrow\sum_{l}b_{l}$;

$j\leftarrow 1$;

while  _$j\leq N_{B}$_ do

$j\leftarrow$ SAMPLE(__) ;

append $e^{-\beta\lambda_{B}B_{j}/N_{B}r}$ to $\mathcal{Q}$;

$j\leftarrow j+1$;

end while

append $\mathcal{T}$ to $\mathcal{X}$;

append $\mathcal{Q}$ to $\mathcal{X}$;

$m\leftarrow m+1$;

end while

return $\mathcal{X}$;

Algorithm 1 Pseudo-code for implementation of an imaginary time composite
simulation using a high-order Trotter formula and a partitioning heuristic to
divide the Hamiltonian terms between the two channels.

The algorithm above is written for illustrative purposes and simply constructs
a composite product formula in the form of an ordered list. This can then be
used as a subroutine for other algorithms in computational physics, such as
those from the field of Quantum Monte-Carlo (QMC). For example, in statistical
physics calculating the partition function $Z$ is of central interest. Given
this returned list, which for clarity, is not a channel, we can approximate
$Z$ as follows:

$\displaystyle Z$ $\displaystyle=\Tr(e^{-\beta
H})=\Tr((e^{-\frac{\beta}{r}A_{1}}...e^{-\frac{\beta}{r}A_{L_{A}}}e^{-\frac{\beta\lambda_{B}}{rN_{B}}B_{1}}...e^{-\frac{\beta\lambda_{B}}{rN_{B}}B_{N_{B}}})^{r})+\epsilon,$
(27)

where $\epsilon$ is bounded by a set of functions implied by Theorem III.3.
Next, we insert the resolution of the identity between each propagator such
that the trace disappears

$\displaystyle
Z=\sum_{j_{1},j_{2},...j_{n}}\bra{j_{n}}e^{-\frac{\beta}{r}A_{1}}\ket{j_{1}}...\bra{j_{L_{A}-1}}e^{-\frac{\beta}{r}A_{L_{A}}}\outerproduct{j_{L_{A}}}{j_{L_{A}}}e^{-\frac{\beta\lambda_{B}}{rN_{B}}B_{1}}\ket{j_{L_{A}+1}}...\bra{j_{n-1}}e^{-\frac{\beta\lambda_{B}}{rN_{B}}B_{N_{B}}}\ket{j_{n}}+\epsilon,$
(28)

given that $n=|\mathcal{X}|r$ or the number of terms in the list multiplied by
the number of iterations. Note the indices are written such that there
$L_{A}r$ Trotter terms and $N_{B}r$ QDrift terms, with the difference that the
subscripts indexing the QDrift terms are dummy indices, as these propagators
are importance sampled at random, as shown in the pseudo-code. Now, with the
above expression, we can either diagonalize these terms and Monte-Carlo sample
this sum, or impose a classical mapping that allows us to apply some of the
well known QMC algorithms. These classical mappings depend on the system in
play, for example, if we are interested in Heisenberg or XXZ chains, it is
possible to map the model onto a higher dimensional Ising model and apply the
method of world lines [20]. The motivation for using this more complicated
propagator list, as opposed to the standard Trotterization, is that a
propagator with a lower error requires a smaller $r$ to achieve the desired
precision, and can therefore reduce the path length of the Monte-Carlo path
integral required to approximate $Z$ above. Our numerics show that this
framework can achieve a substantially lower error in the right conditions via
attaining a smaller $C$ in multiple cost plots in Section V.4.2, therefore,
highlighting the potential applications of this algorithm to improve QMC
simulations.

In this example, it was chosen to show arbitrary inner-order and outer-order
1, given that the goal is simplicity. This can be simply generalized to higher
outer-orders by also requiring an additional subroutine that symmetrizes
$\mathcal{T}$ and $\mathcal{Q}$, and adjusts their time-slices accordingly. We
also remark that the classical function PARTITION() is presented here somewhat
as a black box. This is done for generality, given that the user is free to
use a method of their choosing. In section V.2, we present an explicit
algorithm for this function, which requires a single additional input
$\omega_{c}$ called the chop threshold, and simply involves a looping if
statement that sorts Hamiltonian summands into $A$ or $B$ depending on whether
their spectral norm is greater or less than $\omega_{c}$.

##  IV Local Composite Channels for Real-Time Evolution

We introduce here an algorithm for real time evolution, which utilizes the
local structure of the Hamiltonian and is composite. It combines intuition
from both the Trotter-QDrift composite channels as discussed in Sec. II.5 and
the Lieb-Robinson decomposition, Sec. II.7. Our main result is Eq. (37), which
bounds the error of a local composite channel.

The basic idea behind the method goes as follows: We decompose the Hamiltonian
evolution into local blocks and build composite channels out of the operators
supported on these local blocks. Each block is thus assigned a local partition
$\\{A,B\\}$ and a sample number $N_{B}$.

Based on prior intuition, we expect this approach to perform well on
disordered models yet with local interactions, such that there is a structure
for the local composite channels to exploit. However, by adding QDrift
sampling to the local framework, with the speed limit $v_{LR}$ in mind we are
“fuzzing” the light-cone in the quantum circuit evolution of the system. For
example, in Figure 2 we show a circuit that illustrates the Trotter time
evolution of a local Hamiltonian where the $H_{i,j}$ terms only connect
neighbouring qubits within and between each disjoint subset $A$ and $B$ of the
1-dimensional lattice.

$\ket{\idx}_{A}$$\ket{4}_{B}$$H_{0,1}$$H_{1,2}$$H_{2,3}$$H_{3,4}$ Figure 2:
Brickwork circuit diagram to highlight the emergent light-cones from nearest
neighbour interactions. The notion of “fuzzing” the light-cone occurs when we
use QDrift to sample terms along these light-cone boundaries in the circuit.

The light-cone analogy is now made clear: A term $H_{0,1}$ influences the most
adjacent qubit in $B$ only once 3 “layers” of gates have been applied
(intermittent gates cannot improve this speed). This behavior is analogous to
the Lieb-Robinson velocity that upper-bounds the propagation of information
flow in the lattice. We note though that in a local composite channel, some of
the gates are being randomly sampled in the QDrift channel, thus the evolution
of the light-cone is not exactly reproduced since not all terms are
necessarily implemented due to sampling. Keep in mind that local composite
channel does include more complicated subsets of gates, some of which evolve
boundary region terms like $H_{3,4}$ backwards in time; the circuit above is
presented simply for conceptual purposes.

Towards bounding the error of a local composite channel, we build on the
result from Ref. 23, which concerns with bounds on time evolution operators,
and reformulate it as a bound on channels. The motivation for the
reformulation of the bound into a channel norm is so that QDrift could be
incorporated into this algorithm. This would later allow us to bound the
overall error. Since QDrift produces mixed states, we need to work with the
density matrix formalism.

######  Theorem 6.

Given a time $t$ and a local Hamiltonian $H=\sum_{X}H_{X}$ that generates a
time evolution unitary $U$ and an evolutionary channel $\mathcal{U}$,
$\exists\;(\mu>0)\;|$ for any disjoint regions $A,B,C$ we have the following
operator norm bound:

$||\mathcal{U}_{A\cup B}\mathcal{U}^{\dagger}_{B}\mathcal{U}_{B\cup
C}-\mathcal{U}_{A\cup B\cup C}||\in\mathcal{O}\left(e^{-\mu
dist(A,C)}\right)\sum_{X:bd(AB,C)}||H_{X}||,$ (29)

where the evolutionary channels contain only Hamiltonians terms supported on
their subscripted sets of the lattice, and X:bd(AB,C) indicates $X\subseteq
A\cup B\cup C$ and $X\not\subseteq A\cup B$ and $X\not\subseteq C$.

###### Proof.

We first restate the main result from Ref. 23 on time evolution operators,

$||U_{A\cup B}U^{\dagger}_{B}U_{B\cup C}-U_{A\cup B\cup
C}||\in\mathcal{O}\left(e^{-\mu dist(A,C)}\right)\sum_{X:bd(AB,C)}||H_{X}||.$
(30)

Now, write $V:=U_{A\cup B}U^{\dagger}_{B}U_{B\cup C}$ and $U:=U_{A\cup B\cup
C}$ and consider the following induced channel infinity norm

$\displaystyle\max_{\rho:\norm{\rho}_{1}\leq 1}||V\rho V^{\dagger}-U\rho
U^{\dagger}||$ $\displaystyle=\max_{\rho:\norm{\rho}_{1}\leq 1}||V\rho
V^{\dagger}-U\rho V^{\dagger}+U\rho V^{\dagger}-U\rho U^{\dagger}||$
$\displaystyle=\max_{\rho:\norm{\rho}_{1}\leq 1}||(V-U)\rho
V^{\dagger}+U\rho(V^{\dagger}-U^{\dagger})||$
$\displaystyle\leq\max_{\rho:\norm{\rho}_{1}\leq 1}||(V-U)\rho
V^{\dagger}||+\max_{\rho:\norm{\rho}_{1}\leq 1}||U\rho(V-U)^{\dagger}||$
$\displaystyle\leq 2||V-U||.$ (31)

In the last line we used the fact that the infinity norm of a density matrix
is upper bounded by 1, and the unitary invariance of the Schatten norms.
Inserting the original definitions for $V$ and $U$ into
$\max_{\rho:\norm{\rho}_{1}\leq 1}||\mathcal{V}-\mathcal{U}||\leq 2||V-U||$,
where $\mathcal{V}$ is just $V\rho V^{\dagger}$, and similarly
$\mathcal{U}=U\rho U^{\dagger}$, we successfully upper bounded the desired
channel norm with the result with the constant factor disappearing due to
$\mathcal{O}$ notation. ∎

The bound (6) we have proven above is still not of the form of that given in
QDrift. QDrift bound utilizes the diamond distance, which is a completely
bounded trace norm or Schatten 1-norm defined as
$\left|\left|\mathcal{V}-\mathcal{U}\right|\right|_{\diamond}=\max_{\rho:\norm{\rho}_{1}\leq
1}||(\mathcal{V}-\mathcal{U})\otimes\openone||_{1}$. While this is an upper
bound on the infinity norm in general, it is not immediately clear how tight
this bound may be as factors of dimensionality of the system may come into
play, which we wish to avoid. This can be made clear by a simple example:
Consider an $N\times N$ identity matrix $\openone_{N}$. This matrix has
$||\openone_{N}||_{1}=N$ while $||\openone_{N}||_{\infty}=1$, where $N$ is a
factor of dimensionality. Therefore, if possible, we would also like to
convert the result (6) to a diamond norm. Building on Ref. 22, we accomplish
this in the following lemma:

######  Lemma 7.

Given a time $t$ and a local Hamiltonian $H=\sum_{X}H_{X}$ that generates a
time evolution unitary $U$ and an evolutionary channel $\mathcal{U}$,
$\exists\;(\mu>0)\;|$ for any disjoint regions $A,B,C$ we have the following
diamond distance bound

$\left|\left|\mathcal{U}_{A\cup B}\mathcal{U}^{\dagger}_{B}\mathcal{U}_{B\cup
C}-\mathcal{U}_{A\cup B\cup
C}\right|\right|_{\diamond}\in\mathcal{O}\left(e^{-\mu
dist(A,C)}\right)\sum_{X:bd(AB,C)}||H_{X}||$ (32)

###### Proof.

We use Equations (11-17) from Ref. 22 replacing their Trotter channel with the
local block channel $\mathcal{V}=\mathcal{U}_{A\cup
B}\mathcal{U}^{\dagger}_{B}\mathcal{U}_{B\cup C}$ and set
$\mathcal{U}=\mathcal{U}_{A\cup B\cup C}$ with $V$ and $U$ representing the
unitaries that induce each channel,

$\displaystyle\left|\left|\mathcal{U}(t)-\mathcal{V}(t)\right|\right|_{\diamond}:=$
$\displaystyle\norm{\left(\mathcal{U}(t)-\mathcal{V}(t)\right)\otimes\openone}_{1}$
$\displaystyle=$ $\displaystyle\max_{\rho:\norm{\rho}_{1}\leq
1}\norm{e^{-iHt}\otimes\openone\rho
e^{iHt}\otimes\openone-V(t)\otimes\openone\rho
V^{\dagger}(t)\otimes\openone}_{1}$ $\displaystyle\leq$
$\displaystyle\max_{\rho:\norm{\rho}_{1}\leq
1}\norm{e^{-iHt}\otimes\openone\rho
e^{iHt}\otimes\openone-e^{-iHt}\otimes\openone\rho
V^{\dagger}(t)\otimes\openone}_{1}$
$\displaystyle+\max_{\rho:\norm{\rho}_{1}\leq
1}\norm{e^{-iHt}\otimes\openone\rho V^{\dagger}(t)-V(t)\otimes\openone\rho
V^{\dagger}(t)\otimes\openone}_{1}$ $\displaystyle=$
$\displaystyle\max_{\rho:\norm{\rho}_{1}\leq
1}\norm{\rho\left(e^{-iHt}-V^{\dagger}(t)\right)\otimes\openone}_{1}+\max_{\rho:\norm{\rho}_{1}\leq
1}\norm{\left(e^{-iHt}-V(t)\right)\otimes\openone\rho}_{1}$
$\displaystyle\leq$ $\displaystyle
2\norm{e^{-iHt}-V(t)}_{\infty}\max_{\rho:\norm{\rho}_{1}\leq
1}\norm{\rho}_{1}$ $\displaystyle=$ $\displaystyle
2\norm{e^{-iHt}-V(t)}_{\infty}$ $\displaystyle=$ $\displaystyle
2\norm{U(t)-V(t)}_{\infty}.$ (33)

In the last line, we inserted the Theorem 6 to complete the proof. We observe
the same bound on our Diamond distance as appeared in the induced infinity
norm with the constant factor 2 once again absorbed into $\mathcal{O}$. ∎

Onward, the idea of the the local-composite algorithm is to start with some
Hamiltonian $H=\sum_{X}H_{X}$ defined on some lattice $\Gamma$ with
$X\subseteq\Gamma$. Begin by breaking up the lattice into two regions: $A\cup
B=\Gamma$. Now we take a smaller region of the lattice $Y\subseteq A$ and use
the result to approximate the time evolution operator in the following way,

$e^{-iHt}\approx e^{-iH_{A}t}(e^{-iH_{Y}t})^{\dagger}e^{-iH_{Y\cup B}t}$ (34)

Note that this breakup is distinct from Trotterization in that we are
decomposing the operator into a product formula of which the operator
exponentials still contain sums of operators. However, the operator sums now
only contain terms in each local “block”. We can re-curse this process to make
$m$ blocks of near-equal size. Next, we decide which simulation method to
employ to approximate and implement the block-evolution operators. In the
previous theorem (32), a bound was proved for only one round of decomposition,
but note that when defining $V$, the proof holds regardless of how many
unitary blockings the operator is composed of. In Ref. 23, it was stated that
the error in Theorem 6 for $D$ rounds of decomposition in a $D$-dimensional
lattice is $\in\mathcal{O}\left(e^{-\mu l}\nicefrac{{DL^{D}}}{{l}}\right)$ for
a lattice with chain length $L$ and blocks with overlap $l$. The overlap for
the block in the above example would be the diameter of the set $Y$, which can
be written $diam(Y)=\max_{y,y^{\prime}\in Y}dist(y,y^{\prime})$; the latter
distance is again the distance between lattice indices.

Before investigating this channel numerically, we provide a bound on the
diamond norm between this local composite channel and the exact evolutionary
channel. Along with the fact that the diamond norm is sub-additive and sub-
multiplicative, we will make use of the following two properties:

$\displaystyle\left|\left|\mathcal{A}^{\circ r}-\mathcal{B}^{\circ
r}\right|\right|_{\diamond}$ $\displaystyle\leq
r\left|\left|\mathcal{A}-\mathcal{B}\right|\right|_{\diamond}$
$\displaystyle\left|\left|\mathcal{A}\mathcal{B}-\mathcal{C}\mathcal{D}\right|\right|_{\diamond}$
$\displaystyle\leq\left|\left|\mathcal{A}-\mathcal{C}\right|\right|_{\diamond}+\left|\left|\mathcal{B}-\mathcal{D}\right|\right|_{\diamond},$

Where $\mathcal{A},\mathcal{B},\mathcal{C}$ and ${\mathcal{D}}$ represent
channels.

These properties were shown in Ref. 31. For our purposes, we require a more
general property that holds for the composition of $k$ channels
$\mathcal{A}_{i}$. We provide a simple proof in the following lemma:

######  Lemma 8.

Given a composition of $k$ quantum channels $\mathcal{A}_{i}$, and another
composition of $k$ quantum channels $\widetilde{\mathcal{A}_{i}}$, the diamond
norm of the difference between these two channels has the following bound

$\left|\left|\bigcirc^{k}_{i=1}\mathcal{A}_{i}-\bigcirc^{k}_{i=1}\mathcal{B}_{i}\right|\right|_{\diamond}\leq\sum_{i=1}^{k}\left|\left|\mathcal{A}_{i}-\mathcal{B}_{i}\right|\right|_{\diamond}.$
(35)

###### Proof.

The proof is achieved by induction and it uses the fact that
$\left|\left|\mathcal{A}_{i}\right|\right|_{\diamond}\leq 1$ for any quantum
channel, which is due to the fact that the spectral norm is unitary invariant
and density matrices have a maximum eigenvalue of 1. Starting from the left
and suppressing the composition operation $\circ$ above, we have

$\displaystyle\left|\left|\mathcal{A}_{k}\mathcal{A}_{k-1}...\mathcal{A}_{1}-\mathcal{B}_{k}\mathcal{B}_{k-1}...\mathcal{B}_{1}\right|\right|_{\diamond}$
$\displaystyle=\left|\left|\mathcal{A}_{k}\mathcal{A}_{k-1}...\mathcal{A}_{1}-\mathcal{B}_{k}\mathcal{A}_{k-1}\mathcal{A}_{k-2}...\mathcal{A}_{1}+\mathcal{B}_{k}\mathcal{A}_{k-1}\mathcal{A}_{k-2}...\mathcal{A}_{1}-\mathcal{B}_{k}\mathcal{B}_{k-1}...\mathcal{B}_{1}\right|\right|_{\diamond}$
$\displaystyle=\left|\left|(\mathcal{A}_{k}-\mathcal{B}_{k})\mathcal{A}_{k-1}...\mathcal{A}_{1}+\mathcal{B}_{k}(\mathcal{A}_{k-1}\mathcal{A}_{k-2}...\mathcal{A}_{1}-\mathcal{B}_{k-1}\mathcal{B}_{k-2}...\mathcal{B}_{1})\right|\right|_{\diamond}$
$\displaystyle\leq\left|\left|\mathcal{A}_{k}-\mathcal{B}_{k}\right|\right|_{\diamond}\left|\left|\mathcal{A}_{k-1}...\mathcal{A}_{1}\right|\right|_{\diamond}+\left|\left|\mathcal{B}_{k}\right|\right|_{\diamond}\left|\left|\mathcal{A}_{k-1}\mathcal{A}_{k-2}...\mathcal{A}_{1}-\mathcal{B}_{k-1}\mathcal{B}_{k-2}...\mathcal{B}_{1}\right|\right|_{\diamond}$
$\displaystyle\leq\left|\left|\mathcal{A}_{k}-\mathcal{B}_{k}\right|\right|_{\diamond}+\left|\left|\mathcal{A}_{k-1}\mathcal{A}_{k-2}...\mathcal{A}_{1}-\mathcal{B}_{k-1}\mathcal{B}_{k-2}...\mathcal{B}_{1}\right|\right|_{\diamond}.$
(36)

Repeating this procedure, we decompose the rightmost term into a sum of two
channel differences by induction, allowing one to prove the original identity
in Equation (35). ∎

Next, we utilize Theorem 6, where we converted the spectral norm result on
local block unitaries to a bound on a local decomposition channel. It is now
made clear why this was done; we require a channel description to describe the
QDrift algorithm and relate their bounds. Using Lemma 8, we proceed with the
initial goal of bounding the diamond norm between the local composite channel
and the exact evolutionary channel. Here, for compactness we adopt the
notation of writing blocked composite local channels as $\mathcal{X}_{Y_{i}}$
and “exact” local channels (prior to hybridization into Trotter and QDrift)
$\mathcal{U}_{Y_{i}}$. Here, $Y_{i}$ indicates a set of lattice points
$Y_{i}\subseteq\Gamma$ such that each $Y_{i}$ is constructed via the algorithm
laid out earlier in this Section. Therefore, only certain Hamiltonian
interaction terms are supported on the lattice subsets $Y_{i}$. Additionally,
let $\mathcal{U}_{Y_{i}}$ contain terms that have also undergone the
appropriate conjugation (given that the intersections of some blocks are
evolved backwards in time). Finally, let each set $Y_{i}\in Y$ such that
$Y=\\{Y_{1},Y_{2},...,Y_{m}\\}$ and $m$ is the number of blockings constructed
by the Lieb-Robinson decomposition. This provides convenient notation for the
following theorem:

######  Theorem 9.

Given a time $t$, local channel iterations $r$, and a local Hamiltonian
$H_{\Gamma}$ on a $D$-dimensional lattice $\Gamma$, let $\mathcal{X}_{Y}(t)$
be a local-composite channel with $m$ local blocks $Y_{i}$, and partition
$\\{A,B\\}_{Y}$ on each block, then the diamond norm of the difference between
the local composite channel and the evolutionary channel has the following
bound:

$\displaystyle\left|\left|\bigcirc_{Y}\left(\mathcal{X}_{Y}(\nicefrac{{t}}{{r}})\right)^{\circ
r}-\mathcal{U}(t)\right|\right|_{\diamond}\leq
m\max_{Y}\left|\left|\left(\mathcal{X}_{Y}(\nicefrac{{t}}{{r}})\right)^{\circ
r}-\mathcal{U}_{Y}(t)\right|\right|_{\diamond}+\epsilon_{LR}$
$\displaystyle\quad=\frac{t^{2}m}{2r}\max_{Y}\left(\left(\sum_{i,j}a_{Y_{i}}a_{Y_{j}}\norm{[A_{Y_{i}},A_{Y_{j}}]}+\sum_{i,j}a_{Y_{i}}b_{Y_{j}}\norm{[A_{Y_{i}},B_{Y_{j}}]}\right)+\frac{4\lambda_{B_{Y}}^{2}t^{2}}{N_{B_{Y}}r}\right)+\mathcal{O}\left(e^{-\mu
l}DL^{D}/l\right)$ (37)

where $A_{Y},B_{Y}$ are once again the terms of $H_{Y}$ in Trotter and QDrift
with their respective spectral norms $a_{Y},b_{Y}$, $N_{B}$ is the number of
QDrift samples and $\lambda$ is the sum of the spectral norms of the terms in
QDrift. It should also be stressed that $r$ is a local channel iteration,
which means that we do not iterate the channel that is a composition of blocks
but rather each individual local composite channel built out of the block
terms. The reason for this is that the error in the local decomposition is
independent of time for time independent $H$ with
$t\in\mathcal{O}\left(1\right)$. We also require
$m\in\mathcal{O}\left((L/l)^{D}\right)$ [23].

###### Proof.

We begin with the norm we wish to bound, and then apply the above identities
to format the norm such that known bounds from [23] and [22] can be applied.

$\displaystyle\left|\left|\bigcirc_{i=1}^{m}\left(\mathcal{X}_{Y}(\nicefrac{{t}}{{r}})\right)^{\circ
r}-\mathcal{U}(t)\right|\right|_{\diamond}$
$\displaystyle=\left|\left|\bigcirc_{i=1}^{m}\left(\mathcal{X}_{Y_{i}}(\nicefrac{{t}}{{r}})\right)^{\circ
r}-\bigcirc_{i=1}^{m}\mathcal{U}_{Y_{i}}(t)+\epsilon_{LR}\right|\right|_{\diamond}$
(38)
$\displaystyle\leq\sum_{i=1}^{m}\left|\left|\left(\mathcal{X}_{Y_{i}}(\nicefrac{{t}}{{r}})\right)^{\circ
r}-\mathcal{U}_{Y_{i}}(t)\right|\right|_{\diamond}+\epsilon_{LR}$ (39)
$\displaystyle\leq
m\max_{Y}\left|\left|\left(\mathcal{X}_{Y}(\nicefrac{{t}}{{r}})\right)^{\circ
r}-\mathcal{U}_{Y}(t)\right|\right|_{\diamond}+\epsilon_{LR}$ (40)

Then by inserting the bound given in [22] Equation 31, we arrive at the
original expression. It then immediately follows that we can bound local
composite channels of any order $2k$ by inserting the result of [22] Equation
58, however, we omit this as it follows trivially from above. ∎

The takeaway from this is that the block-local composite channel is bounded by
an error that is very similar to that of the general composite channel asides
from the addition of an extra term that is asymptotically exponentially small.
However, while we expect to see improvement from this localization scheme,
getting more out of this bound does not seem feasible due to the amount of
freedom in this algorithm as almost every parameter is dependent on the
“worst” subset of the set $Y$.

##  V Numerical Simulations

In this section we highlight the details of configuring our numerical
investigations by defining Hamiltonians of interest. We then discuss
partitioning strategies and choice of an error measure before presenting
numerical results. The choice of error measure usually does not require a
dedicated subsection. However, when considering composite algorithms, there is
a subtle feature in that not all algorithms may be treated on equal footing by
the same error measure, which can lead to inconsistencies. For our proposes,
the trace distance is the most useful measure of choice. Rather than computing
upper-bounds, the numerical results we provide here are exact calculations of
the cost $C$, or the precise counts on the number of quantum gates with the
form $e^{-i\frac{t}{r}\bigotimes_{j}\sigma_{j}^{\nu}}$ required to meet some
desired precision $\epsilon$. We write gates in this form to highlight that
the Hamiltonian summands in our simulations are always written as a tensor
product of Pauli operators, allowing for a nice parallel to the well known
rotation gates in the context of quantum computing. For imaginary time,
although this form still holds in our numerics, the same analogy does not
hold, as we are simply constructing $e^{-\frac{\beta}{r}H_{j}}$ on a classical
computer which can be of somewhat arbitrary form.

In order to calculate $C$, we have constructed a library to compile any
desired product formula simulation, given a list of Hamiltonian terms,
partition, QDrift samples, simulation time and desired precision. These
simulator objects are handled by external functions that can partition the
simulators, calculate errors and exact costs, or approximate simulation cost
via Monte-Carlo methods and more. The library is built on NumPy, but also
contains conversion functions to load Hamiltonians generated by quantum
chemistry packages OpenFermion [24] and PySCF [32] to simulate systems of
interest. It also contains methods for geometrically local simulations to
compute blockings using the required set logic.

###  V.1 Hamiltonians of Interest

Within this section, we introduce the Hamiltonians in our numerical
investigation of composite algorithms, as well as briefly outline methods used
for their generation.

####  V.1.1 Electronic Structure Hamiltonians in Second Quantization

The electronic structure problem is perhaps one of the most famous classically
intractable problems that has vast applications in quantum chemistry. In order
to write down these Hamiltonians, it is first necessary to introduce Fermionic
creation and annihilation operators. Fermions are particles with half-integer
spin that obey Fermi-Dirac statistics, meaning they obey the following anti-
commutation relations:

$\displaystyle\\{a_{m},a_{n}^{\dagger}\\}:$
$\displaystyle=a_{m}a_{n}^{\dagger}+a_{n}^{\dagger}a_{m}=\delta_{mn}$ (41)
$\displaystyle\\{a_{m},a_{n}\\}$
$\displaystyle=\\{a_{m}^{\dagger},a_{n}^{\dagger}\\}=0$ (42)

We work in Fock space where the subscripts of operators indicate the
excitation number or the atomic orbital of an electron. Molecular electronic
structure Hamiltonians then take the following form:

$H=\sum_{mn}h_{mn}a_{m}a_{n}^{\dagger}+\sum_{pqrs}h_{pqrs}a_{p}a_{q}a_{r}^{\dagger}a_{s}^{\dagger}$
(43)

Here the first term represents single excitations and the second term keeps
track of double excitations or “hopping” amongst orbitals. The coefficients
$h$ are molecular integrals that depend on the basis of choice to describe the
molecule [26]. Further, more can be done to aid in the implementation of this
problem on a quantum computer. The Jordan-Wigner transformation provides a
one-to-one mapping between the fermionic and spin operators. This will allow
us to write down the time evolution operator in terms of the universal
rotation gates (and CNOT gates). To understand the transformation, first
observe:

$\displaystyle a^{\dagger}$ $\displaystyle=\begin{bmatrix}0&0\\\
1&0\end{bmatrix}=\frac{\sigma^{x}-i\sigma^{y}}{2}:=\sigma^{-}$ (44)
$\displaystyle a$ $\displaystyle=\begin{bmatrix}0&1\\\
0&0\end{bmatrix}=\frac{\sigma^{x}+i\sigma^{y}}{2}:=\sigma^{+}$ (45)

Now to build in the desired commutation relations and generalize this to a
Hilbert space for $N$ qubits, or the tensor (Kronecker) product of $N$
2-dimensional Hilbert spaces
$\mathscr{H}=\bigotimes_{i=1}^{N}\mathscr{H}_{i}$:

$\displaystyle a_{n}^{\dagger}$ $\displaystyle\Leftrightarrow\openone^{\otimes
n-1}\otimes\sigma^{-}\otimes(\sigma^{z})^{\otimes N-n-1}$ (46) $\displaystyle
a_{n}$ $\displaystyle\Leftrightarrow\openone^{\otimes
n-1}\otimes\sigma^{+}\otimes(\sigma^{z})^{\otimes N-n-1}$ (47)

In order to build these Hamiltonians numerically, we use the OpenFermion [24]
and PySCF [32] packages for quantum chemistry. OpenFermion is a library that
allows for the easy manipulation of fermionic operators that arise in quantum
chemistry, as well as it interfaces with a variety of electronic structure
packages that perform molecular integrals in the basis of choice to generate
Hamiltonians in the form of Equation 43. Further, OpenFermion also has the
Jordan Wigner transform built in, allowing one to construct this Hamiltonian
in the Pauli basis. PySCF was our electronic structure package of choice to
compute molecular integrals.

Given the form of Equation 43, we observe that our Hilbert space needs to be
truncated. An active space calculation does exactly this; the Hamiltonian is
written in a space such that only so many orbitals are “active” or such that
an electron can be excited to occupy active orbitals. We generate all of our
electronic structure Hamiltonians in the minimal basis where we use a the
number of qubits equal to the total period of the molecule. For our numerical
investigation, we provide a function to generate chains of hydrogen atoms
given a very simple input; bond length and number of atoms. The function uses
PySCF to compute the molecular integrals, and then uses the data to build the
Hamiltonian in an active space implied by the minimal basis, and using a
minimal spin configuration.

####  V.1.2 Jellium Uniform Electron Gas

Jellium is a model of a uniform electron gas that captures the interactions
between delocalized electrons in a solid with uniformly distributed positive
potentials serving as Nuclei. It is not only a system of interest in Materials
Science, but also as a benchmark system in quantum simulation. More compact
representations of this Hamiltonian have been proposed as a candidate for
experimental simulation on near-term hardware [7]. The system Hamiltonian has
a closed form representation and does not require any additional molecular
integrals to construct:

$H=\frac{1}{2}\sum_{p,\sigma}k^{2}_{p}a^{\dagger}_{p,\sigma}a_{p,\sigma}-\frac{4\pi}{\Omega}\sum_{p\neq
q,j,\sigma}\left(\zeta_{j}\frac{e^{ik_{q-p}\cdot
R_{j}}}{k^{2}_{p-q}}\right)a^{\dagger}_{p,\sigma}a_{q,\sigma}+\frac{2\pi}{\Omega}\sum_{(p,\sigma)\neq(q,\sigma^{\prime}),\nu\neq
0}\frac{a^{\dagger}_{p,\sigma}a^{\dagger}_{q,\sigma^{\prime}}a_{q+\nu,\sigma^{\prime}}a_{p-\nu,\sigma}}{k_{\nu}^{2}}$
(48)

where the $j$th nuclei has position $R_{j}$ and atomic number $\zeta_{j}$, and
$k_{\nu}=\frac{2\pi\nu}{\Omega^{\nicefrac{{1}}{{3}}}}$ with cell volume
$\Omega$ and $\sigma$ containing both up and down spins. For the derivation of
this Hamiltonian see Appendix B of [7]. Conveniently, OpenFermion also
provides simple functions to quickly generate this Hamiltonian, and we do so
in the momentum plane wave basis with periodic boundary conditions. We elect
not to use the more compact plane wave dual basis representation presented in
[7], due to the fact that we are using this Hamiltonian as a benchmark, rather
than studying the outputs of the simulation. For the composite simulation,
Jellium provides many Hamiltonian terms and a very sharply peaked distribution
(see Figure for a system of size equal to that of the spin models we study.
Given that system size is more of a limiting factor than term number in our
numerical study, this presents an opportunity to see how a composite channel
performs on a system with greater $L$. To limit the system size we also use a
spinless model, and then perform the Jordan-Wigner transformation on the
second quantized Hamiltonian to represent our Hamiltonian as a sum of Pauli
operators. This Hamiltonian is constructed with the necessary transformations
using OpenFermion [24].

(a)

(b)

Figure 3: Jellium Spectral Norm Distribution: Semi-log plots of the sorted
normalized spectral norms versus Hamiltonian index for 5 and 7 site Jellium
models in figures (a) and (b) respectively. The plots show the increases in
number of terms as well as how the distributions become increasingly sharply
peaked. In red we provide a potential choice of $\omega_{c}$ for the
partitioning heuristic.

####  V.1.3 Graph Hamiltonian Model

The Hamiltonian we explore here involves a spin chain imposed on a lattice
$\Gamma\in\mathbb{Z}^{D}$ with a graph distance metric
$dist(\mathbf{u},\mathbf{v})=|\mathbf{u}-\mathbf{v}|_{1}$ where $\mathbf{u}$
and $\mathbf{v}$ are vectorized coordinates on the graph with dimension equal
to $D$. For our investigation we only examine lattices with $D=1$, given that
for a fixed number of sites, this gives the most sharply peaked spectral
distribution than any other $D$:

$H=\sum_{i>j}e^{-dist(i,j)}\alpha_{ij}\sigma^{x}_{i}\sigma^{x}_{j}+\sum_{k}\beta_{k}\sigma_{k}^{z}$
(49)

This system is similar to the quantum transverse field Ising model but with
interactions that fall of exponentially with graph distance. The coefficients
$\alpha_{ij}$ and $\beta_{k}$ are site dependant coupling constants that allow
for the introduction of more disorder and/or structure in the Hamiltonian. To
add some disorder to the model, we draw these coefficients pseudo-randomly
from a Gaussian distribution with mean 0 and variance 1.

####  V.1.4 Heisenberg Model

The Heisenberg model describes a quantum spin system in a magnetic field with
nearest-neighbour interactions. The Hamiltonian takes the following form:

$H=\sum_{j}\left(J_{x}\sigma_{j}^{x}\sigma_{j+1}^{x}+J_{y}\sigma_{j}^{y}\sigma_{j+1}^{y}+J_{z}\sigma_{j}^{z}\sigma_{j+1}^{z}\right)+\sum_{i}B_{z}\sigma_{i}^{z}$
(50)

Here $B_{z}$ is the strength of the magnetic field in the $z$ direction and
$J_{\\{x,y,z\\}}$ are coupling constants. Given the intuition of our composite
channel, we expect this model to take advantage of our algorithm when the
coupling constants largely differ in magnitude such that partitioning into
Trotter and QDrift takes advantage of more Hamiltonian structure. Furthermore,
introducing site-dependant coupling constants or writing down a highly
disordered spin system could further add structure that the algorithm can take
advantage of. A Hamiltonian of this nature would look something like the
following:

$H=\sum_{j}\left(J_{x}^{(j)}\sigma_{j}^{x}\sigma_{j+1}^{x}+J_{y}^{(j)}\sigma_{j}^{y}\sigma_{j+1}^{y}+J_{z}^{(j)}\sigma_{j}^{z}\sigma_{j+1}^{z}\right)+\sum_{i}B_{z}\sigma_{i}^{z}$
(51)

This Hamiltonian appears somewhat contrived for extracting the performance of
the local composite channel. However, this system closely resembles the
Edwards-Anderson model of a spin glass, a system of interest in condensed
matter physics. In an attempt to create a sharp distribution, we simply sample
$J_{\nu}^{(j)}$ from an exponential distribution with a scale parameter of
0.1.

###  V.2 Partitioning Schemes

The main difficulty with deploying the composite simulation framework concerns
finding a good partitioning. In introducing composite simulations, Hagan and
Wiebe suggested partitioning schemes derived on the basis of optimizing
analytic cost functions both in deterministic and probabilistic settings [22].
Here, we take a different approach involving the exact calculation of the
simulation error, and an optimization routine that, given convergence, finds
the optimal partitioning and gate count with respect to a chosen error
tolerance and simulation time. This approach is used to answer the question
regarding the best savings one can hope to achieve when deploying composite
methods to simulate a specific Hamiltonian. We are not proposing this as a
pre-processing routine, as it has complexity greater than that of the
simulation itself, which is trivial as the optimization involves solving the
simulation problem recursively. In addition, we also arrive at simple
heuristics that can be used to partition certain Hamiltonians with little
overhead, which we do propose as a strategy for using a composite approach.

Chop is the partition that we introduce in this work. The idea is based on the
heuristic of placing a few terms with larger spectral norms into Trotter-
Suzuki channels and numerous small terms into QDrift, assuming that the
Hamiltonian presents this structure. We start by sorting the terms by their
spectral weights and introducing a “chop threshold”
$\omega_{c}\in[0,\max_{i}h_{i}]$. This scale will determine the partition such
that if a term has spectral norm $h_{i}\geq\omega_{c}$ then $H_{i}\rightarrow
A$ if $h_{i}<\omega_{c}$ then $H_{i}\rightarrow B$. Now we can express the
error tolerance $\epsilon$ as a function of channel iterations $r$, with
partitioning chop $\omega_{c}$ and a sample number $N_{b}$ that will be chosen
in an optimization routine, and for a fixed initial state and time:

$\norm{(\mathcal{X}^{2k})^{\circ
r}(\rho,t/r,\omega_{c},N_{b})-\mathcal{U}{(\rho,t)}}_{1}=\epsilon(r).$ (52)

By fixing an error tolerance for $\epsilon$, the exact cost of the simulation
becomes a black box function with no closed form expression:

$C_{comp}=f(\epsilon_{thresh},r,\omega_{c},N_{b}).$ (53)

This is the cost function we wish to minimize. However, we cannot do that by
conventional methods such as with direct gradient descent. Also, with no
strong intuition for a choice of $N_{B}$, if we wish to optimize this
parameter we have to deal with integer optimization as well. The iterations
$r$, however, while an integer, does not require optimization, but rather
emits a search problem. If we allow an optimizer to pick initial random values
for $N_{b}$ and $\omega_{c}$ from a fixed interval, then we must find the
value $r$ required to meet the error threshold $\epsilon_{thresh}$, which will
ultimately be determined by the optimizer’s choice of the other two
parameters. To complete this, we perform an exponential search on $r$ until we
find some $r$ where $\epsilon(r)\geq\epsilon_{thresh}$ and set this as an
upper bound on $r$. We then perform a binary search to find the smallest value
of $r$ required to meet this condition and count the number of gates in the
channel. This is a very expensive function given that we are precisely
building the composite channel, applying it to the density matrix initial
state in the problem, and counting the gates applied in each iteration of the
search. The expensive nature emerges due to the sheer number of matrix
multiplications required in performing this task, not in the search for $r$,
which is nearly optimal. Note the importance of using the trace distance in
this approach as it guarantees monotonicity of $\epsilon(r)$, which makes the
search possible. This is not so in the framework of sampling the quantum
infidelity, as finding the cost here would require other statistical methods
(see Appendix B.

Now a glaring question left unanswered is the choice of an optimizer. We
implement the Gradient-Boosted Regression Trees (GBRT) algorithm included in
Sci-Kit Optimize [33]. This algorithm is specifically-designed to handle the
optimization of very expensive functions. It is also convenient for our
purposes given that it can handle both integer and real optimization
parameters simultaneously. At a high level, the algorithm works by using a
series of decision trees with an associated loss function. The decision trees
perform regression to fit the input function and are iteratively generated
based on the minimization of the loss function via gradient descent. This
optimizer and cost function 53 can then be easily generalized to the local
composite channels where now we have an $N_{b}$ and $\omega_{c}$ for each
blocking. As the number of local blocks grows, the optimization routine will
need to take a larger number of input parameters in this prescription.
However, the size of the system becomes classically intractable long before we
would consider using this many local blocks, so this is far from a concern.

In some cases, models may exhibit a partitioning that is somewhat canonical
and can lead to excellent performance of composite methods. This occurs when
we have a Hamiltonian that fits naturally into the intuition behind the
algorithm, such that we have a set $A$ containing large terms with small
commutators and a set $B$ with small terms that are known not to commute in
general. We are, therefore, proposing to use the chop partition but by
choosing the chop threshold $\omega_{c}$ based on physical intuition regarding
the Hamiltonian, rather than some expensive optimization routine. A perfect
example of such a system is a Heisenberg model with weak coupling. In this
case, looking at Equation 50, we would set the chop threshold
$\omega_{c}=\max\\{J_{x},J_{y},J_{z}\\}$, which implies we simulate the
interactions with QDrift
$\\{J_{\nu}\sigma_{j}^{\nu}\sigma_{j+1}^{\nu}\\}_{\nu=x,y,z}\rightarrow B$,
and simulate the site energy terms with Trotter-Suzuki
$\\{B_{z}\sigma_{j}^{z}\\}\rightarrow A$. In this way, the terms in the set
$A$ all commute with each other, whereas the terms in the set $B$ are
guaranteed to have a small spectral norm. We bring numerical evidence that
this provides computational advantages in the sections below. In general, any
system with perturbative interactions may benefit from this framework, given
that the commutators within the system are small, as they will avail this
canonical partitioning. In cases where the partitioning is not as obvious, as
is the case with H3 and Jellium, we can achieve similar advantages by choosing
$\omega_{c}=\max\frac{d\norm{H_{j}}}{dj}\textbf{s.t.}\norm{B}\geq\norm{A}$,
meaning we sweep an ordered list of the Hamiltonian spectral norms and track
the largest difference between terms, chopping the list where this occurs,
given $j\geq\frac{L}{2}$. The final condition is just to ensure the majority
of terms are simulated by QDrift. We also use this strategy throughout V.4 and
show advantages.

###  V.3 Error Measures

In this investigation, there is some arbitrariness in the error measure one
can choose in order to quantify the performance of a simulation channel. In
order to evaluate the resources required by an algorithm, one must evaluate
the number of gates required to meet a certain $\epsilon$, which is calculated
by said error measure. In the literature, this $\epsilon$ is often quantified
by the diamond distance utilized in previous sections. However, while
analytically convenient, for any reasonably-sized system, computing this
quantity becomes computationally expensive. While it is possible to evaluate
it efficiently, this requires finding the solution of a semi-definite program,
which is much less efficient than using some other error measures. In
addition, since we are not constrained to analytically solvable expressions or
closed form equations with our numerical methods, we can optimize this cost in
terms of some partition scheme. This is the idea behind the optimal chop
partition, and doing so requires frequent computations of $\epsilon$. With
this in mind, the error of our algorithm should be a quantity that we can
compute in a reasonable amount of time while also being a fair error measure.
The criteria for “fairness” comes within the error measures treating each
algorithm that comprises the composite channel on an equal footing. For
example, if we are to optimize the partitioning with respect to the gate cost
(which is dependent on $\epsilon$), then if Trotter is more performant with
respect to QDrift in one error measure than in another, then our composite
optimizer will favor Trotter, which will be reflected in the partition. As a
result the total cost of the composite channel will be skewed by the error
measure used.

We consider the infidelity and trace distance as possible measures of
$\epsilon$. For definitions and an additional discussion regarding the scaling
and complexity of computing these quantities, see Appendix B. Analytics are
required in order to answer the question of which measure might provide a more
fair comparison. More specifically, we can ask the question of how the error
measures will scale with respect to the total simulation time $t$, and then
test these results numerically. In terms of the infidelity, we provide the
following Theorems:

######  Proposition V.1 (QDrift Infidelity Time Scaling).

Given a QDrift channel $\mathcal{Q}{(\rho,t)}$ and the standard evolutionary
channel $\mathcal{U}(\rho,\nicefrac{{t}}{{N}})$, for a density matrix $\rho$,
time $t/N$, then the infidelity between the outputs of the channels
$1-\mathscr{F}(\mathcal{Q}{(\rho,t)},\mathcal{U}(\rho,t/N))\in\mathcal{O}\left(t^{2}\right)$.

######  Proposition V.2 (Trotter-Suzuki Infidelity Time Scaling).

Given a Trotter-Suzuki channel $\mathcal{T}^{2k}{(\rho,t)}$ and the standard
evolutionary channel $\mathcal{U}{(\rho,t)}$, for a density matrix pure state
$\rho$ and time $t$, then the infidelity between the outputs of the channels
$1-\mathscr{F}(\mathcal{T}^{2k}{(\rho,t)}\mathcal{U}(\rho,t))\in\mathcal{O}\left((t^{2k+1})^{2}\right)$.

The proofs of these theorems are included in Appendix B. Here, we obtain the
non-trivial result that the infidelity is squared when considering Trotter-
Suzuki formulas. This would lead to our optimized chop algorithm to heavily
favor this channel over QDrift, and for this reason, we consider it an
”unfair” error measure. On the other hand, if we consider how the trace
distance scales with simulation time in both algorithms, we obtain the
following Theorems:

######  Proposition V.3 (QDrift Trace Distance Time Scaling).

Given a QDrift channel $\mathcal{Q}{(\rho,t)}$ and the standard evolutionary
channel $\mathcal{U}(\rho,\nicefrac{{t}}{{N}})$, for an arbitrary density
matrix $\rho$ and time $t/N$, then the trace distance between the outputs of
the channels
$\mathscr{T}(\mathcal{Q}{(\rho,t)},\mathcal{U}(\rho,t/N))\in\mathcal{O}\left(t^{2}\right)$.

######  Proposition V.4 (Trotter-Suzuki Trace Distance Time Scaling).

Given a Trotter-Suzuki channel $\mathcal{T}^{2k}{(\rho,t)}$ and the standard
evolutionary channel $\mathcal{U}{(\rho,t)}$, for an arbitrary density matrix
$\rho$ and time $t$, then the trace distance between the outputs of the
channels
$\mathscr{T}(\mathcal{T}^{2k}{(\rho,t)},\mathcal{U}(\rho,t))\in\mathcal{O}\left(t^{2k+1}\right)$.

Here, we see that no such squaring occurs, and the expected time-scaling is
obtained. For this reason, we compute the entire density matrix and $\epsilon$
using the trace distance in all of our numerical simulations. Proofs of the
above theorems, as well as further discussions can again be found in Appendix
B.

###  V.4 Performance Results

In this section, we first numerically analyze the real time quantum algorithm
given by Hagan and Wiebe [22] and then show equivalent numerical calculations
of the imaginary time classical case. To accomplish this, we provide cost
plots in which we provide the minimum $C$, or the number of rotation gates to
achieve a desired simulation accuracy $\epsilon$ (calculated by the trace
distance) for each point in time $t$ or $\beta$. To reiterate, here we exactly
compute entire evolution channels with a random initial state $\rho$ sampled
from the unit hyper-sphere, and directly apply and count gates. We conclude
this section with a brief discussion about numerical studies for the local
composite simulation algorithms. In these plots, we study variants of the
composite channel and display results with the aforementioned notation with
the addition of a tilde over the channel if the partition and $N_{B}$ have
been optimized with GBRT. For example, a composite channel with inner-order 2
and outer-order 1 with an optimized partition and number of QDrift samples
$N_{B}$ is written like so $\widetilde{\mathcal{X}}^{2,1}$.

Throughout the section, we normalize $\norm{H}=1$ and run simulations for
times $t\in(0,\frac{3\pi}{2}]$ so as to ensure the system undergoes non-
trivial dynamics without overlapping the phase. This is done due to the fact
that Trotter formulas have a periodic error for $\norm{H}t\geq 2\pi$, and
running simulations in this range would lead to the optimizer finding the
“good points”, where the error happens to be small, which would provide a very
low cost simulation and a sharp drop in the cost trend. We also report the
cost advantages achieved on crossover points, which are values of $t$ such
that $C_{QD}(t)=C_{TS}(t)$. We denote the composite channel advantage as
$\xi\coloneqq
C_{QD}(t^{\prime})/C_{\mathcal{X}}(t^{\prime})=C_{TS}(t^{\prime})/C_{\mathcal{X}}(t^{\prime})$.
As we are unable to compute these times exactly we use interpolation methods
to report values of $\xi$. This is motivated by the fact that analytics
suggest this to be the point of greatest advantage for a composite channel 22.
This is intuitive, especially given higher-order Trotter-Suzuki formulas which
are known to asymptotically outperform QDrift for large $t$, whilst QDrift is
dominant in the small $t$ limit, suggesting a region where their strengths can
be combined.

Hamiltonian | $\xi$ | Num. Terms | Time
---|---|---|---
Hydrogen-3 | 2.3 | 62 | Real
5 Site Jellium | 9.2 | 56 | Real
6 Site Jellium | 18.8 | 94 | Real
7 Site Jellium | 10.4 | 197 | Real
7 Spin Graph | 4.1 | 49 | Real
8 Spin Graph | 3.9 | 64 | Real
8 Spin Heisenberg | 3.1 | 29 | Imag.
Hydrogen-3 | 2.3 | 62 | Imag.
6 Site Jellium | 18.8 | 94 | Imag.
Table 1: Summary of gate cost improvements observed (contingent on
optimization convergence). We observe that savings tend to somewhat improve as
the number of terms increases (within the same model), with the exception of
Jellium 7 where optimizer struggles with partitioning due to the number of
terms. This is evident in the lack of monotonicity of
$C(\widetilde{\mathcal{X}^{1}}$ in Figure 6. The most significant savings are
seen for the Jellium models. Even in cases where the number of terms are
comparable to other models, larger advantages are persistent in Jellium. This
is further establishes the spectral norm distribution as one of the most
important indicators of performance in the composite framework.

####  V.4.1 Real-Time Composite Quantum Simulations

Beginning with Hydrogen chains, our results are highlighted in Figure 4. This
plot reveals two interesting features: for long-time simulations, heuristics
can be found that essentially match the performance of an optimized formula
simply by inspecting the distribution of the norms of individual summands
$\norm{H_{j}}$, and with the optimized routine, we find a significant
improvement at the crossover point with $\xi=2.3$. These plots begin flat for
most of the simulation channels, which for the most part, indicates that one
application of the channel achieves the desired $\epsilon$ for multiple
sequential simulations at small times. This is expected, and is especially
common with Trotter-Suzuki formulas, given that with one iteration they apply
at least $L$ gates depending on the order, while QDrift provides the option of
sampling single gates.

Figure 4: H3 Cost Plot Simulation (Real-Time): Log-log gate cost plot for the
H3 Hamiltonian generated with OpenFermion using three-dimensional Gaussians in
a minimal basis. The bond distance is chosen to be that which minimizes the
energy surface of H2, which is $\approx 0.8$ angstrom. We achieve a crossover
advantage of $\xi=2.3$, as well as remaining constant factor advantages at
long times.

Figure 4 is interesting given that both the heuristic and optimized channels
provide significant advantages at the crossover point, as well as the
heuristic partition seems to match the asymptotic performance of the optimized
channel. This demonstrates that optimization subroutines are not required to
gain advantages in this framework. Furthermore, the second inner-order
composite channel also shows a consistent advantage over the second-order
Trotter channel.

To gain insight into effective choices of partition and QDrift samples, we
also plot the optimized $N_{B}$ values and the ratio of Trotter terms to total
terms $|A|/|H|$ with time in Figure 5. Here, we find choices that somewhat
agree with our prior intuition. For short times, places almost all terms into
QDrift, and slowly increases $N_{B}$. As $t$ increases, more terms are placed
into Trotter with the partition bouncing around in the regime where Trotter
and QDrift have similar costs, which is also expected. The composite channel
$\widetilde{\mathcal{X}}^{2,1}$ essentially places all terms into the Trotter
simulator, given the favorable asymptotic performance of higher order Trotter
formulas over QDrift, while $\widetilde{\mathcal{X}}^{1}$ finds a balance
between the two at long times, likely due to their equivalent $t$ scaling. The
most interesting behavior is that of $N_{B}$ at mid to long times. Here,
$N_{B}$ peaks near the crossover point and then falls off as Trotter
$t$-scaling becomes dominant in $\widetilde{\mathcal{X}}^{2,1}$. However, for
the $\widetilde{\mathcal{X}}^{1}$ channel, $N_{B}$ experiences somewhat of a
revival after the peak, and stabilizes at 15, which is about 24% of the terms.
We use this percentage to motivate future heuristic choices of $N_{B}$ in our
investigation of Jellium in Figure 6, which turns out to work quite well.

(a)

(b)

Figure 5: Optimized H3 Simulation Parameters: Semi-log plots of parameters
obtained by the GBRT optimization routine for the
$\widetilde{\mathcal{X}^{1}}$ and $\widetilde{\mathcal{X}}^{2,1}$ real time
channels. These parameter choices correspond to the H3 simulation in Figure 4.
In (a) we plot the cardinality of the $A$ set over the total number of terms,
as a function of time. These values are dictated by GBRT optimized value of
$\omega_{c}$. In (b) we present the equivalent plot with $N_{B}$.

When it comes to the simulation of Jellium, we find some of the most
significant performance improvements within this section, including an order
of magnitude cost difference at the crossover point (see Figure 6).
Specifically, in the case of 6-site Jellium, the Trotter and QDrift cost at
the crossover point is approximately 100 gates, versus the composite channel,
which achieves the same precision $\epsilon$ with only about 7 gates. Here, it
is also shown that one can find an adequate partition that leads to advantages
at longer times without the need for any optimization. This is also the only
model whereby the optimization routine struggles to find optimal partitions in
the neighbourhood of the crossover point. This leads to the Jellium 7 model
inheriting a smaller $\xi$ than what is likely achievable.

(a)

(b)

(c)

Figure 6: Jellium Simulation Cost plots (Real-time):  Log-log cost plots of
quantum simulations of Jellium with 5, 6, and 7 sites in (a), (b), and (c)
respectively. In (a) and (b) we have in red, a chop heuristic where no
optimization overhead is used. The distribution of Hamiltonian terms is
chopped immediately before $\max\frac{d\norm{H_{j}}}{dj}$, and approximately
$\frac{1}{5}L$ terms are sampled. This heuristic works quite well, although it
is outperformed by the optimized version, especially at short times. In (a) we
achieve $\xi=9.2$. In (b) we achieve an impressive advantage of $\xi=18.8$,
the largest of all our real time results. In (c), we perform a similar
analysis of the 7 site model with some higher order composite channels, but
find that the optimizer has increased difficulty with larger numbers of
Hamiltonian terms. Here $\xi=10.4$, but through inspecting some neighbouring
points of the crossover region, it surely has the potentially to be much
larger.

The final system we investigate in this section is that of the graph toy model
with exponentially decaying interactions, which is a beyond nearest neighbour
model. In Figure 7, we study this model for chains of length 7 and 8, and find
essentially identical behaviour. When moving from 7 to 8 spins, we only add 15
more terms to the Hamiltonian, which is clearly not enough to see any
significant advantages. In fact, the crossover advantage is slightly smaller
for the bigger model, but this could also be due to the optimizer not fully
converging.

(a)

(b)

Figure 7: Toy Spin Graph Model Cost Plot: In (a) the we have a cost plot for
the 7 spin model where we obtain a crossover advantage of $\xi=4.1$, which is
fairly significant. The plot has multiple regions where different composite
channels are optimal. In (b) we have the 8 spin model where we establish
$\xi=3.9$. This advantage, as well as the channel performance is almost
identical to (a).

####  V.4.2 Imaginary-Time Composite Channels

This section contains cost plots of simulations of the same aforementioned
Hamiltonians, but with imaginary time propagators, where we present cost plots
of the most interesting Hamiltonians from the previous section.

For simulations of the Heisenberg model, we find similar advantages to those
in real time. In Figure 8, we see that our proposed heuristic leads to an
advantage over Trotter-Suzuki and QDrift in the regions of interest. What is
different about this plot is that the optimizer finds the same $N_{B}$ and
partitioning for all short times. This is an artefact of both the Hamiltonian
and how the optimizer is programmed. Since all the splitting (single site)
terms have equal spectral norm, the optimizer is placed in an all or nothing
scenario, as choosing $\omega_{c}<1$ immediately places all terms into QDrift.
Given that after receiving the cost, the program then solves a search problem
to find the minimal $r$ to achieve $\epsilon$ precision, it is rare that it
finds the ideal conditions to build a pure QDrift channel with $r=1$. However,
significant savings are still achieved at the crossover point, and over
$\mathcal{T}^{1}$ and $\mathcal{T}^{2}$ at large $\beta$. Once again, the
validity of heuristic partitions are shown, specifically in the 1st order
composite channel, which exactly matches its optimized version at large
$\beta$.

Figure 8: 8 Spin Heisenberg Model Cost Plot: In this imaginary time simulation
we establish $\xi=3.1$, as well as maintain advantages at large $\beta$. We
also find that our chosen heuristics are essentially optimal at large $\beta$,
with the green and dark blue lines overlapping.

For Hydrogen chains, we obtain strikingly similar results and compared with
those in real time, as seen in Figure 9. We once again obtain a significant
crossover advantage, as well as constant factor advantages at large $\beta$,
or low temperature. Heuristics are also shown to continue to hold in their
effectiveness, in this case, from the crossover point and onward.

Figure 9: H3 Simulation Cost Plot (Imaginary-Time): The same parameters are
used to build the 3 atom Hydrogen chain as we done in real time. We achieve
incredibly similar results, and again recover the real time result of
$\xi=2.3$ that was achieved in Figure 4.

For Jellium, we choose to investigate the system size with the best-behaved
optimizer, as well as the largest $\xi$, which occurs for the 6-site model. In
imaginary time, we once again reproduce a significant advantage, shown in
Figure 10. As in the case with H3, this plot is quite similar to the real time
case in Figure 6. However, here at large $\beta$ the composite channel seems
to do better in imaginary time given that even the first order composite
channel (with optimization) outperforms the second order Trotter channel. This
happens in the final point of the plot where $\mathcal{T}^{2}$ is no longer in
the“flat-regime”. While this is very interesting, it is unclear analytically
why this occurs, and we would likely not expect this trend to continue
asymptotically.

Figure 10: 6 Site Jellium Simulation Cost Plot (Imaginary-Time): Here we
recover the crossover advantage of $\xi=18.8$ from the real time simulation in
Figure 6, which is also the largest advantage achieved in our imaginary time
simulations. We additionally achieve advantages over second order Trotter at
large $\beta$.

Overall, this section nicely complements some of the analytics in Section III
both by reinforcing the fact that composite quantum channels allow for similar
advantages in both real and imaginary time, as well as through calculation of
exact constant factor advantages. In other words, this section provides
convincing evidence on the applicability of composite formulas to classical
imaginary-time Monte-Carlo algorithms.

####  V.4.3 Local Composite Quantum Channels

Distinct from the previous two algorithms, with the local composite channels
we do not immediately expect to see significant simulation advantages for the
small systems we can compute. Recall, this algorithm makes use of the Lieb-
Robinson velocity $v_{LR}$ that limits the propagation of information and thus
correlations in a local lattice model. In our above numerics, our lattices
contain $\leq 8$ sites, meaning even for small $v_{LR}$, the lattice can still
become quickly entangled. In this section, it is important to understand where
the observed advantages are originating, whether they are from the local
decomposition, or something else. For example, the results in prior sections
already suggest that composite channels can outperform Trotter and QDrift
channels in certain regimes. If a block decomposition is introduced, we pay a
small gate cost to break the simulation into subsets (as gates on the boundary
are applied more than once), but the advantages of the composite simulation
are almost guaranteed to outweigh this cost. Thus we wish to find a regime in
which we can perform calculations of exact costs with local composite channels
that explicitly gain advantages via block decomposition. Otherwise, we will
observe essentially the same behaviour as before, but with slightly smaller
constant factors. There are two ways to go about achieving this; one is to add
more sites to the model, which quickly becomes computationally intractable
with standard methods. The second strategy is to decrease the coupling between
sites in the lattice, which naturally decreases $v_{LR}$. This is the strategy
we utilize. In reference to Equation 50 we perform our cost plots on
Heisenberg models with 8 sites with $B_{z}=1$, and $J_{\nu}^{(j)}$ sampled
from an exponential distribution with a scale parameter (serving as a coupling
constant) of 0.00005. To allow for a fair simulation, we then choose
$\epsilon=0.000001$, such that statistically, $~{}98\%$ of terms will be
greater than $\epsilon$, which can be seen from a simple integration of the
PDF. Results of this simulation are shown in Figure 11.

Figure 11: 8 Spin Heisenberg Model Simulated by localized and standard
channels: This cost plot compares Trotter to its localized versions, with the
subscript $l$ indicating the overlap of the boundary region in the block-local
simulation. Channels with no $l$ subscript are the standard simulations from
prior sections. As expected, we observe that the composite channel is most
efficient, however, given that the standard Trotter algorithm outperforms the
localized Trotter channel, we can conclude that we are not in the regime where
locality is providing advantages. The same heuristics were used for the
composite channel as in Figure 8, with $N_{B}=(4,1,4)$ on lattice subsets
$(A,Y,B)$.

Given the gap between $\mathcal{T}^{1}$ and $\mathcal{T}_{l=2}$ in this very
weak coupling regime, it is unclear whether our methods (exact gate counts)
provide the means for investigating advantages gained by locality. In Ref.
[23], via computations of bounds, numerics did not show advantages (in the
form of T-gates) until approximately 100 sites were included in a more
strongly-coupled model with $J_{\nu}^{(j)}\in[-1,1]$ sampled i.i.d., so our
results with far fewer spins are not unexpected. However, we still
theoretically expect to see advantages in the limit where system sizes are
large, and we can take advantage of the Lieb-Robinson bound.

##  VI Discussion and Summary

The main contributions of this work are the extension of the composite channel
to imaginary time and to local Hamiltonians, and the construction of a
composite simulation library for numerically evaluating algorithmic
performance. The imaginary-time bounds provided advocate for the potential use
of QDrift and Composite QDrift-product formulas in classical simulations where
they can potentially improve the efficiency/accuracy of Quantum Monte-Carlo
simulations. This is possible via shortening the path length of discrete time
path integrals in these calculations given that our algorithms can approximate
$e^{-\beta H}$ with fewer propagators than Trotter-Suzuki formulas in most
cases. Our results also highlight that depending on the choice of $\beta$, one
may also perform an imaginary time simulation using QDrift, particularly for
high temperature physics. It is an open question as to whether these methods
can be used to improve the overall accuracy before one runs into the sign
problem. Nevertheless, these methods are interesting as one can use physical
information about the Hamiltonian to choose a partitioning that possibly
grants a computational advantage.

Furthermore, the extension of the composite channel to local Hamiltonians
provides an interesting approach to lattice simulations in which there is a
priori knowledge about the lattice structure. The diamond distance bound on
this algorithm is significant in that the cost of establishing local composite
blocks seemingly minimal in terms of asymptotics, making it a potentially-
useful algorithm. However, it seems more powerful numerical methods may be
required to explore the regimes in which this approach shines, given the size
of the systems likely required. An interesting future trade-off to be explored
here is the size of the blocking $l$, and the size of the commutator sums in
Theorem 9. Furthermore, one can investigate optimizations of partitions and
QDrift samples on lattice subsets. Despite our library containing these
methods, we did not employ them as they can lead to misleading figures. Here,
one must tread cautiously, as adding more blocks with parameterized weights
and QDrift samples can lead to an optimizer finding lower costs, similar to
how adding hidden units to a neural network can lead to the over-fitting of a
data set in machine learning. In this scenario, it may be possible to find
localized composite channels that outperform the standard channels, simply
because the increase in parameters allowed us to find a more tailored
partitioning.

One of the difficulties with the techniques introduced in Ref. 22 is finding
good parameters to yield cost improvements. Our numeric results seem to
suggest that even heuristic approaches for parameters, such as $N_{B}$, can
yield at least factors of 2 reductions for small Hamiltonians, as seen in Fig.
4. The difficulties of multiple parameter optimization becomes more of an
issue when moving to the local simulation framework. With $m$ blocks there are
at least $2m$ parameters that need to be determined with a chop partitioning
scheme. However, simple heuristics still lead to performance advantages over
localized Trotter-Suzuki formulas in Figure 11. The fact that said heuristics
provide advantages with only an 8 spin model suggests that cost savings should
scale favorably with increased lattice sizes.

In Section V.3, a series of short but significant proofs were also provided to
remind the reader to be cautious in examining error within hybridized
algorithms. In our case, the infidelity treated QDrift differently than
Trotter by a square, and using it would have skewed the partitions and
provided inaccurate costs on the composite channel.

Finally, the library built for this project is easy to implement. It is also
flexible, and it reliably evaluates the exact number of operator of
exponentials $e^{-iH_{j}t}$ required to execute on a quantum computer (or
$e^{-iH_{j}\beta}$ for classical) to time evolve a system within a given error
tolerance. To reiterate, the reason this is important is because the cost of
these algorithms is highly dependent on the Hamiltonian, specifically in the
number and size of the terms, and their commutators. Additionally, the costs
and errors discussed in this paper are derived from analytic upper bounds on
the diamond distance and induced Schatten $1\rightarrow 1$ norm, which are
often bounded with repeated applications of the triangle inequality and using
$\max$ operations, and can therefore be loose. These exact numerics can
capture the true costs and errors of these algorithms, and advantages can
motivate further studies into these approaches. This library is available upon
reasonable request to the authors.

##  VII Acknowledgements

M.P. was supported through NSERC discovery grants, N.W. and M.H. were
supported using grants from US DOE National Quantum Information Science
Research Centers, Co-design Center for Quantum Advantage (C2QA) under contract
number DE-SC0012704 as well as support from Google Inc. D.S. acknowledges
support from an NSERC Discovery Grant and the Canada Research Chair program.
JC acknowledges support from the Natural Sciences and Engineering Research
Council (NSERC), the Shared Hierarchical Academic Research Computing Network
(SHARCNET), Compute Canada, and the Canadian Institute for Advanced Research
(CIFAR) AI chair program. Resources used in preparing this research were
provided, in part, by the Province of Ontario, the Government of Canada
through CIFAR, and companies sponsoring the Vector Institute
www.vectorinstitute.ai/#partners.

## References

  * [1] Richard P Feynman.  Simulating physics with computers.  In Feynman and computation, pages 133–153. CRC Press, 2018.
  * [2] Seth Lloyd.  Universal quantum simulators.  Science, 273(5278):1073–1078, 1996.
  * [3] A Yu Kitaev.  Quantum measurements and the abelian stabilizer problem.  arXiv preprint quant-ph/9511026, 1995.
  * [4] Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2002.
  * [5] Markus Reiher, Nathan Wiebe, Krysta M Svore, Dave Wecker, and Matthias Troyer.  Elucidating reaction mechanisms on quantum computers.  Proceedings of the national academy of sciences, 114(29):7555–7560, 2017.
  * [6] Bela Bauer, Sergey Bravyi, Mario Motta, and Garnet Kin-Lic Chan.  Quantum algorithms for quantum chemistry and quantum materials science.  Chemical Reviews, 120(22):12685–12717, 2020.
  * [7] Ryan Babbush, Nathan Wiebe, Jarrod McClean, James McClain, Hartmut Neven, and Garnet Kin-Lic Chan.  Low-depth quantum simulation of materials.  Physical Review X, 8(1):011044, 2018.
  * [8] Stephen P Jordan, Keith SM Lee, and John Preskill.  Quantum algorithms for quantum field theories.  Science, 336(6085):1130–1133, 2012.
  * [9] Dominic W Berry, Graeme Ahokas, Richard Cleve, and Barry C Sanders.  Efficient quantum algorithms for simulating sparse hamiltonians.  Communications in Mathematical Physics, 270:359–371, 2007.
  * [10] Andrew M Childs, Yuan Su, Minh C Tran, Nathan Wiebe, and Shuchen Zhu.  Theory of trotter error with commutator scaling.  Physical Review X, 11(1):011020, 2021.
  * [11] Earl Campbell.  Random compiler for fast hamiltonian simulation.  Physical review letters, 123(7):070503, 2019.
  * [12] Kianna Wan, Mario Berta, and Earl T Campbell.  Randomized quantum algorithm for statistical phase estimation.  Physical Review Letters, 129(3):030503, 2022.
  * [13] Yingkai Ouyang, David R White, and Earl T Campbell.  Compilation by stochastic hamiltonian sparsification.  Quantum, 4:235, 2020.
  * [14] Dominic W Berry, Andrew M Childs, Yuan Su, Xin Wang, and Nathan Wiebe.  Time-dependent hamiltonian simulation with $l1$-norm scaling.  Quantum, 4:254, 2020.
  * [15] Andrew M Childs and Nathan Wiebe.  Hamiltonian simulation using linear combinations of unitary operations.  arXiv preprint arXiv:1202.5822, 2012.
  * [16] Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D Somma.  Simulating hamiltonian dynamics with a truncated taylor series.  Physical review letters, 114(9):090502, 2015.
  * [17] Guang Hao Low and Isaac L Chuang.  Hamiltonian simulation by qubitization.  Quantum, 3:163, 2019.
  * [18] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204, 2019.
  * [19] John M Martyn, Zane M Rossi, Andrew K Tan, and Isaac L Chuang.  Grand unification of quantum algorithms.  PRX Quantum, 2(4):040203, 2021.
  * [20] Masuo Suzuki.  Quantum monte carlo methods—recent developments.  Physica A: Statistical Mechanics and its Applications, 194(1-4):432–449, 1993.
  * [21] WMC Foulkes, Lubos Mitas, RJ Needs, and Guna Rajagopal.  Quantum monte carlo simulations of solids.  Reviews of Modern Physics, 73(1):33, 2001.
  * [22] Matthew Hagan and Nathan Wiebe.  Composite quantum simulations.  arXiv preprint arXiv:2206.06409, 2022.
  * [23] Jeongwan Haah, Matthew B Hastings, Robin Kothari, and Guang Hao Low.  Quantum algorithm for simulating real time evolution of lattice hamiltonians.  SIAM Journal on Computing, (0):FOCS18–250, 2021.
  * [24] Jarrod R McClean, Nicholas C Rubin, Kevin J Sung, Ian D Kivlichan, Xavier Bonet-Monroig, Yudong Cao, Chengyu Dai, E Schuyler Fried, Craig Gidney, Brendan Gimby, et al.  Openfermion: the electronic structure package for quantum computers.  Quantum Science and Technology, 5(3):034014, 2020.
  * [25] Mario Motta, Chong Sun, Adrian TK Tan, Matthew J O’Rourke, Erika Ye, Austin J Minnich, Fernando GSL Brandão, and Garnet Kin-Lic Chan.  Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution.  Nature Physics, 16(2):205–210, 2020.
  * [26] James D Whitfield, Jacob Biamonte, and Alán Aspuru-Guzik.  Simulation of electronic structure hamiltonians using quantum computers.  Molecular Physics, 109(5):735–750, 2011.
  * [27] Masuo Suzuki.  Fractal decomposition of exponential operators with applications to many-body theories and monte carlo simulations.  Physics Letters A, 146(6):319–323, 1990.
  * [28] Naomichi Hatano and Masuo Suzuki.  Finding exponential product formulas of higher orders.  In Quantum Annealing and Other Optimization Methods, pages 37–68. Springer Berlin Heidelberg, nov 2005.
  * [29] Bruno Nachtergaele and Robert Sims.  Lieb-robinson bounds and the exponential clustering theorem.  Communications in mathematical physics, 265:119–130, 2006.
  * [30] Bruno Nachtergaele and Robert Sims.  Much ado about something: why lieb-robinson bounds are useful.  arXiv preprint arXiv:1102.0835, 2011.
  * [31] John Watrous.  The theory of quantum information.  Cambridge university press, 2018.
  * [32] Qiming Sun, Timothy C Berkelbach, Nick S Blunt, George H Booth, Sheng Guo, Zhendong Li, Junzi Liu, James D McClain, Elvira R Sayfutyarova, Sandeep Sharma, et al.  Pyscf: the python-based simulations of chemistry framework.  Wiley Interdisciplinary Reviews: Computational Molecular Science, 8(1):e1340, 2018.
  * [33] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  the Journal of machine Learning research, 12:2825–2830, 2011.

##  Appendix A Proofs of Imaginary Time Bounds

###  A.1 Imaginary time QDrift

See III.1

###### Proof.

Following the supplementary material of [11], we will make use of the
Liouvillian super-operator formalism for the imaginary time channel. In
imaginary time we have:

$\mathcal{L}(\rho)=\\{H,\rho\\}=H\rho+\rho H,$ (54)

which generates the following channel:

$\frac{e^{-\beta\mathcal{L(\rho)}}}{\Tr
e^{-\beta\mathcal{L(\rho)}}}=\frac{e^{-\beta H}\rho e^{-\beta
H}}{\Tr(e^{-\beta H}\rho e^{-\beta H})}.$ (55)

Throughout this section, whenever $\mathcal{L}$ appears, it is implied that we
mean $\mathcal{L}(\rho)$, we use the former for more compact notation. The
difference between $\mathcal{L}$ defined here and that in [11] is that the
commutator becomes an anti-commutator, and the $i$ vanishes due to the nature
of the imaginary time channel. Equivalence can be seen via the Taylor series
expansion of $e^{-\beta\mathcal{L}(\rho)}$, where $\mathcal{L}^{n}$ means $n$
applications of the super-operator:
$\\{H,\\{H,...,\\{H,\\{H,\rho\\}\\}...\\}\\}$ From [11], this super-operator
yields the following properties:

$\norm{\mathcal{L}}_{1\rightarrow 1}\leq 2\norm{H}\leq 2\lambda.$ (56)

Similarly, if we define $\mathcal{L}_{j}$ that generate unitaries under each
$H_{j}$ as in [11], and given that the Hamiltonian terms are normalized by
$h_{j}$ the following also holds:

$\norm{\mathcal{L}_{j}}_{1\rightarrow 1}\leq 2\norm{H_{j}}=2.$ (57)

We are now ready to proceed with the proof. To avoid cumbersome notation we
will write $\mathcal{L}(\rho)$ simply as $\mathcal{L}$. Let us first examine
the induced 1-norm of the difference between exact channel $\mathcal{U}_{N}$
that acts for imaginary time $\nicefrac{{\beta}}{{N}}$ and the QDrift channel
$\mathcal{Q}$, and the generalize this for $N$ compositions of the QDrift
channel representing $N$ samples:

$\displaystyle\norm{\frac{\mathcal{U}_{N}}{\Tr\mathcal{U}_{N}}-\frac{\mathcal{Q}}{\Tr\mathcal{Q}}}_{1\rightarrow
1}$ $\displaystyle=\norm{\frac{e^{-\nicefrac{{\beta}}{{N}}\mathcal{L}}}{\Tr
e^{-\nicefrac{{\beta}}{{N}}\mathcal{L}}}-\frac{\sum_{i}\frac{h_{i}}{\lambda}e^{-\nicefrac{{\beta\lambda}}{{N}}\mathcal{L}_{i}}}{\Tr(\sum_{j}\frac{h_{j}}{\lambda}e^{-\nicefrac{{\beta\lambda}}{{N}}\mathcal{L}_{j}})}}_{1\rightarrow
1}$ (58)
$\displaystyle=\norm{\frac{e^{-\nicefrac{{\beta}}{{N}}\mathcal{L}}}{1+\Tr(\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{N^{n}n!})}-\frac{\sum_{i}\frac{h_{i}}{\lambda}e^{-\nicefrac{{\beta\lambda}}{{N}}\mathcal{L}_{i}}}{1+\Tr(\sum_{j}\frac{h_{j}}{\lambda}{\sum_{m=1}^{\infty}\frac{(-\beta\mathcal{L}_{j})^{m}}{N^{m}m!}})}}_{1\rightarrow
1}$ (59)