# Variational Hamiltonian Simulation for Translational Invariant Systems via
Classical Pre-Processing

Refik Mansuroglu1, Timo Eckstein1,2, Ludwig Nützel1, Samuel A. Wilkinson1,
Michael J. Hartmann1,2<EMAIL_ADDRESS><EMAIL_ADDRESS>1Department of Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg
(FAU), Staudtstraße 7, 91058 Erlangen  2Max Planck Institute for the Science
of Light, Staudtstraße 2, 91058 Erlangen, Germany

###### Abstract

The simulation of time evolution of large quantum systems is a classically
challenging and in general intractable task, making it a promising application
for quantum computation. A Trotter-Suzuki approximation yields an
implementation thereof, where a higher approximation accuracy can be traded
for an increased gate count. In this work, we introduce a variational
algorithm which uses solutions of classical optimizations to predict efficient
quantum circuits for time evolution of translationally invariant quantum
systems. Our strategy can improve upon the Trotter-Suzuki accuracy by several
orders of magnitude. It translates into a reduction in gate count and hence
gain in overall fidelity at the same algorithmic accuracy. This is important
in NISQ-applications where the fidelity of the output state decays
exponentially with the number of gates. The performance advantage of our
classical assisted strategy can be extended to open boundaries with
translational symmetry in the bulk. We can extrapolate our method to beyond
classically simulatable system sizes, maintaining its total fidelity advantage
over a Trotter-Suzuki approximation making it an interesting candidate for
beyond classical time evolution.

Keywords: Variational Algorithms, Hamiltonian Simulation, Time Evolution,
Trotter-Suzuki Decomposition, Noisy Intermediate Scale Quantum, Translational
Invariance

##  1 Introduction

Quantum simulation of real time evolution represents a natural application of
quantum computation. In particular, non-equilibrium dynamics is known to lead
to highly entangled states [1, 2] and is therefore a strong candidate for
quantum advantage, i.e. an application where quantum computers can compute
quantities of interest that are not accessible to classical architectures [3].
Given a time-independent Hamiltonian $H$ of a closed, physical system of $N$
qubits, one can define the time evolution operator as the unitary

$\displaystyle U_{\text{ex}}=\exp{\left(-iH\tau\right)},$ (1)

where here and throughout the paper we set $\hbar=1$. This operator propagates
a quantum state $\ket{\psi(\tau_{0})}$ along a time interval $\tau$ to
$\ket{\psi(\tau_{0}+\tau)}$. The dimension of the underlying Hilbert space
grows exponentially in $N$, which makes simulation of large quantum systems
intractable for classical computation, in general. A quantum computer,
however, can directly operate on the $N$ qubits with quantum gates. When $H$
contains several non-commuting terms, we are unlikely to have a quantum gate
of the form of (1) at hand. Thus, a decomposition into experimentally
accessible quantum gates is needed. One such decomposition can be achieved by
a Trotter-Suzuki, or product formula, approximation [4, 5, 6, 7].

Besides product formula methods, there have been other approaches to implement
algorithms for Hamiltonian simulation, such as by Quantum Signal Processing
[8], Taylor series truncation [9, 10] and linear combination of unitaries
[11]. Recently, the gate complexity in Trotter-Suzuki algorithms has been
shown to scale comparably or even better than these so called post-Trotter
methods [12].

It has been shown [13] that for a permissible Trotter error $\epsilon$ of a
$q^{th}$ order product formula, the gate count $G$ is bounded by

$\displaystyle G\leq N5^{2q}\frac{(NT)^{1+1/2q}}{\epsilon^{2q}}.$ (2)

The Trotter error of $q^{th}$ order with Trotter number $m$ scales as
$\epsilon=\mathcal{O}(\frac{T^{2q+1}}{m^{2q}})$. For a nearest-neighbor
Hamiltonian, for instance, this bound can be further reduced to scale almost
linear in system size and simulation time in the worst case [12]

$\displaystyle
G=\mathcal{O}\left((NT)^{1+\mathcal{O}\left(\frac{1}{q}\right)}\right).$ (3)

For the simulation of Hamiltonians with local interactions [14] the gate
complexity of higher order Trotter-Suzuki formulas meets a lower bound (for
$q\to\infty$ in (2)). This lower bound is realised with product formulas [12]
as well as decompositions using Lieb-Robinson bounds [14].

Since quantum simulation on today’s Noisy Intermediate Scale Quantum (NISQ)
devices [15] suffers from gate imperfections, variational methods have been
put forward to reduce the total gate count and circuit depth [16, 17, 18, 19].
While a Trotter-Suzuki decomposition already yields a good approximation to a
Hamiltonian-generated unitary time evolution, it has been shown that
alternative methods can improve the total fidelity in specific cases by
trading off the improvement of approximation accuracy against the gate count.
These alternative methods include randomization [20, 21, 22, 23, 24, 25],
tensor network techniques [26], variational methods [27, 28, 29, 30, 31, 32,
33] and also a non-variational optimization of product formulas if there are
only few non-commuting terms in the Hamiltonian [34]. While the optimization
potential on the asymptotic behavior of the gate complexity for the simulation
of local interaction Hamiltonians is moderate, constant factor improvements
(in $N$) can already result in quantum advantage under NISQ conditions. This
is due to the fact that the fidelity for a desired output state decreases
exponentially in the gate count for finite fidelity gates.

Apart from the dependence on $N$ and $T$ as in (2), there is a constant factor
$5^{2q}$ which grows with the order of the Trotter-Suzuki formula $q$ [13].
This factor makes higher order product formulas unfavorable for NISQ devices.
On the other hand, finding an optimal gate sequence with variational methods
can improve on that factor by keeping $q$ low and thus represent a good NISQ
solution where simulation times and system sizes are moderate such that
asymptotic behavior is not dominating the gate number. As product formulas
represent a universal, i.e. model-independent, method for Hamiltonian
simulation, it is not surprising that one can optimize it for a specific
Hamiltonian.

Variational methods make use of this optimization potential in different ways.
They might exploit hardware-efficiency for a specific ansatz state [29, 28,
30], propose an efficient diagonalization of the Hamiltonian within the
available gate-set [31, 32] or focus on low-energy initial states [33].

In this work, we tackle two difficulties that arise in these variational
methods. First, employing a quantum optimization routine requires multiple
iterations with high sampling cost for each. This is a manifest problem in
NISQ algorithms [35, 36, 37, 38]. Our method allows to move the optimization
to a classical pre-processor for translationally invariant Hamiltonians.
Second, most of recent variational methods depend on a specific initial state
for which the dynamics are calculated [28, 29, 30, 31]. In contrast, a time-
universal simulation algorithm needs to be applicable for any initial state.
In particular, the same circuit can be applied multiple times on a desired
initial state and all the intermediate states. The method we present is state
agnostic and therefore not limited to applications on states with bound
entanglement as various classical methods [39].

##  2 Main Results

We present an initial state agnostic variational algorithm for the simulation
of translational invariant Hamiltonians where the heavy sampling effort of the
optimization is avoided via classical pre-processing.

For a translationally invariant Hamiltonian, our generalization scheme allows
to obtain the quantum circuit of a large system from the optimal solution of a
small, classically solvable system. The upscaling to larger systems can be
performed qubit-by-qubit if the variational parameters are chosen to be
translationally invariant.

For open boundary conditions, which break translational invariance, we propose
a variant of the method which uses two classical optimization routines, one
which optimizes the open boundary system and one which yields optimal
parameters for a periodic boundary system. The upscaled system then combines
the first set of parameters for the qubits near the boundary and the second
set for qubits within the bulk. Two dimensional upscaled open boundary systems
potentially require classical optimization of more than one small open
boundary system, but otherwise follow the same procedure. This scenario is
favored by modern NISQ devices which would require a number of SWAP operations
to implement periodic boundaries on two dimensional architectures. Also, open
boundary systems are expected to pose larger challenges for classical
simulation methods since no translation symmetry can be exploited in classical
algorithms. The process of gluing together optimized systems is discussed in
detail in section 4 and illustrated in Figure 1.

The presented methods are numerically exemplified on a one- and two-
dimensional Transverse Field Ising Model. On small systems for which we can
compare the quantum circuit with exact time evolution, our variational
approach improves upon the Trotter-Suzuki decomposition in terms of accuracy
by up to 3 orders of magnitude in the cost function and at the same time
reduces the gate count by a factor of 2. We are then able to suggest quantum
simulation circuits for larger systems using the optimal parameters from
classical pre-processing, which outperform the Trotter-Suzuki ansatz in a
similar way.

Under NISQ conditions, i.e. with imperfect gate fidelities, such improvements
turn out to be decisive as to whether an algorithm can be executed with a
Trotter-Suzuki ansatz. These scenarios appear when the simulation time $\tau$
for fixed circuit depth is large, see Figures 2 and 4.

In numerical experiments with the one-dimensional Transverse Field Ising
Model, we found for example that with 50 qubits, a single gate error of
$0.1\%$ and $\tau=2.5$ our algorithm obtains a NISQ preparation infidelity of
ca. $29\%$, being only a $7\%$ increase from the minimal infidelity due to
gate errors. Compared to approx. $57\%$ infidelity of the best Trotter-Suzuki
ansatz. For this numerical experiment, the optimal parameters for the time
evolution on a 6 qubit chain have been used.

##  3 Variational Optimization

We consider a $p$-local quantum many-body Hamiltonian $H$ which consists of
$N\cdot A$ terms

$\displaystyle H=\sum_{n=1}^{N}\sum_{a=1}^{A}c_{a}H_{n,a},$ (4)

where $a\in\\{1,...,A\\}$ counts different interaction types and
$n\in\\{1,...,N\\}$ denotes the first qubit on which $H_{n,a}$ is acting. For
a $k-$local term ($k\leq p$), for instance $H_{n,a}$ acts on qubits
$n,...,n+k-1$ modulo periodic boundary conditions.

Figure 1: a. Upscaling procedure for a one dimensional system with periodic
boundaries. The classical variational optimization is conducted on the gray
small system (here 4 qubits). To extend this small system to a larger system
qubit-by-qubit, for each interaction layer, the variational parameter of the
small system’s two-qubit gates is reused. To include the new qubit (blue), the
small periodic system is cut at one point and two-qubit gates connecting the
open points of the cut system and a new qubit are added (compare also circuit
picture in the middle). This extension strategy is motivated by the
translational invariance of the system. b. Upscaling procedure for a one
dimensional system with open boundaries. For open boundaries the translational
invariance that was used previously as motivation of the upscaling strategy is
no longer fulfilled. However, in this case two small systems, one with open
and one with periodic boundaries, can be used. The open boundary system (gray)
is cut in half and in its middle an arbitrary number of qubits using the
optimal variational parameters of the periodic boundary system (blue) are
inserted applying the same two-qubit gate gluing procedure as before. This
algorithm is illustrated in the middle in the quantum circuit picture.

The $\\{H_{n,a}\\}_{(n,a)}$ of the Hamiltonian are chosen such that each
$H_{n,a}$ is the generator of a quantum gate that can straight-forwardly be
implemented on the considered hardware (depending on the $p$-locality of $H$
and on the quantum architecture, further decompositions into experimentally
accessible gates might be involved).

The unitary time evolution generated by $H$ can be approximated by a first-
order Trotter decomposition [4]

$\displaystyle U_{\text{ex}}$
$\displaystyle=\exp{\left(-i\tau\sum_{n,a}c_{a}H_{n,a}\right)}=\left[\prod_{n,a}^{\leftarrow}\exp{\left(\frac{-i\tau}{m}c_{a}H_{n,a}\right)}\right]^{m}+\mathcal{O}\left(\frac{\tau^{2}}{m}\right),$
(5)

where $m$ denotes the Trotter number and the arrow denotes an ordering which
starts with the lowest pair $(n,a)$ on the right. The ordering inferred on the
tuple $(n,a)$ is not unique and only needs to be fixed when constructing the
corresponding circuit. The leading order error of this approximation is given
by commutators between the Hamiltonian terms [12]. A $q^{\text{th}}$ order
Trotter-Suzuki decomposition while using more gates improves the error to
leading order $\mathcal{O}(\frac{\tau^{q+1}}{m^{q}})$ [4, 5]. We will stick to
first order product formulas as a benchmark in the main text of this paper,
since higher order approaches use exponentially more gates in $q$ and hence
yield worse total fidelities compared to increasing the Trotter number $m$. We
show in A that also smaller Trotter numbers for higher order Trotter sequences
do not save enough gates to yield a better simulation.

Our method is independent of the particular variational quantum circuit ansatz
$U_{\text{var}}$, which is a gate sequence of parametrized unitaries

$\displaystyle U_{\text{var}}(\vec{\theta})=U_{1}(\theta_{1})\cdot
U_{2}(\theta_{2})\cdots$ (6)

We choose the first order Trotter-Suzuki approximation from (5) as a starting
point for the circuit ansatz, as this ensures an accuracy that is at least as
good as Trotter-Suzuki approximation. In order to exhaust the potential of the
ansatz from (5), we will substitute the coefficients $\frac{\tau}{m}c_{a}$ for
each of the Hamiltonians $H_{n,a}$ at each of the $m$ sequences by a
variational parameter $\theta_{r,a}$ with $r\in\\{1,...,m\\}$ labeling the
layers and $a\in\\{1,...,A\\}$ the different gates. The resulting ansatz reads

$\displaystyle
U_{\text{var}}(\vec{\theta})=\prod_{r=1}^{m}\left[\prod_{n,a}^{\leftarrow}\exp{\left(-\frac{i}{2}\theta_{r,a}H_{n,a}\right)}\right].$
(7)

Similar to the Trotter-Suzuki approximation, the variational ansatz admits an
approximation order $m$ which controls the number of parameters. There are
thus $A\cdot m$ different parameters to tune. Similar ansätze as in Eq. (7)
are known as variational Hamiltonian ansatz [37] for adiabatic ground state
preparation, where they are used with a different cost function than in our
case.

The optimization objective is to determine a set of optimal parameters
$\theta_{r,a}$ for which the exact time evolution $U_{\text{ex}}$ is
approximated best by $U_{\text{var}}$. The optimal parameters are found by
minimizing a cost function measuring the distance of unitaries. We choose a
measure of distance induced by the Frobenius norm

$\displaystyle C(U_{\text{var}},\vec{\theta})$

where $\tensor{U}{{}^{a}_{b}}$ denotes the elements of a matrix representation
of the operator $U$. The additional factor $\frac{1}{2^{N}}$ eliminates the
scaling of $C(\vec{\theta})$ in system size $N$. The exact operator of time
evolution $U_{\text{ex}}$ serves as a benchmark which is agnostic of an
initial state. In a similar fashion, we can consider
$C(U_{\text{trot}},\theta)$ as a measure of distance between the Trotter-
Suzuki ansatz and exact time evolution.

The cost function (LABEL:eq:Frob_cost) is a measure of how well time evolution
is approximated on average in the full Hilbert space. As we are considering
finite-dimensional Hilbert spaces, any error of observable dynamics can be
upper bounded by a distance measure

$\displaystyle|\braket{\psi}{U_{\text{ex}}^{\dagger}\,\mathcal{O}\,U_{\text{ex}}}{\psi}-\braket{\psi}{U_{\text{var}}^{\dagger}(\theta)\,\mathcal{O}\,U_{\text{var}}(\theta)}{\psi}|\leq
2||\mathcal{O}||_{\infty}\,||U_{\text{var}}(\theta)-U_{\text{ex}}||_{\infty}=2\kappa||\mathcal{O}||_{\infty}\,||U_{\text{var}}(\theta)-U_{\text{ex}}||_{F},$
(8)

where $||.||_{\infty}$ denotes the operator norm. The constant $\kappa$, in
general, depends on the distribution of the error throughout the Hilbert
space. For a concentrated error, we have $\kappa=1$ while for the pathological
case that only one state is simulated badly while every other state has
perfect fidelity, the upper bound gets loose by a factor of $\kappa=2^{n}$.
For the Trotter algorithm it is known that most of the time, the errors will
be concentrated [25] making the cost function (LABEL:eq:Frob_cost) a good
estimate for the performance of the algorithm.

As an example, the optimization of one- and two-dimensional Transverse Field
Ising Models on rectangular lattices is discussed in section 5, but for now
let us assume, we have found a set of parameters such that $C(\vec{\theta})$
as defined in (LABEL:eq:Frob_cost) is minimal. Note that (LABEL:eq:Frob_cost)
involves the full time evolution operator which is not accessible on a quantum
computer, but has to be calculated on a classical computer. Consequently,
there is no quantum advantage, if $U_{\text{var}}$ is deployed to a quantum
computer. However, a very meaningful, NISQ-friendly application arises, when
the found optimal parameters of the trained smaller system are applied to a
larger unseen system, which is explained in the following.

##  4 Upscaling System Size

Hamiltonians of many-body quantum systems often enjoy translational
symmetries, such that larger systems look similar to smaller systems, locally.
In particular, $H_{n,a}$ of (4) is the same operator for all $n$ except that
it acts on different qubits. Our method exploits this symmetry and generalizes
solutions of variational optimizations to larger quantum systems.

To begin with, the upscaling algorithm is explained for fully translational
invariant systems. These necessarily have periodic boundaries. Since some
modern quantum hardware, however, naturally comes with open boundary
conditions for the implementable gates which disrupt translational invariance,
we also extend the method by systems, where the translational symmetry is only
broken by open boundaries.

###  4.1 Translationally Invariant Systems

Translationally symmetric systems can be scaled up by successive application
of the gates from the small system to the new qubits. Figure 1 a, for
instance, shows a one-dimensional spin chain with periodic boundary
conditions. For nearest-neighbor interactions, the interaction between the
first two qubits is the same as the interaction between any other neighboring
pair of qubits. From this translational invariance, ansätze for larger quantum
systems are generated by copying the quantum gates of the small circuit. More
precisely, we apply the quantum gates, which acted on qubits $N$ and 1 before,
after extension to qubits $N$ and $N+1$ as well as 1 and $N+1$, respectively.
In principle, this qubit-by-qubit upscaling scheme can be applied repetitively
to extend an initially few qubit time evolution circuit to a many qubit one.
The final extended variational circuit can then be easily inferred from (7) to
be

$\displaystyle
U_{\text{var}}(\vec{\theta})=\prod_{r=1}^{m}\left[\prod_{\begin{subarray}{c}n=1,...,N+K\\\
a=1,...,A\end{subarray}}^{\leftarrow}\exp{\left(-\frac{i}{2}\theta_{r,a}H_{n,a}\right)}\right].$
(9)

where $K$ denotes the total number of added qubits. Note that (9) encodes more
than just tensorating copies of the systems, but restructures gates such that
entanglement between the former $N$ qubits and the new $K$ qubits is created
(cf. Figure 1 a).

Very importantly, the method itself is insensitive of the precise model as
long as it is translational invariant. As before, we will stick to the
Trotter-Suzuki inspired variational ansatz and only need to define the
Hamiltonian of the extended system $H^{(N+1)}$ and then use the gate sequence
from the copied addends in the Hamiltonian as before. To get this, we subtract
the old boundary term and insert the new terms for the new qubit

$\displaystyle H^{(N+1)}=H^{(N)}$
$\displaystyle-\sum_{a}\sum_{k=1}^{p-1}H_{N-(k-1)\rightarrow
p-k,a}+\sum_{a}\sum_{k=1}^{p}H_{N-(k-2)\rightarrow p-k,a}$ (10)

where we used the notation $H_{n_{1}\rightarrow n_{2},a}$ to emphasize that
this term acts on the qubits $n_{1},n_{1}+1,...,n_{2}$. This procedure is also
sketched for $p=2$ in Figure 1 a). Ultimately, this upscaling can be repeated
to end up with the Hamiltonian for $N+K$ qubits.

This process might seem redundant at first. However, if we recall the
variational method to improve on trotterized time evolution, we have just
derived a generalization scheme which allows us to use the parameters
$\vec{\theta}$ found for the $N$ qubit system to suggest a decomposition of
quantum gates for the simulation of the $N+K$ qubit system. This is done by,
again, breaking the periodic boundary conditions and gluing the gates
corresponding to the $K$ new qubits but now with the variational parameters
instead of the Trotter parameters $\frac{\tau}{m}c_{a}$.

It becomes apparent that the gluing method is identical to Trotterization on a
coarser level. This can be understood best, when looking at the reverse
process using the Trotter-Suzuki parameters $\frac{\tau}{m}c_{a}$. Given the
Hamiltonian $H^{(N+K)}$ of the large system, $\exp(-i\tau H^{(N+K)})$ is split
into the product $N+K$ identical gates with the only difference of acting on
different qubits as described in (10). Two subsequent Trotterizations, first
into $H^{(N)}$ and $H^{(K)}$ and afterwards into 2-qubit gates yield the same
result as a Trotterization directly into 2-qubit gates. As a result, if the
variational parameters are reset to the Trotter-Suzuki parameters, there is no
difference in a Trotter-Suzuki decomposition of $\exp(-i\tau H^{(N+K)})$ and
gluing $K$ qubits to the system $H^{(N)}$.

Since the leading error of Trotter approximations is given by commutators
between Hamiltonian terms, this error scales linearly in the qubit number for
local Hamiltonians [12]. Due to the similar structure of our ansatz in
equation (9), one may expect comparable scaling for our variational method.
This suggests that the upscaling procedure we employ only causes a degrading
of approximation accuracy that scales linearly in the system size. This
expectation is confirmed by our numerical experiments presented below.

###  4.2 Open Boundary Conditions

On some quantum hardware, such as superconducting circuit devices, the
implementation of periodic boundary conditions on two-dimensional lattice
models requires additional SWAP gates. Hence, the simulation of physical
systems with open boundary conditions is better suited for such devices,
particularly under NISQ conditions where each additional operation
exponentially decreases the success probability. To adapt the upscaling scheme
to systems with open boundary conditions, additionally to a small system with
periodic boundaries, a small system of size $N_{o}$ with open boundaries is
variationally optimized. The corresponding Hamiltonian has the same form as
(4) but without terms reaching beyond the periodic boundary. An intuition for
why this should work is that if the depth of the circuit is short enough, the
bulk qubits of the open system are affected by approximately the same
operations as for a periodic boundary system. This is why we are allowed to
glue an arbitrary periodic system into the middle of an open boundary system.

Time evolution can be approximated by the operators $U_{\text{trot}}{}_{,o}$
and $U_{\text{var}}{}_{,o}(\vec{\theta})$ in an analogue way as before. The
upscaled system is constructed from both the small periodic and non-periodic
one, by cutting the small open boundary system of size $N_{o}$ in half and
inserting an arbitrary number $K$ of periodic boundary qubits in the middle
(see Figure 1 b.). The gluing of the gates for each time step works as
described before. The final circuit for time evolution of a single time step
($m=1$) for the glued system with $N_{o}+K$ qubits with open boundaries reads

$\displaystyle U_{\text{var},o}^{(1,N_{o}+K)}\approx\,$ $\displaystyle
U_{\text{var},o}^{(1,\frac{N_{o}}{2})}(\vec{\theta}_{o})\,U_{\text{var}}(\vec{\theta})\,U_{\text{var},o}^{(K+\frac{N_{o}}{2}+1,K+N_{o})}(\vec{\theta}_{o})$
(11)

where we used $U_{\text{var}}(\theta)$ from (9) for the time evolution of the
bulk. $U_{\text{var},o}^{(1,X)}$ denotes circuit simulating the open boundary
dynamics on qubit 1 to $X$. Furthermore, we split the optimized parameter set
$\theta_{o}$ of the open boundary $N_{o}$ qubit simulation in half in order to
use them for the time evolution of the left (or respectively right) boundary.
For the sake of simplicity, we assume $N_{o}$ to be even. This way, the three
unitaries in (11) correspond to the three segments of Figure 1 b that are
being glued. The boundary parameters $\theta_{o}$ and the bulk parameters
$\theta$ result from two separate optimizations. We are thus making an
approximation neglecting open boundary effects on the bulk qubits. We will see
in numerical experiments that for a one- and two-dimensional Transverse Field
Ising Model it suffices to optimize only few dedicated boundary gates separate
from the periodic bulk.

##  5 Transverse Field Ising Model

To numerically test the performance of the variational gluing method in
comparison to Trotter-Suzuki approaches, the Transverse Field Ising Model
(TFIM) in one and two dimensional is studied. For this, we simulate the exact
time evolution $U_{\text{ex}}$ generated by the Hamiltonian

$\displaystyle H_{TFIM}=J_{z}\sum_{\langle
i,j\rangle}Z_{i}Z_{j}+h_{x}\sum_{i}X_{i},$ (12)

with the Pauli operators $X,Z$ and the coupling constants $J_{z}$ and $h_{x}$.
We denoted a Pauli $Z$ operator acting on the $i^{\text{th}}$ qubit by $Z_{i}$
(and analogously for $X$). Furthermore, the notation $\langle i,j\rangle$ only
counts nearest neighbor terms and we implicitly understood periodic boundary
conditions. The Transverse Field Ising Hamiltonian relates to our previous
notation via $H_{n,a=1}=Z_{n}Z_{n+1},H_{n,a=2}=X_{n}$ and
$c_{1}=J_{z},c_{2}=h_{x}$. To demonstrate our method on another model, we show
additional results for an $XY$ model with a transverse field in C.

The cost values for the Trotter-Suzuki $U_{\text{trot}}$ and variational
sequence $U_{\text{var}}$ are compared for different system sizes, as well as
for multiple repetitions of the circuit (i.e. multiple time steps). To
emphasize that near term NISQ hardware can benefit from our classical pre-
training algorithm for time evolution, simulation results of 2-point spin
correlations including gate infidelities are presented.

We simulated and optimized systems with a total number of up to 12 qubits in
the discussed fashion. For the study of larger systems limitations of memory
and time resources forced us to approximate the cost function
(LABEL:eq:Frob_cost) by sampling over a set of random input vectors
$\mathcal{V}$

$\displaystyle C_{var}(\vec{\theta})$
$\displaystyle\approx\frac{1}{|\mathcal{V}|}\sum_{\psi\in\mathcal{V}}\left(1-|\braket{\psi}{U_{\text{ex}}^{\dagger}U_{\text{var}}(\vec{\theta})}{\psi}|\right),$
(13)

where the absolute value $|.|$ ensures the equality of the two operators up to
a multiplicative global phase. This way, we are able to present results for
system sizes up to 24 qubits. For the circuit simulations we used the python
library Cirq [40] and for optimization we used the Adam gradient descent
method [41]. In all the optimizations, we have used a maximal number of
traning steps of $10^{4}$, a learning rate of $10^{-4}$ and the
hyperparameters $\beta_{1}=0.9,\beta_{2}=0.99$ to control the decay rates of
the gradient and squared gradient average.

There is an ambiguity in the choice of the ordering of the gate sequences
which correspond to $Z_{i}Z_{j}$ and $X_{i}$. This leads to different
operators when doing Trotterization. Here, we fixed the ordering for the sake
of an efficient implementation to

$\displaystyle U_{\text{var}}(\vec{\theta})$
$\displaystyle=\prod_{r=1}^{m}\Bigg{[}\prod_{\begin{subarray}{c}n=1\\\ n\text{
odd}\end{subarray}}^{N}e^{-i\theta_{r,1}Z_{n}Z_{n+1}}\prod_{\begin{subarray}{c}n=1\\\
n\text{
even}\end{subarray}}^{N}e^{-i\theta_{r,1}Z_{n}Z_{n+1}}\prod_{n=1}^{N}e^{-i\theta_{r,2}X_{n}}\Bigg{]}.$
(14)

Each of these gates can be implemented directly on current superconducting
hardware via tunable couplers [42, 43, 44, 45, 46]. In this numeric example,
four cases are studied. These are one- and two-dimensional systems with either
periodic or open boundary conditions. In the case of two-dimensional periodic
lattices, the gluing can be performed in either of the two dimensions by
adding a row or column of qubits analogously to the discussion in section 4.
Demonstration of scaling up in system size on two-dimensional systems with two
open boundaries is computationally prohibitive, as the smallest fully open,
two-dimensional and non-trivially glued system would be of size $5\times 7$
(cf. B). We therefore employ open boundaries in only one direction, while
retaining periodic boundary conditions for the orthogonal direction. Optimal
parameters can be classically found but no measure of accuracy can be
evaluated any longer as we exceed our limit of numerical feasibility with a
total size of 35 qubits.

###  5.1 Single Time Step

To begin with, our variational circuit is compared to Trotter-Suzuki
approaches with equal or higher gate count and depth, where higher gate count
is usually traded against a better approximation. Afterwards, the optimal
parameters are reused in larger systems and also compared to the Trotter-
Suzuki benchmarks. Figures 2 a and d. show that the variational ansatz is able
to beat Trotterizations with a higher Trotter number along all numerically
feasible system sizes (while being less interesting as it requires more gates,
higher order Trotter-Suzuki approximations can also be beaten, see A).
Notably, the cost function values for upscaled systems hardly differ from the
systems on which the optimization is performed. For the one-dimensional
system, optimization was performed on a 6-qubit chain. An upscaling to 24
qubits only raises the cost value by a factor of 3. Analogously, for a two-
dimensional system, the optimal parameters found for a $3\times 3$ grid are
extended to a $4\times 6$ grid by increasing the cost value by around the same
factor. The cost values for the two-dimensional models are generally higher
since each qubit has more interactions partners.

A linear fit in system size gives an estimate of how the improvement develops
for system sizes which can no longer be classically simulated. In all studied
cases, we observe an improvement factor in cost of the variational method
compared to Trotter-Suzuki of more than 1000 in the one-dimensional and more
than 100 in the two-dimensional case while using only half of the number of
gates. The linear fits indicate that these improvement factors are constant
over several larger system sizes.

NISQ gates have finite infidelities, leading to large compound errors when
gate counts are high. It is therefore important to be as gate efficient as
possible when running Hamiltonian simulation on NISQ hardware. To take into
account gate infidelities, we quantify the total infidelity of a simulated
NISQ experiment by

$\displaystyle
1-\mathcal{F}_{NISQ}=1-\braket{\mathcal{F}_{approx}}\mathcal{F}_{gates},$ (15)

where $<\mathcal{F}_{approx}>$ denotes the averaged fidelity of the output
state (variational or Trotter) when the hardware runs without error and
$\mathcal{F}_{gates}$ is the average fidelity of the implementations of the
employed gates. $<\mathcal{F}_{approx}>$ is inferred from the cost function
used for optimization, where (LABEL:eq:Frob_cost) takes an average over the
whole Hilbert space while (13) only takes an average over a randomly chosen
subset. For the gate fidelity $\mathcal{F}_{gates}$, we consider a simple
error model assuming the same fidelity $p_{g}$ for each individual gate. The
circuit fidelity is given as the product of the individual gate fidelities for
every gate in the circuit, i.e. $\mathcal{F}_{gates}=p_{g}^{G}$.

Figure 2:  System size upscaling of variational and Trotter sequence for a
single time step $\tau$. In the simulations of several system sizes and for
one- and two-dimensional systems with either periodic or (half-)open
boundaries, $J_{z}=1.0$ and $h_{x}=0.25$ were used. The data points show
numerical evaluations of the cost function for the periodic boundary case. The
solid and dashed lines in a. and d. represent a linear fit in the number of
qubits for both periodic and (half-)open boundaries. a, d. Comparison of cost
function for a time step of $\tau=0.3$. The optimization was performed on 6
qubits for the one-dimensional and on a $3\times 3$ grid for the two-
dimensional systems. b, e. Variation of the single time step $\tau$ versus
overall infidelity from algorithmic and gate error ($p_{g}=99.9\%$) for $N=15$
qubits ($3\times 5$ for $d=2$). The noise floors generated by the sole gate
infidelity (independent of $\tau$) are depicted by a shaded region in blue
($m=3$) and red ($m=6$). For small time steps, the gate infidelity dominates
the error, while for large time steps, the algorithmic error takes over. For a
one-dimensional system, this turning point lies around $\tau=2.5$ and
$\tau=1.5$ for two-dimensional systems. The turning points of the Trotter-
Suzuki decompositions are earlier. c, f. Scaling of the total infidelity in
system size (analogue to a. and d.) for $\tau=2.5$ ($\tau=1.0$ for $d=2$).
Here, again, the optimal parameters found from a 6-qubit chain ($3\times 3$
grid for $d=2$) are scaled up to a total size of 24 qubits and fitted as in a.
and b..

The total number of gates $G$, in general, is determined by the Trotter number
$m$. In particular, for the Trotter-Suzuki decomposition of a Transverse Field
Ising Model with periodic boundaries, depth $D$ and gate count $G$ of a single
time step read

$\displaystyle D=(d+1)\,m,\qquad\qquad G=(d+1)\,mN-bdm=(N-b)D+bm,$ (16)

where $d$ denotes the spatial dimension and $b$ is 0 for periodic or 1 for
open boundary conditions. We assume an optimistic gate fidelity of
$p_{g}=99.9\%$ [47], so that our results should remain relevant for future
years.

The improvement in $\mathcal{F}_{var}$ and hence in $\mathcal{F}_{NISQ}$
strongly depends on which single time step size $\tau$ is chosen. For the
Trotter-Suzuki approximation, a small ratio of $\tau$ and $m$ is favored to
keep the error small (cf. (4)). This can be done by increasing $m$ and paying
the price of using more gates. Figures 2 b, e show that for small $\tau$ the
gate infidelity dominates the error, whereas for large $\tau$, the algorithmic
error takes over. For the one-dimensional case, Trotter-Suzuki error becomes
important at around $\tau=1.5$ for $m=3$ and at $\tau=2.0$ for $m=6$. The
variational algorithm has the latest turning point, where the algorithmic
error becomes dominant, around $\tau=2.5$. This turning point is earlier
(around $\tau\approx 1.5$) in the two-dimensional model.

If we choose the single time step at this turning point, we find that the
variational sequence increases the lower infidelity bound due to gate
infidelities only marginally, while the Trotter-Suzuki decompositions
introduce a larger additional error.

Figure 3:  Improvement of variational vs. Trotter algorithm for a sector of
the $J_{z}-h_{x}$ parameter space and multiple repetitions. The parameters are
scanned through the values $J_{z}\in[-1,1],h_{x}\in[0,1]$ and $\tau=0.3$ is
fixed. Solid lines denote a moving mean through the simulated data (data
points). Using parameters optimized on a 6 qubit chain both, the 6 qubit (a. -
c.) and 15 qubit (d. - f.) systems, show similar scaling behaviour in the
simulation time. A single time step shows an improvement of cost between
$10^{3}$ and $10^{5}$ throughout the scanned values for $J_{z}$ and $h_{x}$.
Even after 100 repetions for most $J_{z}$ and $h_{x}$, the variational
algorithm beats the Trotter-Suzuki algorithm. Only where the Trotter-Suzuki
decomposition is nearly exact (when $h_{x}$ or $J_{z}$ are almost 0), Trotter-
Suzuki overtakes before 100 repetitions. In the worst case scenario, the
variational algorithm gets overtaken at around 50 time steps. Figure 4 shows
that after 50 repetitions, however, the circuit infidelity is typically
dominated by gate imperfections. The scaling in simulation time is
demonstrated for the cases of open boundaries on a one-dimensional system with
15 qubits (g.), as well as a 2D a system with periodic (h.) and half-open
boundaries (i.) of size $3\times 5$. Here, the interactions $J_{z}=1.0$ and
$h_{x}=0.25$ are chosen as an instance of the average scenario.

For example, with 50 Qubits and periodic boundaries in one dimension the NISQ
infidelity $1-\mathcal{F}_{NISQ}\approx 29\%$ is increased by only $7\%$ due
to the variational approximation infidelity, whereas
$1-\mathcal{F}_{NISQ}\approx 57\%$ for the best ($m=6$) Trotter-Suzuki ansatz.
Both, a higher gate count and approximation infidelity lead to this
deterioration. The standard deviation $\sigma$ of the linear
$\mathcal{F}_{approx}$ extrapolations in Fig 2 c. and f. are small compared to
$\mathcal{F}_{approx}$: $\sigma=0.19\%\cdot\mathcal{F}_{approx}$ (1D) and
$\sigma=0.03\%\cdot\mathcal{F}_{approx}$ (2D). Therefore, the fits in Figure 2
yield a faithful extrapolation in system size.

###  5.2 Multiple Time Steps

Since we determine $U_{\text{var}}$ via a state agnostic optimization, it can
be applied repeatedly to the time evolved state. In this way larger simulation
times $T$ can be covered. Figure 3 show this for a 6- (a. - c.) and a 15-qubit
system (d. - f.) both using the parameters optimized on 6 qubits. The
variational sequence is optimized for a single time step, and improvements in
the cost function decrease for later times. Despite this, the variational
algorithm remains better than the Trotter-Suzuki benchmark even after one
hundred repetitions for most values of the Hamiltonian coefficients $J_{z}$
and $h_{x}$. However, when one of these coefficients is small the Trotter
error decreases, which can lead to Trotter-Suzuki overtaking the variational
algorithm in certain regions of parameter space. Indeed, the Trotter-Suzuki
approximation becomes exact in the limit that either $J_{z}$ or $h_{x}$ goes
to zero. In extreme cases, the accuracy of the Trotter-Suzuki sequence exceeds
that of the variational algorithm after 50 repetitions. However, as we shall
show later, at 50 repetitions the total error is dominated by gate
infidelities rather than approximation error.

We obtain a similar behavior of the cost function for the case of $d=2$ and
open boundary conditions. This indicates that the open boundary variation of
our method can be applied without significant change of performance. In this
numerical demonstration, our variational algorithm exhibits even a better
approximation fidelity for half-open boundaries than for periodic boundaries,
and hence a fully translational invariant system. This further stresses the
feasibility of our method for NISQ-inspired open boundary quantum systems.

###  5.3 Precision of Time-Evolved Observables

When computing the quantum dynamics of a system, one is typically interested
in the evolution of physical observables. From (8), it is clear that the cost
values only yield an upper bound for the error in predicting observables.
Although the variational algorithm has a smaller cost value, Trotter-Suzuki
could thus, in principle, describe the dynamics of a specific observable for a
specific initial state with greater accuracy. The approximation error in the
expectation value of a time-evolved observable hence depends on the precise
observable and on the chosen initial state. In this section, we show results
for exemplary observable dynamics and an initial state.

Note that the gluing procedure is not sensitive to entanglement of the input
state, since the time evolution operator is agnostic of the initial state.
States of high entanglement (or strongly entangling dynamics) are particularly
interesting, as the Trotter error of localized physics has been found to be
independent of simulation time [48].

In Figure 4, we illustrate the relative error on two-point correlation of two
neighboring spins on a 15 qubit system ($3\times 5$ for $d=2$) far away from
the boundary, $\mathcal{O}=\braket{\psi(t)}{Z_{8}Z_{9}}{\psi(t)}$
($\mathcal{O}=\braket{\psi(t)}{Z_{(2,3)}Z_{(2,4)}}{\psi(t)}$ for $d=2$). We
find that for the initial state $\ket{\psi(0)}=\ket{+}^{\otimes N}$, the two-
point correlation $\mathcal{O}$ predicted by $U_{\text{var}}$ lies closer to
the exact dynamics than any of the Trotter-Suzuki benchmarks.

Taking a product state as initial state is particularly interesting, as it is
accessible in experiment, yet also yields an entangled state after a few
repetitions. Although the $m=3$ variational circuit and the $m=6$ Trotter-
Suzuki benchmark achieve similar approximation error, the latter contains
twice as many gates and thus suffers from a squared gate fidelity.

These imperfections are illustrated in Figures 4 c. - d., where we show the
overall infidelity of (15), for up to 20 repetitions of the single time step
$\tau=1.0$. Besides the variational sequence beating both Trotter-Suzuki
benchmarks, the overall infidelity is minimized even more when optimizing for
a larger time step $2\tau$ using $m=6$ layers. Despite a significant noise
floor set by the gate infidelities, we achieve an improvement of infidelity by
a factor between 3 ($d=2$) and 4 ($d=1$) for a single repetition. Since the
best Trotter-Suzuki decomposition starts with an infidelity of $44.5\%$, this
shows that the simulation time $2\tau=5.0$ cannot yield satisfactory results
in the sketched NISQ scenario, while the variational sequence does (infidelity
of $18.8\%$). Since also twice as many parameters are involved, it is expected
that the variational sequence optimized on longer times performs better. This
comes at the cost of a more difficult optimization task in the classical pre-
processing step. The trade-off between more time spent once in classical pre-
processing against a permanently better quantum circuit is very favourable on
NISQ devices. Generalizing this to longer simulation times $k\tau>5.0$ and
more parameters $m>6$, however, does not further improve the long time
performance and thus is omitted in Figure 4.

All the numerical experiments support the maxim to keep the number of
repetitions low while using as large time steps $\tau$ as possible. To do
this, an increase of the Trotter number for the variational single time step
$m$ might be worth evaluating (cf. Figure 4). For small $\tau$ (as $\tau=0.3$
in Figures 2 and 3, for instance), raising the Trotter number $m$ does not
significantly improve the cost values. Therefore, the data points for $m>3$
are omitted.

Figure 4:  a. - b. Error of the two-point correlation
$\mathcal{O}=<\psi|Z_{8}Z_{9}|\psi>$ between qubit 8 and 9 of a 15 qubit
system ($3\times 5$ for $d=2$) for multiple time steps starting with the
initial state $\psi=\ket{+}^{\otimes 15}$. The interaction coefficients and
the single time step read $J_{z},=1.0$, $h_{x}=0.25$ and $\tau=1.0$. Solid
lines denote a moving mean through the simulated data (data points). Although
the variational sequence quickly converges to the same error as Trotter with
$m=6$, it will be less affected by gate infidelities. c. - d. Infidelity
including imperfect gates for multiple time steps with single time step
$\tau=2.5$ ($\tau=1.0$ for $d=2$). The infidelities have been normalized to
the infidelity value of the variational sequence optimized on a single time
step with $m=3$ layers. To further reduce the infidelity for longer times, it
pays off to optimize a single time step of $2\tau$ including $m=6$ layers and
hence with double the number of parameters. Both variational algorithms beat
the Trotter-Suzuki benchmarks for up to 8 repetitions. Afterwards every
approach is dominated by noise.

Also remarkably, the (half-)open boundary systems yield even smaller
infidelity values than their periodic counterparts. Besides supporting the
statement that open boundaries can be glued to get NISQ-friendly systems, this
phenomenon can be explained, since in open boundary systems fewer gates
contribute to both approximation and gate error.

##  6 Conclusion

In this work, we demonstrated how classical pre-processing can be leveraged to
determine a variational representation of time evolution operator of a
translation symmetric quantum system that reliably exceeds comparable
Trotterization in approximation accuracy. The main idea for our variational
algorithm consists of classically optimizing a variational time evolution
ansatz circuit on a small tractable system and then reusing these optimal
parameters for larger, classically intractable systems by gluing gate
sequences together.

These results are particularly interesting due to their applicability to NISQ
devices. Most promisingly, it can avoid any sampling cost and limitation
associated with variational optimization on quantum hardware. Further, the
significant increase in approximation fidelity compared to Trotter-Suzuki with
the same gate count for practically relevant system sizes can be taken
advantage of to realize the same approximation fidelity at a notably reduced
gate count. We found that, for one- and two-dimensional example problems, the
observed advantages persist and are in fact enhanced for systems with open
boundaries, despite the loss of translational symmetry on the edges. Such open
boundary systems are important as they are easier to realize on NISQ hardware.

Our results open up a highly versatile method for variational quantum
simulation. Since we have shown that it is sometimes sufficient to maximize
the simulation quality locally, possible limitations of this gluing scheme are
worth investigating further. We also expect to gain further insights into the
scope and performance of our method from investigations of the relations
between the range of correlations in simulated quantum states and the size of
the system considered in classical optimization. Since the optimization is
however state agnostic, depending only on the Hamiltonian of the system, it is
expected to perform well in simulating the time evolution of highly-entangled
states with long-ranged correlations.

Another way of optimizing quantum simulation algorithms can be achieved by a
quantum optimization. The cost function of (13) can be measured by comparing
initial and output state of a circuit which first performs the variational
ansatz circuit which propagates forward in time. Afterwards the state is
propagated back to its initial state by a high order Trotter-Suzuki
approximation. This way, the variational circuit can be optimized to invert
the Trotter-Suzuki decomposition with fewer gates and later scaled up in
system size. Further investigation is needed to determine which system sizes
are suitable for optimization.

Randomization techniques have attracted much attention recently and have been
shown to compress Trotter circuits [21, 24]. It remains to be investigated
whether randomized ansätze can be used for variational algorithms like the one
presented in this work.

We thank Petr Zapletal for valuable feedback and discussions and Naeimeh
Mohseni for feedback on the manuscript. This work received funding from the
European Union’s Horizon 2020 research and innovation program under Grant
Agreement No. 828826 “Quromorphic.”, from the German Federal Ministry of
Education and Research via the funding program quantum technologies - from
basic research to the market under contract number 13N15577 “MANIQU”. It is
also part of the Munich Quantum Valley, which is supported by the Bavarian
state government with funds from the Hightech Agenda Bayern Plus. TE
acknowledges funding from the International Max Planck Research School Physics
of Light.

## Data Availability Statement

The data that support the findings of this study are available upon reasonable
request from the authors.

##  Appendix A Higher Order Trotter-Suzuki Decomposition

Figure 5:  Comparision of approximation and preperation fidelity of
variational algorithm and higher order $q$ Trotter-Suzuki simulation for up to
100 repetitions of time evolution under a Transverse Field Ising Hamiltonian
with $J_{z}=1.0,h_{x}=0.25,\tau=1.0$ of a $3\times 6$ qubit system. In a, b.
the Trotter number is $m=1$ and in c, d. $m=2$. The variational algorithm
always has $m=3$. Solid lines denote a moving mean through the simulated data
(data points). a, c. The variational algorithm beats the $4^{\text{th}}$ order
product formula for $m=1$ while using only $30\%$ of the gates. The
$6^{\text{th}}$ and $8^{\text{th}}$ order product formula always beat the
variational algorithm in approximation accuracy but need more gates by a
factor of 50 and 250. b, d. Taking into account a finite gate fidelity of
$p_{g}=99.9\%$ together with the algorithmic fidelity $\mathcal{F}_{approx}$
yields the total infidelity of the algorithm. Already after a single time
step, the gate fidelity becomes dominant, such that the variational algorithm
outruns all higher order product formulas.

A more accurate decomposition as in (5) can be achieved by higher order
product formulas. A second order Trotter-Suzuki decomposition is achieved by
only reordering the operators in a symmetric way. For a single Trotter step
($m=1$) this reads

$\displaystyle U_{\text{ex}}$
$\displaystyle=\left[\prod_{n,a}^{\leftarrow}\exp{\left(\frac{-i\tau}{2}c_{a}H_{n,a}\right)}\prod_{n,a}^{\rightarrow}\exp{\left(\frac{-i\tau}{2}c_{a}H_{n,a}\right)}\right]+\mathcal{O}\left(\tau^{3}\right)=:S_{2}(\tau),$
(17)

Furthermore, one can achieve arbitrarily accurate approximations of an
exponential $e^{-i\tau H}$ by applying the $q^{\text{th}}$ order Trotter-
Suzuki formula which is defined recursively by [4]

$\displaystyle
S_{q}(\tau)=[S_{q-2}(\nu_{q}\tau)]^{2}S_{q-2}((1-4\nu_{q})\tau)[S_{q-2}(\nu_{q}\tau)]^{2},\quad\text{
for }q>2$ (18)

with $\nu_{q}=\left(4-4^{\frac{1}{q-1}}\right)^{-1}$. In addition to (18), one
can also pick a Trotter number $m>1$ to raise the accuracy further. Employing
a higher $m$ is achieved by substituting $\tau\mapsto\frac{\tau}{m}$ in (18),
but repeating the sequence $m$ times, thus using more gates by a factor of
$m$.

We will stick to the $m\in\\{1,2\\}$ case here and show that the use of higher
order formulas as a benchmark for our method is practically not relevant for
NISQ-focused variational optimisation, since the number of required gates,
which scales with $5^{q}$ [13], will dominate the error for finite gate
infidelities. We provide results for $q\in\\{2,4,6,8\\}$.

The overall circuit fidelity is given by the product of the imperfect gate
fidelity $\mathcal{F}_{gate}$ and the fidelity due to the approximate
algorithm

$\displaystyle\mathcal{F}_{var}(\ket{\psi})=\braket{\psi}{U_{\text{ex}}^{\dagger}U_{\text{var}}}{\psi}$
(19)

for an arbitrary initial state $\ket{\psi}$. Analogously, one can define the
algorithmic fidelity $\mathcal{F}_{trot}$ of a Trotter-Suzuki algorithm
described by the product formula in (18). We combine the notations
$\mathcal{F}_{var}$ and $\mathcal{F}_{trot}$ into $\mathcal{F}_{approx}$.
Figure 5 shows the total infidelity
$1-\mathcal{F}_{gate}\cdot\mathcal{F}_{approx}$ for up to 100 time steps. One
can see that while second order Trotter-Suzuki would admit fewer gates than
the variational algorithm ($m=3$), it suffers from greater algorithmic errors.
Higher order Trotter-Suzuki algorithms are already dominated by gate
infidelities after a single time step. The variational algorithm hence outruns
all of the benchmark algorithms in the presented simulation problem. For this
reason, no higher order $q$ Trotter-Suzuki formulas are discussed in the main
text. Other product formulas [49, 50] behave very similar to Trotter in this
respect. As Trotter exhaust theoretical complexity lower bounds for higher
orders [12], we only discuss higher order Trotter formulas here.

##  Appendix B Smallest Two-dimensional Transverse Field Ising Model with
Fully Open Boundary Conditions

Figure 6:  Smallest possible, non-trivially glued configuration of a fully
open two-dimensional qubit grid.  The building blocks, on which the
optimization is performed, are shown on the left and the fully glued grid on
the right and bold lines depict open boundaries. The minimal size of a half-
open grid is $3\times 4$ (orange) and of an open grid $4\times 4$ (green). The
minimal extension of the bulk then consists of three qubits (blue) thus
leaving an open $5\times 7$ grid.

We construct the smallest non-trivially glued fully open boundary system in
two dimensions. It turns out that, although optimal parameters to use on fully
open two-dimensional systems can be found, no measure of accuracy can be
numerically evaluated. In Figure 6 every colored square corresponds to a qubit
and thick lines represent open boundaries. As every qubit contributes to the
Hamiltonian with a number of $A$ interactions, the squares can be thought of
as a set of $A\cdot m$ optimal parameters and the figure depicts how
parameters are copied when scaling up. The blue squares all share the same
parameters for each gate.

Optimized open boundary systems have to be cut in half to insert bulk qubits
inbetween (cf. section 4). In two dimensions however, this means the system
needs to be cut along both dimensions. Then first, a half-open system is
inserted to faithfully keep the open boundary upright (orange blocks in Figure
6). Finally, completely periodic qubits can be filled in the very bulk of the
system (blue blocks), but also the green blocks have added at least one column
of qubits to keep everything rectangular. The additional green column copies
parameters from the inner most qubit of one half of the boundary (for one
column it does not matter whether it copies the parameter from the left or the
right). This way, the inner most qubit of each half of the fully open system
has to be considered equivalent to a qubit from a half-open system. If we
would not choose at least 4 qubits along this open boundary, at least one
boundary qubit would coincide with a bulk qubit which is contradictory.

If there is a periodic boundary along one dimension nothing has to be cut in
half, however, the optimization should treat every qubit such that it has two
distinct neighbors along the periodic dimension. This can only be achieved
with at least three qubits, as otherwise left and right neighbors would
coincide. As a result, the smallest periodic boundary block is of size
$3\times 3$, the smallest half-open block is $3\times 4$ qubits (3 along the
periodic boundary and 4 along the open boundary) and the smallest fully open
block is of size $4\times 4$.

Thus, the smallest glued open boundary system in two dimensions is of size
$5\times 7$ consisting of three optimized blocks, a half-open one of size
$3\times(2+2)$, a fully open block of size $5\times(2+2)$ (the optimized block
is of size $(2+2)\times(2+2)$) and a row of 3 qubits in the middle, see Figure
6. The optimization of the single elements can all be performed on a classical
computer, but the cost function for a $5\times 7$ qubit system can no longer
be evaluated, as we exceed our limit of numerical feasibility. The results for
a half-open system of size $3\times(4+K)$ can be evaluated and hence represent
a numerical example of gluing open boundaries to two-dimensional spin lattice
models.

##  Appendix C The XY Model with a Transverse Field

Next to the study of the Transverse Field Ising Model, in this appendix the
$XY$ model with a transverse field (TFXY) is discussed using similar quality
measures as in section 5. Coming from the Transverse Field Ising Model, we
just add a $YY$ interaction to neighboring qubits and get an $XY$ model with
an external field. The Hamiltonian reads

$\displaystyle H_{TFXY}=J_{y}\sum_{\langle
i,j\rangle}Y_{i}Y_{j}+J_{z}\sum_{\langle
i,j\rangle}Z_{i}Z_{j}+h_{x}\sum_{i}X_{i}.$ (20)

Since this model includes more non-commuting terms, we expect the Trotter
error to be larger, making it another interesting model for NISQ simulation.
This claim is supported by Figure 7 which compares cost values of variational
and Trotter-Suzuki sequence scaled in system size and in simulation time, as
well as errors of observable dynamics.

Figure 7:  System size and simulation time scaling of variational and Trotter
sequence for the XY Model. In the simulations of several system sizes and for
one- and two-dimensional systems with periodic boundaries, the interaction
coefficients $J_{y}=0.5,J_{z}=1.0,h_{x}=0.25$ and a single time step of
$\tau=0.3$ ($\tau=0.6$ for Variational $m=6$) were used. As in before, solid
lines in a. and b. are linear fits and moving means through the simulated data
(data points) in c - f.. a. – b. Comparison of cost function for a single time
step $\tau$ and different system sizes. The optimization was performed on 6
qubits for the one-dimensional and on a $3\times 3$ grid for the two-
dimensional system. The data points represent numerical evaluations of the
cost function whereas the solid lines represent a linear fit in the number of
qubits. c. - d. Cost values for multiple time steps on a one-dimensional
system with 15 qubits and a two-dimensional system of size $3\times 5$. For
the XY model, the choice of $\tau=0.6$ and $m=6$ yields a better result for
scaling in time. e. - f. Relative error of the two-point correlation
$\mathcal{O}=\braket{Z_{8}Z_{9}}$ for the same systems as in c. - d.. The
initial state is chosen to be $\psi=\ket{+}^{\otimes 15}$. Although the cost
value is generally smaller in the one-dimensional case (cf. a. - d.), the
Trotter-Suzuki sequence approximates the chosen observable better in the two-
dimensional case, reducing the improvement by optimization. In c. - f. solid
lines denote a moving mean through the simulated data (data points).

Figures 7 a. - b. show that the first order Trotter-Suzuki approximation has a
larger cost value compared to the TFXY Model with the same interaction
coefficients as in section 5 except adding a $YY$ interaction with
$J_{y}=0.5$. As before, the improvement from optimization of a single time
step is larger in one-dimensional than two-dimensional systems and also the
scaling in system size, resembles the behavior of the Transverse Field Ising
Model well. The extrapolation of the cost value in systems size keeps up an
improvement in cost value of a factor of 1000 (100 for $d=2$) compared to
Trotter-Suzuki with Trotter number $m=3$. For a one-dimensional system with 50
qubits, for instance, the Trotter-Suzuki decomposition cost value is
extrapolated to around $0.01$ for $m=3$ and $5\cdot 10^{-3}$ for $m=6$. The
variational sequence has a cost value of $10^{-5}$. For a two-dimensional grid
of 50 qubits, in total, the cost values are of the order of $10^{-1}$ (Trotter
$m=3$), $10^{-2}$ (Trotter $m=6$) and $10^{-3}$ (Variational $m=3$).

The scaling in simulation time is shown in Figures 7 c. - f.. The general
trend of the plots again resembles the case of the Transverse Field Ising
Model which indicates a certain model independence of the variational method.
The slight increase in approximation error in the $XY$ model however also
reduces the number of circuit repetitions which can be performed until the
variational sequence falls back to accuracies of Trotter-Suzuki. For an
optimization on $m=3$ layers, this happens after around 40 repetitions in the
one-dimensional, and 20 repetitions in the two-dimensional case. To push this
limit to longer times, an optimization on double the time step $2\tau=0.6$ and
double the layers $m=6$ can be performed. The improvement gained for larger
time steps is traded off against a more complex classical optimization
algorithm, here (cf. Figure 4 and discussion).

The simulation of observable dynamics shall complete our excursion to the $XY$
model. In Figures 7 e. - f., the relative error of a two-point correlation
along 50 repetitions is shown. The initial state $\psi=\ket{+}^{\otimes 15}$
and observable $\braket{Z_{8}Z_{9}}$ are the same as in the analogue
discussion for the Transverse Field Ising Model. What is different is the
single time step $\tau=0.3$ and that we also compare the variational sequence
optimized on $m=6$ layers. Although the cost values of the two different
variational sequences show a clear hierarchy for multiple time steps, there is
no clear winner identifiable when it comes to the approximation of the
observable dynamics. For the two-dimensional system, the $m=6$ Trotter-Suzuki
sequence sometimes yields better approximations to the chosen observable,
whereas it is beaten in all repetitions in the one-dimensional case. This
emphasizes the volatility of observable errors, although also here an
advantageous behavior of the variational algorithms is visible. Note that
still $m=6$ Trotter-Suzuki uses double the number of gates than the other
sequences which needs to be taken into account in the case of imperfect gate
fidelities.

0.50.2 0.50.2

## References

  * [1] S. Trotzky, Y-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert, and I. Bloch.  Probing the relaxation towards equilibrium in an isolated strongly correlated one-dimensional bose gas.  Nature Physics, 8(4):325–330, Feb 2012.
  * [2] J. Eisert, M. Friesdorf, and C. Gogolin.  Quantum many-body systems out of equilibrium.  Nature Physics, 11(2):124–130, Feb 2015.
  * [3] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph Bardin, et al.  Quantum supremacy using a programmable superconducting processor.  Nature, 574:505–510, 2019.
  * [4] Masuo Suzuki.  Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics and statistical physics.  Journal of Mathematical Physics, 26(4):601–612, 1985.
  * [5] Masuo Suzuki.  General theory of fractal path integrals with applications to many‐body theories and statistical physics.  Journal of Mathematical Physics, 32(2):400–407, 1991.
  * [6] Adam Smith, M. S. Kim, Frank Pollmann, and Johannes Knolle.  Simulating quantum many-body dynamics on a current digital quantum computer.  npj Quantum Information, 5(1), Nov 2019.
  * [7] Ian D. Kivlichan, Craig Gidney, Dominic W. Berry, Nathan Wiebe, Jarrod McClean, Wei Sun, Zhang Jiang, Nicholas Rubin, Austin Fowler, Alán Aspuru-Guzik, Hartmut Neven, and Ryan Babbush.  Improved Fault-Tolerant Quantum Simulation of Condensed-Phase Correlated Electrons via Trotterization.  Quantum, 4:296, July 2020.
  * [8] Guang Hao Low and Isaac L. Chuang.  Optimal hamiltonian simulation by quantum signal processing.  Phys. Rev. Lett., 118:010501, Jan 2017.
  * [9] 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), Mar 2015.
  * [10] Andrew M. Childs, Dmitri Maslov, Yunseong Nam, Neil J. Ross, and Yuan Su.  Toward the first quantum simulation with quantum speedup.  Proceedings of the National Academy of Sciences, 115(38):9456–9461, Sep 2018.
  * [11] Andrew M. Childs and Nathan Wiebe.  Hamiltonian simulation using linear combinations of unitary operations.  Rinton Press, 12(11&12), Nov 2012.
  * [12] Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe, and Shuchen Zhu.  Theory of trotter error with commutator scaling.  Phys. Rev. X, 11:011020, Feb 2021.
  * [13] Dominic W. Berry, Graeme Ahokas, Richard Cleve, and Barry C. Sanders.  Efficient quantum algorithms for simulating sparse hamiltonians.  Communications in Mathematical Physics, 270(2):359–371, Dec 2006\.
  * [14] 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, Jan 2021.
  * [15] John Preskill.  Quantum Computing in the NISQ era and beyond.  Quantum, 2:79, August 2018.
  * [16] Arthur G. Rattew, Shaohan Hu, Marco Pistoia, Richard Chen, and Steve Wood.  A domain-agnostic, noise-resistant, hardware-efficient evolutionary variational quantum eigensolver, 2020.
  * [17] Ying Li and Simon C. Benjamin.  Efficient variational quantum simulator incorporating active error minimization.  Phys. Rev. X, 7:021050, Jun 2017.
  * [18] Cristina Cîrstoiu, Zoë Holmes, Joseph Iosue, Lukasz Cincio, Patrick J. Coles, and Andrew Sornborger.  Variational fast forwarding for quantum simulation beyond the coherence time.  npj Quantum Information, 6(1), Sep 2020.
  * [19] Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik.  The theory of variational hybrid quantum-classical algorithms.  New Journal of Physics, 18(2):023023, Feb 2016.
  * [20] Earl Campbell.  Shorter gate sequences for quantum computing by mixing unitaries.  Phys. Rev. A, 95:042306, Apr 2017.
  * [21] Earl Campbell.  Random compiler for fast hamiltonian simulation.  Phys. Rev. Lett., 123:070503, Aug 2019.
  * [22] Yingkai Ouyang, David R. White, and Earl T. Campbell.  Compilation by stochastic Hamiltonian sparsification.  Quantum, 4:235, February 2020.
  * [23] Paul K. Faehrmann, Mark Steudtner, Richard Kueng, Maria Kieferova, and Jens Eisert.  Randomizing multi-product formulas for improved hamiltonian simulation, 2021.
  * [24] Andrew M. Childs, Aaron Ostrander, and Yuan Su.  Faster quantum simulation by randomization.  Quantum, 3:182, September 2019.
  * [25] Chi-Fang Chen, Hsin-Yuan Huang, Richard Kueng, and Joel A. Tropp.  Concentration for random product formulas.  PRX Quantum, 2(4), Oct 2021.
  * [26] F. Barratt, James Dborin, Matthias Bal, Vid Stojevic, Frank Pollmann, and A. G. Green.  Parallel quantum simulation of large systems on small nisq computers.  npj Quantum Information, 7(1):79, May 2021.
  * [27] Benjamin D. M. Jones, David R. White, George O. O’Brien, John A. Clark, and Earl T. Campbell.  Optimising Trotter-Suzuki decompositions for quantum simulation using evolutionary strategies.  In Proceedings of the Genetic and Evolutionary Computation Conference, GECCO ’19, page 1223–1231, New York, NY, USA, 2019. Association for Computing Machinery.
  * [28] Stefano Barison, Filippo Vicentini, and Giuseppe Carleo.  An efficient quantum algorithm for the time evolution of parameterized circuits.  Quantum, 5:512, Jul 2021.
  * [29] Marcello Benedetti, Mattia Fiorentini, and Michael Lubasch.  Hardware-efficient variational quantum algorithms for time evolution.  Physical Review Research, 3(3), Jul 2021.
  * [30] Jonathan Wei Zhong Lau, Tobias Haug, Leong Chuan Kwek, and Kishor Bharti.  NISQ algorithm for hamiltonian simulation via truncated taylor series, 2021.
  * [31] Joe Gibbs, Kaitlin Gili, Zoë Holmes, Benjamin Commeau, Andrew Arrasmith, Lukasz Cincio, Patrick J. Coles, and Andrew Sornborger.  Long-time simulations with high fidelity on quantum hardware, 2021.
  * [32] Benjamin Commeau, M. Cerezo, Zoë Holmes, Lukasz Cincio, Patrick J. Coles, and Andrew Sornborger.  Variational hamiltonian diagonalization for dynamical quantum simulation, 2020.
  * [33] Kentaro Heya, Ken M Nakanishi, Kosuke Mitarai, and Keisuke Fujii.  Subspace variational quantum simulator, 2019.
  * [34] Thomas Barthel and Yikang Zhang.  Optimized Lie–Trotter–Suzuki decompositions for two and three non-commuting terms.  Annals of Physics, 418:168165, 2020.
  * [35] Matthew P Harrigan, Kevin J Sung, Matthew Neeley, Kevin J Satzinger, Frank Arute, Kunal Arya, Juan Atalaya, Joseph C Bardin, Rami Barends, Sergio Boixo, et al.  Quantum approximate optimization of non-planar graph problems on a planar superconducting processor.  Nature Physics, 17(3):332–336, 2021.
  * [36] Kishor Bharti, Alba Cervera-Lierta, Thi Ha Kyaw, Tobias Haug, Sumner Alperin-Lea, Abhinav Anand, Matthias Degroote, Hermanni Heimonen, Jakob S. Kottmann, Tim Menke, Wai-Keong Mok, Sukin Sim, Leong-Chuan Kwek, and Alán Aspuru-Guzik.  Noisy intermediate-scale quantum algorithms.  Rev. Mod. Phys., 94:015004, Feb 2022.
  * [37] M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles.  Variational quantum algorithms.  Nature Reviews Physics, 3(9):625–644, aug 2021.
  * [38] Jules Tilly, Hongxiang Chen, Shuxiang Cao, Dario Picozzi, Kanav Setia, Ying Li, Edward Grant, Leonard Wossnig, Ivan Rungger, George H. Booth, and Jonathan Tennyson.  The variational quantum eigensolver: a review of methods and best practices, 2021.
  * [39] Román Orús.  Tensor networks for complex quantum systems.  Nature Reviews Physics, 1(9):538–550, aug 2019.
  * [40] Cirq Developers.  Cirq, March 2021.  See full list of authors on Github: https://github.com/quantumlib/Cirq/graphs/ contributors.
  * [41] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  CoRR, abs/1412.6980, 2015.
  * [42] Y. Salathé, M. Mondal, M. Oppliger, J. Heinsoo, P. Kurpiers, A. Potočnik, A. Mezzacapo, U. Las Heras, L. Lamata, E. Solano, S. Filipp, and A. Wallraff.  Digital quantum simulation of spin models with circuit quantum electrodynamics.  Phys. Rev. X, 5:021027, Jun 2015.
  * [43] Michele C. Collodo, Johannes Herrmann, Nathan Lacroix, Christian Kraglund Andersen, Ants Remm, Stefania Lazar, Jean-Claude Besse, Theo Walter, Andreas Wallraff, and Christopher Eichler.  Implementation of conditional phase gates based on tunable $zz$ interactions.  Phys. Rev. Lett., 125:240502, Dec 2020.
  * [44] B. Foxen, C. Neill, A. Dunsworth, P. Roushan, B. Chiaro, et al.  Demonstrating a continuous set of two-qubit gates for near-term quantum algorithms.  Phys. Rev. Lett., 125:120504, Sep 2020.
  * [45] Mahdi Sameti and Michael J. Hartmann.  Floquet engineering in superconducting circuits: From arbitrary spin-spin interactions to the kitaev honeycomb model.  Phys. Rev. A, 99:012333, Jan 2019.
  * [46] Nathan Lacroix, Christoph Hellings, Christian Kraglund Andersen, Agustin Di Paolo, Ants Remm, Stefania Lazar, Sebastian Krinner, Graham J. Norris, Mihai Gabureac, Johannes Heinsoo, Alexandre Blais, Christopher Eichler, and Andreas Wallraff.  Improving the performance of deep quantum optimization algorithms with continuous gate sets.  PRX Quantum, 1:110304, Oct 2020.
  * [47] Morten Kjaergaard, Mollie E. Schwartz, Jochen Braumüller, Philip Krantz, Joel I.-J. Wang, Simon Gustavsson, and William D. Oliver.  Superconducting qubits: Current state of play.  Annual Review of Condensed Matter Physics, 11(1):369–395, 2020\.
  * [48] Markus Heyl, Philipp Hauke, and Peter Zoller.  Quantum localization bounds Trotter errors in digital quantum simulation.  Science Advances, 5(4), 2019.
  * [49] Haruo Yoshida.  Construction of higher order symplectic integrators.  Physics Letters A, 150(5):262–268, 1990.
  * [50] Naomichi Hatano and Masuo Suzuki.  Finding Exponential Product Formulas of Higher Orders, pages 37–68.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2005.