# Exponential optimization of quantum state preparation via adiabatic
thermalization

Davide Cugini  Dipartimento di Fisica, Università di Pavia, via Bassi 6,
27100, Pavia, Italy  Davide Nigro  Dipartimento di Fisica, Università di
Pavia, via Bassi 6, 27100, Pavia, Italy  Mattia Bruno  Dipartimento di Fisica,
Università di Milano–Bicocca, and INFN, sezione di Milano–Bicocca, Piazza
della Scienza 3, I-20126 Milano, Italy  Dario Gerace  Dipartimento di Fisica,
Università di Pavia, via Bassi 6, 27100, Pavia, Italy

###### Abstract

The preparation of a given quantum state on a quantum computing register is a
typically demanding operation, requiring a number of elementary gates that
scales exponentially with the size of the problem. Using the adiabatic theorem
for state preparation, whose error decreases exponentially as a function of
the thermalization time, we report an explicit analytic expression for the
dependence of the characteristic time on the Hamiltonian used in the adiabatic
evolution. Exploiting this knowledge, we then design a preconditioning term
that modifies the adiabatic preparation, thus reducing its characteristic time
and hence giving an exponential advantage in state preparation. We prove the
efficiency of our method with extensive numerical experiments on prototypical
spin-models, which gives a promising strategy to perform quantum simulations
of manybody models via Trotter evolution on near-term quantum processors.

Introduction.– We are currently living the so-called Noisy Intermediate-Scale
Quantum (NISQ) era [1], in which actual devices are made of a limited number
of noisy qubits available in a number of different technological platforms.
Restricting to the realm of digital quantum processors satisfying the
DiVincenzo criteria [2], typically 30 all-to-all connected qubits are
available in trapped ions quantum processors [3] and more than 100 in
superconducting quantum circuits with next-neighbor coupling [4]. The physics
of manybody quantum systems may thus be efficiently simulated through the
direct application of such devices [5, 6]. A proven advantage in quantum
simulations is yet to be presented, at time of writing. In fact, in principle
one expects NISQ hardware to overcome known issues in simulating quantum
systems on classical devices, such as the sign problem for fermionic systems
[7, 8]. Nowadays, a substantial limitation resides in the severely restricted
hardware resources; nevertheless first proof-of-concept quantum simulations of
chemical properties have been reported [9, 10, 11], as well as the time-
dynamics of manybody Hamiltonians [12, 13, 14] and quantum field theories [15,
16, 17, 18, 19].

In quantum simulations on NISQ devices, one of the most pressing targets is
the efficient generation of a desired quantum state of a given model
Hamiltonian on a quantum register, which is essential, e.g., to directly
compute expectation values of the observables under study. This is the focus
of the present work. We define $H_{1}$ as the target Hamiltonian, properly
mapped on a quantum computer, and denote with $\\{|\Omega_{n}(1)\rangle\\}$
the set of its eigenstates. Using a sequence of single- and two-qubit gates
[20], a unitary operator $U_{n}$ acting on a register of $L$ qubits such that
$|\Omega_{n}(1)\rangle=U_{n}|0\rangle^{\otimes L}$ may always be found, and
the number of gates required to exactly implement $U_{n}$ scales exponentially
with $L$, i.e., with the size of the system. Currently, the noisy character of
state-of-the art qubits and quantum gates available in NISQ devices limits the
number of operations that can be reliably executed on actual hardware. Hence,
significant effort has been devoted to propose efficient algorithms allowing
to obtain a well approximated $U_{n}$ within limited quantum computing
resources. For state preparation, variational quantum eigensolvers have been
proposed as a valid solution, consisting of hybrid quantum/classical
algorithms in which quantum circuits with parametrized operations are
classically optimized to achieve the desired state [21, 22, 23, 24]. Despite
their usefulness for simple models, a general strategy to define such circuits
for arbitrary systems is still not available; moreover, they may involve
demanding classical optimizations in high-dimensional spaces [25]. An
alternative approach has been recently proposed [26], based on imaginary time
evolution, which might prove more effective than variational eigensolvers on
NISQ devices, despite its intrinsic limitations in the preparation of quantum
states with long correlation lengths.

Here we take a different route, by applying the Adiabatic Thermalization (AT)
[27] to show that any eigenstate of a target Hamiltonian can be generated, in
principle, with high fidelity and limited quantum computing resources. The key
idea is to start from a Hamiltonian $H_{0}$, different from $H_{1}$, whose
eigenstates are analytically known and trivial to prepare on a quantum
register. Then, its $n$-th eigenstate $|\Omega_{n}(0)\rangle$ is evolved in
time from $t=0$ to $t=\tau$, according to the auxiliary Hamiltonian

$\displaystyle\small{H\left(s\right)=\left(1-f\left(s\right)\right)H_{0}+f\left(s\right)H_{1}}\,,\quad
s=t/\tau\,,$ (1)

which is parametrized in terms of a dimensionless variable $s\in[0,1]$.
Essentially, $H\left(s\right)$ interpolates between $H_{0}$ and $H_{1}$
through the time-dependent holomorphic function $f(s)$ such that $f(0)=0$ and
$f(1)=1$. In what follows we denote the time evolution operator associated
with $H(s)$ in the time interval $[0,\tau]$ as $U_{\tau}$. Its application to
the n-th eigenstate of $H_{0}$ produces the final state
$U_{\tau}\ket{\Omega_{n}(0)}$, which is an approximation of the actual
eigenstate $\ket{\Omega_{n}(1)}$ of $H_{1}$, with a discrepancy quantified by
the infidelity

$\mathcal{I}(\tau)\equiv
1-\absolutevalue{\bra{\Omega_{n}(1)}U_{\tau}\ket{\Omega_{n}(0)}}^{2}\,.$ (2)

The present work is based on the observation that an appropriate choice of
$H_{0}$, or equivalently $\Delta=H_{1}-H_{0}$, may lead to computational speed
up. Starting from known analytical results on the exponential convergence of
this procedure [29], we hereby focus on the analytic evaluation of the error
in quantum state preparation as a function of the thermalization time, which
we then use to formulate a novel strategy to select $H_{0}$ such that the
residual thermalization effect is further reduced when computing the
approximate quantum state. To support our formal result, we show the
effectiveness of this approach by explicitly analysing the Heisenberg model in
either one or two dimensions, and demonstrating by numerical simulations the
suppression of residual thermalization effects at a fixed quantum circuit
depth. Given the known limitations of the AT to prepare singly targeted
eigenstates [28], we hereby extend our protocol also to the preparation of
eigenspaces spanned by more than one eigenstate. These results are
considerably relevant for prospective quantum simulations of manybody models
on NISQ hardware with unprecedented system size, possibly towards quantum
advantage.

Exponential scaling.– Let $\sigma_{0}(s)$ be a subset of the full spectrum
$\sigma(s)$ of $H(s)$, and $\mathcal{H}$ the Hilbert space of the system
configurations. We address the problem of preparing a state belonging to the
subspace $V(s)\subset\mathcal{H}$ spanned by those eigenstates with
eigenvalues in $\sigma_{0}(s)$. Let $P(s)$ be the projector on $V(s)$, our aim
is then to efficiently evolve the subspace $P(0)\mathcal{H}$ into
$P(1)\mathcal{H}$ through AT. In what follows, the thermalization error is
quantified as the operator norm

$\epsilon_{\mathrm{AT}}(\tau)\equiv\norm{\left(\mathds{1}-P(1)\right)U_{\tau}P(0)}\,,$
(3)

which is the generalization of Eq. (2) to the case of
$\mathrm{dim}\left(V(s)\right)>1$, where the norm of an operator $A$ acting on
the Hilbert space $\mathcal{H}$ is defined as

$\norm{A}\equiv\mathrm{sup}_{\ket{\Psi}\in\mathcal{H}}\frac{\norm{A\ket{\Psi}}}{\norm{\ket{\Psi}}}\,.$
(4)

Our starting point is the result presented in Ref. [29], which guarantees the
exponential falloff of $\epsilon_{\mathrm{AT}}$ as a function of $\tau$,
provided that $H(s)$ is holomorphic, the gap $d(s)$ between $\sigma_{0}(s)$
and $\sigma_{1}(s)=\sigma(s)\char 92\relax\sigma_{0}(s)$ is positive, and the
bandwidth $D(s)$ of $\sigma_{0}(s)$ is finite for any $s\in[0,1]$, as
schematically represented in Fig. 1. The existence of a finite and positive
$C,g(s)\in\mathbb{R}$ such that, for $\tau$ large enough,

$\begin{split}\epsilon_{\mathrm{AT}}(\tau)&\leq
C\int_{0}^{1}ds\,\mathrm{exp}\left(-\left\lfloor\tau/g(s)\right\rfloor\right)\,,\end{split}$
(5)

in which $\left\lfloor\cdot\right\rfloor$ is the floor operation, has been
proved in Ref. [29].

Figure 1: Schematic illustration of the adiabatic evolution of the spectral
features related to a time-dependent Hamiltonian $H(s)$, in which $\sigma(s)$
is the full spectrum, partitioned into the subsets of eigenvalues
$\sigma_{0}(s)$ and $\sigma_{1}(s)$. Here, $D(s)$ denotes the bandwidth of
$\sigma_{0}(s)$, while $d(s)=\mathrm{min}\left[d_{+}(s),d_{-}(s)\right]$
defines the minimal gap between the two subsets of eigenvalues.

For finite-dimensional Hilbert spaces, the parameter $D$ is inherently finite,
which guarantees the exponential falloff of the error just by the conditions
that $H(s)$ is holomorphic and $d(s)>0$. By defining the characteristic decay
time as

$\tilde{g}\geq\mathrm{sup}_{s\in[0,1]}g(s)\,,$ (6)

the exponential bound

$\epsilon_{\mathrm{AT}}\leq
C\,\mathrm{exp}\left(-\left\lfloor\tau/\tilde{g}\right\rfloor\right)$ (7)

is straightforwardly obtained. Clearly, any reduction of $\tilde{g}$
exponentially reduces the error bound in Eq. (7). The purpose of this work is
then to find an analytic expression for $g(s)$ that will ultimately allow us
to define an efficient optimization strategy. In particular, we hereby report
our extension of the result of Ref. [29] which relates $g(s)$ to the metrics
$\\{\norm{\Delta},D(s),d(s)\\}$, therefore providing a concrete way to control
$\tilde{g}$. For the purposes of the present manuscript, it is relevant to
notice that these metrics can be efficiently estimated by classical means,
e.g., via the Lanczos Algorithm [31, 32]. Given the functions

$W(d,D)=\left[2^{8}\left(1+2\frac{D}{\pi
d}\right)^{2}\left(\frac{\pi^{2}}{3}+\frac{4}{d}\right)\right]^{7/3}\,,$ (8)

and

$\varrho\left(s,\tfrac{d}{\norm{\Delta}}\right)=\mathrm{sup}_{\rho\in\mathbf{C}}\left[|\rho|:\left|f(s+\rho)-f(s)\right|<\tfrac{d}{4\norm{\Delta}}\right]$
(9)

in the case of $H(s)$ defined according to Eq. (1), we obtain that the ratio

$g(s)=4\frac{W\left(d(s),D(s)\right)}{\varrho\left(s,\frac{d(s)}{\norm{\Delta}}\right)}$
(10)

satisfies Eq. (5) (our proofs will be reported elsewhere). Since $g(s)$ is a
non-increasing function of $d(s)$ and $1/D(s)$, we fix $\tilde{g}$ to

$\tilde{g}=4\frac{W(d_{\mathrm{min}},D_{\mathrm{max}})}{\mathrm{inf}_{s\in[0,1]}\varrho(s,\frac{d_{\mathrm{min}}}{\norm{\Delta}})}\,,$
(11)

where

$d_{\mathrm{min}}=\mathrm{inf}_{s\in[0,1]}d(s)\,,\quad
D_{\mathrm{max}}=\mathrm{sup}_{s\in[0,1]}D(s)\,.$ (12)

For the commonly employed linear holomorphic function $f(s)=s$, the expression
above reduces to

$\tilde{g}=16\frac{\norm{\Delta}}{d_{\mathrm{min}}}W(d_{\mathrm{min}},D_{\mathrm{max}}))\,.$
(13)

Given a target Hamoltonian $H_{1}$, the naive choice for $H_{0}$ is the
diagonal part of $H_{1}$, such that its eigenstates are manifest. In what
follows we present a universal method to select a more convenient auxiliary
Hamiltonian, $H_{0}$, in the context of digital quantum simulations, where
$U_{\tau}$ is constructed using the Trotter-Suzuki approximation [33]. Our
proposal may not lead to the absolute best choice for specific problems, but
it guarantees an exponential advantage, in the quantum state preparation of a
target Hamoltonian $H_{1}$, over the naively picked $H_{0}$. The most generic
starting Hamiltonian can be expressed in terms of its matrix elements with
respect to the computational basis as

$\bra{i}H_{0}\ket{j}\equiv\delta_{ij}\bra{i}H_{1}\ket{j}+\bra{i}M\ket{j}\,,$
(14)

in which $M$ denotes a parameterized Hermitian operator (called the
preconditioner) to be optimized in order to minimize $\tilde{g}$. We impose
three conditions on $M$: (i) $H_{0}$ should be diagonal, such that its
eigenstates are the elements of the computational basis and it is then trivial
to initialize them on a quantum computing register; (ii) the additional
computational cost for the implementation of
$\mathcal{U}_{M}(\theta)=\mathrm{exp}\left(-i\theta M\right)$ at each Trotter
step (where $\theta\in\mathbf{R}$) should not frustrate the cost reduction
stemming from the minimization of $\tilde{g}$; finally, (iii) the number of
real parameters needed to describe $M$ should not increase exponentially with
the number of qubits, $L$. The first condition is satisfied by any diagonal
operator. We hereby propose the extremely “cheap” form

$\bra{i}M\ket{i}=\Sigma_{j=0}^{L-1}\alpha_{j}\,\mathrm{bin}_{j}\left(i\right)\,,$
(15)

in which $\ket{i}$ is the i-th element of the computational basis, and
$\mathrm{bin}_{j}\left(i\right)$ is the j-th digit of $i$ once expressed in
binary notation. The key idea behind the choice of Eq. (15) will be reported
elsewhere. Most importantly, this choice allows to fulfill both conditions
(ii) and (iii). In fact, the number of free parameters is $L$ in this case,
thus increasing linearly with the system size. This result is essential to
reduce the quantum computational cost of this procedure. In particular, the
corresponding time evolution operator is

$\mathcal{U}_{M}(\theta)=\bigotimes_{j=0}^{L-1}\mathrm{P}_{j}(\alpha_{j}\theta)\,,$
(16)

in which $\mathrm{P}_{i}$ is the single-qubit phase gate, applied to the i-th
qubit. Notice that the quantum circuit depth of $\mathcal{U}_{M}(\theta)$
corresponds to a single qubit gate, independently from the size of the system.

Figure 2: Ground state preparation for the 9 sites 1D and 2D Heisenberg models
with periodic boundary conditions, schematically illustrated in panels (a) and
(b), respectively, where the error dependence is plotted as a function of
thermalization time, $\tau$, for $J^{x}=5J^{z}$.

Essentially, we propose to interpolate between $H_{0}$ (Eqs. (14)-(15)) and
the targeted $H_{1}$ by AT, after classically minimizing $\tilde{g}$ through a
scanning of the free parameters $\\{\alpha_{j}\\}$ in $H_{0}$. This simple
protocol will guarantee an exponentially faster quantum state preparation.

Figure 3: Characteristic times $g_{\mathrm{fit}}$ obtained from fitting the
numerical curves $\epsilon_{\mathrm{AT}}(\tau)$ (as the ones in Fig. 2) with a
single exponential (see text), for different values of the ratio $J^{x}/J^{z}$
of the 1D Heisenberg model. Standard deviations of the fit are plotted as
error bars.

Numerical results.– Here we apply the AT to the preparation of the ground
state for the Heisenberg model in either one- or two-dimensional (1D, 2D)
square lattices with periodic boundary conditions. The purpose is to
numerically confirm the previous theoretical result in case models of
practical relevance for quantum simulations. The corresponding quantum
circuits are implemented and executed on the Qiskit qasm simulator [34].

To focus only on $\epsilon_{\mathrm{AT}}$, we use sufficiently short Trotter
steps, to suppress the systematic error coming from the Trotter-Suzuki
approximation. In the following, the state preparation is performed by
applying three different approaches, all with the same adiabatic function
$f(s)=s$. The first is the standard AT, where $M$ is set to zero in Eq. (14).
Since the calculation of $\norm{\Delta}$ is computationally cheap and
independent from the chosen target subspace $V(s)$, we consider a second
approach where the parameters $\alpha_{j}$ are chosen to minimize
$\norm{\Delta}$ in Eq. (13). Notice that this does not guarantee that
$\tilde{g}$ is minimized. In fact, we finally consider a third approach
consisting in the full classical evaluation of $\tilde{g}$ in Eq. (13)
involving the estimation $d_{\mathrm{min}}$ and $D_{\mathrm{max}}$. Despite
yielding a significantly better scaling for our case studies, for more
complicated models it may be less practical, since this approach requires the
knowledge of the spectrum at several intermediate values of $s$.

The XZ Heisenberg model is a prototypical manybody model describing systems of
spin-1/2 particles on a lattice, defined by the interacting Hamiltonian

$H_{1}=-\frac{1}{2}\Sigma_{\langle
i,j\rangle}\left(J^{z}Z_{i}Z_{j}+J^{x}X_{i}X_{j}\right),$ (17)

where the sum runs over next-neighboring sites, while $X_{i}$ and $Z_{i}$ are
the Pauli matrices acting on the i-th qubit, conventionally defined as

$X=\begin{bmatrix}0&1\\\ 1&0\\\ \end{bmatrix}\,,\quad Y=\begin{bmatrix}0&-i\\\
i&0\\\ \end{bmatrix}\,,\quad Z=\begin{bmatrix}1&0\\\ 0&-1\\\ \end{bmatrix}\,.$
(18)

In this work we consider real and positive coupling constants, $J^{z}$ and
$J^{x}$. The AT process is then described by the operator

$H(s)=-\frac{1}{2}\Sigma_{\langle
i,j\rangle}\left(J^{z}Z_{i}Z_{j}+sJ^{x}X_{i}X_{j}\right)+(1-s)M\,.$ (19)

At $s=0$ the Hamiltonian is diagonal and, in the unpreconditioned case
($M=0$), has a degenerate ground state. The corresponding eigenstates,
$\ket{0}^{\otimes L}$ and $\ket{1}^{\otimes L}$, are orthonormal and invariant
under translations. It is possible to preserve the translational invariance of
the AT process by imposing the condition $\alpha_{i}=\alpha$ ($\forall$ $i$)
in Eq. (16). This choice, albeit not unique, allows to reduce the number of
parameters to a single one. The effect of the preconditioner ($M$) on the
physical system is equivalent to the introduction of an external uniform
magnetic field along the $z$ axis, with field intensity $\alpha$. The results
for the ground state preparation of this model are reported in Fig. 2, for
either the 1D or the 2D models. The outcome is quite striking: independently
on the dimensionality of the model or the parameters regime, the AT provides a
way of determining the ground state of the XZ Heisenberg model with an error
that scales almost exponentially with the thermalization time. In addition, by
comparing the three different approaches above, we notice a significant
reduction in the required thermalization time to achieve a given error, with
the third (most optimal) approach allowing to obtain the shortest $\tau$. The
performance differences between the different approaches become more and more
evident on increasing $J^{x}/J^{z}$, as shown in Fig. 3. Here, the
characteristic time obtained by fitting
“$\ln{\\{\epsilon_{\mathrm{AT}}(\tau)\\}}$” with the function
“$C_{\mathrm{fit}}-\tau/g_{\mathrm{fit}}$” is reported as a function of
$J^{x}/J^{z}$. In particular, the $g_{\mathrm{fit}}$ values obtained with the
three approaches are nearly indistinguishable when $J^{x}\ll J^{z}$, while
they exhibit a clear advantage after applying our preconditioning strategy for
$J^{x}\gg J^{z}$, i.e. when the target Hamiltonian is significantly different
from the starting one. The peak in the plot occurring around
$J^{x}/J^{z}\simeq 1$ corresponds to the critical point where, for an infinite
number of spins, the system would exhibit a phase transition. In that region
the preparation of the ground state is still a rather expensive task and our
proposal is less effective. Similar results have been obtained for the Ising
model in a transverse magnetic field, in both one and two dimensions, and will
be reported elsewhere.

Conclusions.– We have applied the theory of Adiabatic Thermalization (AT) to
the framework of quantum simulations on NISQ devices, specifically focusing on
the efficient preparation of a given quantum state. Starting from the seminal
result of Nenciu [29], we have reported an exponential boundary for the AT
error (analytical proof reported elsewhere), with an explicit expression for
its characteristic time in terms of quantities that can be classically
estimated. Using this improved analytic control, we propose a universal
preconditioning term, denoted as $M$ and tailored for NISQ computing, which
can be effectively optimized by using the aforementioned formulation. Upon
incorporation into the Hamiltonian model to be simulated on quantum hardware,
$M$ allows to exponentially suppress the AT error. Finally, we have tested our
protocol by numerical simulations performed for the representative Heisenberg
and Ising models (only the former were explicitly shown), in either one or two
dimensions. Our results demonstrate an exponential improvement, particularly
notable under more strenuous computational scenarios. Despite performing
numerical tests on simple spin models, we expect similar behaviours to occur
for more complicated systems. We expect that our work will significantly
enhance the efficiency of quantum state preparation, thereby boosting the
prospective achievement of quantum utility in state-of-art digital quantum
devices.

Acknowledgments. This research was supported by the Italian Ministry of
Research (MUR) through NRR project CN00000013 - National Research Center on
“HPC, Big Data and Quantum Computing” and PNRR project PE0000023 - National
Quantum Science Technology Institute (NQSTI). The research of M.B. for this
work was partially funded through the MUR program for young researchers “Rita
Levi Montalcini”. The authors warmly acknowledge G. Guarnieri, F. Scala, and
F. Tacchino for several scientific discussions.

## References

  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * DiVincenzo [2000] D. P. DiVincenzo, The physical implementation of quantum computation, Fort. der Phys. 48, 771 (2000).
  * Moses  _et al._ [2023] S. A. Moses _et al._ , A race-track trapped-ion quantum processor, Phys. Rev. X 13, 041052 (2023).
  * Kim  _et al._ [2023] Y. Kim _et al._ , Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * Tacchino  _et al._ [2020] 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 (2020).
  * Daley  _et al._ [2022] A. J. Daley _et al._ , Practical quantum advantage in quantum simulation, Nature 607, 667 (2022).
  * 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).
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Computational complexity and fundamental limitations to fermionic quantum monte carlo simulations, Physical review letters 94, 170201 (2005).
  * Kandala [2017] A. a. Kandala, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242 (2017).
  * Kandala  _et al._ [2019] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Error mitigation extends the computational reach of a noisy quantum processor, Nature 567, 491 (2019).
  * Motta and Rice [2022] M. Motta and J. E. Rice, Emerging quantum computing algorithms for quantum chemistry, WIREs Computational Molecular Science 12, e1580 (2022).
  * Chiesa  _et al._ [2019] A. Chiesa, F. Tacchino, M. Grossi, P. Santini, I. Tavernelli, D. Gerace, and S. Carretta, Quantum hardware simulating four-dimensional inelastic neutron scattering, Nature Physics 15, 455 (2019).
  * Arute  _et al._ [2020] F. Arute _et al._ , Observation of separated dynamics of charge and spin in the fermi-hubbard model (2020), arXiv:2010.07965 [quant-ph] .
  * Klymko  _et al._ [2022] K. Klymko _et al._ , Real-time evolution for ultracompact hamiltonian eigenstates on quantum hardware, PRX Quantum 3, 020323 (2022).
  * Martinez  _et al._ [2016] E. A. Martinez _et al._ , Real-time dynamics of lattice gauge theories with a few-qubit quantum computer, Nature 534, 516 (2016).
  * Klco  _et al._ [2018] N. Klco _et al._ , Quantum-classical computation of schwinger model dynamics using quantum computers, Physical Review A 98, 032331 (2018).
  * de Jong  _et al._ [2022] W. A. de Jong _et al._ , Quantum simulation of nonequilibrium dynamics and thermalization in the schwinger model, Physical Review D 106, 054508 (2022).
  * Chakraborty  _et al._ [2022] B. Chakraborty _et al._ , Classically emulated digital quantum simulation of the schwinger model with a topological term via adiabatic state preparation, Physical Review D 105, 094503 (2022).
  * Nguyen  _et al._ [2022] N. H. Nguyen _et al._ , Digital quantum simulation of the schwinger model and symmetry protection with trapped ions, PRX Quantum 3, 020324 (2022).
  * Nielsen and Chuang [2000] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, 2000).
  * Peruzzo  _et al._ [2014] A. Peruzzo _et al._ , A variational eigenvalue solver on a photonic quantum processor, Nature communications 5, 4213 (2014).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18, 023023 (2016).
  * Moll  _et al._ [2018] N. Moll _et al._ , Quantum optimization using variational algorithms on near-term quantum devices, Quantum Science and Technology 3, 030503 (2018).
  * Anand and oders [2022] A. Anand and oders, A quantum computing view on unitary coupled cluster theory, Chemical Society Reviews 51, 1659 (2022).
  * McClean  _et al._ [2018] J. R. McClean _et al._ , Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 4812 (2018).
  * Motta  _et al._ [2020] M. Motta _et al._ , Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nature Physics 16, 205 (2020).
  * Born and Fock [1928] M. Born and V. Fock, Beweis des adiabatensatzes, Zeitschrift für Physik 51, 165 (1928).
  * Jansen  _et al._ [2007] S. Jansen, M.-B. Ruskai, and R. Seiler, Bounds for the adiabatic approximation with applications to quantum computation, Journal of Mathematical Physics 48, 102111 (2007).
  * Nenciu [1993] G. Nenciu, Linear adiabatic theory. exponential estimates, Comm. Math. Phys. 152, 479 (1993).
  * Note [1] While this result was previously published [29], we report an alternative demonstration $\to$ proof in the Supplementary Information, for completeness.
  * Lanczos [1950] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, (1950).
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Trotter [1959] H. F. Trotter, On the product of semi-groups of operators, Proceedings of the American Mathematical Society 10, 545 (1959).
  * Qiskit contributors [2023] Qiskit contributors, Qiskit: An open-source framework for quantum computing (2023).


# Experimental error mitigation using linear rescaling for variational quantum
eigensolving with up to 20 qubits

Eliott Rosenberg<EMAIL_ADDRESS>Department of Physics, Cornell University,
Ithaca, NY 14853, USA  Paul Ginsparg  Department of Physics, Cornell
University, Ithaca, NY 14853, USA  Peter L. McMahon<EMAIL_ADDRESS>School of Applied and Engineering Physics, Cornell University, Ithaca, NY
14853, USA

###### Abstract

Quantum computers have the potential to help solve a range of physics and
chemistry problems, but noise in quantum hardware currently limits our ability
to obtain accurate results from the execution of quantum-simulation
algorithms. Various methods have been proposed to mitigate the impact of noise
on variational algorithms, including several that model the noise as damping
expectation values of observables. In this work, we benchmark various methods,
including a new method proposed here. We compare their performance in
estimating the ground-state energies of several instances of the 1D mixed-
field Ising model using the variational-quantum-eigensolver algorithm with up
to 20 qubits on two of IBM’s quantum computers. We find that several error-
mitigation techniques allow us to recover energies to within 10% of the true
values for circuits containing up to about 25 ansatz layers, where each layer
consists of CNOT gates between all neighboring qubits and Y-rotations on all
qubits.

##  I Introduction

Quantum computation promises to revolutionize physics and chemistry. Our
ability to make quantitative predictions about many-body and strongly coupled
quantum mechanical systems, such as atoms with many electrons, or systems of
quarks and gluons, has, with classical computers, been severely hampered by
the exponential memory requirements of storing and operating on many-body
wavefunctions. Quantum computers overcome this difficulty by representing the
wavefunction in the state of the qubits themselves. Noise, however, has
prevented current quantum computers from surpassing the capabilities of
classical computers in almost all applications, including quantum simulation.

In the long term, fault-tolerant quantum computation may be achieved using
error-correction techniques that employ many physical qubits to encode each
logical qubit. In the near term, there is much effort aimed at less costly
methods to mitigate or compensate for the effects of noise in current quantum
hardware. Some of these error-mitigation techniques involve artificially
increasing the noise and extrapolating the output to zero noise [1, 2, 3, 4,
5, 6]. Others involve representing the ideal circuit in terms of noisy ones
(the so-called “quasi-probability method”) [2, 3]. Still others involve using
replica copies of the system to effectively take powers of the density matrix
and suppress errors [7, 8], while others attempt to learn the effect of the
noise by training on classically simulable circuits [6, 9, 10, 11].

Several proposed error-mitigation methods have successfully modeled the noise
as damping the expectation values of observables towards their expectation
values in maximally mixed states, i.e., $\langle\hat{O}\rangle_{\rm
noisy}\approx C\langle\hat{O}\rangle_{\rm
exact}+(1-C)\mbox{Tr}[\hat{O}]/2^{n}$, where $0<C<1$, $\hat{O}$ is an
observable, and $n$ is the number of qubits [6, 9, 10, 12, 13, 14, 15, 16,
11]. In this work, we consider traceless observables, so the last term is not
present (although its inclusion would not present any significant
difficulties). Knowing the damping factor, $C$, the exact expectation value
can be recovered from the noisy one. Various methods for estimating the
damping factor are proposed in the works cited above. Here we propose another
method and compare the performance of these methods to each other and to zero
noise extrapolation when applied to the test problem of measuring the ground
state energy of a 1D mixed field Ising spin chain. We find that several
methods offer significant improvements over no mitigation, enabling circuits
of up to about 25 ansatz layers to yield expectation values within 10% of the
correct expectation value. Which method works best is problem-dependant, but
for our test problem, we find that estimating the damping factor from the
perturbative regime of the Hamiltonian works particularly well.

##  II Test problem

Our goal is to benchmark error mitigation methods, in order to get a sense for
when excessive circuit depth impedes obtaining meaningful results. To
benchmark these methods, we consider the problem of finding ground states of
the 20-qubit 1D mixed-field Ising model.

The mixed-field Ising model is defined by the following Hamiltonian:

$H=-J\sum_{i}Z_{i}Z_{i+1}-h_{x}\sum_{i}X_{i}-h_{z}\sum_{i}Z_{i},$ (1)

where we impose cyclic boundary conditions on the spins. If $h_{z}=0$, this
reduces to the transverse-field Ising model, which is exactly solvable using a
Jordan-Wigner transformation when the number of qubits is large. In this case,
a term arising from the cyclic boundary conditions can be neglected [17]. It
is, of course, also solvable in the limits $h^{2}_{x}/(J^{2}+h_{z}^{2})\ll 1$
and $J^{2}/h_{x}^{2}\ll 1$. For intermediate values of the field strengths,
however, it is non-integrable and quantum chaotic [18, 19]. Approximate
methods have been developed to study this model outside of its integrable
regime [20]. For the system size studied here (20 qubits), such approximate
methods are not necessary since the Hamiltonian can be exactly diagonalized on
a classical computer. For a discussion of the perturbative regimes of this
Hamiltonian, see Appendix A. Fig. A.1 compares the exact ground state energies
of this model to those computable in second-order perturbation theory from the
analytically solvable limits. Throughout this work, we set $J=1$, which sets
the overall scale of $H$.

We use an ansatz to parameterize approximate ground states of this
Hamiltonian. In particular, we consider a variant of the Alternating Layered
Ansatz (ALT) [21, 22], shown in Fig. C.2 and described in Appendix C. We
consider both the full ALT ansatz and a constrained version of it, with even
cyclic permutation symmetry imposed. The permutation symmetry reduces the
number of parameters from $n(l+1)$ to $2(l+1)$, where $n$ is the number of
qubits and $l$ is the number of ansatz layers. Fig. C.3 shows that this ansatz
is able to capture the ground states of these Hamiltonians as the number of
layers is increased. In fact, we do not need to make $l$ very large to capture
the ground state, but we do so anyway to study error mitigation.

For future applications of the Variational Quantum Eigensolver (VQE), in which
the ansatz circuits are wide and deep enough that they are not classically
simulable, the optimization step will need to be performed using measurements
from the quantum computer. If the effect of the noise is indeed to damp
expectation values of observables, and the parameter dependence of the damping
factor does not noticeably shift the location of the minimum, then it is not
necessary to apply error mitigation during the optimization routine. In this
case, correcting the damping is more relevant after the optimization, when we
evaluate the optimized circuit. We demonstrate that our symmetrized ansatz
circuits can be optimized for the mixed-field Ising model on 20 qubits using
unmitigated measurements from quantum computers (Figs. 1 and 2). Having shown
this, and because our primary goal is to benchmark error-mitigation
techniques, we proceed to optimize the remaining circuits classically. In all
Figures other than 1 and 2, the results were obtained with circuits that were
optimized classically.

##  III Error mitigation techniques

We consider the following questions: Given optimized circuits that approximate
ground states of the Hamiltonian (1), how can we best estimate the damping
factor that characterizes the suppression of the expectation value of the
Hamiltonian when evaluated on a quantum computer? How deep can we make our
ansatz using this method? The remainder of this section contrasts techniques
for estimating the damping factor. We do not include the method of Ref. [15],
which requires simultaneous measurements of all of the qubits, because
readout-error mitigation on a large number of qubits is a significant
challenge beyond the scope of this work. We also do not include the motion-
reversal method of Ref. [11] because, for circuits of the depths considered
here, repeating them and their inverses several times leads to very small
fidelities, to which it is impossible to accurately fit an exponential using a
reasonable number of measurements. Further, this method relies on the
assumption that the circuit and its inverse lead to similar decays in
fidelity, which may not be true. We also do not include the quasi-probability
method because it requires gateset tomography. We benchmark the following
methods. “Extrapolation in circuit depth” is a new method proposed here,
whereas the others have been proposed in previous works. We also tested two
methods that did not perform as well, which are described in Appendix D.

###  III.1 Zero noise extrapolation

In our implementation of zero noise extrapolation (ZNE) [1, 2, 3, 4, 5, 6], we
artificially increase the noise by replacing CNOT gates with odd numbers of
CNOT gates. For example, if we want to increase the noise by a factor of 3, we
replace each CNOT gate by 3 CNOT gates. We can scale the noise by factors
other than odd integers by randomly picking CNOT gates to repeat an odd number
of times. We then use an exponential fit to extrapolate the measured energy to
zero noise. We use three points (one of which is the unmodified circuit) to do
this extrapolation.

###  III.2 Calibrating using the perturbative regime

Recently, the Google team observed that, at least for the specific system they
studied, the damping factor was independent of the coupling strength [12].
They proposed, therefore, that one can measure the damping factor in a
perturbative regime of the Hamiltonian, and then use the same damping factor
outside of the perturbative regime to correct the results from the quantum
computer.

This method is, of course, limited to Hamiltonians that have a perturbative
regime, which doesn’t include all Hamiltonians of scientific interest. For
example, instances of the Sachdev-Kitaev-Ye (SYK) Hamiltonian do not have a
perturbative regime, but simulations of the model may give insights into
holography and quantum gravity (e.g., [23, 24, 25, 26, 27]). Nevertheless,
many Hamiltonians of interest do have a perturbative regime, so this technique
may still prove useful.

The mixed-field Ising model has several different perturbative regimes,
discussed in Appendix A. For the purposes of this error-mitigation technique,
we consider only perturbation in small $h_{x}$; we fix $h_{z}=0.1$ and
consider $0<h_{x}\leq 0.5$ to be the perturbative regime and $h_{x}=1.5$ to be
outside of the perturbative regime. While we might obtain a better estimate of
the damping factor by averaging measurements over multiple perturbative
regimes, approaching the desired field strength from both sides, multiple
perturbative regimes is not a generic feature of Hamiltonians so we do not
employ this trick.

###  III.3 Extrapolation in circuit depth

In contrast to Ref. [12], but in agreement with the observation that
exponentials work well for ZNE [3], we expect that the damping factor should
decay exponentially in the number of ansatz layers and generally observe this
to be the case. Further, circuits with small depth may be classically
simulable because only qubits in the backwards light cone of the measured
qubits need to be included. For 1D ansätze with nearest-neighbor 2-qubit
gates, such as the ansatz considered here, the number of qubits needed for
exact classical simulation is $\mbox{min}(m+2l,n),$ where $l$ is the number of
ansatz layers, $m$ is the number of qubits on which the observable acts,
assumed to be adjacent, and $n$ is the total number of qubits. If 36 qubits
are classically simulable, then for the Hamiltonian considered here (Eq. 1),
with $m=2$, we can simulate up to 17 ansatz layers for arbitrarily large $n$.
More generally, one needs to consider the size of the backwards light cone as
a function of the number of ansatz layers to determine the maximum number of
layers that can be simulated classically. This method will not work as well or
at all for observables that act on many qubits in the same term, such as the
SYK Hamiltonian in either the Jordan-Wigner or Bravyi-Kitaev [28]
representation. To test this mitigation method for our system, we fit an
exponential to the damping factor for up to 15 ansatz layers and then use the
fit to predict the damping factor for up to 50 layers.

###  III.4 Calibrating by omitting 1-qubit gates

Another recent paper [14] proposed estimating the damping factor by
eliminating all of the 1-qubit gates in the circuit. The remaining 2-qubit
gates, which are all controlled-NOTs, leave the state in $|00\ldots 0\rangle$,
and the damping factor is estimated as the fidelity of measuring all zeros.
For our system, eliminating the single-qubit gates is achieved by setting
$\vec{\theta}=0$. In addition to estimating the damping factor using the
fidelity, as suggested by Ref. [14], we also compare the measured energy to
the energy in the ideal state to estimate the damping factor.

This method can be regarded as a special case of a method in which the ansatz
is evaluated at classically solvable points. For our ansatz, if we pick the
Y-rotations to be $\theta_{i}=n_{i}\pi/2$, for integers $n_{i}$, then the
circuit is within the Clifford group and hence classically simulable. Using
Clifford circuits to estimate the damping factor has been proposed previously
[9]. For generic choices of the $n_{i}$, however, the expectation value of the
Hamiltonian becomes vanishingly small — and hence impossible to measure with a
finite number of shots — as the number of ansatz layers becomes large. The
fidelity also cannot be measured without adding more gates to the circuit
because the final state is not generically a computational basis state, so
estimating the predicted damping factor from the fidelity is not possible.
Therefore, these circuits cannot be used to measure the damping factor using
the method described in the previous paragraph. We might instead consider
picking all of the $n_{i}$ to be even, so that the state remains a
computational basis state. In generic computational basis states, at large
$n$, the Hamiltonian will still have a small expectation value, although not
as small as in a generic state. Nevertheless, we could estimate the damping
factor by measuring the fidelity of the final state in such a circuit. In the
present work, we focus on the case $\vec{\theta}=0$, which produces the ground
state of the free ($h_{x}=0$) Hamiltonian, a state close to the states of
interest. We leave the study of other Clifford circuits to future work.

In Ref. [14], this method is combined with ZNE. We also benchmark the
combination of this method with ZNE. There are two ways of combining these
methods. One could apply ZNE to the original circuit and to the calibration
circuit (the one with the single-qubit gates removed) separately and then
divide the former result by the predicted damping factor from the latter
result (we call this method “ZNE first”), or we could divide the measured
energy from each noise-scaled target circuit by the estimated damping factor
from the corresponding noise-scaled calibration circuit and then extrapolate
the quotient to zero noise (we call this method “ZNE last”). We benchmark both
methods.

##  IV Experimental techniques

Independent of the error-mitigation technique under consideration, we always
do the following to reduce noise:

  1. 1.

We omit gates outside of the backwards light cone of the measured observable.
This technique has previously been called a light-cone filter and has been
found to be beneficial [13].

  2. 2.

In order to implement the previous technique to its fullest advantage, and
also to reduce the number of qubits measured simultaneously, we measure each
Hamiltonian term separately.

  3. 3.

The qubits in our ansatz are connected in a loop, as are the qubits on the
quantum computer. For a loop of $n$ qubits, there are $2n$ ways to map the
circuit qubits to the physical qubits that preserve this connectivity (cyclic
permutations and a reflection). We average over four of these mappings. Such
an average has been called qubit assignment averaging [12]. Instead of picking
these mappings randomly, we pick configurations for which the damping factor,
as predicted by multiplying fidelities, is the largest (closest to one). This
ensures that we avoid particularly noisy qubits or gates to the extent
possible while also averaging out systematic biases.

  4. 4.

Additionally, we implement randomized compiling [29, 30, 31, 14], which helps
convert coherent errors into depolarizing-channel errors. In this method,
every CNOT gate is randomly replaced by one of 16 dressed CNOT gates, to which
it is equivalent in the absence of noise. The 16 dressed gates are described
succinctly in Figure 1 of [14]. We apply this method when benchmarking error
mitigation methods (i.e. Figs. 3-6) but not when running the VQE experiments
(Figs. 1-2).

  5. 5.

Finally, we apply readout error mitigation, assuming that readout errors are
uncorrelated and that the error rates are as measured by IBM and reported in
the backend properties. For details, see Appendices E-F.

##  V Results and Discussion

First, we demonstrate that the training phase of VQE can converge even when we
don’t apply error-mitigation techniques to correct the damping. Then we
compare the performance of the various error-mitigation techniques at
estimating the damping factor for tuned circuits.

###  V.1 VQE using unmitigated quantum circuit evaluations

Figure 1: VQE for the 20-qubit Ising Hamiltonian, with energy evaluations
performed on ibmq_sydney, with $h_{x}=1.5$ and $h_{z}=0.1$. Unmitigated
energies are fed into the SPSA algorithm. We impose cyclic permutation
symmetry on the ansatz to reduce the number of parameters. Note that SPSA
involves computing $E(\vec{\theta}_{i}\pm c_{i}\Delta_{i})$, not
$E(\vec{\theta}_{i})$. In the upper panels, the dark-blue dots are the
unmitigated energies measured on the quantum computer. These are damped
compared to the dark-yellow dots, which show exact classical evaluations of
the same circuits, from which we compute the observed damping factor (blue
dots in lower panels). The purple horizontal line in each upper panel shows
the energy of the classically optimized ansatz circuit, and the magenta dashed
line shows the energy of the exact ground state. The horizontal black dashed
line indicates the first excited state energy. In all trials except the
4-ansatz-layer trial, the optimizer succeeds in finding a state whose energy
is below the first excited state energy (we would pick the point with the
lowest energy, not necessarily the last point).

Figure 2: Same as Fig. 1 except on ibmq_toronto. The discontinuities,
particularly apparent in the 4-layer experiment, coincide with daily
recalibrations of the device. For the 3- and 4-layer trials the optimizer
succeeds in obtaining states with energies below the first excited state
energy.

Consider the training phase of VQE. In order to improve the convergence of our
optimizer, we impose a cyclic permutation symmetry on the ansatz, so that each
row of $Y$-rotations is described by two angles instead of $n$. This reduces
the number of parameters from $n(l+1)$ to $2(l+1)$, where $n$ is the number of
qubits and $l$ is the number of ansatz layers. We begin at a random point in
parameter space and use the quantum computer to measure the energy. We feed
the measured energy into the Simultaneous Perturbation Stochastic
Approximation (SPSA) optimizer, running on a classical computer, which picks
the updated parameters [32]. The first 50 function evaluations are used to
calibrate an SPSA hyperparameter as in [33], following which the optimizer
iterates through the SPSA algorithm. At each iteration, SPSA uses two
evaluations of the energy, at $\vec{\theta}=\vec{\theta}_{i}\pm
c_{i}\Delta_{i}$, and we submit both of these evaluations together in the same
job for improved speed.

Figs. 1 and 2 show the results of this VQE optimization for a 20-qubit mixed-
field 1D Ising model on ibmq_sydney and ibmq_toronto, respectively, with 3–6
ansatz layers. We see that the optimizer is often able to converge to a state
close to the true ground state, that is, below the energy of the first excited
state. The lowest energy state achieved during our VQE runs has an energy of
$-33.54$. For comparison, the true ground state energy is $-33.90$, and the
first excited state energy is $-32.60$. We might worry that this convergence
would be spoiled by a parameter-dependent damping factor. Indeed, in some of
our runs, the damping factor can be seen to increase as the optimizer
converges to the ground state, indicating its parameter dependence.

Although the parameter dependence of the damping factor does not prevent
convergence in our system, we cannot guarantee that it never will. If the
damping factor shifts the locations of the minima, techniques to minimize this
effect will be necessary. These could include additional experimental
techniques to symmetrize over measurement outcomes, such as exchanging the
roles of “1” and “0”, or applying an error mitigation technique during the
optimization routine. Many of the error-mitigation techniques discussed here,
however, assume either that the damping factor is parameter independent or
that the circuit has been optimized. Of the methods that we benchmark, only
ZNE could be used during the optimization phase.

One might naturally guess that some of the parameter dependence of the damping
factor comes from asymmetric readout errors. We show in Appendix F that,
indeed, the asymmetric readout errors lead to a predictable additive shift of
expectation values, and hence a parameter-dependent damping factor (Eq. 7).
However, because the shift is constant, its effect is to shift the overall
energy by a constant amount, leaving the location in parameter space of the
minimum unchanged. For the mixed-field Ising Hamiltonian, this overall shift
is approximately $-n(h_{x}+h_{z})(e_{1}-e_{0})\approx-1.6$ for $n=20$,
$h_{x}=1.5$, $h_{z}=0.1$, and $(e_{1}-e_{0})\approx 0.05.$ This is a larger
relative shift when the energy is closer to zero. It should therefore have the
effect of making the damping factor closer to unity at negative energies
further from the minimum (closer to zero). We see instead the opposite effect,
so this is not the source of the parameter-dependent damping factor.

###  V.2 Comparison of mitigation techniques

Having established that, in at least some cases, error mitigation is not
needed to tune the ansatz parameters, we now compare the performance of the
error-mitigation techniques, applied to optimized circuits. For the methods
that work by estimating a damping factor, we plot the observed damping factor
on the same axes as those predicted by the various methods (Fig. 3).
Additionally, for all of the methods, we plot the relative error of the
reconstructed energy as a function of the number of ansatz layers (Fig. 4).
Finally, we show to how many layers each method allows us to extend our ansatz
circuits while successfully mitigating the energy to within 10% of its true
value (Figs. 5-6).

Figure 3: Measured and predicted damping factors for ground states of the
20-qubit mixed field Ising model with $h_{x}=1.5$ and $h_{z}=0.1$. The methods
used to predict the damping factor are described in Secs. III.2, and III.3,
and III.4, respectively. For the first method, we take
$h_{x}=\\{0.1,0.2,0.3,0.4,0.5\\}$ to represent the perturbative regime and
average the measured damping factors for those five circuits. ZNE is absent
from this plot because it does not predict a damping factor. When computing
energies, we apply readout error mitigation using the reported readout error
rates, as discussed in Appendix F. We see that the “from pert” method
correctly predicts the damping factor for circuits up to about 25 layers and
that all of the methods fail by 50 layers. Error bars indicate statistical
errors only (see Appendix B).

Figure 4: Relative error in the measured ground state energy of the 20-qubit
mixed field Ising model with $h_{x}=1.5$ and $h_{z}=0.1$, mitigated using each
of the methods studied in this paper.

##  VI Discussion and Conclusion

We have used the variational-quantum-eigensolver algorithm to solve a quantum
spin model with up to 20 qubits (spins), and have compared various error-
mitigation techniques that can be applied to the VQE output to correct the
measured ground-state energy. We see that, for studying ground states of the
mixed-field Ising model on 20 qubits on ibmq_toronto and ibmq_sydney, all of
these methods provide improvements over no mitigation. Zero noise
extrapolation fails for large circuit depths because the measured energy, even
with no added noise, is damped very close to zero, and it is difficult to
extrapolate these measurements to zero noise. For the run on ibmq_sydney, ZNE
first fails at 9 layers (Figs.5-6). Some of the fits used in our
implementation of ZNE are shown in Fig. 7 and illustrate the difficulty of
extrapolating to zero noise when the initial circuit is deep.

Using the Hamiltonian’s perturbative regime to estimate the damping factor
worked relatively well, perhaps best out of the methods that we tested. In
Fig. 3, we see that it correctly predicts the damping factor on both devices
most of the time even though it, too, seems to fail when the circuit extends
to 50 ansatz layers. The reader may have noticed that the error bars in Fig. 4
are larger for this method than for the others. This is because the variance
of the damping factor over the five Hamiltonians in the perturbative regime
contributes to the statistical uncertainty of the predicted damping factor
(see Fig. 8). This could be reduced by picking more points in the perturbative
regime. However, the fact that the spread of the measured damping factors is
larger than one would expect from their individual error bars indicates that
this method is unlikely to predict the correct damping factor exactly.
Further, this method has a significant overhead, requiring the optimization of
the ansatz circuits for each of the perturbative regime Hamiltonians in
addition to the Hamiltonian of interest. Its applicability is also limited
because, unlike ZNE, it can only be applied to optimized circuits and only to
systems that have a perturbative regime.

As seen in Fig. 3, the damping factor did decay approximately exponentially in
the number of ansatz layers. Therefore, the exponential fit to the points with
at most 15 layers worked reasonably well out to about 25 layers, which is
better than ZNE performed on ibmq_sydney but not as well as the “from pert”
method worked. Extrapolating in the number of ansatz layers also only works
for optimized circuits, but it does not require the Hamiltonian to have a
perturbative regime. It does, however, require the Hamiltonian and the ansatz
to be local in order for the shallow-depth ansatz circuits to be exactly
classically simulable.

The method of removing the single-qubit gates, like ZNE, does not require the
ansatz to be optimized and does not place any restrictions on the Hamiltonian
or the ansatz. However, it does assume that the damping factor is independent
of the ansatz parameters, which we have seen is not true. Consequently, this
method, whether combined with ZNE or not, does not perform as well as the
other methods (c.f. Fig. 3).

In summary, we have seen that error-mitigation methods that work by estimating
a damping factor can outperform zero noise extrapolation and extend feasible
circuit depths to about 25 ansatz layers. However, these methods are also more
specialized, best suited to specific types of Hamiltonians and ansätze and
poorly suited to the optimization phase, during which one could still apply
ZNE. The Hamiltonian studied in this work was particularly simple and did not
require error mitigation during the optimization phase. However, that will not
always be the case.

By carefully studying error mitigation techniques and the circuit depths that
they enable on current devices, we are progressing toward the day when quantum
comptuers will be a tool for learning new things about physics and chemistry
systems.

Figure 5: Effectiveness of error mitigation. “3” means that the reconstructed
energy is within 10% of the correct value. “2” means statistically
indeterminate at 1 standard deviation. “1” means statistically indeterminate
at 2 standard deviations. “0” means that the reconstructed energy deviated by
more than 10% from the true energy. In this figure, the methods involving
calibration circuits do not include separate readout error mitigation.

Figure 6: Same as Fig. 5, but applying readout error mitigation (Appendix F)
to all energy evaluations, including in calibration circuits.

(a) 1 ansatz layer, without readout error mitigation

(b) 1 ansatz layer, with readout error mitigation

(c) 6 ansatz layers, without readout error mitigation

(d) 6 ansatz layers, with readout error mitigation

(e) 30 ansatz layers, without readout error mitigation

(f) 40 ansatz layers, without readout error mitigation

Figure 7: Some of the exponential fits used for zero noise extrapolation on
ibmq_sydney. It is not clear that readout error mitigation, as we perform it,
is beneficial for this method. In fact, for 30 and 40 ansatz layers, some of
the readout-mitigated energies were negative, causing the exponential fits to
fail. In future work, we intend to to explore better ways of performing
readout error mitigation. It is also clear why this method fails at large
circuit depth. The measured energy is very small even before we artifically
increase the noise, and it is therefore difficult to extrapolate back to zero.

Figure 8: Some of the extrapolations used when estimating the damping factor
from the perturbative regime. The solid horizontal line indicates the
predicted damping factor (the mean of the damping factors for $h_{x}\leq
0.5$), and the dashed horizontal lines indicate the 1$\sigma$ uncertainty of
the mean. Readout error mitigation has already been applied in these plots.

## Data and code availability

The code used to perform the numerical simulations, the experimental runs on
IBM’s hardware, and the data analysis (including error mitigation) is
available at https://github.com/mcmahon-lab/error_mitigation_vqe. The data
used to produce the figures may be made available upon request.

## Acknowledgements

We thank Thomas Hartman for helpful discussions and Mandar Sohoni for
providing feedback on a draft of this paper. Additionally, we thank Eun-Ah Kim
and three anonymous referees for useful feedback. We also thank Sarah Kaiser
for finding a bug in an earlier (unpublished) version of the accompanying
code. We gratefully acknowledge financial support from US DOE grant DE-
SC0020397 and access to IBM quantum-computer time via the AFRL IBM Q Network
Hub, facilitated by Laura Wessing and Jon Maggiolino. The views expressed are
those of the authors, and do not reflect the official policy or position of
IBM or the IBM Quantum team. We used Qiskit [34] to interact with IBM’s
devices. PLM acknowledges membership in the CIFAR Quantum Information Science
Program as an Azrieli Global Scholar.

## References

  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Physical Review X 7, 021050 (2017).
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error mitigation for short-depth quantum circuits, Physical Review Letters 119, 180509 (2017).
  * Endo  _et al._ [2018] S. Endo, S. C. Benjamin, and Y. Li, Practical quantum error mitigation for near-future applications, Physical Review X 8, 031027 (2018).
  * Kandala  _et al._ [2019] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Error mitigation extends the computational reach of a noisy quantum processor, Nature 567, 491 (2019).
  * Otten and Gray [2019] M. Otten and S. K. Gray, Recovering noise-free quantum observables, Physical Review A 99, 012338 (2019).
  * He  _et al._ [2020] A. He, B. Nachman, W. A. de Jong, and C. W. Bauer, Zero-noise extrapolation for quantum-gate error mitigation with identity insertions, Physical Review A 102, 012426 (2020).
  * Koczor [2021] B. Koczor, Exponential error suppression for near-term quantum devices, Physical Review X 11, 10.1103/physrevx.11.031057 (2021).
  * Huggins  _et al._ [2021] W. J. Huggins, S. McArdle, T. E. O’Brien, J. Lee, N. C. Rubin, S. Boixo, K. B. Whaley, R. Babbush, and J. R. McClean, Virtual distillation for quantum error mitigation (2021), arXiv:2011.07064 [quant-ph] .
  * Czarnik  _et al._ [2021] P. Czarnik, A. Arrasmith, P. J. Coles, and L. Cincio, Error mitigation with Clifford quantum-circuit data (2021), arXiv:2005.10189 [quant-ph] .
  * Montanaro and Stanisic [2021] A. Montanaro and S. Stanisic, Error mitigation by training with fermionic linear optics (2021), arXiv:2102.02120 [quant-ph] .
  * Shaw [2021] A. Shaw, Classical-quantum noise mitigation for NISQ hardware (2021), arXiv:2105.08701 [quant-ph] .
  * Google AI Quantum and collaborators [2020] Google AI Quantum and collaborators, Observation of separated dynamics of charge and spin in the Fermi-Hubbard model (2020), arXiv:2010.07965 [quant-ph] .
  * Mi  _et al._ [2021] X. Mi, P. Roushan, C. Quintana, S. Mandrà, J. Marshall, C. Neill, F. Arute, K. Arya, J. Atalaya, R. Babbush, and et al., Information scrambling in quantum circuits, Science 10.1126/science.abg5029 (2021).
  * Urbanek  _et al._ [2021] M. Urbanek, B. Nachman, V. R. Pascuzzi, A. He, C. W. Bauer, and W. A. de Jong, Mitigating depolarizing noise on quantum computers with noise-estimation circuits (2021), arXiv:2103.08591 [quant-ph] .
  * Vovrosh  _et al._ [2021] J. Vovrosh, K. E. Khosla, S. Greenaway, C. Self, M. S. Kim, and J. Knolle, Simple mitigation of global depolarizing errors in quantum simulations, Physical Review E 104, 10.1103/physreve.104.035309 (2021).
  * Ville  _et al._ [2021] J.-L. Ville, A. Morvan, A. Hashim, R. K. Naik, B. Mitchell, J.-M. Kreikebaum, K. P. O’Brien, J. J. Wallman, I. Hincks, J. Emerson, E. Smith, E. Younis, C. Iancu, D. I. Santiago, and I. Siddiqi, Leveraging randomized compiling for the QITE algorithm (2021), arXiv:2104.08785 [quant-ph] .
  * Pfeuty [1970] P. Pfeuty, The one-dimensional Ising model with a transverse field, Annals of Physics 57, 79 (1970).
  * Kim and Huse [2013] H. Kim and D. A. Huse, Ballistic spreading of entanglement in a diffusive nonintegrable system, Physical Review Letters 111, 127205 (2013).
  * Kim  _et al._ [2014] H. Kim, T. N. Ikeda, and D. A. Huse, Testing whether all eigenstates obey the eigenstate thermalization hypothesis, Physical Review E 90, 052105 (2014).
  * Wurtz and Polkovnikov [2020] J. Wurtz and A. Polkovnikov, Emergent conservation laws and nonthermal states in the mixed-field Ising model, Physical Review B 101, 195138 (2020).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nature Communications 12, 10.1038/s41467-021-21728-w (2021).
  * Nakaji and Yamamoto [2021] K. Nakaji and N. Yamamoto, Expressibility of the alternating layered ansatz for quantum computation, Quantum 5, 434 (2021).
  * Babbush  _et al._ [2019] R. Babbush, D. W. Berry, and H. Neven, Quantum simulation of the Sachdev-Ye-Kitaev model by asymmetric qubitization, Physical Review A 99, 040301 (2019).
  * Martyn and Swingle [2019] J. Martyn and B. Swingle, Product spectrum ansatz and the simplicity of thermal states, Physical Review A 100, 032107 (2019).
  * Luo  _et al._ [2019] Z. Luo, Y.-Z. You, J. Li, C.-M. Jian, D. Lu, C. Xu, B. Zeng, and R. Laflamme, Quantum simulation of the non-fermi-liquid state of Sachdev-Ye-Kitaev model, npj Quantum Information 5, 53 (2019).
  * Brown  _et al._ [2021] A. R. Brown, H. Gharibyan, S. Leichenauer, H. W. Lin, S. Nezami, G. Salton, L. Susskind, B. Swingle, and M. Walter, Quantum gravity in the lab: Teleportation by size and traversable wormholes (2021), arXiv:1911.06314 [quant-ph] .
  * Nezami  _et al._ [2021] S. Nezami, H. W. Lin, A. R. Brown, H. Gharibyan, S. Leichenauer, G. Salton, L. Susskind, B. Swingle, and M. Walter, Quantum gravity in the lab: Teleportation by size and traversable wormholes, Part II (2021), arXiv:2102.01064 [quant-ph] .
  * Bravyi and Kitaev [2002] S. B. Bravyi and A. Y. Kitaev, Fermionic quantum computation, Annals of Physics 298, 210–226 (2002).
  * Wallman and Emerson [2016] J. J. Wallman and J. Emerson, Noise tailoring for scalable quantum computation via randomized compiling, Physical Review A 94, 10.1103/physreva.94.052325 (2016).
  * Cai and Benjamin [2019] Z. Cai and S. C. Benjamin, Constructing smaller Pauli twirling sets for arbitrary error channels, Scientific Reports 9, 10.1038/s41598-019-46722-7 (2019).
  * Cai  _et al._ [2020] Z. Cai, X. Xu, and S. C. Benjamin, Mitigating coherent noise using Pauli conjugation, npj Quantum Information 6, 10.1038/s41534-019-0233-0 (2020).
  * [32] We use Qiskit’s implementation of SPSA, which is documented at https://qiskit.org/documentation/stubs/qiskit.aqua.components.optimizers.SPSA.html.
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242–246 (2017).
  * Abraham  _et al._ [2019] H. Abraham _et al._ , Qiskit: An open-source framework for quantum computing (2019).
  * Asfaw  _et al._ [2020] A. Asfaw _et al._ ,  _Learn Quantum Computation Using Qiskit_ (2020).
  * [36] See https://quantum-computing.ibm.com/lab/docs/iql/manage/systems/properties.
  * Funcke  _et al._ [2020] L. Funcke, T. Hartung, K. Jansen, S. Kühn, P. Stornati, and X. Wang, Measurement error mitigation in quantum computers through classical bit-flip correction (2020), arXiv:2007.03663 [quant-ph] .

##  Appendix A Perturbation theory for the mixed-field Ising model

The mixed-field Ising Hamiltonian (Eq. 1) is analytically solvable in several
limits, as discussed in the main text. In the small-$h_{x}$ limit, it is
diagonal in the computational basis, and the ground state (assuming $J>0$ and
$h_{z}>0$) approaches $|00\ldots 00\rangle$. The ground state energy, to
second order in $h_{x}$ is

$E_{\rm gs}=-n\left(h_{z}+J+\frac{h_{x}^{2}}{2h_{z}+4J}+O(h_{x}^{3})\right).$
(2)

Further, in the large $h_{x}$ limit, the ground state approaches
$|++\ldots++\rangle$, where
$|+\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle).$ In this limit, the ground
state energy is

$E_{\rm
gs}=-n\left(h_{x}+\frac{2h_{z}^{2}+J^{2}}{4h_{x}}+O\left(h_{x}^{-2}\right)\right).$
(3)

Finally, in the small $h_{z}$ limit, the Hamiltonian approaches the
transverse-field Ising model, which, in the large-$n$ limit, is exactly
solvable via a Jordan-Wigner transformation [17]. Deriving an analytical
expression for the ground state energy in the small $h_{z}$ limit is beyond
the scope of this work. Instead, we numerically compute the energy to second-
order in perturbation theory in this limit (i.e. using the exact numerical
eigenvectors from $h_{z}=0$).

The perturbative results are compared to the exact ground state energies in
Fig. A.1.

Figure A.1: Exact ground state energy of the mixed-field Ising model (Eq. 1)
and predictions from second-order perturbation theory, for $h_{z}=0.1$.
Relevant to this work, we see that $h_{x}\leq 0.5$ is in the small $h_{x}$
perturbative regime, while $h_{x}=1.5$ is outside of it (at least to second
order). For the small-$h_{z}$ limit, we cheat slightly and use exact
diagonalization to compute the eigenstates and eigenenergies when $h_{z}=0$.
We thus include the cyclic boundary term that is neglected in analytic
treatments (cf., Eq. 2.4 in [17], and the following discussion), so the
small-$h_{z}$ line should be regarded as slightly better than one could do
analytically, but approaching the analytical solution as the number of spins
increases.

##  Appendix B A note about error bars

Throughout this paper, error bars indicate statistical errors only. There are,
of course, a variety of systematic errors as well. However, the present goal
is to correct these systematic errors, so it is more appropriate for error
bars to only show statistical errors, which result from having a finite number
of counts. We consider the variation among the qubit assignments used in the
qubit assignment averaging to be a systematic error, as the same choices of
qubit assignments are used in the calibration circuits as in the circuit that
we wish to correct. However, we consider the variation over the different
“perturbative” points in the method introduced in Sec. III.2 to be statistical
error, since a larger sample of such points can always be taken.

##  Appendix C Variational ansatz circuit

As described in the main text, our ansatz consists of CNOTs between adjacent
qubits and Y-rotations on all qubits (Fig. C.2). CNOT is a basis gate on IBM’s
quantum computers. However, the Y-rotations and the change of basis required
to measure $X_{i}$ must also be written in terms of basis gates. We use the
following basis-gate decompositions:

$\begin{split}R_{y}(\theta)&=R_{z}(\pi)\sqrt{X}R_{z}(\pi+\theta)\sqrt{X}\\\
HR_{y}(\theta)&=R_{z}(\pi)\sqrt{X}R_{z}(3\pi/2-\theta)\sqrt{X}R_{z}(-\pi).\end{split}$
(4)

ibmq_toronto and ibmq_sydney include 12- and 20-qubit loops, the topologies
studied here, so we have not needed to apply any gates between non-neighboring
qubits.

Figure C.2: The Alternating Layered Ansatz that we use throughout this work,
shown for 6 qubits and 3 layers. When we impose partial cyclic permutation
symmetry, we demand that, within a layer of $y$ rotations, all of the
rotations on even qubits have the same angle, and all of the rotations on odd
qubits have the same angle. Figure C.3: Ansatz performance in representing
ground states of the mixed-field Ising Hamiltonian on 20 qubits, with
$h_{x}=1.5$ and $h_{z}=0.1$. The symmetric ansatz is able to represent the
ground states nearly as well as the full ansatz with the same number of
layers. Further, relatively few layers are needed to capture the ground state.
However, we experiment with many more layers in this paper so as to benchmark
error mitigation techniques.

##  Appendix D Poorly performing methods

In addition to the methods discussed in the main text, we also benchmark two
methods that do not perform as well. In one of these methods, which we call
“multiplying fidelities,” we multiply $(1-e_{i})$, where $e_{i}$ is the
reported error rate for gate $i$, over all of the gates appearing in the
backwards light cone of the observed qubits. This is then taken as an estimate
of the overall damping factor. In Fig. D.4, it is clear that this method
predicts more suppression than actually exists. The other poorly performing
method is to use the noise model in Qiskit Aer, which includes local
depolarizing channels on all of the gates, thermal relaxation errors, and
readout errors. This method is not scalable since it requires simulating the
circuit, but we benchmark it for comparison. It consistently predicts less
suppression than exists.

Figure D.4: The measured damping factor (black “x”s) versus the number of
ansatz layers, for classically optimized circuits that are then run on
ibmq_sydney or ibmq_toronto. The circuits are optimized to the ground state of
the 12-qubit mixed-field Ising Hamiltonian with $h_{x}=1.5$ and $h_{z}=0.1$.
The blue line shows the damping factor that is predicted from multiplying the
fidelities of the gates in the backwards light cone of the measured qubits.
The orange line shows the prediction from a Qiskit Aer noise model with local
depolarizing errors, thermal relaxation errors, and readout errors. In both
cases, the error rates are obtained from the device’s calibration. For deep
circuits, multiplying the gate fidelities consistently predicts too much
suppression, whereas the Qisit Aer model with local gate errors predicts too
little suppression.

##  Appendix E Readout error mitigation

In the standard technique for readout error mitigation on $n$ measured qubits,
one prepares the qubits in each of the $2^{n}$ computational basis states and
measures the probabilities of reading out the different outcomes (e.g. [35]).
This results in a $2^{n}\times 2^{n}$ matrix, which is inverted to perform the
mitigation. When $n$ is large, this method becomes intractable and one must
make simplifying assumptions, such as ignoring correlations among the qubits
or among sets of qubits. For the Mixed-Field Ising Hamiltonian (Eq. 1), we
only need to measure at most two qubits simultaneously, so the standard
readout error mitigation method is easily tractable. Nevertheless, we assume
uncorrelated readout errors and use the methods developed in Appendix F to
perform the mitigation. We use the readout error rates measured by IBM during
their approximately daily calibrations [36]. We have found that this works
about as well as doing the full readout error mitigation using our own
calibrations.

##  Appendix F Uncorrelated readout errors

In this appendix, we show that the effect of uncorrelated readout errors is to
reduce the expectation value of a Pauli operator by the factor $(1-2e_{m})$,
for each qubit measured, where $e_{m}=(e_{0}+e_{1})/2$ is the average readout
error, and offset it by a small additive amount. Some similar results appear
in [37].

Suppose that we want to measure a Pauli string of length $N$ by applying
readout gates to all of the qubits. If the readout gates were noiseless, we
would measure outcome $q$ with probability $f(q)$. The noiseless expectation
value of the Pauli operator is then $\langle P\rangle=\sum_{q}P(q)f(q)$, where
$P(q)=\pm 1$ is the bit parity of $q$. Now, suppose that the measurement gates
add noise such that the probability of a 0 getting measured as a 1 is $e_{0}$
and the probability of a 1 getting measured as a 0 is $e_{1}$. Let
$\tilde{f}(q)$ be the probability distribution of the outcomes, including this
readout error. Let $\tilde{\langle P\rangle}=\sum_{q}P(q)\tilde{f}(q)$.

The main result that we will establish is that

$\tilde{\langle
P\rangle}=\sum_{q}P(q)f(q)(1-2e_{0})^{n_{0}(q)}(1-2e_{1})^{n_{1}(q)},$ (5)

where $n_{0}(q)$ is the number of bits in $q$ that are 0, and $n_{1}(q)$ is
the number of bits in $q$ that are 1. If the two error rates are equal
($e_{0}=e_{1}=e_{m}$), Eq. 5 becomes

$\langle\tilde{P}\rangle=\langle P\rangle(1-2e_{m})^{N},\hskip
10.0pte_{0}=e_{1}=e_{m}.$ (6)

However, we can also simplify Eq. 5 without assuming the error rates are
equal. For $N=1$, $\langle P\rangle=f(0)-f(1)$, so Eq. 5 becomes

$\langle\tilde{P}\rangle=\langle P\rangle(1-e_{0}-e_{1})+e_{1}-e_{0},\hskip
10.0ptN=1.$ (7)

For larger $N$, $\langle P\rangle$ does not determine $f(q)$. However, the
largest and smallest $\langle\tilde{P}\rangle$ for a given $\langle P\rangle$
will always result from letting $f(q)$ be nonzero for at most two choices of
$n_{1}(q)$, one with $n_{1}$ even and the other with $n_{1}$ odd. Define these
values as $n_{1}^{+}$ and $n_{1}^{-}$, respectively. Define
$n_{0}^{\pm}=N-n_{1}^{\pm}$ and $f_{\pm}$ to be the corresponding
probabilities. Then $\langle P\rangle=f_{+}-f_{-}$, and so Eq. 5 becomes

$\begin{split}\langle\tilde{P}\rangle&=\left(\frac{1+\langle
P\rangle}{2}\right)\left(1-2e_{0}\right)^{n_{0}^{+}}\left(1-2e_{1}\right)^{n_{1}^{+}}-\left(\frac{1-\langle
P\rangle}{2}\right)\left(1-2e_{0}\right)^{n_{0}^{-}}\left(1-2e_{1}\right)^{n_{1}^{-}}\end{split}$
(8)

Without loss of generality, we assume that $e_{1}\geq e_{0}$. We also assume
that $e_{1}<0.5$. Then $\langle\tilde{P}\rangle$ is maximized by picking
$n_{1}^{+}=0$ and $n_{1}^{-}$ to be either $N$ or $N-1$, whichever is odd.
Similarly, the smallest $\langle\tilde{P}\rangle$ is obtained by picking
$n_{1}^{+}=N$ or $N-1$, whichever is even, and $n_{1}^{-}=1$. For the case of
$N=2$, this becomes

$(1-2e_{1})\left[\langle
P\rangle(1-e_{1}-e_{0})-(e_{1}-e_{0})\right]\leq\langle\tilde{P}\rangle\leq(1-2e_{0})\left[(1-e_{0}-e_{1})\langle
P\rangle+e_{1}-e_{0}\right],\hskip 10.0ptN=2.$ (9)

The midpoint between these two bounds gives an estimate of
$\langle\tilde{P}\rangle$ when $N=2$:

$\langle\tilde{P}\rangle\approx\langle P\rangle(1-e_{0}-e_{1})^{2}\langle
P\rangle+(e_{1}-e_{0})^{2},\hskip 10.0ptN=2.$ (10)

Note that the $y$-intercept is suppressed compared to the $N=1$ case (Eq. 7).
Indeed, it becomes more suppressed as we go to higher $N$. Numerically, we
find that, for $N\gtrsim 3$, for randomly chosen $f(q)$,

$\langle\tilde{P}\rangle\approx(1-e_{0}-e_{1})^{N}\langle P\rangle,\hskip
10.0ptN\gtrsim 3$ (11)

where this becomes a better approximation as $N$ increases. We present
numerical evidence for this result in Figs. F.5–F.6.

###  F.1 Proof of Eq. 5

Consider two bit strings, both of length $N$, $q$ and $q^{\prime}$. Define the
following quantities:

$\begin{split}k_{0}(q,q^{\prime})&=\sum_{i=0}^{N-1}\mathbbm{1}_{q^{\prime}_{i}=0}\mathbbm{1}_{q_{i}=1},\hskip
50.0ptu_{0}(q,q^{\prime})=\sum_{i=0}^{N-1}\mathbbm{1}_{q^{\prime}_{i}=0}\mathbbm{1}_{q_{i}=0},\\\
k_{1}(q,q^{\prime})&=\sum_{i=0}^{N-1}\mathbbm{1}_{q^{\prime}_{i}=1}\mathbbm{1}_{q_{i}=0},\hskip
50.0ptu_{1}(q,q^{\prime})=\sum_{i=0}^{N-1}\mathbbm{1}_{q^{\prime}_{i}=1}\mathbbm{1}_{q_{i}=1}.\end{split}$
(12)

Then, upon adding readout error, the probability $\tilde{f}(q)$ of measuring
bit string $q$ is

$\tilde{f}(q)=f(q)(1-e_{0})^{n_{0}(q)}(1-e_{1})^{n_{1}(q)}+\sum_{q^{\prime}\neq
q}f(q^{\prime})e_{0}^{k_{0}(q,q^{\prime})}e_{1}^{k_{1}(q,q^{\prime})}(1-e_{0})^{u_{0}(q,q^{\prime})}(1-e_{1})^{u_{1}(q,q^{\prime})}.$
(13)

The expectation value of bit parity in the distribution $\tilde{f}$ is

$\begin{split}\langle\tilde{P}\rangle=&\sum_{q}P(q)f(q)(1-e_{0})^{n_{0}(q)}(1-e_{1})^{n_{1}(q)}\\\
&+\sum_{q}\sum_{q^{\prime}\neq
q}f(q^{\prime})P(q^{\prime})(-1)^{k_{0}(q,q^{\prime})+k_{1}(q,q^{\prime})}e_{0}^{k_{0}(q,q^{\prime})}e_{1}^{k_{1}(q,q^{\prime})}(1-e_{0})^{n_{0}(q^{\prime})-k_{0}(q,q^{\prime})}(1-e_{1})^{n_{1}(q^{\prime})-k_{1}(q,q^{\prime})}\\\
=&\sum_{q}P(q)f(q)(1-e_{0})^{n_{0}(q)}(1-e_{1})^{n_{1}(q)}\\\
&+\sum_{q^{\prime}}\sum_{k_{0}=0}^{n_{0}(q^{\prime})}\sum_{k_{1}=0}^{n_{1}(q^{\prime})}\begin{pmatrix}n_{0}(q^{\prime})\\\
k_{0}\end{pmatrix}\begin{pmatrix}n_{1}(q^{\prime})\\\
k_{1}\end{pmatrix}f(q^{\prime})P(q^{\prime})e_{0}^{k_{0}}e_{1}^{k_{1}}(1-e_{0})^{n_{0}(q^{\prime})-k_{0}}(1-e_{1})^{n_{1}(q^{\prime})-k_{1}}(-1)^{k_{0}+k_{1}}\mathbbm{1}_{k_{0}+k_{1}>0}\\\
=&\sum_{q^{\prime}}\sum_{k_{0}=0}^{n_{0}(q^{\prime})}\sum_{k_{1}=0}^{n_{1}(q^{\prime})}\begin{pmatrix}n_{0}(q^{\prime})\\\
k_{0}\end{pmatrix}\begin{pmatrix}n_{1}(q^{\prime})\\\
k_{1}\end{pmatrix}f(q^{\prime})P(q^{\prime})e_{0}^{k_{0}}e_{1}^{k_{1}}(1-e_{0})^{n_{0}(q^{\prime})-k_{0}}(1-e_{1})^{n_{1}(q^{\prime})-k_{1}}(-1)^{k_{0}+k_{1}}\\\
=&\sum_{q}P(q)f(q)(1-2e_{0})^{n_{0}(q)}(1-2e_{1})^{n_{1}(q)},\end{split}$ (14)

which is Eq. 5.

Figure F.5: In this Figure and in Fig. F.6, we generate random probability
distributions $f(q)$ that give the desired expected values of bit parity
$\langle P\rangle$, which are equally spaced between $-1$ and $1$.
$\langle\tilde{P}\rangle$ is then calculated using Eq. 5. The resulting data
are compared to Eq. 10, and the residuals $\hat{R}$ are plotted. In this
Figure, two qubits are measured, and the effect of the readout error is
bounded by Eq. 9 (blue lines) and approximated by 10 (red line). We set
$e_{1}=0.1$ and $e_{0}=0.05$.

Figure F.6: Same as Fig. F.5 but reading out 3, 4, 10, and 11 qubits,
respectively. Eq. 11 (red line) is a good approximation and becomes better as
the number of measured qubits increases.