11institutetext:  Leibniz-Institut für Polymerforschung Dresden, Institut
Theorie der Polymere, 01069 Dresden, Germany
Technische Universität Dresden, Institut für Theoretische Physik, 01069
Dresden, Germany
Institut für Theoretische Physik II: Weiche Materie, Heinrich-Heine-
Universität Düsseldorf, Düsseldorf, 40225, Germany

#  Escape Dynamics in an Anisotropically Driven Brownian Magneto-System

I. Abdoli 11 J.-U. Sommer 11 2 2 H. Löwen 33 A. Sharma 11 2 2112233

###### Abstract

Thermally activated escape of a Brownian particle over a potential barrier is
well understood within Kramers theory. When subjected to an external magnetic
field, the Lorentz force slows down the escape dynamics via a rescaling of the
diffusion coefficient without affecting the exponential dependence on the
barrier height. Here, we study the escape dynamics of a charged Brownian
particle from a two-dimensional truncated harmonic potential under the
influence of Lorentz force due to an external magnetic field. The particle is
driven anisotropically by subjecting it to noises with different strengths
along different spatial directions. We show that the escape time can largely
be tuned by the anisotropic driving. While the escape process becomes
anisotropic due to the two different noises, the spatial symmetry is restored
in the limit of large magnetic fields. This is attributed to the Lorentz force
induced coupling between the spatial degrees of freedom which makes the
difference between two noises irrelevant at high magnetic fields. The
theoretical predictions are verified by Brownian dynamics simulations. In
principle, our predictions can be tested by experiments with a Brownian
gyrator in the presence of a magnetic field.

##  1 Introduction

Escape of a particle from a metastable potential well is one of the most
celebrated problems in statistical physics [1]. The escape problem was first
conceived and studied by Kramers to model the rate of chemical reactions [2].
In the classical escape problem, which includes only thermal fluctuations, the
escape rate decreases exponentially with the barrier height.

The classical trap model has been generalized to include fluctuations other
than thermal fluctuations [3, 4, 5, 6, 7]. Such a scenario arises naturally in
active gels [8, 9], in which embedded tracer particles are subjected to both
thermal fluctuations and motor-induced, athermal fluctuations [10], and for
active particles [11, 12, 13, 14, 15, 16, 17]. While one models the thermal
fluctuations in the usual fashion as white Gaussian noise, the other noise is
often described as a colored noise to model the effect of active, athermal
fluctuations. The presence of an additional noise can dramatically modify the
escape dynamics of a particle [10, 18, 19, 20].

Here we explore the escape dynamics of a single Brownian particle from a
potential well subjected to two noises with different strengths along
different spatial directions, i.e., acting on the different Cartesian
components of the momentum vector. Such a system constitutes a primitive
Brownian engine composed of a gyrating Brownian particle [21]. This kind of
anisotropic driving can be experimentally realized by applying a strongly
fluctuating electric field to a charged Brownian particle in one direction
mimicking the role of an additional temperature [22]. Consequently, the
effective temperature along this direction takes a value which is higher than
the temperature of the bath. We show that breaking the spatial symmetry via
two different noises results in anisotropic escape dynamics; more particles
escape the potential well along the axis with larger noise strength which we
refer to as the hot axis. This gives rise to a large tunability of the escape
time by the difference between noise strengths or effective temperature
difference. However, when the particle is further subjected to Lorentz force
due to an external magnetic field, we show that the spatial symmetry is
restored in the limit of large magnetic fields; the escape is driven by an
effective noise with the average strength of the two noises. This is
attributed to the Lorentz force induced coupling between the spatial degrees
of freedom which makes the difference between noise strengths irrelevant at
high magnetic fields.

In a spatially isotropic Brownian magneto-system, the Lorentz force slows down
the escape dynamics via a trivial rescaling of the diffusion coefficient [23,
24, 25] without affecting the exponential dependence on the barrier height. In
contrast, in an anisotropically driven system, a magnetic field affects the
dynamics in a qualitatively different way [26]. While curving the trajectory
of a particle, there also occurs energy transfer in form of heat from the hot
source with a higher effective temperature to the cold source with lower
effective temperature, which is mediated by the magnetic field [27]. As a
consequence, at large magnetic fields, the two spatial degrees of freedom
become identical regardless of the effective temperature difference. In
general, it is difficult to study the escape problem in a multidimensional
potential well. Here we use an asymptotic method to compute the mean escape
time for a large barrier height, obtaining an explicit analytical prediction.

##  2 Anisotropically driven Brownian magneto-system

Figure 1: A single charged Brownian particle, subjected to an external
magnetic field $B$ and trapped in an isotropic harmonic potential
$V(x,y)=k(x^{2}+y^{2})/2$ with the parameter $k$. The particle is
simultaneously subjected to two noises with different strengths along its $x$
and $y$ degrees of freedom. The particle escapes more along the $y$ axis due
to the broken spatial symmetry by the two noises, which is shown in reddish
color. We show that escape distribution can be tuned and symmetrized by
varying the magnitude of the applied magnetic field.

We consider a Brownian magneto-system, which is made of a single diffusing
particle of the charge $q$ and mass $m$ in the presence of a constant magnetic
field $B$ in the $z$ direction. The $x$ and $y$ positions of the particle,
which we indicate in the vectorial form $\bm{r}=(x,y)^{\top}$, are coupled to
noises with different strengths mapping to effective temperatures $T_{x}$ and
$T_{y}$, respectively (see Fig. 1). The transpose is indicated by $\top$.
Since the Lorentz force does not affect the motion of the particle in the
direction of the applied magnetic field, we effectively investigate a two-
dimensional system in the $xy$ plane.

The dynamics of the particle with velocity $\bm{v}=(v_{x},v_{y})^{\top}$ which
is trapped in a potential in the form of
$\mathop{}\\!\mathrm{V}(\bm{r})=\frac{1}{2}\bm{r}^{\top}\cdot\hat{\mathrm{U}}\cdot\bm{r}$
can be described by the following underdamped Langevin equation:

$m\dot{\bm{q}}(t)=-\hat{\mathrm{F}}\bm{q}(t)+\bm{\xi}(t),$ (1)

where $\bm{q}(t)=(x(t),y(t),v_{x}(t),v_{y}(t))^{\top}$ and
$\bm{\xi}(t)=(0,0,\xi_{x}(t),\xi_{y}(t))^{\top}$ is Gaussian white noise with
zero mean and time correlation
$\langle\bm{\xi}(t)\bm{\xi}^{\top}(t^{\prime})\rangle=2\gamma\hat{\mathrm{T}}\delta(t-t^{\prime})$
where $\gamma$ is the friction coefficient. Throughout the paper we set the
Boltzmann constant $k_{B}$ to unity. For an isotropic potential
$\hat{\mathrm{U}}=k\hat{\mathbb{I}}$ with $k$ being the stiffness of the
potential. We consider a reference system at temperature $T$. The two
effective temperatures along Cartesian coordinates are measured in units of
$T$. Length and time are measured in units of $\sqrt{T/k}$ and $\gamma/k$,
respectively. Here
$\hat{\mathrm{T}}=\mathop{}\\!\bm{\mathrm{diag}}(0,0,T_{x},T_{y})$ is a
diagonal matrix and the matrix $\hat{\mathrm{F}}$ is defined as

$\hat{\mathrm{F}}=\left(\begin{array}[]{cc}\hat{\mathrm{0}}&-m\hat{\mathbb{I}}\\\
\hat{\mathrm{U}}&\hat{\mathrm{G}}\\\
\end{array}\right),\,\,\,\text{where}\,\,\,\,\hat{\mathrm{G}}=\gamma\left(\begin{array}[]{cc}1&-\kappa\\\
\kappa&1\\\ \end{array}\right).$ (2)

Here $\hat{\mathbb{I}}$ is the identity matrix and $\kappa=qB/\gamma$ is the
diffusive Hall parameter which quantifies the strength of the Lorentz force
relative to the frictional force. Note that the hat over the symbols indicates
the matrices and the vectors are shown by bold symbols.

We are interested in the overdamped dynamics of the particle. Usually this is
done by setting the inertia term in Eq. (1) to zero. However, in the presence
of a magnetic field, this yields an incorrect description of the overdamped
dynamics [28, 29] . A careful small-mass limit of the Langevin equation
reveals that the overdamped dynamics, though diffusive in nature, are
described an unusual odd-diffusion tensor that has both even (diagonal) and
odd (off-diagonal) elements. We note that the variance of the particle
position is determined only by the even part of the diffusion tensor [24, 26].
The odd part of the tensor, which has its origin in the broken time reversal
symmetry due to the magnetic field, gives rise to additional Lorentz
(rotational) fluxes which are missed if one takes the route of setting the
inertial term to zero. The Lorentz fluxes can play a crucial role in the
dynamical properties of Brownian systems [24, 26].

##  3 Escape from an isotropic trap

We consider a trapped particle in an isotropic potential
$\mathop{}\\!\mathrm{V}(\bm{r})$ and postulate that the particle escapes the
trap when it reaches the boundary, truncated at $r=a$, where $r=|\bm{r}|$ is
the distance from the origin, as shown in Fig. 1. An exact calculation of the
mean escape time is possible in a single temperature magneto-system, as
derived in details in the Supplemental Information (SI) given in [25], taking
advantage of the spatial isotropy. However, in an anisotropically driven
magneto-system, the broken spatial symmetry prevents such an approach. We
therefore use an asymptotic method to investigate the escape dynamics in our
Brownian magneto-system. We assume that the barrier height $\Delta
E=\frac{1}{2}ka^{2}$, is sufficiently large that the particle leaks out slowly
across the trap and settles into a quasistationary state: the escape is a
Poisson process with the inverse rate of mean escape time $\langle
t_{esc}\rangle$. The quasistationary probability density is given by
$P(\bm{r},t)\sim\rho_{ss}(\bm{r})e^{-t/\langle t_{esc}\rangle}$ where
$\rho_{ss}(\bm{r})$ is the steady-state probability density obtained in the
limit of $a\rightarrow\infty$ which we have obtained in previous works [26,
27] and in the SI in [25].

The total outgoing flux at $r=a$ is given by
$J(t)=-\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\int_{0}^{2\pi}\int_{0}^{a}P(r,\theta,t)r\mathop{}\\!\mathrm{d}r\mathop{}\\!\mathrm{d}\theta$
where $P(r,\theta,t)$ is the quasistationary probability distribution in the
polar coordinates. The outgoing flux can be alternatively defined as the
probability $a\int_{0}^{2\pi}\rho_{ss}(a,\theta)e^{-t/\langle
t_{esc}\rangle}\mathop{}\\!\mathrm{d}\theta$ to be on the boundary at time
$t$, times the velocity of the fluctuation path leading to the boundary [30,
10]. Using the two equivalent definitions of the outgoing flux the mean escape
time is equal to the inverse of the stationary-state probability density on
the boundary besides a prefactor. The prefactor can be determined from the
exact analytical result for a spatially isotropic magneto-system, as shown in
the SI in [25]. By doing so, we obtain the mean escape time (see the SI in
[25] for details), which can be written as

$\langle
t_{esc}\rangle^{-1}\approx\frac{2}{1+\kappa^{2}}\sqrt{\frac{T\delta_{E}^{2}}{\bar{T}(1-\delta_{\kappa}^{2})}}e^{-\frac{\delta_{E}}{1-\delta_{\kappa}^{2}}}\mathrm{I}_{0}\left(-\frac{\delta_{\kappa}\delta_{E}}{1-\delta_{\kappa}^{2}}\right),$
(3)

where $\mathrm{I}_{0}(x)$ is the modified Bessel function of the first kind of
the order zero, $\bar{T}=(T_{x}+T_{y})/2$ is the effective temperature
average, $\delta_{E}=\Delta E/\bar{T}$ is the scaled barrier height, and

$\delta_{\kappa}=\frac{\Delta T}{2\bar{T}\sqrt{1+\kappa^{2}}},$ (4)

is the scaled effective temperature difference, given as $\Delta
T=T_{y}-T_{x}$, which is an explicit function of the diffusive Hall parameter.

Figure 2: The mean escape time as a function $\kappa$ for different barrier
heights $\delta_{E}$ with the effective temperature difference $\Delta T=2.0$
and effective temperature average $\bar{T}=3.0$. Solid lines show the
theoretical prediction from Eq. (3) and the symbols depict the simulation
results. The black, dashed line is a plot of the prefactor $1+\kappa^{2}$. In
the limit of large magnetic fields, the overall trend of the mean escape time
is determined by this prefactor. The inset shows the distribution of the
escape time for $\kappa=4.0$ and $\delta_{E}=9.0$. The solid line is the
Poisson distribution $P(t)=\langle t_{esc}\rangle^{-1}\exp(-t/\langle
t_{esc}\rangle)$ and the symbols depict the simulation results. The
theoretical predictions overestimate the escape time for small barrier
heights. Figure 3: The mean escape time as function of the scaled effective
temperature difference $\delta_{\kappa}$ for different barrier heights. Solid
lines show the results from the theoretical prediction in Eq. (3) and the
symbols depict the simulation results. While a small $\delta_{\kappa}$
corresponds to either large magnetic fields or small effective temperature
difference, $\delta_{\kappa}\rightarrow 1$ corresponds to either small
magnetic fields or large effective temperature difference. Simulation data is
obtained for $\Delta T=1.0,2.0,3.0,4.0,5.0$ with $\bar{T}=3.0$. The diffusive
Hall parameter varies from $0.0$ to $10.0$. The dashed lines show the
saturation to the system with temperature $\bar{T}$ corresponding to either
the effectively zero-temperature difference or large-magnetic-field limit. The
inset shows the angular distribution for different values of the diffusive
Hall parameter, $\kappa=0.0,0.5,2.0,10.0$ from left to right, respectively.
The effective temperature difference is $\Delta T=2.0$ and the effective
temperature average is $\bar{T}=3.0$. In the absence of a magnetic field the
particle mostly escapes along the $y$ axis due to the broken spatial symmetry
via the two noises with different strengths, mapped to effective temperatures.

In a spatially isotropic system the mean escape time in Eq. (3) reduces to the
familiar expression from Kramers theory $\langle
t_{esc}\rangle^{-1}\approx(1+\kappa^{2})^{-1}2\delta_{E}\exp(-\delta_{E})$, as
shown in the SI in [25]. In the general case of two different noises, the
escape time also depends on the effective temperature difference via the
parameter $\delta_{\kappa}$. For a fixed $\kappa$, increasing the effective
temperature difference accentuates the spatial anisotropy of the escape
dynamics. However, in the limit of large magnetic fields
$\delta_{\kappa}\rightarrow 0$, independent of the effective temperature
difference, in which case one again obtains the Kramers escape rate for a
spatially isotropic system with temperature $\bar{T}$. The $\kappa$-governed
crossover from anisotropic to isotropic escape can be qualitatively understood
as following. Lorentz force curves the trajectory of a moving particle.
However, since the two spatial degrees of freedom are subjected to noises with
different strengths, the curving of the trajectory also results in energy
transfer in the form of heat from the hot source to the cold one.
Consequently, in the limit of large magnetic fields, the two spatial degrees
of freedom become identical regardless of the effective temperature
difference.

To validate our theoretical predictions we perform Brownian dynamics
simulations using the Langevin equations of motion. It has been shown that the
overdamped Langevin equation for Brownian motion in the presence of a magnetic
field can give rise to unphysical values for velocity-dependent variables such
as fluxes. Therefore, we use the underdamped Langevin equation, given in Eq.
(1), with a sufficiently small mass.  The change in the velocity due to the
Gaussian noise is of the order $\sqrt{dt}/m$. By choosing the integration time
step $dt<m^{2}$ we ensure that the velocity increments are small. Throughout
the paper, we fix the mass to $m=0.001$ and the integration time step to
$dt=5\times 10^{-7}\tau$, where $\tau$ serves as a natural time scale which we
set to unity. The temperature $T$ of the reference system and the unit of the
length are set to unity in the simulations. We consider a total number of
particles of $1000$.

Figure 2 shows the mean escape time as a function of the diffusive Hall
parameter $\kappa$ for different values of the scaled barrier height
$\delta_{E}$. The data are obtained from simulations of a system with the
effective temperature average $\bar{T}=3.0$ and effective temperature
difference $\Delta T=2.0$. Our theoretical predictions are in good agreement
with the simulation results. In the SI given in [25] we show that the large
barrier approximation of the exact expression for escape time in a single
temperature system is highly accurate compared to the known Kramers results
when the expression is scaled by an empirical factor $1.2$. We use the same
factor in our asymptotic approach and observe that our theoretical predictions
are in good agreement with the simulation results.

While in the limit of large magnetic fields the overall trend of the curves
are determined by the prefactor $1+\kappa^{2}$, as shown in Fig. 2 by dashed
line, it deviates from this prefactor for small magnetic fields. In other
words, while the anisotropic driving of the system effects the escape dynamics
for small magnetic fields, the spatial symmetry is restored in the limit of
large magnetic fields.

Figure 3 shows the scaled mean escape time as a function of the parameter
$\delta_{\kappa}$. The Brownian dynamics simulation has been done for systems
with effective temperature differences varying from $\Delta T=1.0$ to $5.0$
and the same effective temperature average $\bar{T}=3.0$. The diffusive Hall
parameter is varying from $\kappa=0.0$ to $10.0$. In the limit of large
magnetic fields or small temperatures, the scaled effective temperature
difference parameter $\delta_{\kappa}\rightarrow 0$. Small magnetic fields or
large effective temperature difference correspond to
$\delta_{\kappa}\rightarrow 1$. The scaled mean escape time decreases with
increasing parameter $\delta_{\kappa}$. The dashed lines show the saturation
to a system with temperature $\bar{T}$ corresponding to effectively zero-
temperature difference or alternatively to large-magnetic-field limit. In the
inset, we show the angular distribution from simulations for different values
of the diffusive Hall parameter. As can be seen, in the absence of a magnetic
field the particle mostly escapes along the $y$ axis due to the broken spatial
symmetry via the two different noises. The escape becomes symmetric in the
limit of a large magnetic field, namely the effective temperature difference
becomes irrelevant.

##  4 Conclusions

We studied the escape dynamics in an anisotropically driven Brownian magneto-
system. We derived an approximate expression for the mean escape time taking
into account two noises with different strength along different spatial
directions, mapped to two effective temperatures. We showed that Lorentz force
induces a coupling between the spatial degrees of freedom which, in the limit
of large magnetic fields, restores the spatial symmetry; the two spatial
degrees of freedom become identical regardless of the effective temperature
difference.

While experimental realisation of the proposed magneto-system might be
particularly difficult in colloidal systems, it is more realistic in a complex
plasma where extremely large magnetic fields can be generated [31]. In such a
set up, not only can one study the escape dynamics, one can also measure the
momenta governed energy flow between the hot source and the cold bath. Such
heat currents have been predicted theoretically, however, the studies employed
toy models in which a particle is coupled to multiple thermostats [32, 33]. A
complex plasma, with anisotropic driving, thus presents a potentially
realizable system to study momenta governed heat currents. The anisotropic
noise strength can also be realized by putting the charged object in contact
with active particles that provide the noise. If the latter are propelling in
an anisotropic way, as for instance on an anisotropic patterned substrate [34]
or in a nematic liquid crystalline cell [35], their effect on the charged
particle would realize an anisotropic noise strength. Other possible
realizations in which a confined particle experiences random kicks in an
anisotropic environment under the additional action of a Lorentz-, Coriolis-
or Magnus force are also conceivable.

In future work we aim to extend our model to the escape through a narrow hole
[36, 37, 38]. An interesting generalization would be to replace the isotropic
potential by an anisotropic potential where the stationary state of the system
is characterized by not only a non-Boltzmann probability distribution, but
also additional Lorentz fluxes [27]. Also, it could be interesting to study
the escape in an anisotropically driven Brownian magneto-system from a
bistable state, which might be realized by trapping the particle using two
optical tweezers [39]. Another study of interest would be a system in the
presence of a fluctuating magnetic field where a new turnover is observed as a
generic signature of the system [40, 41]. Finally, our analysis might be
applicable to other systems such as chiral colloidal microswimmers in
parabolic potentials [42], active Janus particles in a complex plasma [43],
and even rotating skyrmions [44, 45].

##  5 Acknowledgments

A. Sharma acknowledges the support by the Deutsche Forschungsgemeinschaft
(DFG) within the project SH 1275/3-1.

## References

  * [1] Hänggi P., Talkner P. Borkovec M. Reviews of Modern Physics621990251.
  * [2] Kramers H. A. Physica71940284.
  * [3] Hänggi P., Marchesoni F. Grigolini P. Zeitschrift für Physik B Condensed Matter561984333.
  * [4] Jung P. Hänggi P. Physical Review Letters61198811.
  * [5] Sharma A., Wittmann R. Brader J. M. Physical Review E952017012115.
  * [6] Scacchi A., Brader J. M. Sharma A. Physical Review E1002019012601.
  * [7] Caprini L., Marini Bettolo Marconi U., Puglisi A. Vulpiani A. The Journal of Chemical Physics1502019024902.
  * [8] Mizuno D., Tardin C., Schmidt C. F. MacKintosh F. C. Science3152007370.
  * [9] Sheinman M., Sharma A., Alvarado J., Koenderink G. MacKintosh F. Physical Review Letters1142015098104.
  * [10] Woillez E., Kafri Y. Gov N. S. Physical Review Letters1242020118002.
  * [11] Burada P. Lindner B. Physical Review E852012032102.
  * [12] Koumakis N., Maggi C. Di Leonardo R. Soft Matter1020145695.
  * [13] Schneider E. Stark H. EPL (Europhysics Letters)127201964003.
  * [14] Wexler D., Gov N., Rasmussen K. Ø. Bel G. Physical Review Research22020013003.
  * [15] Debnath T., Chaudhury P., Mukherjee T., Mondal D. Ghosh P. K. The Journal of Chemical Physics1552021194102.
  * [16] Schwarzendahl F. J. Löwen H. EPL (Europhysics Letters)134202148005.
  * [17] Zanovello L., Faccioli P., Franosch T. Caraglio M. The Journal of Chemical Physics1552021084901.
  * [18] Militaru A., Innerbichler M., Frimmer M., Tebbenjohanns F., Novotny L. Dellago C. Nature Communications1220212446.
  * [19] Caprini L., Cecconi F. Marini Bettolo Marconi U. The Journal of Chemical Physics2021.
  * [20] Gera T. Sebastian K. The Journal of Chemical Physics1552021014902.
  * [21] Filliger R. Reimann P. Physical Review Letters992007230602.
  * [22] Argun A., Soni J., Dabelow L., Bo S., Pesce G., Eichhorn R. Volpe G. Physical Review E962017052106.
  * [23] Filliger R. Reimann P. EPL (Europhysics Letters)77200730008.
  * [24] Abdoli I., Vuijk H. D., Sommer J.-U., Brader J. M. Sharma A. Physical Review E1012020012120.
  * [25] In the Supplemental Information we show that the large-barrier approximation for a system with temperature $\bar{T}$ gives $\langle t_{esc}\rangle_{eq}\approx(1+\kappa^{2})\gamma e^{\delta_{E}}/2k\delta_{E}$ where the temperature difference is zero and the escape occurs due to the equilibrium fluctuations.
  * [26] Abdoli I., Kalz E., Vuijk H. D., Wittmann R., Sommer J.-U., Brader J. M. Sharma A. New Journal of Physics222020093057.
  * [27] Abdoli I., Wittmann R., Brader J. M., Sommer J.-U., Löwen H. Sharma A. arXiv preprint arXiv:2110.052842021.
  * [28] Chun H.-M., Durang X. Noh J. D. Physical Review E972018032117.
  * [29] Vuijk H. D., Brader J. M. Sharma A. Journal of Statistical Mechanics: Theory and Experiment20192019063203.
  * [30] Risken H. Fokker-Planck Equation (Springer) 1996 pp. 63–95.
  * [31] Kählert H., Carstensen J., Bonitz M., Löwen H., Greiner F. Piel A. Physical Review Letters1092012155003.
  * [32] Murashita Y. Esposito M. Physical Review E942016062148.
  * [33] Lee J. S., Park J.-M., Chun H.-M., Um J. Park H. Physical Review E1012020052132.
  * [34] Volpe G., Buttinoni I., Vogt D., Kümmerer H.-J. Bechinger C. Soft Matter720118810.
  * [35] Zhou S., Tovkach O., Golovaty D., Sokolov A., Aranson I. S. Lavrentovich O. D. New Journal of Physics192017055006.
  * [36] Schuss Z., Singer A. Holcman D. Proceedings of the National Academy of Sciences104200716098.
  * [37] Bressloff P. C. Stochastic Processes in Cell Biology Vol. 41 (Springer) 2014.
  * [38] Nayak S., Debnath T., Das S., Debnath D. Ghosh P. K. The Journal of Physical Chemistry C124202018747.
  * [39] Ciliberto S. Physical Review X72017021051.
  * [40] Baura A., Ray S. Bag B. C. The Journal of Chemical Physics1382013244110.
  * [41] Mondal S., Baura A., Das S. Bag B. C. Physica A: Statistical Mechanics and its Applications502201858.
  * [42] Jahanshahi S., Löwen H. Ten Hagen B. Physical Review E952017022606.
  * [43] Nosenko V., Luoni F., Kaouk A., Rubin-Zuzic M. Thomas H. Physical Review Research22020033226.
  * [44] Zhang S., Wang W., Burn D., Peng H., Berger H., Bauer A., Pfleiderer C., Van Der Laan G. Hesjedal T. Nature Communications920182115.
  * [45] Brown B. L., Täuber U. C. Pleimling M. Physical Review B972018020405.


# A blueprint for a Digital-Analog Variational Quantum Eigensolver using
Rydberg atom arrays

Antoine Michel<EMAIL_ADDRESS>Electricité de France, EDF Recherche et
Développement, Département Matériaux et Mécanique des Composants, Les
Renardières, F-77250 Moret sur Loing, France  Université Paris-Saclay,
Institut d’Optique Graduate School, CNRS, Laboratoire Charles Fabry, F-91127
Palaiseau Cedex, France  Sebastian Grijalva  Loïc Henriet  PASQAL, 7 rue
Léonard de Vinci, F-91300 Massy, France  Christophe Domain  Electricité de
France, EDF Recherche et Développement, Département Matériaux et Mécanique des
Composants, Les Renardières, F-77250 Moret sur Loing, France  Antoine Browaeys
Université Paris-Saclay, Institut d’Optique Graduate School, CNRS, Laboratoire
Charles Fabry, F-91127 Palaiseau Cedex, France

(August 27, 2024)

###### Abstract

We address the task of estimating the ground-state energy of Hamiltonians
coming from chemistry. We study numerically the behavior of a digital-analog
variational quantum eigensolver for the H2, LiH and BeH2 molecules, and we
observe that one can estimate the energy to a few percent points of error
leveraging on learning the atom register positions with respect to selected
features of the molecular Hamiltonian and then an iterative pulse shaping
optimization, where each step performs a derandomization energy estimation.

##  I Introduction

Quantum simulation holds the promises to solve outstanding questions in many-
body physics, in particular finding the ground state of strongly interacting
quantum systems [1, 2]. The determination of the ground state energies of
complex molecules, one of the main tasks in quantum chemistry, is therefore an
example of application where quantum simulation could be of interest. Along
this line, proof-of-principle demonstrations were obtained using photons [3,
4], ions [5, 6, 7] or quantum circuits [8]. The last two examples used an
hybrid approach were a classical computer optimizes in an iterative way the
results obtained by a quantum device that was operating in a digital mode,
i.e. as a series of one and two-qubit gates.

Rydberg quantum simulators are another example of promising quantum simulation
platforms thanks to their potential for scaling the number of qubits and their
programmability [9]. They rely on individual atoms trapped in arrays of
optical tweezers that can interact when promoted to Rydberg states. The
platform naturally implements spin Hamiltonians. Analog quantum simulation
with hundreds of atoms has now been achieved [10, 11, 12].

One appealing feature of this platform is the ability to place the atoms in
arbitrary position in two and three dimensions, thus allowing large
flexibility in their connectivity. Another feature is their ability to prepare
different initial product states as heuristic trials before the unitary
evolution (whether it is by a set of digital gates or the action of an analog
Hamiltonian evolution). However, this freedom in register preparation has a
significant time cost that adds to the repetition clock rate [13].

Figure 1: Digital-Analog VQE on Rydberg atoms studied in this paper. The
computer icons indicate where classical processing enhances and informs the
cycle. First, a target molecular Hamiltonian $\hat{H}_{\text{T}}$ is passed to
computer _(a)_. This computer extracts a subset of terms that define a target
interaction matrix $V^{\text{T}}$. Computer _(b)_ then optimizes atom
positions (with interaction matrix $V^{\text{R}}$) to approximate
$V^{\text{T}}$. With a chosen register, computer _(a)_ tests virtually some
product states that can be experimentally implemented, warm-starting the
algorithm. Meanwhile, computer _(c)_ takes the Pauli strings
$\\{\hat{P}^{(s)}\\}$ that constitute $\hat{H}_{\text{T}}$ and outputs a
_derandomized_ set of Pauli measurement basis $\\{\hat{R}(\hat{P}^{(m)})\\}$
that will be used at the end of the circuit where the prepared state is
measured in the $Z$-axis. A parameterized pulse acts on the initial state, and
the readout data is used by computer _(d)_ to estimate the expectation value
of the target Hamiltonian $\langle\hat{H}_{\text{T}}\rangle$. Finally,
computer _(d)_ calculates new pulse parameters to update the quantum
evolution. When a desired precision or some stopping criteria is reached, the
best energy value is returned.

Neutral atom devices are naturally suited for _analog quantum algorithms_ ,
where the analog blocks are represented by control pulses that drive the
system (or subsets of it). Given a prepared state, the parameterized pulses
can be adjusted to variationally improve on a given score of the state.
Methods for the optimization of parameters have been the subject of intense
exploration in recent years [14, 15, 2, 16, 17, 18, 19, 20, 21, 22].
Additionally, the information and “cost functions” from the prepared quantum
system are obtained by repeatedly measuring the state in the computational
basis, which constitutes an operational overhead. Recent results [23, 24, 25,
26, 11, 27] on protocols for the estimation of quantum observables are
available and have helped establishing efficient measurement procedures based
on generalized random measurements and a series of post-processing steps that
are performed on a classical computer and that alleviate the measurement
overhead. The types of randomized measurements that we shall describe in this
paper require local rotations on the qubits of the register, thus constituting
another “digital” layer, from a quantum circuit perspective. In fact,
_digital-analog_ algorithms [28], benefit from the fact that analog operations
can be performed with much higher fidelities than when using digital gates,
while local single-qubit gates can be added explicitly in crucial steps of the
process (state preparation and measurement).

In this paper, we explore the implementation of a digital-analog VQE algorithm
in a Rydberg quantum simulator. We account for typical constraints of the
platform: the local action is restrained to the initial state preparation and
measurement, with Hamiltonian time-evolution acting on the entire system as a
“global” gate. We study numerically this version of a VQE for the
$\text{H}_{\text{2}}$ molecule using common ansatze, followed by a more
efficient protocol for larger molecules. We discuss the embedding of the
Hamiltonian in the atom register, the way in which the optimization of the
pulse sequence can be performed and the necessity of including an efficient
estimation of energies (namely, we explore the effect of a _derandomization_
estimation [23]) at each iteration step. We apply this numerically to the
examples of LiH and BeH2. The manuscript is organized as follows: In section
II we recall how the Variational Quantum Eigensolver (VQE) estimates the
energy of the ground state of a molecular-based Hamiltonian. We then describe
the basic ingredients of Rydberg Atom Quantum Processors and the Hamiltonians
that they implement. We end the section by explaining the optimization cycle
of variational quantum algorithms on these devices. In section III we describe
the strategies for implementation of the VQE, going from a direct application
of a Unitary Couple Cluster Ansatz, to the Quantum Alternating Operator ansatz
and finally to a more hardware-oriented approach that combines elements of
register preparation, pulse optimization and observable estimation. This is
followed in section IV by numerical results of the error in energy obtained as
a function of the number of repetitions of the experiment, an informative
measure of the performance of hybrid classical-quantum implementations.

##  II Analog Variational Quantum Eigensolver with Rydberg atoms

The Variational Quantum Eigensolver (or VQE) is a hybrid quantum-classical
algorithm designed to find the lowest eigenvalue of a given Hamiltonian [29].
We describe below the origin of the Hamiltonians that we consider and how VQE
can be studied with a Rydberg Quantum Processor.

###  II.1 Hamiltonians from Quantum Computational Chemistry

We first recall the method used to express the electronic Hamiltonian as a
spin model (see e.g. [7]). We start from the Born-Oppenheimer approximation of
the Hamiltonian of the system, which considers the nuclei of the molecules as
classical point charges:

$\hat{H}=-\sum_{i}\frac{\nabla_{i}^{2}}{2}-\sum_{i,I}\frac{\mathcal{Z}_{I}}{|\mathbf{r}_{i}-\mathbf{R}_{I}|}+\frac{1}{2}\sum_{i\neq
j}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}$ (1)

(in atomic units) where $\nabla_{i}$ is the kinetic energy term for the $i$-th
electron, $\mathcal{Z}_{I}$ is the charge of the $I$-th nucleus, and
$\mathbf{r}$, $\mathbf{R}$ denote the distance of the $i$-th electron and the
$I$-th nucleus with respect to the center of mass, respectively. We aim to
obtain the ground state energy of (1).

One needs to define a basis set in which to represent the electronic
wavefunctions. We shall concentrate on the Slater-type orbital approximation
for the basis set, with three Gaussian functions, STO-3G. This minimal basis
set $\\{\phi_{i}(\mathbf{x}_{i})\\}$ (where
$\mathbf{x}_{i}=(\mathbf{r}_{i},\sigma_{i})$ encodes the $i$-th electron’s
spatial and spin coordinates) includes the necessary orbitals to represent the
valence shell of an atom. Moreover, the wavefunctions need to be anti-
symmetric under the exchange of electrons. This can be achieved through
_second quantization_ , where one defines anticommuting fermionic
creation/annihilation operators
$\\{\hat{a}_{p}^{\dagger}\\},\\{\hat{a}_{p}\\}$ and rewrites the initial
Slater determinant form of the wavefunction as
$|\Psi\rangle=\prod_{p}(\hat{a}_{p}^{\dagger})^{\phi_{p}}|\mathrm{vacuum}\rangle$,
representing the occupation of each molecular orbital.

The fermionic operators are used to rewrite (1) as:

$\hat{H}=\sum_{p,q}h_{pq}\hat{a}_{p}^{\dagger}\hat{a}_{q}+\frac{1}{2}\sum_{p,q,r,s}h_{pqrs}\hat{a}_{p}^{\dagger}\hat{a}_{q}^{\dagger}\hat{a}_{r}\hat{a}_{s}.$
(2)

where the coefficients $h_{pq}$ and $h_{pqrs}$ encode the spatial and spin
configuration of each of the electrons and depend on the inter-nuclear and
inter-electron distances $\mathbf{R},\mathbf{r}$:

$\displaystyle h_{pq}$ $\displaystyle=\int
d\mathbf{x}\phi_{p}^{*}(\mathbf{x})\left(-\frac{\nabla^{2}}{2}-\sum_{i}\frac{\mathcal{Z}_{i}}{|\mathbf{R}_{i}-\mathbf{r}|}\right)\phi_{q}(\mathbf{x})$
(3) $\displaystyle h_{pqrs}$ $\displaystyle=\int
d\mathbf{x}_{1}d\mathbf{x}_{2}\frac{\phi_{p}^{*}(\mathbf{x}_{1})\phi_{p}^{*}(\mathbf{x}_{2})\phi_{r}(\mathbf{x}_{1})\phi_{s}(\mathbf{x}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}.$

Next, we map the fermionic operators acting on Fock states of $n$ orbitals to
a Hilbert space of operators acting on spin states of $N$ qubits. This
corresponds to the quantum processors’ effective interaction Hamiltonians,
quantum gates and measurement basis. Useful maps of this kind include the
Jordan-Wigner (JW) [30] or the Bravyi-Kitaev (BK) [31] transformations. The
obtained Hamiltonian is a sum of tensor products of single-qubit Pauli
matrices:

$\displaystyle\hat{H}_{T}=\sum_{s=1}^{S}{\bf
c}_{s}\bigg{(}\bigotimes_{j=1}^{N}\hat{P}_{j}^{(s)}\bigg{)}$ (4)

where $\hat{P}_{j}\in\\{\mathbb{1},X,Y,Z\\}$, $S$ is the number of Pauli
strings in the Hamiltonian and $N$ the number of qubits.

###  II.2 Rydberg Atom Quantum Processor

Rydberg atom arrays are now well-established quantum simulation platforms [13,
9]. Briefly, atoms are trapped in optical tweezers, each containing exactly
one atom. The tweezers may be arranged in any 1D, 2D or 3D geometrical
configurations. The register can be rebuilt after each computational cycle. To
perform quantum processing, we use the fact that the platform implements spin-
like Hamiltonians, where the interactions originate from strong dipole-dipole
couplings between atoms laser-excited to Rydberg states.

Depending on the choice of atomic levels, the atoms experience different
effective interactions. In the case of the Ising mode, $\ket{0}$ is a “ground”
state prepared by optical pumping [9] and $\ket{1}$ is a Rydberg state of the
atom. The Hamiltonian term for this interaction is:

$\displaystyle\hat{H}_{\text{Ising}}$
$\displaystyle=\sum_{i>j}\frac{C_{6}}{r_{ij}^{6}}\hat{n}_{i}\hat{n}_{j},$ (5)

with $\hat{n}_{i}=|1\rangle_{i}\langle 1|=(\mathbb{1}_{i}+Z_{i})/2$ the
projector on the Rydberg state and $r_{ij}$ the distance between atoms. Here
and below, $X_{i},Y_{i}$ and $Z_{i}$ indicate the local Pauli operators.

If instead the two states chosen are two dipole-coupled Rydberg states (for
example $\ket{0}=\ket{nS}$ and $\ket{1}=\ket{nP}$ for large $n$), the
interaction is resonant and realizes a so-called “XY” or “flip-flop” term:

$\displaystyle\hat{H}_{\text{XY}}$ $\displaystyle=\sum_{i\neq
j}\frac{C_{3}}{r_{ij}^{3}}(X_{i}X_{j}+Y_{i}Y_{j}),$ (6)

where $C_{3}$ depends on the chosen Rydberg orbitals and their orientation
with respect to the interatomic axis. It corresponds to a coherent exchange of
neighboring spin states $\ket{10}$ to $\ket{01}$.

In addition, we can include time-dependent terms on the Hamiltonian, by means
of a laser pulse (Ising mode) or a microwave field (XY mode) targeting the
transition between the ground and excited states. This is represented by the
following “drive” terms:

$\displaystyle\hat{H}_{\text{drive}}=\frac{\hbar}{2}\sum_{i=1}^{N}\Omega_{i}(t)X_{i}-\hbar\sum_{i=1}^{N}\delta_{i}(t)\hat{n}_{i}\
.$ (7)

Here, $\Omega(t)$ is Rabi frequency and $\delta(t)$ the detuning of the field
with respect to the resonant transition frequency. The addressing can be
either global or local. In the procedure used in this work, the local
addressing is restricted to the initial state preparation and the register
readout stages.

###  II.3 Variational Algorithms on a Rydberg atoms device

In the analog VQE algorithm, we seek to estimate the energy of the ground
state of a qubit Hamiltonian called the _target_ Hamiltonian,
$\hat{H}_{\text{T}}$, by using an iterative method. The _resource_ Hamiltonian
is the one realized by the hardware, and can be configured with different
types of interactions ($\hat{H}_{\mathrm{inter}}$) (5, 6) and driving fields
($\hat{H}_{\mathrm{drive}}$):

$\hat{H}_{\text{R}}=\hat{H}_{\rm inter}+\hat{H}_{\rm drive}.$ (8)

Experimentally, the transition from the ground to the excited state is
typically generated by a two-photon process, from which an approximate two
level system is extracted, driven by an effective Rabi frequency $\Omega$ and
detuning $\delta$ during the quantum processing stage. We use their values as
parameters in our analog presentation of a VQE algorithm: The first step is to
prepare the register of $N$ atoms with a geometry that determines the
interaction terms $\hat{H}_{\rm inter}$ and then to initialize the system in a
state $\ket{\psi_{0}}$. Then, a pulse sequence is applied to evolve the system
under the resource Hamiltonian $\hat{H}_{\text{R}}(\Omega(t),\delta(t))$ whose
corresponding time-ordered unitary evolution operator is
$U(t)=\mathcal{T}\exp\big{(}-i\int_{0}^{t}\hat{H}_{\text{R}}(\Omega(\tau),\delta(\tau))d\tau\big{)}$.
The final prepared state is:

$\ket{\psi(\Omega,\delta,t)}=U(t)\ket{\psi_{0}}.$ (9)

The energy of a prepared state will be calculated with respect to the target
Hamiltonian:

$E(\Omega,\delta,t)=\bra{\psi(\Omega,\delta,t)}\hat{H}_{\text{T}}\ket{\psi(\Omega,\delta,t)}.$
(10)

After each cycle, a classical optimizer adjusts the parameters
$\Omega\rightarrow\Omega^{\prime}$, $\delta\rightarrow\delta^{\prime}$ and
$t\rightarrow t^{\prime}$ and we repeat the evolution of the initial quantum
state $\ket{\psi_{0}}$ with the new parameter set
$U(\Omega^{\prime},\delta^{\prime},t^{\prime})$. We aim to obtain for each
iteration $E(\Omega^{\prime},\delta^{\prime},t^{\prime})\leq
E(\Omega,\delta,t)$ 111This classical optimization problem can be addressed
for example by obtaining the gradient of the energy function.. After several
iterations of this loop, the variational scheme attempts to prepare a state
whose energy is a good approximation of the ground state energy of
$\hat{H}_{\text{T}}$ [2].

##  III Description of the Protocols

In this section, we describe two analog variational quantum algorithms for the
estimation of the ground state energy and apply them to quantum chemistry
problems. The protocols differ mainly by the choice of ansatz: one is the
_Unitary Coupled Cluster_ (UCC) ansatz [32], while the other is an adaptation
of a _hardware-efficient ansatz_ [8], based on repeating alternating values of
amplitude, frequency or phase of the applied pulses. We verify numerically the
performance of these two types of ansatz in a Rydberg-based Quantum Processor
(QP). Next we discuss a protocol for larger molecules tailored after the
hardware capabilities. We begin by considering the prototypical example of the
H2 molecule.

###  III.1 UCC ansatz on an analog quantum processor: application on H2

Numerous implementations of the VQE algorithm rely on the use of digital
gates. Recent experimental implementations for the H2, LiH and BeH2 molecules
have been realized in [8, 7], with superconducting and trapped ions devices
respectively. For the analog version of this algorithm on H2, we consider the
target Hamiltonian and the ansatz as in [7]. The Jordan-Wigner and Bravyi-
Kitaev transformations lead to two different spin Hamiltonians of this
molecule:

$\displaystyle\begin{split}\hat{H}_{\text{JW}}=&\,{\bf c_{0}}\mathbb{1}+{\bf
c_{1}}(Z_{0}+Z_{1})+{\bf c_{2}}(Z_{2}+Z_{3})+\\\ &{\bf c_{3}}Z_{3}Z_{2}+{\bf
c_{4}}Z_{2}Z_{0}+{\bf c_{5}}(Z_{2}Z_{0}+Z_{3}Z_{1})+\\\ &{\bf
c_{6}}(Z_{2}Z_{1}+Z_{3}Z_{0})+{\bf c_{7}}(X_{3}Y_{2}Y_{1}X_{0}+\\\
&Y_{3}X_{2}X_{1}Y_{0}-X_{3}X_{2}Y_{1}Y_{0}+Y_{3}Y_{2}X_{1}X_{0})\end{split}$
(11)

and

$\displaystyle\begin{split}\hat{H}_{\text{BK}}=&\,{\bf f_{0}}\mathbb{1}+{\bf
f_{1}}Z_{0}+{\bf f_{2}}Z_{1}+{\bf f_{3}}Z_{2}+{\bf f_{4}}Z_{1}Z_{0}+\\\ &{\bf
f_{5}}Z_{2}Z_{0}+{\bf f_{6}}Z_{3}Z_{1}+{\bf f_{7}}X_{2}Z_{1}X_{0}+{\bf
f_{8}}Y_{2}Z_{1}Y_{0}+\\\ &{\bf f_{9}}Z_{2}Z_{1}Z_{0}+{\bf
f_{10}}Z_{3}Z_{2}Z_{0}+{\bf f_{11}}Z_{3}Z_{2}Z_{1}+\\\ &{\bf
f_{12}}Z_{3}X_{2}Z_{1}X_{0}+{\bf f_{13}}Z_{3}Y_{2}Z_{1}Y_{0}+{\bf
f_{14}}Z_{3}Z_{2}Z_{1}Z_{0}\end{split}$ (12)

where the coefficients $\\{\bf c_{i}\\}$ and $\\{\bf f_{j}\\}$ are calculated
from (3). Since in (12) qubits $1$ and $3$ are only affected by the operators
$\mathbb{1}$ and $Z$ one can actually work with the following two-qubit
effective Hamiltonian [33]:

$\displaystyle\begin{split}\hat{H}_{\text{BK}}^{\text{(eff)}}=&\,{\bf
g_{0}}\mathbb{1}+{\bf g_{1}}Z_{0}+{\bf g_{2}}Z_{1}+{\bf g_{3}}Z_{0}Z_{1}+\\\
&{\bf g_{4}}X_{0}X_{1}+{\bf g_{5}}Y_{0}Y_{1}.\end{split}$ (13)

Usually, a good ansatz
$|\psi(\boldsymbol{\theta})\rangle=U(\boldsymbol{\theta})|\psi_{0}\rangle$
requires a balance between hardware constraints and symmetries in target
Hamiltonian. However, using the ‘knobs’ available on the hardware is often not
efficient, and one thus needs additional guidance to reach the states we are
looking for in a potentially very large Hilbert space. In this sense, the
well-established Unitary Coupled Cluster (UCC) ansatz allows one to perform an
unitary operation while keeping advantages of coupled cluster ansatz from
chemistry [34].

In most cases, implementing the UCC ansatz in a quantum processor requires
constructing a digital quantum circuit with full local addressing. An example
where global addressing is sufficient is the H2 molecule. The initial guess of
the molecular wave function is a product state obtained from the classical
Hartree-Fock calculation performed to determine the coefficients $\\{\bf
c_{i}\\}$ and $\\{\bf f_{j}\\}$. Considering only relevant single and double
excitations in the unitary coupled-cluster operator (UCC-SD) yields the
following one-parameter unitary:

$U_{\text{UCC-
SD}}(\theta)=e^{\theta(c_{2}^{\dagger}c_{3}^{\dagger}c_{1}c_{0}-c_{0}^{\dagger}c_{1}^{\dagger}c_{3}c_{2})}$
(14)

where the minimal set of orbitals are represented by the fermionic
annihilation and creation operators $c,c^{\dagger}$ [7]. A Jordan-Wigner
transformation on these operators leads to the UCC ansatz
$\ket{\psi(\theta)}=\exp(-i\theta X_{3}X_{2}X_{1}Y_{0})\ket{0011}$, where
$\ket{0011}$ is the Hartree-Fock state. In the case of the (effective) Bravyi-
Kitaev transform (13), we obtain the simpler UCC ansatz
$\ket{\psi(\theta)}=\exp(-i\theta X_{1}Y_{0})\ket{01}$.

Since the evolution Hamiltonian commutes with the XY Hamiltonian (6), one can
use the latter ansatz and attempt to drive the Rydberg QP in the XY mode,
using $\Omega=0$ and non-zero local detunings, leaving the rest of the
parameters to be set by variational optimization:

$\displaystyle\begin{split}\ket{\psi(\delta_{0},\delta_{1},t)}&=\exp\Big{(}-it(\delta_{0}Z_{0}+\delta_{1}Z_{1}+\hat{H}_{\text{XY}})\Big{)}\ket{01}\\\
&=a(\delta_{0},\delta_{1},t)\ket{01}+b(\delta_{0},\delta_{1},t)\ket{10},\end{split}$
(15)

which coincides with the subspace reached with the UCC ansatz:

$\exp(-i\theta X_{1}Y_{0})\ket{01}=a(\theta)\ket{01}+b(\theta)\ket{10}.$ (16)

A numerical implementation of this protocol is shown in Fig. 2, where the
classical optimization was performed with a differential evolution algorithm
[35]. We observe that the ground-state energy can be obtained with an error
smaller than $5\%$ using less than $36500$ shots for each point.

Such examples of a UCC ansatz implementable with an analog approach, often
rely on finding symmetries between target and resource Hamiltonians [25].
Nevertheless, this kind of protocol remains impractical for larger molecules
because of the increasingly higher number of qubits and Pauli strings in the
Hamiltonian. In order to use the analog approach for larger encodings, we
explore other approaches below.

Figure 2: Numerical implementation of an analog VQE algorithm using a UCC
ansatz. (a) A zoom on the smallest inter-atomic distance ($0.2$ Å) shows the
evolution of the optimization with respect to the number of shots. The
differential evolution was set to perform at most 4 iterations (red squares).
The red scale shows the errorbar over 20 realizations. It takes approximately
$3.5$ hours of runtime for a QP operating at 3 Hz to achieve $\epsilon=10\%$
(light grey scale) of error and $4$ hours to achieve $\epsilon=5\%$ (dark grey
scale). (b) Relative error in percentage (red circles) between the mean VQE
result and the numerically computed lowest eigenvalue of the target
Hamiltonian (in STO-3G basis) over 10 realizations (gray crosses). The
expected error is below $\varepsilon=5\%$ (grey area). The inset depicts the
same result on an energy scale and compares it with the exact solution (green
squares). The result is drawn as a function of hydrogen inter-atomic distance.

###  III.2 Alternating pulses

We now describe an alternating operator approach, based on the QAOA algorithm
[36]. Let $|\psi_{0}\rangle$ be the state composed of all qubits in the ground
state. The whole sequence is composed by alternating constant (global) pulses,
corresponding to two non-commuting Hamiltonians $\hat{H}_{a},\hat{H}_{b}$:

$\displaystyle\hat{H}_{a}$
$\displaystyle=\frac{\hbar}{2}\sum_{i=1}^{N}\Big{(}\Omega X_{i}-\delta
Z_{i}\Big{)}+\hat{H}_{\text{inter}}$ (17) $\displaystyle\hat{H}_{b}$
$\displaystyle=\frac{\hbar}{2}\sum_{i=1}^{N}\Omega
X_{i}+\hat{H}_{\text{inter}}.$ (18)

These Hamiltonians define evolution operators $U_{a}(t)$ and $U_{b}(t)$,
during a certain time $t$ (see (9)). The ansatz of $L$ layers is written as:

$\ket{\psi(\mathbf{t}_{a},\mathbf{t}_{b})}=\prod_{\ell=1}^{L}U_{a}(t^{\ell}_{a})U_{b}(t^{\ell}_{b})\ket{\psi_{0}},$
(19)

where the arrays of parameters
$\mathbf{t}_{k}=(t^{1}_{k},\ldots,t^{\ell}_{k},\ldots,t^{L}_{k})$,
$k\in\\{a,b\\}$, fix the duration of each pulse in the layer, as described in
[27]. As another example, a different choice of parameters was used in [11],
considering a single Hamiltonian:

$\displaystyle\begin{split}\hat{H}=\frac{\hbar}{2}&\sum_{i=1}^{N}\Big{(}\Omega(t)e^{i\phi(t)}\ket{0}_{i}\bra{1}+\text{h.c.}\Big{)}+\hat{H}_{\text{inter}}\end{split}$
(20)

with different time $\mathbf{t}=(t^{1},\ldots,t^{\ell},\ldots,t^{L})$ and
phase $\boldsymbol{\phi}=(\phi^{1},\ldots,\phi^{\ell},\ldots,\phi^{L})$ arrays
defining $L$ segments of the pulse. The corresponding ansatz is then:

$\ket{\psi(\mathbf{t},\boldsymbol{\phi})}=\prod_{\ell=1}^{L}U(t^{\ell},\phi^{\ell})\ket{\psi_{0}}.$
(21)

The two approaches can be implemented in existing experimental setups,
especially when the target Hamiltonian is equal to to the resource Hamiltonian
(such as the case of the Maximal Independent Set problem with Unit Disks,
which is native to the Rydberg atoms setting). However, these methods struggle
to minimize the molecular target Hamiltonian energies within a limited number
of iterations and measurement repetitions. The alternating pulse ansatz
assumes an initial register configuration and initial guesses for the
durations of the pulses in each layer, two tasks that are the subject of
active research. The expectation is that a properly chosen register and an
optimized pulse will drive the system to a low-energy state. In Fig. 3, we
compare numerically the performance of the two alternating pulse ansatze
discussed above ([27], [11]). We also included the procedure described in Sec.
III.3, which addresses the embedding of the problem in the register and an
estimation protocol for the observables. Comparing the required number of
shots for these approaches highlights the necessity of including an efficient
estimation protocol for the observables.

Figure 3: Evolution of the ground-state energy as a function of the
accumulated number of shots for the H2 molecule at a fixed inter-nuclear
distance. We have averaged numerically 200 realizations of VQE with 3 and 8
layers using the two alternating pulse ansatz, [27] (straight line), and [11]
(dashed line). In the alternating operator approach, each Pauli string mean
value is performed with 1000 shots, which for the H2 Hamiltonian represents
$1.5\times 10^{4}$ shots before obtaining the first energy data point.
Achieving energy errors below $5\%$ (gray area) requires at least
$\mathcal{O}(10^{6})$ shots in total. For the optimized procedure (adding more
control over atom positions, pulse shaping and derandomization estimation),
the same energy error typically requires $\mathcal{O}(10^{5})$ shots.

###  III.3 Optimized Register and Iteratively Parameterized Pulses

In this section, we present a more refined approach to deal with larger
systems, aiming at exploiting the capabilities already available in Rydberg
simulators. To exemplify the procedure, we consider in the following the Ising
mode with the resource Hamiltonian (5).

####  III.3.1 Atom register and initial state

Even though we only consider global pulses for the processing stage, there
still remains freedom in the choice of the positions of the atoms. This
determines the strength of pairwise interactions and defines a connectivity
graph whose edges correspond to the atoms that experience a blockade effect
[13] (a different graph structure can be defined for the XY mode (6)).

In order to find suitable atomic positions, the coordinates are optimized in
the plane so that the associated interaction energy matrix resembles as much
as possible the information contained in the target Hamiltonian. Since the
latter contains general Pauli strings, we consider a subset of terms whose
coefficients can be expressed in terms of the coordinates of the atoms222A
broader series of techniques for embedding the problem information into the
atom register has been considered in [37, 38]. A simple choice consists in
selecting the terms that can be directly compared with the Ising-like
interaction of the atoms: Let the matrix $V^{\text{T}}$ be given by the
positive coefficients of the terms with only two $Z$ operators in the
molecular target Hamiltonian and $V^{\text{R}}$ (our “register” matrix) the
resulting values of interaction strength $C_{6}/r_{ij}^{6}$ for each pair
$i,j$ of atom positions in the register. This defines a score function
$\sum_{i,j}(V^{\text{T}}_{ij}-V^{\text{R}}_{ij})^{2}$ that we minimize
numerically by varying the atom coordinates.

The set of atomic positions that arises from this minimization will be our
optimized register. Its geometry will be used to simulate the target
Hamiltonian, but has no intrinsic chemical meaning. The information that is
taken from the Hamiltonian can be chosen from other subsets of the Pauli
strings (e.g. terms with 3 or more $Z$ operators) and different
interpretations of how the coefficients constitute a register matrix. A
different resource Hamiltonian, such as one with XY interactions, would imply
a different choice of subset. In Fig. 4, we summarize graphically the
procedure for the case of the $\text{H}_{\text{2}}$ molecule with the Jordan-
Wigner transformation. It turns our that we obtain at a geometry very similar
to the one heuristically picked for the alternating pulse ansatz.

Figure 4: Protocol to optimize the positions of the atoms in a register based
on a target Hamiltonian. Blue: we begin with a register of randomly placed
atoms and all the Ising interaction terms are entered in a $N\times N$ matrix.
Yellow: all positive coefficients before the Pauli strings with only two $Z$
operators in the target Hamiltonian are combined in another (target) matrix
$N\times N$. The coordinates of the atoms are optimized to minimize the
distance between the two matrices. We then obtain a new register on which we
will apply the VQE sequence.

####  III.3.2 Optimization of the parameterized pulse sequence

We constructed a variation of the so-called ctrl-VQE protocol [18] for the
case of a global pulse on the register, in which the number of parameters
increases at every optimization iteration, while the total time
$t_{\text{tot}}$ remains fixed:

Consider a set of Rabi frequencies $\\{\Omega_{i}\\}_{i=1}^{K}$ and detunings
$\\{\delta_{i}\\}_{i=1}^{K}$ defined discretely over a set of time labels
$0<t_{1}<\ldots<t_{K}=t_{\text{tot}}$. Then, at iteration $k$, a new time
label $0<t_{k}<t_{\text{tot}}$ is generated at random, lying between two
previous time labels, $t_{i-1}<t_{k}<t_{i}$. To avoid labels too close to each
other, we will accept $t_{k}$ if the intervals $|t_{i-1}-t_{k}|,|t_{k}-t_{i}|$
are large enough compared to the response time of the waveform generator of
the machine (in the order of a few ns). The corresponding Rabi frequency
$\Omega_{i}$ and detuning $\delta_{i}$ from the parent interval
$[t_{i-1},t_{i}]$ are then split into two independent parameters
$\Omega_{i}^{\prime},\Omega^{\prime}_{k}$ and
$\delta^{\prime}_{i},\delta^{\prime}_{k}$ whose initial values are set equal
to their parent parameters (see Fig. 5). Finally, the new set of parameters
$\\{\Omega_{1},\ldots,\Omega^{\prime}_{i},\Omega^{\prime}_{k},\ldots,\Omega_{K}\\}$
(likewise for $\\{\delta_{i}\\}_{i=1}^{K})$ is optimized starting from the
previous iteration values. This algorithm acts therefore as a pulse shaping
process. From time $t_{i-1}$ to $t_{i}$ the acting Hamiltonian is:

$\hat{H}_{i}=\frac{\hbar}{2}\Big{(}\Omega_{i}\sum_{j=1}^{N}X_{j}-\delta_{i}\sum_{j=1}^{N}Z_{j}\Big{)}+\hat{H}_{\text{inter}}$
(22)

and our ansatz, for $K$ iterations, becomes:

$\ket{\psi(\boldsymbol{\Omega},\boldsymbol{\delta})}=\mathcal{T}\prod_{i=1}^{K}\exp\Big{[}-i\int_{t_{i}}^{t_{i+1}}\hat{H}_{i}(\tau)d\tau\Big{]}\ket{\psi_{0}}$
(23)

Note that while the interval involves a constant Hamiltonian, we include a
time-dependent integration at each interval, to indicate that the waveforms
that compose the pulse can be adapted to hardware conditions (e.g. by being
interpolated, or by adapting the shape with an envelope function).

Figure 5: Iterative splitting and optimization of the pulse parameters: (a)
Choose at random a time $t_{k}$, which will fall in the interval
$(t_{i-1},t_{i})$, and accept it if $|t_{i-1}-t_{k}|$ and $|t_{k}-t_{i}|$
conform to the device response time. Split the corresponding $\Omega_{i}$ into
two parameters $\Omega^{\prime}_{i},\Omega^{\prime}_{k}$ with initial value
equal to $\Omega_{i}$. Do the same to the set $\\{\delta_{i}\\}$. (b) Optimize
the new set of parameters to lower the energy of the prepared state (10).

####  III.3.3 Energy estimation by derandomization

In our algorithm implementation, we take as a figure of performance of the
variational optimization run the _total_ number of shots required to achieve a
given error threshold $\varepsilon$ for the energy. A bounded number of
processing cycles is required to remain within a realistic time lapse for the
entire implementation process. Rather than measuring several times each of the
Pauli observables in the Hamiltonian, we use an estimation protocol
(derandomization [23]) based on fixing local Pauli measurements from an
originally random set. This allows to efficiently predict the energy of the
prepared state $\ket{\psi}$, $\langle\hat{H}_{\text{T}}(\psi)\rangle$, at each
loop of the optimization of the parameters.

More specifically, the derandomization algorithm starts with an initial
measurement set of $M$ random Pauli strings $\\{\hat{P}^{(m)}\\}_{m=1}^{M}$. A
greedy algorithm improves the overall expected performance of the measurement
set, effectively “derandomizing” the operators of each random Pauli string in
sequence. The improvement is quantified by the average of the _confidence
bound_ , which ensures that the empirical average 333 A Pauli string $A$
_hits_ $B$, if by changing some operators in $A$ to $\mathbb{1}$, we form $B$
(for example $ZX\mathbb{1}$ hits $\mathbb{1}X\mathbb{1}$ and
$Z\mathbb{1}\mathbb{1}$). The empirical average is obtained using those Pauli
measurement basis $\\{\hat{P}^{(m)}\\}$ that hit an observable
$\hat{P}^{(s)}$, with the relevant measured bits expressed as $\pm$:
$\omega_{s}=\frac{1}{N_{\mathrm{h}}}\sum^{M}_{\begin{subarray}{c}m:\\\
\hat{P}^{(m)}\text{hits}\hat{P}^{(s)}\end{subarray}}\Bigg{(}\prod^{N}_{\begin{subarray}{c}j:\hat{P}^{(s)}_{j}\neq\mathbb{1}\end{subarray}}\mathbf{b}^{(m)}_{j}\Bigg{)},$
where $N_{\text{h}}$ counts how many Pauli strings in the set
$\\{\hat{P}^{(m)}\\}$ hit $\hat{P}^{(s)}$, and
$\mathbf{b}^{(m)}=\mathbf{b}^{(m)}_{1}\cdots\mathbf{b}^{(m)}_{N}$ is the
bitstring measured with the basis $\hat{P}^{(m)}$ $\omega_{s}$ corresponding
to the $s$-th term of $\hat{H}_{\text{T}}$ is within a desired accuracy
$|\omega_{s}-\langle\hat{P}^{(s)}\rangle|/|\langle\hat{P}^{(s)}\rangle|<\epsilon$
and with a high probability. The total energy is finally estimated as
$\langle\hat{H}_{\text{T}}(\psi)\rangle\approx\sum_{s=1}^{S}\omega_{s}$.

While the pulses that prepare the state are global, the measurement itself
requires the implementation of local rotations on the qubits. This can be
achieved experimentally by using a toolbox such as the one described in [39],
thus emphasizing the digital-analog interplay that is now within reach for
next-generation neutral atom devices.

##  IV Numerical Results

###  IV.1 Application on LiH and BeH2 molecules

Figure 6: Numerical results of the VQE algorithm with our digital-analog
protocol (a) for BeH2 at an intermolecular distance of $0.5$ Å where we
increased the number of shots beyond 1 day of experiment and observed the
expected improvement over several days of calculations. The light grey shade
and the dark grey shade indicate respectively $\varepsilon=10\%$ and
$\varepsilon=5\%$ error benchmark. The red line shows the improvement mean
value over 100 run. (b) Result of BeH2 molecule with an encoding of 6 qubits
and 165 Pauli strings and (c) result of LiH molecule result with an encoding
of 6 qubits and 118 Pauli strings. The insets show the register geometry at
specific inter-nuclear distances. For the case of LiH a single heuristic
choice was used, while for BeH2 an optimized geometry was prepared at each
inter-nuclear distance, minimizing the distance between selected terms of the
target Hamiltonian and the interaction energies of the atoms in the register
(see Sec. III.3.1). Blue squares: mean value for several simulations. Green
line: result from exact diagonalization. The gray shade indicates an
$\varepsilon=5\%$ error benchmark. The total number of shots for each
optimization result (red crosses) is set to 350000, corresponding roughly to a
day of processing in a Rydberg QP.

We have applied the method described in section III.3 to the LiH and BeH2
molecules. Using the Qiskit [40] framework combined with Pyquante [41], we
calculate the one and two-body integrals of (3), encoding the problem into 6
qubits using the Bravyi-Kitaev method. The Hamiltonians contain 118 and 165
Pauli strings respectively. To design the pulse sequence and include realistic
device constraints into the simulations we used the open source package Pulser
[42]. The Powell algorithm [43] was used for the classical optimization of the
pulse values with 20 function evaluations for each iteration. The two initial
Rabi frequency and detuning are chosen randomly in the interval $[0,2\times
2\pi]$MHz for each optimization procedure. During the optimization, Rabi
frequencies are bounded to this interval to remain within experimentally
accessible values [10], while the interval for the detuning was taken as
$[-2\times 2\pi,2\times 2\pi]$MHz. We ran the algorithm five times for four
different inter-nuclear distances $R$ yielding the results shown in Fig. 6.
The algorithm converges with small errors in most cases, but we notice the
impact of the initial parameters on the obtained energies. For instance at
$R=1.5$ Å for the BeH2 molecule, the obtained energy values are up to $0.4$
Hartree apart.

To optimize the register geometry, we took the coordinates as parameters,
starting from random positions and minimized the score function described in
Sec. III.3.1. We optimized the atom register based on the Nelder-Mead method
[44], with a few thousand function evaluations, and we also compared to
several heuristic choices obtained by a term-by-term comparison with the
interaction matrix (a well-performing choice of positions is shown for LiH in
Fig. 6)

We define our error as

$\varepsilon=|{E_{\text{exact}}-E_{\text{estimated}}|/|E_{\text{exact}}}|,$
(24)

where $E_{\text{exact}}=\langle\hat{H}_{\text{T}}\rangle$ is the exact
diagonalization solution with respect to the target Hamiltonian and
$E_{\text{estimated}}$ is the energy calculated with the optimized geometry,
the optimized pulse sequence, and the derandomization estimation. The
optimized configurations, together with optimized pulse parameters and
energies estimated at each iteration with derandomized measurements give rise
to energy errors typically below the $\varepsilon=5\%$ threshold in less than
350000 shots.

The implementation of the derandomization algorithm allows us to choose the
number of measurements that we wish to take (our budget), for a given target
accuracy of estimation, which we set to correspond to our $\varepsilon=5\%$
benchmark. The resulting “derandomized” Pauli measurements included typically
close to 20 different Pauli strings, calculated from the minimization of the
average confidence bound that ensures an empirical average within the chosen
$\varepsilon$. Since some derandomized Pauli strings have more operators in
common with the terms in the target Hamiltonian (they “hit” more target
observables), we adjusted the measurement repetitions to be spent
proportionally more in them, which improved statistics. We also verified that
the obtained accuracy improves upon increasing the allowed number of shots,
although we don’t expect a full convergence, given the incomplete information
used to define $\hat{H}_{\text{T}}$.

###  IV.2 Roadmap for more complex molecules

We discuss in this section some observations about the presented protocol for
larger molecules and more complex basis sets, where the number of terms in the
Hamiltonian and the required qubits to encode it grows quickly. Currently
available neutral-atom devices can load hundreds of traps [45], but the
available space on the register will eventually become a resource limitation.
In Fig. 7 we show the embedding results for H2O and CH4 in different basis
sets, where thousands of terms would need to be measured. Our simple
restriction to $Z$-terms captures limited features of the Hamiltonian, mostly
concentrating atoms where the largest values need to be reproduced. In fact,
as the system size grows, we observe that most of the atoms in the register
act as a “background” for these clusters. Choosing different terms from the
Target Hamiltonian will bring forward other features, highlighting the
opportunities of using learning methods to find more performing atom
positions. We have not addressed here the possibilities offered by three-
dimensional registers [46], which allow for more complex embeddings and have
been already studied for graph-combinatorial problems [47], although they can
be straightforwardly included in the protocol.

Figure 7: Example embeddings for the H2O and CH4 molecules. Using a Jordan-
Wigner encoding, for two different basis, our H2O Hamiltonian (left column)
consisted of 14 qubits and 2110 terms (STO-3G basis), of which 595 were chosen
as features for the reference Hamiltonian. Choosing a 6-31G* basis gives 36
qubits, with 83003 terms and 5594 of them relevant. For CH4 (right column),
the STO-3G basis gives 18 atoms and 6892 Pauli terms, of which 1359 terms were
used for the embedding. The 6-31G* basis requires 44 qubits, 297075 terms and
we have picked an embedding with 11772 terms.

Two aspects of the optimization that rely on classical computation can be
further refined: the choice of initial state and the selection of parameters:
initializing the optimization with a product state from a Hartree-Fock state
approximation can help exploring a lower energy set of output states. In Fig.
8 we have scanned through all product states of an 8-qubit system to select
those who benefit the most from the first step of the pulse-optimization
protocol presented above. The best choice of initial product state can then be
used for the Rydberg QP implementation444Preparing the initial product state
requires for example masking atoms with the help of a spatial light
modulator.. Comparing the energy
$\langle\psi_{0}|U(\theta)^{\dagger}\hat{H}_{\text{T}}U(\theta)|\psi_{0}\rangle$
of a candidate initial product state $|\psi_{0}\rangle$ evolving under a
constant pulse parameterized by $\theta$ can be performed for example using
tensor-network techniques over an HPC backend [48, 49]. On the other hand, the
number of parameters can be chosen at will, and do not depend on the number of
qubits. More advanced control techniques can be applied here, and there is
large choice of techniques and numerical tools in the subject of quantum
optimal control. We recall that the optimization of the pulse sequence is not
dependent on the embedding itself – it is a global property of the system
evolution, where nearby atoms constitute blockade regions that characterize
the final state.

Figure 8: Relative error $\varepsilon$ after the first step of the pulse
optimization III.3.2, using different initial product states. While the
default $\ket{\psi_{0}}=\ket{0}^{\otimes N}$ (black line) barely improves its
$\varepsilon$ by the optimization step, other product states (blue, dashed
line and red, dashed-dotted line) achieve a very low error. We have used an
10-qubit Jordan-Wigner encoding of the LiH Molecular Hamiltonian under the
STO-3G basis (276 terms), and averaged several instances for each product
state.

##  V Discussion

In this work, we have numerically studied a digital-analog quantum algorithm
in the context of quantum chemistry using an ideal Rydberg quantum processor
as a hardware. We have considered small molecules with resource Hamiltonians
of two qubits for the H2 molecule and six qubits for the LiH and BeH2
molecules to demonstrate the applicability of our methods. Our purpose was to
describe the construction of such an algorithm, discussing the cost of each
stage in terms of the number of measurement repetitions. Our numerical results
should be viewed as a first benchmark and should trigger further explorations.
Besides, it provides a roadmap for the improvement of Rydberg quantum
processor, in particular in terms of cycle time.

By considering the symmetries of H2 Hamiltonian, we show how the UCC method
efficiently and accurately approximates the ground state energy. However,
finding a two-body Hamiltonian which commutes with a more general molecular
target Hamiltonian is a hard problem. We therefore proposed another protocol
for larger molecules: we optimized the geometries of the atomic array, pulse
sequences and included an estimation method (derandomization) for the energy
measurement. We targeted $5\%$ of accuracy compared to the exact
diagonalization method for Hamiltonian with 6 qubits and more than a hundred
Pauli strings.

We observed that the geometry of the array has a significant impact on the
result: In the case of LiH, where the target matrix $V^{T}$ does not provide
much information due to the few terms with only two $Z$ operators, the
optimized positions underperform with respect to a careful choice of
positions, although scaling the heuristics that gives rise to such a geometry
for molecules with large number of qubits in their encoding is impractical.
This calls for the design of more advanced embedding algorithms. Indeed, for
the case of BeH2, the register optimization achieves rather small energy
errors, especially for larger distances.

Previous studies [14] have quantified the demanding resource requirements for
practical VQE applications. After several iterations of pulse optimization
with energy estimation via derandomization, the error on the average energy of
the final prepared state descends to $\varepsilon<5\%$, and expected to be
obtained within a day of measurement in a typical current-day Rydberg QP.

In the numerical implementation we have used out-of-the-box optimizers with
limitations for the numerical task at hand. Other possibilities include the
use of an interpolated waveform for each set of parameters and shaping the
pulse using bayesian optimization routines as explored in [38] for the study
of combinatorial graph problems. Note that experimentally, one could clone
several times the atom layout in spatially separated regions of the register
(at least for a small number of qubits), multiplying the obtained number of
bitstrings. To achieve close to $1\%$ relative error, we expect that at least
a week of Rydberg QP runtime would be necessary (see extrapolation shown in
Fig. 6). The capacity of a circuit ansatz to construct a desired quantum state
while keeping a small depth and number of parameters is studied by its
_expressibility_ [50, 51]. In the case of analog systems, this is an emerging
topic of research [52], with the goal of ensuring that a given ansatz could
potentially lead to a good approximation of the ground state and achieve
chemical accuracy.

The impact of experimental errors in a real-life implementation will also lead
to performance reductions. SPAM (State Preparation And Measurement) errors are
typically the largest source of discrepancy for the neutral atom devices [53],
but the energy errors observed in numerical simulations remain low as long as
the failure rates are small, given that the variational nature of the
algorithm shows robustness to several types of errors [54]. Recently [55], VQE
was experimentally implemented in a superconducting quantum processor for H2,
LiH and F2 with 4, 6 and 12 qubits, respectively, using the UCC ansatz and a
different flavor of measurement protocol [56]. The readouts went through an
error mitigation post-processing routine, which showed that these techniques
can greatly compensate the noise effects from their quantum processor, with an
error reduction of up to two orders of magnitude and leading to chemical
accuracy in some circumstances. We expect that such mitigation can be added to
the protocols considered in this paper and will help experimental
implementations in Rydberg QP.

To tackle quantum simulation algorithms for energy estimation, more
developments in both quantum and classical parts of the hybrid algorithm are
needed. Reaching chemical accuracy for molecules with a few tens of qubits
remains an open challenge that can now begin to be explored in experimental
devices. This will provide evidence to generate new and fundamental insights
to understand under what conditions a computational advantage can be achieved.

###### Acknowledgements.

We thank Thierry Lahaye and Louis Vignoli for discussions and reading of the
manuscript. This work was supported by the European Union’s Horizon 2020
research and innovation program under grant agreement No. 817482 (PASQuanS),
and the European Research Council (Advanced grant No. 101018511-ATARAXIA). It
was also supported by EDF R&D, the Research and Development Division of
Electricité de France under the ANRT contract N°2020/0011.

##  Appendix A LiH and BeH${}_{\text{2}}$ Hamiltonians

In this section, examples of complete Hamiltonians of molecules LiH (for an
inter-atomic distance of $1.5$ $\mathrm{\mathring{A}}$) and BeH${}_{\text{2}}$
(for an inter-atomic distance of $1.17$ $\mathrm{\mathring{A}}$) obtained with
the method described in Sec.II.1 are shown.

$\begin{split}H_{\mathrm{LiH}}=&{\bf-0.19975}+{\bf
0.05393}Z_{0}{\bf-0.12836}Z_{1}{\bf-0.31773}Z_{0}Z_{1}{\bf-0.31773}Z_{3}+{\bf
0.0605}Z_{1}Z_{3}\\\ &+{\bf 0.11409}Z_{0}Z_{1}Z_{3}+{\bf 0.05362}Z_{4}+{\bf
0.11434}Z_{2}Z_{4}{\bf-0.03787}Z_{2}Z_{3}Z_{4}+{\bf
0.05362}Z_{1}Z_{2}Z_{3}Z_{4}\\\ &+{\bf
0.0836}Z_{0}Z_{1}Z_{2}Z_{3}Z_{4}{\bf-0.03787}Z_{5}+{\bf
0.05666}Z_{1}Z_{5}+{\bf 0.11434}Z_{0}Z_{1}Z_{5}+{\bf 0.0836}Z_{3}Z_{5}\\\
&+{\bf 0.05666}Z_{2}Z_{3}Z_{5}{\bf-0.12836}Z_{4}Z_{5}+{\bf
0.0847}Z_{1}Z_{4}Z_{5}+{\bf 0.0605}Z_{0}Z_{1}Z_{4}Z_{5}+{\bf
0.05393}Z_{3}Z_{4}Z_{5}\\\ &+{\bf 0.12357}Z_{2}Z_{3}Z_{4}Z_{5}+{\bf
0.01522}X_{1}{\bf-0.01522}Z_{0}X_{1}+{\bf
0.01089}X_{1}Z_{3}{\bf-0.01089}Z_{0}X_{1}Z_{3}\\\ &+{\bf
0.00436}X_{1}Z_{2}Z_{3}Z_{4}{\bf-0.00436}Z_{0}X_{1}Z_{2}Z_{3}Z_{4}+{\bf
0.01273}X_{1}Z_{5}{\bf-0.01273}Z_{0}X_{1}Z_{5}{\bf-0.00901}X_{1}Z_{4}Z_{5}\\\
&+{\bf 0.00901}Z_{0}X_{1}Z_{4}Z_{5}+{\bf
0.00448}X_{0}X_{2}{\bf-0.00479}X_{0}Z_{1}X_{2}{\bf-0.03512}X_{0}Z_{1}X_{2}Z_{3}{\bf-0.03512}Y_{0}Y_{2}Z_{4}\\\
&{\bf-0.00479}Y_{0}Y_{2}Z_{3}Z_{4}+{\bf
0.00448}Y_{0}Z_{1}Y_{2}Z_{3}Z_{4}{\bf-0.03306}X_{0}Z_{1}X_{2}Z_{5}+{\bf
0.00237}Y_{0}Y_{2}Z_{3}Z_{5}+{\bf 0.00237}X_{0}Z_{1}X_{2}Z_{4}Z_{5}\\\
&{\bf-0.03306}Y_{0}Y_{2}Z_{3}Z_{4}Z_{5}-{\bf(4\times
10^{-5})}X_{0}X_{1}X_{2}{\bf-0.00277}Y_{0}Y_{1}X_{2}+{\bf
0.01054}X_{0}X_{1}X_{2}Z_{3}+{\bf 0.01054}X_{0}Y_{1}Y_{2}Z_{4}\\\ &+{\bf
0.00277}Y_{0}X_{1}Y_{2}Z_{3}Z_{4}-{\bf(4\times
10^{-5})}X_{0}Y_{1}Y_{2}Z_{3}Z_{4}+{\bf
0.01173}X_{0}X_{1}X_{2}Z_{5}{\bf-0.00154}X_{0}Y_{1}Y_{2}Z_{3}Z_{5}\\\
&{\bf-0.00154}X_{0}X_{1}X_{2}Z_{4}Z_{5}+{\bf
0.01173}X_{0}Y_{1}Y_{2}Z_{3}Z_{4}Z_{5}+{\bf
0.01522}X_{3}X_{4}{\bf-0.00901}Z_{1}X_{3}X_{4}+{\bf
0.01089}Z_{0}Z_{1}X_{3}X_{4}\\\ &+{\bf
0.00436}Y_{3}Y_{4}{\bf-0.01273}Z_{2}Y_{3}Y_{4}+{\bf
0.00436}X_{3}X_{4}Z_{5}{\bf-0.01273}Z_{2}X_{3}X_{4}Z_{5}+{\bf
0.01522}Y_{3}Y_{4}Z_{5}{\bf-0.00901}Z_{1}Y_{3}Y_{4}Z_{5}\\\ &+{\bf
0.01089}Z_{0}Z_{1}Y_{3}Y_{4}Z_{5}+{\bf
0.00658}X_{1}X_{3}X_{4}{\bf-0.00658}Z_{0}X_{1}X_{3}X_{4}+{\bf
0.00658}X_{1}Y_{3}Y_{4}Z_{5}{\bf-0.00658}Z_{0}X_{1}Y_{3}Y_{4}Z_{5}\\\
&{\bf-0.00776}X_{0}Z_{1}X_{2}X_{3}X_{4}+{\bf 0.00776}Y_{0}Y_{2}Y_{3}Y_{4}+{\bf
0.00776}Y_{0}Y_{2}X_{3}X_{4}Z_{5}{\bf-0.00776}X_{0}Z_{1}X_{2}Y_{3}Y_{4}Z_{5}\\\
&+{\bf
0.00211}X_{0}X_{1}X_{2}X_{3}X_{4}{\bf-0.00211}X_{0}Y_{1}Y_{2}Y_{3}Y_{4}{\bf-0.00211}X_{0}Y_{1}Y_{2}X_{3}X_{4}Z_{5}+{\bf
0.00211}X_{0}X_{1}X_{2}Y_{3}Y_{4}Z_{5}+{\bf 0.00004}X_{5}\\\
&{\bf-0.00154}Z_{1}X_{5}+{\bf 0.01054}Z_{0}Z_{1}X_{5}+{\bf
0.00277}Z_{3}X_{5}{\bf-0.01173}Z_{2}Z_{3}X_{5}+{\bf 0.00004}Z_{4}X_{5}+{\bf
0.00154}Z_{1}Z_{4}X_{5}\\\
&{\bf-0.01054}Z_{0}Z_{1}Z_{4}X_{5}{\bf-0.00277}Z_{3}Z_{4}X_{5}+{\bf
0.01173}Z_{2}Z_{3}Z_{4}X_{5}+{\bf
0.00211}X_{1}X_{5}{\bf-0.00211}Z_{0}X_{1}X_{5}\\\
&{\bf-0.00211}X_{1}Z_{4}X_{5}+{\bf
0.00211}Z_{0}X_{1}Z_{4}X_{5}{\bf-0.00837}X_{0}Z_{1}X_{2}X_{5}+{\bf
0.00837}Y_{0}Y_{2}Z_{3}X_{5}+{\bf 0.00837}X_{0}Z_{1}X_{2}Z_{4}X_{5}\\\
&{\bf-0.00837}Y_{0}Y_{2}Z_{3}Z_{4}X_{5}+{\bf
0.00303}X_{0}X_{1}X_{2}X_{5}{\bf-0.00303}X_{0}Y_{1}Y_{2}Z_{3}X_{5}{\bf-0.00303}X_{0}X_{1}X_{2}Z_{4}X_{5}\\\
&+{\bf 0.00303}X_{0}Y_{1}Y_{2}Z_{3}Z_{4}X_{5}+{\bf
0.00448}X_{3}X_{4}X_{5}+{\bf
0.03306}Z_{2}X_{3}X_{4}X_{5}{\bf-0.00479}Y_{3}Y_{4}X_{5}+{\bf
0.00237}Z_{1}Y_{3}Y_{4}X_{5}\\\ &{\bf-0.03512}Z_{0}Z_{1}Y_{3}Y_{4}X_{5}+{\bf
0.00448}Y_{3}X_{4}Y_{5}+{\bf 0.03306}Z_{2}Y_{3}X_{4}Y_{5}+{\bf
0.00479}X_{3}Y_{4}Y_{5}{\bf-0.00237}Z_{1}X_{3}Y_{4}Y_{5}\\\ &+{\bf
0.03512}Z_{0}Z_{1}X_{3}Y_{4}Y_{5}{\bf-0.00776}X_{1}Y_{3}Y_{4}X_{5}+{\bf
0.00776}Z_{0}X_{1}Y_{3}Y_{4}X_{5}+{\bf
0.00776}X_{1}X_{3}Y_{4}Y_{5}{\bf-0.00776}Z_{0}X_{1}X_{3}Y_{4}Y_{5}\\\
&{\bf-0.03074}Y_{0}Y_{2}X_{3}X_{4}X_{5}+{\bf
0.03074}X_{0}Z_{1}X_{2}Y_{3}Y_{4}X_{5}{\bf-0.03074}Y_{0}Y_{2}Y_{3}X_{4}Y_{5}{\bf-0.03074}X_{0}Z_{1}X_{2}X_{3}Y_{4}Y_{5}\\\
&+{\bf
0.00837}X_{0}Y_{1}Y_{2}X_{3}X_{4}X_{5}{\bf-0.00837}X_{0}X_{1}X_{2}Y_{3}Y_{4}X_{5}+{\bf
0.00837}X_{0}Y_{1}Y_{2}Y_{3}X_{4}Y_{5}+{\bf
0.00837}X_{0}X_{1}X_{2}X_{3}Y_{4}Y_{5}\end{split}$ (25)
$\begin{split}H_{\text{BeH}_{\text{2}}}=&{\bf-1.90305}{\bf-0.48894}Z_{0}+{\bf
0.14357}Z_{1}{\bf-0.18803}Z_{0}Z_{1}+{\bf 0.12314}Z_{2}+{\bf
0.18326}Z_{0}Z_{2}\\\ &+{\bf 0.10964}Z_{1}Z_{2}+{\bf
0.18222}Z_{0}Z_{1}Z_{2}{\bf-0.48894}Z_{3}+{\bf 0.1288}Z_{0}Z_{3}+{\bf
0.1136}Z_{0}Z_{1}Z_{3}\\\ &+{\bf 0.11249}Z_{2}Z_{3}+{\bf
0.11746}Z_{1}Z_{2}Z_{3}+{\bf 0.14357}Z_{4}{\bf-0.18803}Z_{3}Z_{4}+{\bf
0.1136}Z_{0}Z_{3}Z_{4}\\\ &+{\bf 0.10602}Z_{0}Z_{1}Z_{3}Z_{4}+{\bf
0.10306}Z_{2}Z_{3}Z_{4}+{\bf 0.10577}Z_{1}Z_{2}Z_{3}Z_{4}+{\bf
0.12314}Z_{5}+{\bf 0.11249}Z_{0}Z_{5}\\\ &+{\bf 0.10306}Z_{0}Z_{1}Z_{5}+{\bf
0.10451}Z_{2}Z_{5}+{\bf 0.10785}Z_{1}Z_{2}Z_{5}+{\bf 0.18326}Z_{3}Z_{5}+{\bf
0.10964}Z_{4}Z_{5}\\\ &+{\bf 0.11746}Z_{0}Z_{4}Z_{5}+{\bf
0.10577}Z_{0}Z_{1}Z_{4}Z_{5}+{\bf 0.10785}Z_{2}Z_{4}Z_{5}+{\bf
0.11352}Z_{1}Z_{2}Z_{4}Z_{5}+{\bf 0.18222}Z_{3}Z_{4}Z_{5}\\\
&{\bf-0.00743}X_{0}X_{1}{\bf-0.00229}Y_{0}Y_{1}{\bf-0.00229}X_{0}X_{1}Z_{2}{\bf-0.00743}Y_{0}Y_{1}Z_{2}{\bf-0.00711}Y_{0}Y_{1}Z_{3}\\\
&{\bf-0.00711}X_{0}X_{1}Z_{2}Z_{3}{\bf-0.00875}Y_{0}Y_{1}Z_{3}Z_{4}{\bf-0.00875}X_{0}X_{1}Z_{2}Z_{3}Z_{4}{\bf-0.00352}Y_{0}Y_{1}Z_{5}{\bf-0.00352}X_{0}X_{1}Z_{2}Z_{5}\\\
&{\bf-0.0072}Y_{0}Y_{1}Z_{4}Z_{5}{\bf-0.0072}X_{0}X_{1}Z_{2}Z_{4}Z_{5}{\bf-0.04165}X_{0}X_{2}+{\bf
0.03769}X_{0}Z_{1}X_{2}{\bf-0.00396}Y_{0}Y_{2}\\\ &+{\bf
0.00839}X_{1}X_{2}{\bf-0.01015}Z_{0}X_{1}X_{2}{\bf-0.01015}Y_{1}Y_{2}+{\bf
0.00839}Z_{0}Y_{1}Y_{2}+{\bf 0.01355}Z_{0}X_{1}X_{2}Z_{3}\\\ &+{\bf
0.01355}Y_{1}Y_{2}Z_{3}+{\bf 0.01082}Z_{0}X_{1}X_{2}Z_{3}Z_{4}+{\bf
0.01082}Y_{1}Y_{2}Z_{3}Z_{4}+{\bf 0.00854}Z_{0}X_{1}X_{2}Z_{5}+{\bf
0.00854}Y_{1}Y_{2}Z_{5}\\\ &+{\bf 0.01408}Z_{0}X_{1}X_{2}Z_{4}Z_{5}+{\bf
0.01408}Y_{1}Y_{2}Z_{4}Z_{5}+{\bf
0.03611}X_{0}X_{3}{\bf-0.03611}X_{0}Z_{1}X_{3}{\bf-0.03611}X_{0}X_{3}Z_{4}\\\
&+{\bf 0.03611}X_{0}Z_{1}X_{3}Z_{4}{\bf-0.02498}X_{1}X_{3}+{\bf
0.02498}Z_{0}X_{1}Z_{2}X_{3}+{\bf
0.02498}X_{1}X_{3}Z_{4}{\bf-0.02498}Z_{0}X_{1}Z_{2}X_{3}Z_{4}\\\
&{\bf-0.03615}X_{2}X_{3}+{\bf 0.03615}Z_{1}X_{2}X_{3}+{\bf
0.03615}X_{2}X_{3}Z_{4}{\bf-0.03615}Z_{1}X_{2}X_{3}Z_{4}{\bf-0.01573}X_{0}X_{1}X_{2}X_{3}\\\
&{\bf-0.01573}Y_{0}X_{1}Y_{2}X_{3}+{\bf 0.01573}X_{0}X_{1}X_{2}X_{3}Z_{4}+{\bf
0.01573}Y_{0}X_{1}Y_{2}X_{3}Z_{4}{\bf-0.02498}X_{0}X_{4}+{\bf
0.02498}X_{0}Z_{1}X_{4}\\\ &+{\bf
0.02498}X_{0}Z_{3}X_{4}Z_{5}{\bf-0.02498}X_{0}Z_{1}Z_{3}X_{4}Z_{5}+{\bf
0.02085}X_{1}X_{4}{\bf-0.02085}Z_{0}X_{1}Z_{2}X_{4}{\bf-0.02085}X_{1}Z_{3}X_{4}Z_{5}\\\
&+{\bf 0.02085}Z_{0}X_{1}Z_{2}Z_{3}X_{4}Z_{5}+{\bf
0.02464}X_{2}X_{4}{\bf-0.02464}Z_{1}X_{2}X_{4}{\bf-0.02464}X_{2}Z_{3}X_{4}Z_{5}+{\bf
0.02464}Z_{1}X_{2}Z_{3}X_{4}Z_{5}\\\ &+{\bf 0.01532}X_{0}X_{1}X_{2}X_{4}+{\bf
0.01532}Y_{0}X_{1}Y_{2}X_{4}{\bf-0.01532}X_{0}X_{1}X_{2}Z_{3}X_{4}Z_{5}{\bf-0.01532}Y_{0}X_{1}Y_{2}Z_{3}X_{4}Z_{5}{\bf-0.00743}X_{3}X_{4}\\\
&{\bf-0.00229}Y_{3}Y_{4}{\bf-0.00711}Z_{0}Y_{3}Y_{4}{\bf-0.00875}Z_{0}Z_{1}Y_{3}Y_{4}{\bf-0.00352}Z_{2}Y_{3}Y_{4}{\bf-0.0072}Z_{1}Z_{2}Y_{3}Y_{4}\\\
&{\bf-0.00229}X_{3}X_{4}Z_{5}{\bf-0.00711}Z_{0}X_{3}X_{4}Z_{5}{\bf-0.00875}Z_{0}Z_{1}X_{3}X_{4}Z_{5}{\bf-0.00352}Z_{2}X_{3}X_{4}Z_{5}{\bf-0.0072}Z_{1}Z_{2}X_{3}X_{4}Z_{5}\\\
&{\bf-0.00743}Y_{3}Y_{4}Z_{5}+{\bf 0.01972}Y_{0}Y_{1}Y_{3}Y_{4}+{\bf
0.01972}X_{0}X_{1}Z_{2}Y_{3}Y_{4}+{\bf 0.01972}Y_{0}Y_{1}X_{3}X_{4}Z_{5}+{\bf
0.01972}X_{0}X_{1}Z_{2}X_{3}X_{4}Z_{5}\\\
&{\bf-0.0173}Z_{0}X_{1}X_{2}Y_{3}Y_{4}{\bf-0.0173}Y_{1}Y_{2}Y_{3}Y_{4}{\bf-0.0173}Z_{0}X_{1}X_{2}X_{3}X_{4}Z_{5}{\bf-0.0173}Y_{1}Y_{2}X_{3}X_{4}Z_{5}{\bf-0.03615}X_{0}X_{5}\\\
&+{\bf 0.03615}X_{0}Z_{1}X_{5}+{\bf
0.03615}X_{0}Z_{4}X_{5}{\bf-0.03615}X_{0}Z_{1}Z_{4}X_{5}+{\bf
0.02464}X_{1}X_{5}{\bf-0.02464}Z_{0}X_{1}Z_{2}X_{5}\\\
&{\bf-0.02464}X_{1}Z_{4}X_{5}+{\bf 0.02464}Z_{0}X_{1}Z_{2}Z_{4}X_{5}+{\bf
0.04177}X_{2}X_{5}{\bf-0.04177}Z_{1}X_{2}X_{5}{\bf-0.04177}X_{2}Z_{4}X_{5}\\\
&+{\bf 0.04177}Z_{1}X_{2}Z_{4}X_{5}+{\bf 0.01232}X_{0}X_{1}X_{2}X_{5}+{\bf
0.01232}Y_{0}X_{1}Y_{2}X_{5}{\bf-0.01232}X_{0}X_{1}X_{2}Z_{4}X_{5}{\bf-0.01232}Y_{0}X_{1}Y_{2}Z_{4}X_{5}\\\
&{\bf-0.04165}X_{3}X_{5}+{\bf
0.03769}X_{3}Z_{4}X_{5}{\bf-0.00396}Y_{3}Y_{5}+{\bf
0.00839}X_{4}X_{5}{\bf-0.01015}Z_{3}X_{4}X_{5}\\\ &+{\bf
0.01355}Z_{0}Z_{3}X_{4}X_{5}+{\bf 0.01082}Z_{0}Z_{1}Z_{3}X_{4}X_{5}+{\bf
0.00854}Z_{2}Z_{3}X_{4}X_{5}+{\bf
0.01408}Z_{1}Z_{2}Z_{3}X_{4}X_{5}{\bf-0.01015}Y_{4}Y_{5}\\\ &+{\bf
0.01355}Z_{0}Y_{4}Y_{5}+{\bf 0.01082}Z_{0}Z_{1}Y_{4}Y_{5}+{\bf
0.00854}Z_{2}Y_{4}Y_{5}+{\bf 0.01408}Z_{1}Z_{2}Y_{4}Y_{5}+{\bf
0.00839}Z_{3}Y_{4}Y_{5}\\\
&{\bf-0.0173}Y_{0}Y_{1}Z_{3}X_{4}X_{5}{\bf-0.0173}X_{0}X_{1}Z_{2}Z_{3}X_{4}X_{5}{\bf-0.0173}Y_{0}Y_{1}Y_{4}Y_{5}{\bf-0.0173}X_{0}X_{1}Z_{2}Y_{4}Y_{5}+{\bf
0.01858}Z_{0}X_{1}X_{2}Z_{3}X_{4}X_{5}\\\ &+{\bf
0.01858}Y_{1}Y_{2}Z_{3}X_{4}X_{5}+{\bf 0.01858}Z_{0}X_{1}X_{2}Y_{4}Y_{5}+{\bf
0.01858}Y_{1}Y_{2}Y_{4}Y_{5}{\bf-0.01573}X_{0}X_{3}X_{4}X_{5}\\\ &+{\bf
0.01573}X_{0}Z_{1}X_{3}X_{4}X_{5}{\bf-0.01573}X_{0}Y_{3}X_{4}Y_{5}+{\bf
0.01573}X_{0}Z_{1}Y_{3}X_{4}Y_{5}+{\bf
0.01532}X_{1}X_{3}X_{4}X_{5}{\bf-0.01532}Z_{0}X_{1}Z_{2}X_{3}X_{4}X_{5}+{\bf
0.01532}X_{1}Y_{3}X_{4}Y_{5}\\\
&{\bf-0.01532}Z_{0}X_{1}Z_{2}Y_{3}X_{4}Y_{5}+{\bf
0.01232}X_{2}X_{3}X_{4}X_{5}{\bf-0.01232}Z_{1}X_{2}X_{3}X_{4}X_{5}+{\bf
0.01232}X_{2}Y_{3}X_{4}Y_{5}\\\ &{\bf-0.01232}Z_{1}X_{2}Y_{3}X_{4}Y_{5}+{\bf
0.01415}X_{0}X_{1}X_{2}X_{3}X_{4}X_{5}+{\bf
0.01415}Y_{0}X_{1}Y_{2}X_{3}X_{4}X_{5}+{\bf
0.01415}X_{0}X_{1}X_{2}Y_{3}X_{4}Y_{5}\\\ &+{\bf
0.01415}Y_{0}X_{1}Y_{2}Y_{3}X_{4}Y_{5}\end{split}$ (26)

## References

  * Georgescu  _et al._ [2014] I. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Reviews of Modern Physics 86, 153 (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).
  * Lanyon  _et al._ [2010] B. P. Lanyon, J. D. Whitfield, G. G. Gillett, M. E. Goggin, M. P. Almeida, I. Kassal, J. D. Biamonte, M. Mohseni, B. J. Powell, M. Barbieri, A. Aspuru-Guzik, and A. G. White, Towards quantum chemistry on a quantum computer, Nature Chemistry 2, 106 (2010).
  * 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, 4213 (2014).
  * Shen  _et al._ [2017] Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung, and K. Kim, Quantum implementation of the unitary coupled cluster for simulating molecular electronic structure, Phys. Rev. A 95, 020501 (2017).
  * Shen  _et al._ [2018] Y. Shen, Y. Lu, K. Zhang, J. Zhang, S. Zhang, J. Huh, and K. Kim, Quantum optical emulation of molecular vibronic spectroscopy using a trapped-ion device, Chem. Sci. 9, 836 (2018).
  * Hempel  _et al._ [2018] C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt, and C. F. Roos, Quantum chemistry calculations on a trapped-ion quantum simulator, Phys. Rev. X 8, 031022 (2018).
  * 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 (2017).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Many-body physics with individually controlled Rydberg atoms, Nature Physics 16, 132 (2020).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Programmable quantum simulation of 2D antiferromagnets with hundreds of Rydberg atoms, Nature 595, 233 (2021), arXiv: 2012.12268.
  * Ebadi  _et al._ [2022] S. Ebadi _et al._ , Quantum optimization of maximum independent set using rydberg atom arrays, Science 376, 1209 (2022).
  * Chen  _et al._ [2023] C. Chen _et al._ , Continuous Symmetry Breaking in a Two-dimensional Rydberg Array, Nature , 1 (2023).
  * Henriet  _et al._ [2020] L. Henriet, L. Beguin, A. Signoles, T. Lahaye, A. Browaeys, G.-O. Reymond, and C. Jurczak, Quantum computing with neutral atoms, Quantum 4, 327 (2020), arXiv: 2006.12326.
  * Wecker  _et al._ [2015] D. Wecker, M. B. Hastings, and M. Troyer, Progress towards practical quantum variational algorithms, Physical Review A 92, 042303 (2015).
  * 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), arXiv: 2012.09265.
  * Barkoutsos  _et al._ [2020] P. K. Barkoutsos, G. Nannicini, A. Robert, I. Tavernelli, and S. Woerner, Improving Variational Quantum Optimization using CVaR, Quantum 4, 256 (2020), arXiv: 1907.04769.
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 4812 (2018).
  * Meitei  _et al._ [2021] O. R. Meitei, B. T. Gard, G. S. Barron, D. P. Pappas, S. E. Economou, E. Barnes, and N. J. Mayhall, Gate-free state preparation for fast variational quantum eigensolver simulations, npj Quantum Information 7, 155 (2021).
  * Wakaura  _et al._ [2021] H. Wakaura, T. Tomono, and S. Yasuda, Evaluation on Genetic Algorithms as an optimizer of Variational Quantum Eigensolver(VQE) method, arXiv:2110.07441 [quant-ph] (2021), arXiv: 2110.07441.
  * Banchi and Crooks [2021] L. Banchi and G. E. Crooks, Measuring Analytic Gradients of General Quantum Evolution with the Stochastic Parameter Shift Rule, Quantum 5, 386 (2021), arXiv: 2005.10299.
  * Gacon  _et al._ [2021] J. Gacon, C. Zoufal, G. Carleo, and S. Woerner, Simultaneous perturbation stochastic approximation of the quantum fisher information, Quantum 5, 567 (2021).
  * Piskor  _et al._ [2022] T. Piskor, J.-M. Reiner, S. Zanker, N. Vogt, M. Marthaler, F. K. Wilhelm, and F. G. Eich, Using gradient-based algorithms to determine ground-state energies on a quantum computer, Physical Review A 105, 062415 (2022).
  * Huang  _et al._ [2021] H.-Y. Huang, R. Kueng, and J. Preskill, Efficient estimation of Pauli observables by derandomization, Physical Review Letters 127, 030503 (2021), arXiv: 2103.07510.
  * Elben  _et al._ [2019] A. Elben, B. Vermersch, C. F. Roos, and P. Zoller, Statistical correlations between locally randomized measurements: A toolbox for probing entanglement in many-body quantum states, Physical Review A 99, 052323 (2019).
  * Kokail  _et al._ [2019] C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, and P. Zoller, Self-Verifying Variational Quantum Simulation of the Lattice Schwinger Model, Nature 569, 355 (2019), arXiv: 1810.03421.
  * Nam  _et al._ [2020] Y. Nam _et al._ , Ground-state energy estimation of the water molecule on a trapped-ion quantum computer, npj Quantum Information 6, 33 (2020).
  * Dalyac  _et al._ [2021] C. Dalyac, L. Henriet, E. Jeandel, W. Lechner, S. Perdrix, M. Porcheron, and M. Veshchezerova, Qualifying quantum approaches for hard industrial optimization problems. A case study in the field of smart-charging of electric vehicles, EPJ Quantum Technology 8, 1 (2021).
  * Parra-Rodriguez  _et al._ [2020] A. Parra-Rodriguez, P. Lougovski, L. Lamata, E. Solano, and M. Sanz, Digital-analog quantum computation, Physical Review A 101, 022305 (2020).
  * Fedorov  _et al._ [2022] D. A. Fedorov, B. Peng, N. Govind, and Y. Alexeev, VQE method: a short survey and recent developments, Materials Theory 6, 2 (2022).
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Über das Paulische Äquivalenzverbot, Zeitschrift für Physik 47, 631 (1928).
  * Bravyi and Kitaev [2002] S. Bravyi and A. Kitaev, Fermionic quantum computation, Annals of Physics 298, 210 (2002), arXiv: quant-ph/0003137.
  * Bartlett  _et al._ [1989] R. J. Bartlett, S. A. Kucharski, and J. Noga, Alternative coupled-cluster ansätze II. The unitary coupled-cluster method, Chemical Physics Letters 155, 133 (1989).
  * O’Malley  _et al._ [2016] P. J. J. O’Malley _et al._ , Scalable quantum simulation of molecular energies, Phys. Rev. X 6, 031007 (2016).
  * Helgaker  _et al._ [2000] T. Helgaker, P. Jørgensen, and J. Olsen, Coupled-cluster theory, in  _Molecular Electronic-Structure Theory_ (John Wiley & Sons, Ltd, 2000) Chap. 13, pp. 648–723.
  * Storn and Price [1997] R. Storn and K. Price, Differential Evolution – A Simple and Efficient Heuristic for global Optimization over Continuous Spaces, Journal of Global Optimization 11, 341 (1997).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A Quantum Approximate Optimization Algorithm, arXiv:1411.4028 (2014).
  * Leclerc  _et al._ [2022] L. Leclerc, L. Ortiz-Guitierrez, S. Grijalva, B. Albrecht, J. R. K. Cline, V. E. Elfving, A. Signoles, L. Henriet, G. Del Bimbo, U. A. Sheikh, M. Shah, L. Andrea, F. Ishtiaq, A. Duarte, S. Mugel, I. Caceres, M. Kurek, R. Orus, A. Seddik, O. Hammammi, H. Isselnane, and D. M’tamon, Financial risk management on a neutral atom quantum processor, arXiv:2212.03223 (2022).
  * Coelho  _et al._ [2022] W. d. S. Coelho, M. D’Arcangelo, and L.-P. Henry, Efficient protocol for solving combinatorial graph problems on neutral-atom quantum processors, arXiv:2207.13030 (2022).
  * Notarnicola  _et al._ [2021] S. Notarnicola, A. Elben, T. Lahaye, A. Browaeys, S. Montangero, and B. Vermersch, A randomized measurement toolbox for Rydberg quantum technologies, arXiv:2112.11046 (2021).
  * Qiskit contributors [2023] Qiskit contributors, Qiskit: An open-source framework for quantum computing (2023).
  * Muller [2022] R. Muller, PyQuante2 (2022).
  * Silvério  _et al._ [2022] H. Silvério, S. Grijalva, C. Dalyac, L. Leclerc, P. J. Karalekas, N. Shammah, M. Beji, L.-P. Henry, and L. Henriet, Pulser: An open-source package for the design of pulse sequences in programmable neutral-atom arrays, Quantum 6, 629 (2022).
  * Powell [1964] M. J. D. Powell, An efficient method for finding the minimum of a function of several variables without calculating derivatives, The Computer Journal 7, 155 (1964).
  * Nelder and Mead [1965] J. A. Nelder and R. Mead, A Simplex Method for Function Minimization, The Computer Journal 7, 308 (1965).
  * Schymik  _et al._ [2022] K.-N. Schymik, B. Ximenez, E. Bloch, D. Dreon, A. Signoles, F. Nogrette, D. Barredo, A. Browaeys, and T. Lahaye, In situ equalization of single-atom loading in large-scale optical tweezer arrays, Physical Review A 106, 022611 (2022).
  * Barredo  _et al._ [2018] D. Barredo, V. Lienhard, S. de Léséleuc, T. Lahaye, and A. Browaeys, Synthetic three-dimensional atomic structures assembled atom by atom, Nature 561, 79 (2018).
  * Dalyac and Henriet [2022] C. Dalyac and L. Henriet, Embedding the mis problem for non-local graphs with bounded degree using 3d arrays of atoms, arXiv:2209.05164 (2022).
  * Bidzhiev  _et al._ [2023] K. Bidzhiev, A. Wennersteen, M. Beji, M. Dagrada, M. D’Arcangelo, S. Grijalva, A.-C. L. Henaff, A. Quelle, and A. S. Naik, Cloud on-demand emulation of quantum dynamics with tensor networks, arXiv:2302.05253 (2023).
  * Rudolph  _et al._ [2022] M. S. Rudolph, J. Miller, J. Chen, A. Acharya, and A. Perdomo-Ortiz, Synergy between quantum circuits and tensor networks: Short-cutting the race to practical quantum advantage, arXiv:2208.13673 (2022).
  * Sim  _et al._ [2019] S. Sim, P. D. Johnson, and A. Aspuru-Guzik, Expressibility and entangling capability of parameterized quantum circuits for hybrid quantum-classical algorithms, Advanced Quantum Technologies 2, 1900070 (2019).
  * Holmes  _et al._ [2022] Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles, Connecting ansatz expressibility to gradient magnitudes and barren plateaus, PRX Quantum 3, 010313 (2022).
  * Tangpanitanon  _et al._ [2020] J. Tangpanitanon, S. Thanasilp, N. Dangniam, M.-A. Lemonde, and D. G. Angelakis, Expressibility and trainability of parametrized analog quantum systems for machine learning applications, Phys. Rev. Res. 2, 043364 (2020).
  * de Léséleuc  _et al._ [2018] S. de Léséleuc, D. Barredo, V. Lienhard, A. Browaeys, and T. Lahaye, Analysis of imperfections in the coherent optical excitation of single atoms to Rydberg states, Physical Review A 97, 053803 (2018).
  * Henriet [2020] L. Henriet, Robustness to spontaneous emission of a variational quantum algorithm, Physical Review A 101, 012335 (2020).
  * Guo  _et al._ [2022] S. Guo _et al._ , Scalable quantum computational chemistry with superconducting qubits, arXiv:2212.08006 (2022).
  * Wu  _et al._ [2023] B. Wu, J. Sun, Q. Huang, and X. Yuan, Overlapped grouping measurement: A unified framework for measuring quantum states, Quantum 7, 896 (2023).