# Adaptive projected variational quantum dynamics

David Linteau<EMAIL_ADDRESS>Institute of Physics, École Polytechnique
Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Center for Quantum
Science and Engineering, École Polytechnique Fédérale de Lausanne (EPFL),
CH-1015 Lausanne, Switzerland  Stefano Barison  Institute of Physics, École
Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Center for Quantum Science and Engineering, École Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  National Centre for
Computational Design and Discovery of Novel Materials MARVEL, EPFL, Lausanne,
Switzerland  Netanel H. Lindner  Physics Department, Technion, 320003 Haifa,
Israel  Giuseppe Carleo  Institute of Physics, École Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Center for Quantum Science and
Engineering, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland  National Centre for Computational Design and Discovery
of Novel Materials MARVEL, EPFL, Lausanne, Switzerland

###### Abstract

We propose an adaptive quantum algorithm to prepare accurate variational time
evolved wave functions. The method is based on the projected Variational
Quantum Dynamics (pVQD) algorithm, that performs a global optimization with
linear scaling in the number of variational parameters. Instead of fixing a
variational ansatz at the beginning of the simulation, the circuit is grown
systematically during the time evolution. Moreover, the adaptive step does not
require auxiliary qubits and the gate search can be performed in parallel on
different quantum devices. We apply the new algorithm, named Adaptive pVQD, to
the simulation of driven spin models and fermionic systems, where it shows an
advantage when compared to both Trotterized circuits and non-adaptive
variational methods. Finally, we use the shallower circuits prepared using the
Adaptive pVQD algorithm to obtain more accurate measurements of physical
properties of quantum systems on hardware.

##  I Introduction

Simulation of static and dynamic properties of quantum systems is a
notoriously hard task for classical computers. While analytical solutions are
available only for specific cases, the amount of time and computing resources
required in general by exact numerical methods is exponential in the system
size, making the calculations quickly unfeasible. While several approximated
many-body numerical techniques have been proposed [1, 2, 3, 4], the accurate
description of important physical and chemical phenomena is a very active
research problem [5, 6, 7, 8].

In recent years, quantum computers have seen significant developments [9, 10,
11], opening potential opportunities for scientific discoveries. Hardware
capabilities continue to advance steadily, and we can already create and
manipulate complex many-body quantum systems [12, 13, 14, 15, 16, 17].
However, large-scale fault-tolerant quantum computers remain far in the
future, and contemporary devices show limitations in connectivity, size, and
coherence times.

Accounting for these constraints, Variational Quantum Algorithms (VQAs) have
emerged as the leading strategy to take advantage of near-term quantum devices
[18, 19, 20, 21]. In this class of algorithms, the solution of a given problem
(e.g. finding the ground state of a physical system) is encoded in a quantum
circuit that depends on some parameters optimized with the aid of a classical
device. VQAs have not only been proposed for quantum simulations but also for
a variety of different applications, such as machine learning [22, 23],
combinatorial optimization [24, 25], quantum error correction [26, 27] and
compilation [28, 29, 30]. Variational schemes have also been introduced in
quantum dynamics [31, 32, 33, 34, 35, 36, 37, 38, 39], as a more efficient
alternative to Trotterization [40, 41, 42, 43, 44]. The accuracy of a
variational quantum simulation is then tied to the ability of a parameterized
circuit to describe time-evolved wave functions. Even if the initial wave
function is well-described by the chosen parameterized circuit, the complexity
of the time-evolved wave functions varies with time and the chosen circuit may
fail to describe them. The choice of the parameterized circuit is therefore
crucial and it remains an open problem in variational simulations of quantum
dynamics.

Adaptive schemes have been proposed in the context of variational ground state
search [45, 46, 47, 48] especially to avoid committing to a particular
parameterized circuit. The key idea is to construct the parameterized circuit
during optimization. By systematically appending specific quantum gates to the
parameterized circuit, adaptive methods have been shown to surpass standard
approaches in the number of operations required and in the accuracy of the
final results. Moreover, adaptive methods provide flexible circuits suited for
dynamics simulations [33, 49]. However, including an adaptive step for
dynamics usually requires measurements of additional quantities, that might be
difficult to perform, or auxiliary qubits.

In this work, we introduce an adaptive variational algorithm for real-time
evolution based on the projected Variational Quantum Dynamics (pVQD) algorithm
[36], denoted Adaptive pVQD. The method inherits all the properties of the
original pVQD algorithm and integrates the adaptive modification of the
parameterized circuit without requiring auxiliary qubits. The structure of
this paper is as follows: in Section II we present the algorithm and describe
how the adaptive routine is performed; in Section III we apply the method to
study a time-dependent and a fermionic system, benchmarking the method against
Trotter evolution and the original pVQD algorithm; Section IV concludes the
paper with some considerations and outlooks on the proposed method.

##  II Method

Figure 1: Flowchart of the time evolution of the Adaptive pVQD algorithm.
Starting with a parameter-free circuit, we discretize the time evolution into
multiple time steps. At each time step we optimize the parameters to
approximate the real time evolution of the quantum system. If the optimization
does not converge to the required accuracy, or the ansatz does not contain any
parameter, then rotations $\\{R_{A_{i}^{*}}\\}$ based on the generators
$\\{A_{i}^{*}\\}$ are appended to the circuit according to the adaptive step
procedure described in Section II.1. The algorithm stops once the final time
$t_{f}$ is reached.

Consider a physical system governed by a Hamiltonian $H$. For clarity of
exposition, we focus on time-independent Hamiltonians. However, this is not a
requirement of the algorithm, as we explicitly show in Section III. To
simulate the dynamics of quantum systems on a quantum computer, we have to
prepare the time-evolved wave function $|\Psi(t)\rangle=U(t)|\psi_{0}\rangle$,
where $|\psi_{0}\rangle=U_{0}|0\rangle^{\otimes N}$ is the initial state, $N$
indicates the number of qubits representing the physical system and $U(t)$ is
the unitary time evolution operator. The Adaptive pVQD algorithm aims to
approximate the state $|\Psi(t)\rangle$ using parameterized states of the form

$|\psi(\bm{\theta},\bm{A})\rangle=U(\bm{\theta},\bm{A})|\psi_{0}\rangle=\prod_{i}e^{-i\theta_{i}A_{i}}|\psi_{0}\rangle,$
(1)

where each real parameter $\theta_{i}\in\bm{\theta}$ is associated to a
Hermitian generator $A_{i}\in\bm{A}$. The parameterized state is therefore
specified by the set of parameters and operators $\\{\bm{\theta},\bm{A}\\}$,
and it can be implemented as a quantum circuit. From now on, we adopt the
notation $|\psi(\bm{\theta})\rangle\equiv|\psi(\bm{\theta},\bm{A})\rangle$ and
$U(\bm{\theta})\equiv U(\bm{\theta},\bm{A})$.

To simulate a physical model until a final time $t_{f}$, we divide the
evolution into small time intervals $\Delta t$. We further assume that the
parameterized state $|\psi(\bm{\theta})\rangle$ is a good approximation of the
time-evolved wave function at time $t$. The wave function at time $t+\Delta t$
can thus be represented by $U_{\text{TS}}(\Delta t)|\psi(\bm{\theta})\rangle$,
where $U_{\text{TS}}(\Delta t)$ is a Trotter-Suzuki decomposition of the time
evolution operator $U(\Delta t)$ [40, 41]. In this manuscript we use a first
order decomposition, but higher orders can be considered. The choice of the
optimal $\Delta t$ is problem dependent and will be discussed in Section III.
We then approximate the evolution step $t\to t+\Delta t$ using a new set of
parameters $\bm{\theta}\to\bm{\theta}+\bm{d\theta}$ that maximizes the overlap
between $U_{\text{TS}}(\Delta t)|\psi(\bm{\theta})\rangle$ and
$|\psi(\bm{\theta}+\bm{d\theta})\rangle$. This can be achieved by minimizing,
with respect to $\bm{d\theta}$, the infidelity

$\mathcal{I}(\bm{d\theta},\Delta t)=1-\mathcal{F}(\bm{d\theta},\Delta t),$ (2)

where the fidelity

$\mathcal{F}(\bm{d\theta},\Delta
t)=|\langle\psi(\bm{\theta}+\bm{d\theta})|U_{\text{TS}}(\Delta
t)|\psi(\bm{\theta})\rangle|^{2}$ (3)

can be measured on a quantum device [36].

At each time step, the initial parameters and operators
$\\{\bm{\theta},\bm{A}\\}$ are those obtained at the previous time step.
Assuming that the set of operators $\bm{A}$ is sufficient to describe the
state at time $t+\Delta t$, we find the parameter shift $\bm{d\theta}^{*}$
that minimizes $\mathcal{I}(\bm{d\theta},\Delta t)$. Details about the
minimization routine can be found in Appendix A. If the minimization routine
is not successful, new gates built using generators
$(A_{0}^{*},A_{1}^{*},\cdots,A_{k}^{*})$ from the operator pool are added to
the parameterized circuit following the adaptive procedure described in
Section II.1. This adaptive procedure is repeated up until the convergence
criteria are met.

The algorithm starts with the initial state $|\psi_{0}\rangle$ represented by
an empty set of operators. As needed, new gates are added through the time
evolution until the chosen final time $t_{f}$. The complete procedure is
illustrated in Fig. 1. We note that the original pVQD scheme [36] can be
recovered by fixing the set of operators $\bm{A}$ through the entire
simulation.

###  II.1 Adaptive step

When the parameterized circuit $|\psi(\bm{\theta})\rangle$ is not expressive
enough to accurately describe the time step evolution by only shifting the
variational parameters, we add new gates to it. This is referred to as the
adaptive step of the algorithm. Given an operator pool, we determine the best
gate to grow the quantum circuit. As first proposed in [45], we look for the
operator whose gate maximizes the derivative of the cost function with respect
to its parameter. This is achieved by iterating over all the operators in the
pool, a step that can be performed in parallel even on different quantum
devices.

For ground state methods, the cost function is the energy of the system, while
the gradient is obtained by measuring the expectation value of the commutator
between the trial operator and the Hamiltonian [45, 50]. We must ensure that
is possible to apply a similar procedure when dynamics is considered. In the
adaptive scheme proposed in [33], this step requires an additional measurement
of the variance of the Hamiltonian with respect to the non-adaptive case. In
our method, the gradient of the fidelity with respect to the shift
$d\theta_{a}$ of parameter $\theta_{a}$ associated with a trial operator
$A_{a}$ has the form

$\frac{\partial\mathcal{F}}{\partial d\theta_{a}}=\bra{\phi(\bm{\theta},\Delta
t)}e^{-id\theta_{a}A_{a}}[P_{0},iA_{a}]e^{id\theta_{a}A_{a}}\ket{\phi(\bm{\theta},\Delta
t)},$ (4)

where we define the projector $P_{0}=\outerproduct{\psi_{0}}{\psi_{0}}$ and
the state $\ket{\phi(\bm{\theta},\Delta
t)}=U^{\dagger}(\bm{\theta})U_{\text{TS}}(\Delta t)\ket{\psi(\bm{\theta})}$
(see Appendix B for the full derivation). To ensure continuity of time
evolution, we initially set $\theta_{a}=0$. We note that measuring the
derivative of the fidelity corresponds to measuring the Hermitian operator
$[P_{0},iA_{a}]$ with respect to the pVQD circuit
$U^{\dagger}(\bm{\theta})U_{\text{TS}}(\Delta t)\ket{\psi(\bm{\theta})}$
modified by the addition of the gate $e^{id\theta_{a}A_{a}}$. However, we
evaluate the derivative using the parameter shift rule [51], as for the
minimization routine (for more details, see Appendix A). This operator search
is still parallelizable on multiple devices and does not require auxiliary
qubits.

The adaptive step has been lately extended and optimized [46, 48, 50], with
new protocols that greatly reduce the computational resources required with
respect to the first proposal. In particular, we adopt the scheme presented in
[48], which increases the depth of the parameterized circuit
$|\psi(\bm{\theta})\rangle$ by 1 at every adaptive step. While the infidelity
defined in Eq. 2 remains above a fixed threshold $\varepsilon$, additional
adaptive steps are performed. For a detailed description, see Appendix C.

###  II.2 Operator pool

The choice of the operator pool is a key ingredient in the success and
efficiency of adaptive variational algorithms. Having a complete pool of
operators is exponentially complex in the size of the physical system,
therefore, one has to make some restrictions in its selection. Many different
strategies have been proposed, such as the creation of a minimally complete
pool [46, 52], the inclusion of symmetries directly in the operator pool [53],
or the extension of a complete pool acting on a subsystem of the studied model
[47].

In the study of the dynamics, we can refer to the Trotterization of the time
evolution operator to select the pool. In particular, we consider local (L)
and non-local (NL) operator pools, respectively, given by

$\displaystyle\mathcal{A}_{\text{L}}$
$\displaystyle=\\{X_{i},Y_{i},Z_{i}\\}_{i=0}^{N-1}$ (5) $\displaystyle\ \ \ \
\cup\\{X_{i}X_{i+1},Y_{i}Y_{i+1},Z_{i}Z_{i+1}\\}_{0\leq i\leq N-2},$
$\displaystyle\mathcal{A}_{\text{NL}}$
$\displaystyle=\\{X_{i},Y_{i},Z_{i},X_{i}X_{j},Y_{i}Y_{j},Z_{i}Z_{j}\\}_{0\leq
i<j\leq N-1},$ (6)

where $X_{i},Y_{i}$ and $Z_{i}$ are the Pauli gates acting on site $i$. Given
that $\mathcal{A}_{\text{L}}\subset\mathcal{A}_{\text{NL}}$, we expect that
$\mathcal{A}_{\text{NL}}$ will generate more flexible parameterized states.
However, not only the choice of $\mathcal{A}_{\text{NL}}$ leads to a
measurement overhead, but the non-local nature of this pool may add long-range
controlled-NOT (CNOT) gates to the circuit, according to the device
connectivity. In Section III, we report the comparison of the two pools in the
study of a fermionic system.

##  III Results

We apply the Adaptive pVQD method to the study of the 1D Heisenberg XYZ model
with an external driving field and the 2D Fermi-Hubbard model. Both have non-
trivial dynamics and open the pVQD method to the study of time-dependent and
fermionic systems. In both cases, open boundary conditions were imposed.

###  III.1 Driven Heisenberg model

Figure 2:  Dynamics of the driven Heisenberg XYZ model studied with the
Adaptive pVQD algorithm with local pool (L), compared to standard Trotter
evolution, pVQD and pVQD with block extensions. The plot shows the results for
an open chain of $L=8$ spins with $J_{x}=1,J_{y}=0.8$ and $J_{z}=0.6$. The top
and middle panel show the measurements of a single spin observable and a
correlator, respectively. The bottom panel shows the number of CNOTs in the
circuit describing the time-evolved wave function. The simulation started in
the antiferromagnetic state $|\psi_{0}\rangle=|01010101\rangle$, and the
infidelity threshold was set to $\varepsilon=10^{-4}$ for all the variational
methods.

Given an open chain of $L$ spins, the driven Heisenberg XYZ Hamiltonian can be
written as:

$\displaystyle
H(t)=\sum_{i=0}^{L-2}(J_{x}X_{i}X_{i+1}+J_{y}Y_{i}Y_{i+1}+J_{z}Z_{i}Z_{i+1})+D(t)$
(7)

where $J_{x},J_{y}$ and $J_{z}$ are coupling parameters and $D(t)$ is the
time-dependent driving term. Many different driving terms can be applied to
the system. Among those we choose

$\displaystyle D(t)$ $\displaystyle=\sum_{i=0}^{L-1}(-1)^{i}\sin(\omega
t)Z_{i}\,,$ (8)

where $\omega$ is the driving frequency.

First, we investigate the performance of the Adaptive pVQD algorithm with a
local pool on a perfect simulator and compare to Trotterized circuits and the
original implementation of pVQD. We consider $J_{x}=1,J_{y}=0.8,J_{z}=0.6$, an
antiferromagnetic initial state $|\psi_{0}\rangle=|0101\rangle$ and a final
evolution time $t_{f}=2$. In the classic version of the pVQD algorithm, we
have to choose an ansatz for the time evolved wave function. We consider a
circuit equivalent to a Trotter step where all the rotations are defined by
variational parameters. The Trotter step circuit implementation for this model
is shown in Appendix E. Both the Trotter and the pVQD full circuits are then
obtained repeating this structure $n_{\text{TS}}$ times. In particular, we fix
$n_{\text{TS}}=10$ for the Trotter circuit and $n_{\text{TS}}=3$ for the pVQD
ansatz.

After running the algorithms, we compare the different circuits obtained and
use them to measure expectation values of single- and two-spin observables.
The results are shown in Fig. 2. The Trotter circuit lags behind variational
methods both in terms of accuracy and resource required. The pVQD method
instead achieves accurate results until $t=1.0$, where the associated circuit
becomes shallower than the one of Adaptive pVQD. This phenomenon suggests that
in that time step the fixed representation power is the main source of error
in the variational calculations.

In order to show the flexibility of the Adaptive pVQD, we implement a naive
modification of the pVQD algorithm, that we indicate as pVQD with block
extensions. In this case, a new step of the Trotterized variational ansatz is
added to the circuit once the optimization procedure does not reach the
desired accuracy. While this approach does improve the performance of the pVQD
algorithm, we remark that it is not general, as it depends on the ansatz
structure we have chosen. Furthermore, we can see from the bottom panel of
Fig. 2 that the Adaptive pVQD method always produces shallower circuits, with
resources tailored to the needs of the specific time step.

Figure 3:  Adaptive pVQD algorithm with local pool compared to standard
Trotter evolution and pVQD for the driven Heisenberg XYZ model. We employ the
same setup indicated in Fig. 2 for multiple systems of size
$L\in\left[3,11\right]$. The top panel shows the integrated exact infidelity
of pVQD and Trotterization over an entire time evolution with final time
$t_{f}=2$ as a function of the system size. The bottom panel shows the circuit
depth at the end of the time evolution.  Figure 4:  Observables measured with
the IBM Manila device for the driven Heisenberg XYZ model on an open chain
with 4 sites, $J_{x}=1,J_{y}=0.8,J_{z}=0.6$ and an antiferromagnetic initial
state $|\psi_{0}\rangle=|0101\rangle$. The Trotter simulation is performed
with a fixed Trotter step size of $dt=0.2$. The Adaptive pVQD circuits
$|\psi(\bm{\theta})\rangle$ were obtained with a noiseless simulation that
used a local operator pool. The shaded areas correspond to 50 noisy
simulations using the noise model of IBM Manila. Each data point and error bar
correspond to the mean and the standard deviation, respectively, of 50
experiments performed on hardware. Zero noise extrapolation was applied to
both noisy simulations and hardware experiments. Idle qubits were also
dynamically decoupled from the active ones.

Then, we extend the study to systems with different sizes. To this end, we
define the integrated exact infidelity

$\Delta_{\mathcal{I}}^{\mathrm{ex}}(t_{f})=\int_{0}^{t_{f}}\left(1-|\langle\Psi(t)|\psi(\bm{\theta})\rangle|^{2}\right)dt$
(9)

with respect to the exact wave function $\ket{\Psi(t)}$ computed on a
classical device. We again fix a final evolution time $t_{f}=2$ and evaluate
$\Delta_{\mathcal{I}}^{\mathrm{ex}}(t_{f})$ for each method for systems of
$L\in\left[3,11\right]$ spins. In particular, we consider a Trotter circuit
with a fixed depth of $n_{\mathrm{TS}}=10$ and one with fixed Trotter step
size $dt=J_{x}t/n_{\mathrm{TS}}=0.05$, the same we use in the Trotter step of
the pVQD algorithm. The results are shown in Fig. 3, together with the circuit
depth at the end of the time evolution.

We note that the depth of the Adaptive pVQD circuits increases with the system
size and converges to the Trotter circuit with fixed depth, while having a
lower integrated exact infidelity. We highlight that Fig. 3 only indicates the
depth of the final circuit. In the case of Adaptive pVQD, this corresponds to
the deepest circuit prepared. The Trotterized circuits with a fixed Trotter
step size yield the lowest values for $\Delta_{\mathcal{I}}^{\mathrm{ex}}$,
but $n_{\mathrm{TS}}=40$ Trotter steps are required to evolve the system to
$t_{f}=2$, resulting in circuits almost one order of magnitude deeper than any
other. We performed multiple pVQD simulations with different variational
ansätze equivalent to $n_{\mathrm{TS}}=1,2,3,8$ Trotter steps. We note that
the integrated exact infidelities of pVQD with $n_{\mathrm{TS}}=1,2,3$ have a
steep transition when the number of gates becomes smaller than the adaptive
circuit. This phenomenon suggests that the ansatz limitation is the main
source of error in the variational calculations, while the adaptive circuit is
able to increase effectively its representation power. On the other hand, the
standard pVQD calculation with $n_{\mathrm{TS}}=8$ never undergoes this
transition. While the integrated exact infidelity is always lower than the
adaptive approach, we have to note that the entire time evolution is performed
with a deeper circuit. Finally, we note a plateau in the depth of the circuit
required by the adaptive algorithm when $L>8$. This is similar to what
observed in [33], where the system size at which the number of gates required
saturates depends on the evolution time.

The adaptive method is able to produce circuits that are orders of magnitude
shallower than Trotterization while keeping the accuracy comparable to it.
Those circuit can be used to improve the measurement of observables at long
times on current quantum devices, which are otherwise limited by the depth of
the Trotterization. For this reason, we first run the Adaptive pVQD algorithm
on the simulator and use the resulting sets of variational parameters to
prepare quantum circuit on the hardware for a system of $L=4$ spins. In Fig.
4, we compare observables measured both on those variational wave functions
and on Trotterized circuits with a fixed Trotter step size of $dt=0.2$.

In this experiment, the final Trotter circuit has 180 CNOTs. This circuit is
beyond what is currently accessible on quantum devices, settling the
expectation value of the correlator close to $0$ for $J_{x}t>0.8$. On the
other hand, the Adaptive pVQD parameterized circuit
$|\psi(\bm{\theta})\rangle$ has 28 CNOTs at the end of the evolution. This
improvement in the number of gates is crucial for the application of error
mitigation techniques, especially at longer times. In particular, zero noise
extrapolation (ZNE [54, 31]) was applied both on the noisy simulations and
hardware experiments. We choose a quadratic fit on values obtained with noise
scaling factors $\left[1,2,3\right]$. Moreover, when running our algorithm on
hardware, we dynamically decouple the idle qubits from the active ones using
the standard procedure available in Qiskit [55]. We expect that more advanced
noise mitigation techniques, such as the one presented in [56], will improve
the results on the Trotter circuit. However, this is also true for the
variational circuit prepared by the Adaptive pVQD.

###  III.2 Fermi-Hubbard model

The Hamiltonian of the Fermi-Hubbard model on a $L_{x}\times L_{y}$
rectangular lattice is given by

$H=-J\sum_{\langle
ij\rangle,\sigma}(c_{i\sigma}^{\dagger}c_{j\sigma}+c_{j\sigma}^{\dagger}c_{i\sigma})+U\sum_{i=0}^{L_{x}L_{y}-1}n_{i\uparrow}n_{i\downarrow},$
(10)

where $c_{i\sigma}^{\dagger}$ ($c_{i\sigma}$) is the creation (annihilation)
fermionic operator of spin $\sigma\in\\{\uparrow,\downarrow\\}$ at site $i$,
$n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ counts the number of fermions
with spin $\sigma$ at site $i$ and $\langle ij\rangle$ denotes nearest
neighbor sites on the lattice. The first term in the Hamiltonian accounts for
the hopping between nearest neighbor lattice sites, while the second term
describes the on-site interactions.

There are several ways to encode fermionic Hamiltonians into qubit operators
[57, 58, 59, 60, 61, 62, 63]. In this work, we consider the Jordan-Wigner
mapping [57] to encode each fermionic mode into a qubit. Since every lattice
site can host two modes ($\uparrow$, $\downarrow$), $N=2L_{x}L_{y}$ qubits are
required to simulate the Fermi-Hubbard model on a $L_{x}\times L_{y}$ grid.
Before performing a fermionic encoding, we eliminate the spin index via
$c_{i\uparrow}\to c_{i}$ and $c_{i\downarrow}\to c_{i+N/2}$ (and analogously
for the number operator $n_{i\sigma}$). We then map each fermionic operator
into a spin operator:

$\displaystyle c_{i}$ $\displaystyle\to Z^{\otimes
i}\otimes\sigma^{+}\otimes\mathbb{I}^{\otimes N-i-1},$ (11) $\displaystyle
c_{i}^{\dagger}$ $\displaystyle\to Z^{\otimes
i}\otimes\sigma^{-}\otimes\mathbb{I}^{\otimes N-i-1},$ (12)

where $\sigma^{\pm}=(X\pm iY)/2$. The local occupation number can then be
identified with the local spin number according to $n_{i}\in\\{0,1\\}\mapsto
Z_{i}\in\\{\uparrow,\downarrow\\}$. More details on the fermionic indexing
convention and implementing a Trotter step can be found in Appendix E.

Given that the mapping requires an ordering of the fermionic modes, operators
that are local in space might generate very long Pauli strings. For example,
considering the snake-like pattern, vertical hopping terms generate strings of
Pauli $Z$ with sizes up to $2L_{x}-2$. This represents a bottleneck in
studying fermionic systems with dimensionality higher than 1 on current
quantum devices. By restricting the operator pool, we investigate the
possibility of describing time-evolved wave functions of the 2D Hubbard model
using only local gates. We perform noiseless simulations of a $2\times 2$
square lattice, comparing local and non-local operator pools. In particular,
we measure the expectation values of a local density operator and a density
correlator and count the number of CNOTs in the circuits. We use a fixed-depth
Trotter simulation and a pVQD with block extension as a benchmark. The results
are shown in Fig. 5.

Figure 5:  Adaptive pVQD schemes for the Fermi-Hubbard model on a $2\times 2$
open square lattice (8 qubits) with $U/J=0.8$. Local (L) and non-local (NL)
operator pools are used to perform noiseless simulations and the results are
compared to a Trotter evolution with $n_{\mathrm{TS}}=5$ Trotter steps and
pVQD with block extensions. The system starts in the half-filled
antiferromagnetic state
$|\psi_{0}\rangle=|n_{0\uparrow}n_{1\uparrow}n_{2\uparrow}\cdots\rangle=|10100101\rangle$.
We fixed the infidelity threshold to $\varepsilon=10^{-4}$. The top and middle
panels show the expectation values of an on-site number density operator and a
number density correlator over time. The bottom panel shows the number of
CNOTs in the circuit describing the time-evolved wave function.

We do not restrict ourselves to specific quantum hardware to keep the
comparison as general as possible. Instead, we count the number of CNOTs in a
circuit by transpiling it into an abstract device with all-to-all connectivity
that is able to perform arbitrary single qubit rotations and CNOTs. The local
and non-local pool variants show different behavior over time in the count of
CNOTs. We note that the non-local variant always requires fewer CNOTs than its
local counterpart. However, some CNOTs are long-range, and their
implementation on an actual device can be challenging on hardware with fixed
topology and limited connectivity. In contrast, the circuit structure produced
by the local pool variant is already suited for current hardware
implementation. More details about the Adaptive pVQD output circuits can be
found in Appendix D. Moreover, the plot highlights another limitation of the
naive pVQD with block extensions approach. Indeed, it always prepare more
expensive circuits than the Adaptive pVQD with non local pool and in the end
it has similar CNOT requirement to the local variant, while being restricted
to use long range gates as required by the Trotter step.

##  IV Conclusions

We presented an adaptive version of pVQD, called Adaptive pVQD, to simulate
the real-time evolution of quantum systems. This algorithm importantly
circumvents the need to choose a fixed ansatz from the beginning of the time
evolution. The parameterized quantum circuits are grown adaptively to be both
problems and hardware-tailored. This is obtained with a measurement overhead
required to determine the best gate among those included in the operator pool.

However, the gate search can be operated in parallel and, in our scheme, does
not involve circuits with auxiliary qubits. This makes the Adaptive pVQD
algorithm more hardware-efficient than standard methods, as exemplified in
this work with the driven Heisenberg model on the IBM quantum hardware.
Finally, we have simulated the dynamics of the 2D Hubbard model with only
local gates, using the adaptive procedure to mitigate one of the bottlenecks
that current quantum devices face in studying fermionic systems. Given the
ease of introduction to the standard pVQD algorithm and its benefits, we
believe that the adaptive procedure described here can be of great use in the
simulation of dynamics both for current and future quantum devices.

## Data availability

The code used to run the simulations is open source and can be found at [64].
It was written in Python using Qiskit [55]. Exact classical simulations were
performed using Qutip [65].

## Acknowledgments

We thank S. Economou for insightful discussions. This research was supported
by the NCCR MARVEL, a National Centre of Competence in Research, funded by the
Swiss National Science Foundation (grant number 205602).

##  Appendix A Minimization routine

Here we present additional details on the minimization routine that we applied
throughout the simulations we presented in the main text. In particular, we
follow a gradient-based approach, with gradient computed using the parameter-
shift rule. Gradient-based and non-gradient-based optimization algorithms for
dynamics were previously used for instance in [36] and [37], for both ideal
and noisy quantum simulations. The parameter shift rule readily applies here
since every Pauli string $A_{i}$ is involutory, i.e. $A_{i}^{2}=\mathbb{I}$
[51]. For a fixed set of operators $\bm{A}$, the gradient of the infidelity
was thus computed via the parameter shift rule:

$\frac{\partial\mathcal{I}}{\partial
d\theta_{i}}=\frac{\mathcal{I}(\bm{\theta}+\bm{d}\theta+s\bm{e}_{i})-\mathcal{I}(\bm{\theta}+\bm{d}\theta-s\bm{e}_{i})}{2\sin
s},$ (13)

where $\bm{e}_{i}$ is the standard unit vector, and we fixed $s=\pi/2$. The
gradient was then fed to Adam [66], implemented with the default
hyperparameters and a learning rate $\alpha=0.005$. The shift parameters
$\bm{d\theta}^{*}$ were consequently obtained using Adam.

Two stopping criteria for the optimizer were used: (1) the
$\ell_{\infty}$-norm of the gradient of the infidelity is below a tolerance
and (2) a maximum number of iterations is reached. Fianlly, as showed in [36],
an optimization threshold independent from $\Delta t$ can be used if
$\mathcal{I}$ is substituted with $\mathcal{I}/\Delta t^{2}$ as cost function.

##  Appendix B Gradient of the Fidelity

In this Appendix, we derive the expression for the gradient of the adaptive
step presented in Eq. 4. Given the quantum circuit $U(\bm{\theta})$ that
prepares the state $|\psi(\bm{\theta})=U(\bm{\theta})|\psi_{0}\rangle$, we
want to add the gate $e^{-id\theta_{a}A_{a}}$ to it, defining the new state
$|\psi(\bm{\theta}+\bm{d\theta})\rangle=U(\bm{\theta})\,e^{-id\theta_{a}A_{a}}|\psi_{0}\rangle$.
To obtain the gradient of the fidelity with respect to this added parameter
$d\theta_{a}$, it is convenient to first rewrite the fidelity given in Eq. 3
as follows

$\displaystyle\mathcal{F}(\bm{d\theta},\Delta t)$
$\displaystyle=|\langle\psi(\bm{\theta}+\bm{d\theta})|U_{\text{TS}}(\Delta
t)|\psi(\bm{\theta})\rangle|^{2}$
$\displaystyle=|\langle\psi_{0}|e^{id\theta_{a}A_{a}}U^{\dagger}(\bm{\theta})U_{\text{TS}}(\Delta
t)U(\bm{\theta})|\psi_{0}\rangle|^{2}$
$\displaystyle=\langle\psi_{0}|e^{id\theta_{a}A_{a}}U^{\dagger}(\bm{\theta})U_{\text{TS}}(\Delta
t)U(\bm{\theta})|\psi_{0}\rangle$
$\displaystyle\quad*\langle\psi_{0}|U^{\dagger}(\bm{\theta})U_{\text{TS}}^{\dagger}(\Delta
t)U(\bm{\theta})e^{-id\theta_{a}A_{a}}|\psi_{0}\rangle$
$\displaystyle=\langle\psi_{0}|U^{\dagger}(\bm{\theta})U_{\text{TS}}^{\dagger}(\Delta
t)U(\bm{\theta})e^{-id\theta_{a}A_{a}}|\psi_{0}\rangle$
$\displaystyle\quad*\langle\psi_{0}|e^{id\theta_{a}A_{a}}U^{\dagger}(\bm{\theta})U_{\text{TS}}(\Delta
t)U(\bm{\theta})|\psi_{0}\rangle$
$\displaystyle=\langle\phi(\bm{\theta},\Delta
t)|e^{-id\theta_{a}A_{a}}P_{0}e^{id\theta_{a}A_{a}}|\phi(\bm{\theta},\Delta
t)\rangle,$ (14)

where we defined $|\phi(\bm{\theta},\Delta
t)\rangle=U^{\dagger}(\bm{\theta})U_{\text{TS}}(\Delta
t)U(\bm{\theta})|\psi_{0}\rangle$ and the projector
$P_{0}=|\psi_{0}\rangle\langle\psi_{0}|$. One can then readily differentiate
with respect to $d\theta_{a}$ to obtain

$\displaystyle\frac{\partial\mathcal{F}}{\partial d\theta_{a}}$
$\displaystyle=\bra{\phi(\bm{\theta},\Delta
t)}e^{-id\theta_{a}A_{a}}[P_{0},iA_{a}]e^{id\theta_{a}A_{a}}\ket{\phi(\bm{\theta},\Delta
t)}$ (15)

which precisely corresponds to Eq. 4.

##  Appendix C Adaptive step implementation

In this Appendix we illustrate the adaptive procedure we have used in our
simulations, based on what was initially proposed in [48]. The overall
procedure can be divided in the following steps:

  1. 1.

Compute the gradient of the fidelity for each operator in the pool. To process
the pool, the gate $e^{-i\theta_{a}A_{a}}$ associated to each trial operator
$A_{a}\in\mathcal{A}$ is appended one at a time to the current parameterized
circuit $\\{\bm{\theta},\bm{A}\\}$, resulting in the trial circuit
$\\{(\bm{\theta},0),(\bm{A},A_{a})\\}$. For the trajectory in parameter space
to remain continuous, the new parameter $\theta_{a}$ is set to 0. The gradient
of the fidelity with respect to the new parameter is computed for each trial
circuit using the parameter shift rule, given explicitly in Eq. 13.

  2. 2.

Pick the operator in the pool that maximizes the gradient. Update the
parameters and operators to $\bm{\theta}\to(\bm{\theta},0)$ and
$\bm{A}\to(\bm{A},A^{*})$, where $A^{*}$ is the operator $A_{a}$ that
maximizes the fidelity gradient.

  3. 3.

Remove the operators in the pool that act on qubit(s) already acted on. Given
that the operator $A^{*}$ obtained in Item 2 acts on the qubits indices
$\bm{\alpha}$, the subset of the operator pool that also acts on at least one
index in $\bm{\alpha}$, namely

$\mathcal{A}_{\bm{\alpha}}=\\{A_{a}|A_{a}\in\mathcal{A}\ \text{acts on}\
\bm{\beta},\bm{\beta}\cup\bm{\alpha}\neq\emptyset\\}$ (16)

should be removed from the current operator pool. Hence the pool can be
updated as follows:
$\mathcal{A}\to\mathcal{A}\setminus\mathcal{A}_{\bm{\alpha}}$.

  4. 4.

Go back to Item 2 until the operator pool is empty.

  5. 5.

Return the new circuit. The new parameterized circuit is characterized by
$\bm{\theta}\to(\bm{\theta},0,\cdots,0)$ and
$\bm{A}\to(\bm{A},A_{0}^{*},A_{1}^{*},\cdots,A_{k}^{*})$, assuming that $k$
new operators were added.

As stated in the main text, this procedure guarantees that the depth of the
parameterized circuit $|\psi(\bm{\theta})\rangle$ is increased by 1 in each
adaptive step [48].

##  Appendix D Adaptive pVQD output circuits

We illustrate in Figs. 6 and 7 examples of parameterized circuits obtained
with the Adaptive pVQD algorithm in simulations shown in the main text. Each
column of operators in the circuits corresponds to an adaptive step.

Figure 6:  Variational circuit obtained at $J_{x}t=2$ in the simulation shown
in Fig. 3, using the Adaptive pVQD algorithm and local operator pool.  Figure
7:  Variational circuit obtained at $Jt=4$ in the simulation shown in Fig. 5,
using the Adaptive pVQD algorithm and local operator pool.

##  Appendix E Trotter step circuit encodings

In this Appendix we provide the circuits we used to implement a single Trotter
step of the driven Heisenberg and the Hubbard models. The Trotter step in the
driven Heisenberg model is implemented with a checkerboard pattern of the two
qubit gates $R_{XX},R_{YY},R_{ZZ}$, with a layer of single qubit $R_{Z}$ at
the end. We show a sketch in Fig. 8.

Figure 8: Implementation of an antiferromagnetic initial state and a Trotter
step for the driven Heisenberg model given in Eq. 7.

To realize the Trotter circuit for the Hubbard model, we first have to
establish an ordering in the latices sites and the modes. We number the sites
using a snake-like pattern and, as indicated in the main text, we eliminate
the spin index via $c_{i\uparrow}\to c_{i}$ and $c_{i\downarrow}\to
c_{i+N/2}$. Under this ordering, the Jordan-Wigner transformation of the
Hamiltonian terms reads

$\displaystyle
c_{i\uparrow}^{\dagger}c_{j\uparrow}+c_{j\uparrow}^{\dagger}c_{i\uparrow}$
$\displaystyle\mapsto\frac{1}{2}\left[X_{i}\prod_{k=i+1}^{j-1}Z_{k}X_{j}+Y_{i}\prod_{k=i+1}^{j-1}Z_{k}Y_{j}\right],$
(17) $\displaystyle
c_{i\downarrow}^{\dagger}c_{j\downarrow}+c_{j\downarrow}^{\dagger}c_{i\downarrow}$
$\displaystyle\mapsto\frac{1}{2}\Bigg{[}X_{i+N/2}\prod_{k=i+1}^{j-1}Z_{k+N/2}X_{j+N/2}\,+$
(18)
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,+Y_{i+N/2}\prod_{k=i+1}^{j-1}Z_{k+N/2}Y_{j+N/2}\Bigg{]},$
$\displaystyle n_{i\uparrow}n_{i\downarrow}$
$\displaystyle\mapsto\frac{1}{4}(\mathbb{I}-Z_{i})(\mathbb{I}-Z_{i+N/2}),$
(19)

where we assumed $j>i$ without loss of generality. Given the mapped
Hamiltonian, the Trotter step can not be implemented using only
$R_{XX},R_{YY},R_{ZZ}$ and $R_{Z}$ gates. Indeed, the non locality of the
mapping requires some multi-qubit rotation with size up to $2L_{x}$. The two
multi-qubit gates are the rotations generated by the Pauli strings $XZZX$ and
$YZZY$, which can be decomposed as shown in [44]. Fig. 9 presents our
implementation.

Figure 9:  (a) Gates used to define a Trotter step. (b) Quantum circuit
encoding the first order Trotter step of the Hubbard model with an half-filled
antiferromagnetic initial state.

## References

  * Sandvik  _et al._ [2010] A. W. Sandvik, A. Avella, and F. Mancini, Computational studies of quantum spin systems, in  _AIP Conference Proceedings_ (AIP, 2010).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Orús [2019] R. Orús, Tensor networks for complex quantum systems, Nature Reviews Physics 1, 538 (2019).
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová , Machine learning and the physical sciences, Reviews of Modern Physics 91, 10.1103/revmodphys.91.045002 (2019).
  * Cady  _et al._ [2008] C. W. Cady, R. H. Crabtree, and G. W. Brudvig, Functional models for the oxygen-evolving complex of photosystem ii, Coordination Chemistry Reviews 252, 444 (2008), the Role of Manganese in Photosystem II.
  * Schimka  _et al._ [2010] L. Schimka, J. Harl, A. Stroppa, A. Grüneis, M. Marsman, F. Mittendorfer, and G. Kresse, Accurate surface and adsorption energies from many-body perturbation theory, Nature Materials 9, 741 (2010).
  * Leggett [2006] A. J. Leggett, What DO we know about high Tc?, Nature Physics 2, 134 (2006).
  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Arute  _et al._ [2019] F. Arute _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Zhong  _et al._ [2020] H.-S. Zhong _et al._ , Quantum computational advantage using photons, Science 370, 1460 (2020).
  * Huang  _et al._ [2022] H.-Y. Huang, M. Broughton, J. Cotler, S. Chen, J. Li, M. Mohseni, H. Neven, R. Babbush, R. Kueng, J. Preskill, and J. R. McClean, Quantum advantage in learning from experiments, Science 376, 1182 (2022).
  * Dolde  _et al._ [2014] F. Dolde, V. Bergholm, Y. Wang, I. Jakobi, B. Naydenov, S. Pezzagna, J. Meijer, F. Jelezko, P. Neumann, T. Schulte-Herbrüggen, J. Biamonte, and J. Wrachtrup, High-fidelity spin entanglement using optimal control, Nature Communications 5, 10.1038/ncomms4371 (2014).
  * Waldherr  _et al._ [2014] G. Waldherr, Y. Wang, S. Zaiser, M. Jamali, T. Schulte-Herbrüggen, H. Abe, T. Ohshima, J. Isoya, J. F. Du, P. Neumann, and J. Wrachtrup, Quantum error correction in a solid-state hybrid spin register, Nature 506, 204 (2014).
  * Nam  _et al._ [2019] Y. Nam _et al._ , Ground-state energy estimation of the water molecule on a trapped ion quantum computer (2019).
  * Wan  _et al._ [2020] Y. Wan, R. Jördens, S. D. Erickson, J. J. Wu, R. Bowler, T. R. Tan, P.-Y. Hou, D. J. Wineland, A. C. Wilson, and D. Leibfried, Ion transport and reordering in a 2d trap array, Advanced Quantum Technologies 3, 2000028 (2020).
  * Hughes  _et al._ [2020] A. C. Hughes, V. M. Schäfer, K. Thirumalai, D. P. Nadlinger, S. R. Woodrow, D. M. Lucas, and C. J. Ballance, Benchmarking a high-fidelity mixed-species entangling gate, Phys. Rev. Lett. 125, 080504 (2020).
  * Kim [2023] Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational quantum algorithms, Nature Reviews Physics 3, 625 (2021).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Noisy intermediate-scale quantum algorithms, Rev. Mod. Phys. 94, 015004 (2022).
  * Yuan  _et al._ [2019] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, Theory of variational quantum simulation, Quantum 3, 191 (2019).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5, 10.1038/ncomms5213 (2014).
  * Biamonte  _et al._ [2017] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Quantum machine learning, Nature 549, 195 (2017).
  * Cong  _et al._ [2019] I. Cong, S. Choi, and M. D. Lukin, Quantum convolutional neural networks, Nature Physics 15, 1273 (2019).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm (2014).
  * Wang  _et al._ [2018] Z. Wang, S. Hadfield, Z. Jiang, and E. G. Rieffel, Quantum approximate optimization algorithm for maxcut: A fermionic view, Phys. Rev. A 97, 022304 (2018).
  * Johnson  _et al._ [2017] P. D. Johnson, J. Romero, J. Olson, Y. Cao, and A. Aspuru-Guzik, Qvector: an algorithm for device-tailored quantum error correction (2017).
  * Xu  _et al._ [2021] X. Xu, S. C. Benjamin, and X. Yuan, Variational circuit compiler for quantum error correction, Phys. Rev. Appl. 15, 034068 (2021).
  * Khatri  _et al._ [2019] S. Khatri, R. LaRose, A. Poremba, L. Cincio, A. T. Sornborger, and P. J. Coles, Quantum-assisted quantum compiling, Quantum 3, 140 (2019).
  * Sharma  _et al._ [2020] K. Sharma, S. Khatri, M. Cerezo, and P. J. Coles, Noise resilience of variational quantum compiling, New Journal of Physics 22, 043006 (2020).
  * Jones and Benjamin [2022] T. Jones and S. C. Benjamin, Robust quantum compilation and circuit optimisation via energy minimisation, Quantum 6, 628 (2022).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Phys. Rev. X 7, 021050 (2017).
  * Cîrstoiu  _et al._ [2020] C. Cîrstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger, Variational fast forwarding for quantum simulation beyond the coherence time, npj Quantum Information 6, 10.1038/s41534-020-00302-0 (2020).
  * Yao  _et al._ [2021] Y.-X. Yao, N. Gomes, F. Zhang, C.-Z. Wang, K.-M. Ho, T. Iadecola, and P. P. Orth, Adaptive variational quantum dynamics simulations, PRX Quantum 2, 030307 (2021).
  * Lin  _et al._ [2021] S.-H. Lin, R. Dilip, A. G. Green, A. Smith, and F. Pollmann, Real- and imaginary-time evolution with compressed quantum circuits, PRX Quantum 2, 010342 (2021).
  * Barratt  _et al._ [2021] F. Barratt, J. Dborin, M. Bal, V. Stojevic, F. Pollmann, and A. G. Green, Parallel quantum simulation of large systems on small NISQ computers, npj Quantum Information 7, 79 (2021).
  * Barison  _et al._ [2021] S. Barison, F. Vicentini, and G. Carleo, An efficient quantum algorithm for the time evolution of parameterized circuits, Quantum 5, 512 (2021).
  * Berthusen  _et al._ [2022] N. F. Berthusen, T. V. Trevisan, T. Iadecola, and P. P. Orth, Quantum dynamics simulations beyond the coherence time on noisy intermediate-scale quantum hardware by variational trotter compression, Physical Review Research 4, 10.1103/physrevresearch.4.023097 (2022).
  * Barison  _et al._ [2022] S. Barison, F. Vicentini, I. Cirac, and G. Carleo, Variational dynamics as a ground-state problem on a quantum computer, Phys. Rev. Res. 4, 043161 (2022).
  * Miessen  _et al._ [2023] A. Miessen, P. J. Ollitrault, F. Tacchino, and I. Tavernelli, Quantum algorithms for quantum dynamics, Nature Computational Science 3, 25 (2023).
  * Trotter [1959] H. F. Trotter, On the product of semi-groups of operators, Proceedings of the American Mathematical Society 10, 545 (1959).
  * Suzuki [1976] M. Suzuki, Generalized trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems, Communications in Mathematical Physics 51, 183 (1976).
  * Abrams and Lloyd [1997] D. S. Abrams and S. Lloyd, Simulation of many-body fermi systems on a universal quantum computer, Phys. Rev. Lett. 79, 2586 (1997).
  * Ortiz  _et al._ [2001] G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Quantum algorithms for fermionic simulations, Phys. Rev. A 64, 022319 (2001).
  * Tacchino  _et al._ [2019] F. Tacchino, A. Chiesa, S. Carretta, and D. Gerace, Quantum computers as universal quantum simulators: State-of-the-art and perspectives, Advanced Quantum Technologies 3, 1900052 (2019).
  * Grimsley  _et al._ [2019] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nature Communications 10, 3007 (2019).
  * Tang  _et al._ [2021] H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, Qubit-ADAPT-VQE: An adaptive algorithm for constructing hardware-efficient ansätze on a quantum processor, PRX Quantum 2, 10.1103/prxquantum.2.020310 (2021).
  * Van Dyke  _et al._ [2022] J. S. Van Dyke, G. S. Barron, N. J. Mayhall, E. Barnes, and S. E. Economou, Scaling adaptive quantum simulation algorithms via operator pool tiling (2022).
  * Anastasiou  _et al._ [2022] P. G. Anastasiou, Y. Chen, N. J. Mayhall, E. Barnes, and S. E. Economou, Tetris-adapt-vqe: An adaptive algorithm that yields shallower, denser circuit ansätze (2022).
  * Gomes  _et al._ [2023] N. Gomes, D. B. Williams-Young, and W. A. de Jong, Computing the many-body green’s function with adaptive variational quantum dynamics (2023).
  * Anastasiou  _et al._ [2023] P. G. Anastasiou, N. J. Mayhall, E. Barnes, and S. E. Economou, How to really measure operator gradients in adapt-vqe (2023), arXiv:2306.03227 [quant-ph] .
  * Mari  _et al._ [2021] A. Mari, T. R. Bromley, and N. Killoran, Estimating the gradient and higher-order derivatives on quantum hardware, Physical Review A 103, 10.1103/physreva.103.012405 (2021).
  * Shkolnikov  _et al._ [2021] V. O. Shkolnikov, N. J. Mayhall, S. E. Economou, and E. Barnes, Avoiding symmetry roadblocks and minimizing the measurement overhead of adaptive variational quantum eigensolvers (2021), arXiv:2109.05340 [quant-ph] .
  * Yordanov  _et al._ [2021] Y. S. Yordanov, V. Armaos, C. H. W. Barnes, and D. R. M. Arvidsson-Shukur, Qubit-excitation-based adaptive variational quantum eigensolver, Communications Physics 4, 228 (2021).
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error mitigation for short-depth quantum circuits, Phys. Rev. Lett. 119, 180509 (2017).
  * tA v  _et al._ [2021] A. tA v _et al._ , Qiskit: An open-source framework for quantum computing (2021).
  * van den Berg  _et al._ [2023] E. van den Berg, Z. K. Minev, A. Kandala, and K. Temme, Probabilistic error cancellation with sparse Pauli–Lindblad models on noisy quantum processors, Nature Physics 10.1038/s41567-023-02042-2 (2023).
  * E.P. Jornan [1993] E. W. E.P. Jornan,  _The Collected Works of Eugene Paul Wigner_ (Springer, Berlin, Heidelberg, 1993).
  * Bravyi and Kitaev [2002] S. B. Bravyi and A. Y. Kitaev, Fermionic quantum computation, Annals of Physics 298, 210 (2002).
  * Verstraete and Cirac [2005] F. Verstraete and J. I. Cirac, Mapping local hamiltonians of fermions to local hamiltonians of spins, Journal of Statistical Mechanics: Theory and Experiment 2005, P09012 (2005).
  * Whitfield  _et al._ [2016] J. D. Whitfield, V. c. v. Havlíček, and M. Troyer, Local spin operators for fermion simulations, Phys. Rev. A 94, 030301 (2016).
  * Setia  _et al._ [2019] K. Setia, S. Bravyi, A. Mezzacapo, and J. D. Whitfield, Superfast encodings for fermionic quantum simulation, Phys. Rev. Res. 1, 033033 (2019).
  * Chen and Xu [2022] Y.-A. Chen and Y. Xu, Equivalence between fermion-to-qubit mappings in two spatial dimensions (2022).
  * Nys and Carleo [2023] J. Nys and G. Carleo, Quantum circuits for solving local fermion-to-qubit mappings, Quantum 7, 930 (2023).
  * Linteau [2022] D. Linteau, adaptive-pvqd, https://github.com/dalin27/adaptive-pvqd (2022).
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, Qutip: An open-source python framework for the dynamics of open quantum systems, Comp. Phys. Comm. 183, 10.1016/j.cpc.2012.02.021 (2012).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2014).