# Tailgating quantum circuits for high-order energy derivatives

Jack Ceroni<EMAIL_ADDRESS>Xanadu, Toronto, ON, M5G 2C8, Canada
Department of Mathematics, University of Toronto, Toronto, ON, M5S 3E1, Canada
Alain Delgado  Xanadu, Toronto, ON, M5G 2C8, Canada  Soran Jahangiri  Xanadu,
Toronto, ON, M5G 2C8, Canada  Juan Miguel Arrazola  Xanadu, Toronto, ON, M5G
2C8, Canada

###### Abstract

To understand the chemical properties of molecules, it is often important to
study derivatives of energies with respect to nuclear coordinates or external
fields. Quantum algorithms for computing energy derivatives have been
proposed, but only limited work has been done to address the specific
challenges that arise in this context, where calculations are more complicated
and involve more stringent requirements on accuracy compared to single-point
energy calculations. In this work, we introduce a technique to improve the
performance of variational quantum circuits calculating energy derivatives.
The method, which we refer to as tailgating, is an adaptive procedure that
selects gates based on their gradient with respect to the expectation value of
Hamiltonian derivatives. These gates are then added at the end of a quantum
circuit originally designed to calculate ground- or excited-state energies. A
distinguishing feature of this approach is that the appended gates do not need
to be optimized: their parameters can be set to zero and varied only for the
purpose of computing energy derivatives, via calculating derivatives with
respect to circuit parameters. We support the validity of this method by
establishing sufficient conditions for a circuit to compute accurate energy
gradients. This is achieved through a connection between energy derivatives
and eigenstates of Taylor approximations of the Hamiltonian. We illustrate the
advantages of the tailgating approach by performing simulations calculating
the vibrational modes of beryllium hydride and water: quantities that depend
on second-order energy derivatives.

##  I Introduction

Quantum algorithms have been studied as a potential avenue to address some of
the challenges of accurately simulating the electronic structure of molecules
McArdle et al. (2020); Aspuru-Guzik et al. (2005); Quantum et al. (2020);
Lanyon et al. (2010). Advanced quantum algorithms for performing these
calculations require a large number of qubits and deep circuits Cao et al.
(2019); O’Brien et al. (2021); Su et al. (2021); Low and Chuang (2019), which
makes them challenging to implement on currently available quantum hardware or
classical simulators. Variational quantum algorithms (VQAs) have been instead
considered as a platform to perform proof-of-principle simulations of
molecules on existing devices. VQAs are a paradigm for designing quantum
algorithms in which the structure of a quantum circuit is fixed and the gate
parameters are optimized to minimize a suitable cost function, typically the
expectation value of an observable with respect to the prepared wavefunction
Cerezo et al. (2021).

In the context of quantum chemistry, one of the most-studied algorithms is the
variational quantum eigensolver (VQE) Kandala et al. (2017); Peruzzo et al.
(2014). Here, a quantum circuit is optimized in order to minimize the
expectation value of a molecular Hamiltonian with respect to the output state
of the circuit, yielding an approximation of the ground-state energy. Similar
algorithms have been proposed to calculate other electronic properties such as
excited-state energies Higgott et al. (2019); Nakanishi et al. (2019); McClean
et al. (2017); Parrish et al. (2019), equilibrium geometries Delgado et al.
(2021), dipole transition amplitudes Ibe et al. (2020), transition states Azad
and Singh (2022), and molecular dynamics Sokolov et al. (2021). It has also
been shown that computing first- and second-order energy derivatives of ground
and excited state energies can be performed using variational methods Mitarai
et al. (2020); O’Brien et al. (2021); Azad and Singh (2022). Second-order
energy derivatives can be used to characterize the vibrational structure of
molecules, which allows for the prediction of chemical properties such as
vibronic spectra Huh et al. (2015), simulating local mode dynamics Sparrow et
al. (2018); Jahangiri et al. (2020) and electron transport Jahangiri et al.
(2021), and computing thermodynamic observables Stober et al. (2022).

Beyond their scope of application, the techniques used for constructing
variational quantum circuits have also become more sophisticated. One popular
approach is the ADAPT-VQE method Grimsley et al. (2019); Tang et al. (2021);
Bilkis et al. (2021). Starting from a pool of gates, this algorithm
iteratively adds gates to a circuit by selecting those with large gradients
with respect to the target cost function. This procedure has the advantage of
being tailored to an individual Hamiltonian and typically leads to shorter-
depth circuits without sacrificing accuracy. ADAPT-VQE has shown good results
when used to compute energy eigenvalues of some fixed Hamiltonian Claudino et
al. (2020), and circuits constructed with adaptive methods have been effective
for molecular geometry optimization procedures in small molecules, where it is
required to calculate first-order energy derivatives Delgado et al. (2021).
However, no prior work has been done to analyze the capability of adaptive
circuits for computing second-order energy derivatives beyond the simple case
of the hydrogen molecule Mitarai et al. (2020).

In this work, we argue that standard adaptive procedures fail to produce
circuits capable of calculating accurate high-order energy derivatives because
they only take into account limited information about the problem Hamiltonian.
We then address this issue by introducing the method of _tailgating_ : a
procedure in which parameterized gates are adaptively added to the end of a
variational circuit, with their “optimized” variational parameters set to
zero. The gates are selected by evaluating their gradient with respect to the
expectation value of derivatives of the Hamiltonian. They do not contribute to
the preparation of the ground state, but do have an effect when calculating
energy derivatives, increasing the accuracy of the calculations. This
procedure is particularly suitable for adaptive circuits, but in principle any
circuit can be tailgated to increase accuracy when calculating energy
derivatives.

We support the validity of this method by establishing sufficient conditions
for a quantum circuit to calculate accurate energy derivatives. This is
achieved through a connection between Taylor approximations of the Hamiltonian
and energy derivatives. Finally, we illustrate the advantages of tailgating
with numerical examples, calculating the energy Hessians of beryllium hydride
and water molecules. We subsequently use the Hessians to compute the
vibrational frequencies of each molecule. We demonstrate that, in both cases,
tailgated circuits provide values that are considerably more accurate than
those obtained from standard adaptive circuits.

##  II Computing Accurate Energy Derivatives

We begin by introducing notation that will be used throughout this work. Let
$U(\theta)$ be a parameterized quantum circuit that prepares an output state
$|\psi(\theta)\rangle=U(\theta)|0\rangle$, where
$\theta=(\theta_{1},\theta_{2},\ldots,\theta_{M})$ are the parameters of the
circuit. We use $H(R)$ to represent a parameterized Hamiltonian, where
$R=(R_{1},R_{2},\ldots,R_{N})$ is a vector of parameters. For concreteness, we
focus on the case where $R$ represents the nuclear coordinates of a molecule.
Let $\theta^{*}(R)$ denote the optimal parameters such that the output state
$|\psi(\theta^{*}(R))\rangle$ closest approximates the ground state
$\ket{\psi_{0}(R)}$. We also describe the corresponding ground-state energy as

$E(R)=\langle\psi_{0}(R)|H(R)|\psi_{0}(R)\rangle.$ (1)

We focus on the ground state for simplicity, but the analysis described
throughout this work applies equally to any eigenstate of the Hamiltonian.
Finally, we define

$\tilde{E}(\theta,R)=\langle\psi(\theta)|H(R)|\psi(\theta)\rangle,$ (2)

and $\tilde{E}(R)=\tilde{E}(\theta^{*}(R),R)$ as the approximation of the
energy eigenvalue resulting from the output state of the circuit. For
simplicity, we use $|\psi(R)\rangle$ as a shorthand for
$|\psi(\theta^{*}(R))\rangle$. We work under the assumption that the family of
Hamiltonians $H(R)$ is non-degenerate, so we can truly speak of a single
ground state.

After identifying and optimizing a circuit $U(\theta)$ such that
$\tilde{E}(R_{0})\approx E(R_{0})$ for a given $R_{0}$, e.g., the equilibrium
geometry of the molecule, our goal is to compute accurate $n$-th order energy
derivatives at $R_{0}$. In other words, we want to ensure that

$\frac{\partial^{n}\tilde{E}(R_{0})}{\partial R_{j_{1}}\cdots
R_{j_{n}}}\approx\frac{\partial^{n}E(R_{0})}{\partial R_{j_{1}}\cdots
R_{j_{n}}},$ (3)

for some choice of $n$ and coordinates $\\{R_{j_{1}},\ \ldots,\ R_{j_{n}}\\}$.

Now that we have introduced the necessary notation, we briefly discuss the
computation of energy derivatives, and the associated difficulties of
performing these calculations. Computing quantities which rely on $n$-th order
energy derivatives, with variational quantum circuits, already suffer from the
issue of high sampling cost. This occurs even in the case of second-order
derivatives, where the number of samples required to compute a single second-
order energy derivative to error $\epsilon$ scales as Mitarai et al. (2020)

$N_{\text{samples}}=\tilde{O}\left(\frac{n^{4}M^{2}}{\epsilon^{\prime
2}}\right).$ (4)

where $n$ is the number of qubits, $M$ is the number of gate parameters in the
variational circuit, and $\epsilon^{\prime}$ is an error parameter defined in
Ref. Mitarai et al. (2020). It depends linearly on $\epsilon$, and implicitly
on $n$ and $M$. In this paper, we do not focus on the important problem of
sample complexity, but rather on a different roadblock related to the choice
of circuit ansatz. When computing energy derivatives, there will generally be
contributions due to derivatives of the state with respect to the nuclear
coordinates. This can be seen for example in the expression for second-order
derivatives Azad and Singh (2022):

$\displaystyle\frac{\partial^{2}E(R)}{\partial R_{i}\partial R_{j}}$
$\displaystyle=\langle\psi_{0}(R)|\frac{\partial^{2}H(R)}{\partial
R_{i}\partial R_{j}}|\psi_{0}(R)\rangle$
$\displaystyle+2\text{Re}\left[\langle\psi_{0}(R)|\frac{\partial
H(R)}{\partial R_{i}}\frac{\partial|\psi_{0}(R)\rangle}{\partial
R_{j}}\right],$ (5)

which contains a contribution from the state derivative
$\partial|\psi_{0}(R)\rangle/\partial R_{j}$. Expressions for derivatives of
this form are derived in Appendix B. This suggests, at a high level, that in
order to evaluate energy derivatives, we need a circuit that can prepare
approximate ground states not only at $R_{0}$, but also in some neighbourhood
around $R_{0}$.

Guaranteeing that a particular circuit possesses this property is not
straightforward. For standard adaptive methods that only select gates based on
the Hamiltonian $H(R_{0})$, we argue that there may be situations in which
certain omitted gates have negligible gradient at $R_{0}$, but are required to
produce accurate energy derivatives. More precisely, these omitted gates,
while having vanishingly small gradients during standard adaptive algorithms
like ADAPT-VQE, have considerable effects on the degree of accuracy to which
the state $|\psi(R)\rangle$ matches the state $|\psi_{0}(R)\rangle$, in some
neighbourhood around $R_{0}$. This change in accuracy can affect state
derivatives. For instance, suppose that

$|\psi(R)\rangle=U(\theta^{*}(R))|0\rangle=|\psi_{0}(R)\rangle+\epsilon(R)|\zeta(R)\rangle,$
(6)

for all $R$ in a neighbourhood around $R_{0}$, and
$\epsilon(R)=|||\psi(R)\rangle-|\psi_{0}(R)\rangle||$. If $\epsilon(R)$ is
sufficiently small for a neighbourhood around $R_{0}$, the circuit $U$
prepares a state $|\psi(R)\rangle$ which is a good approximation for
$|\psi_{0}(R)\rangle$. However, it could be the case that the magnitude of
$\partial\epsilon(R)|\zeta(R)\rangle/\partial R_{j}$ is large, implying that
first derivatives of $|\psi(R)\rangle$ and $|\psi_{0}(R)\rangle$ differ
significantly. Intuitively, even though $|\psi(R)\rangle$ is close to
$|\psi_{0}(R)\rangle$, the error between the two states could change rapidly
as a function of $R$. This could then lead to large differences in second-
order energy derivatives, as is indicated by Eq. (II).

Evidence for this phenomenon occurring in circuits prepared by adaptive
methods is illustrated in Fig. 1. Here, we plot the fidelity
$|\langle\psi(R_{0}+\delta e_{H})|\psi_{0}(R_{0}+\delta e_{H})\rangle|^{2}$,
where $|\psi_{0}(R)\rangle$ is the ground state of the BeH2 molecule at $R$
and $\delta e_{H}$ is a vector which stretches the two hydrogen atoms a
distance $\delta$ in opposite directions. We adaptively prepare a circuit $U$
with respect to the Hamiltonian $H(R_{0})$ using a pool of single and double
excitation gates Arrazola et al. (2021a). We then perform VQE (to a pre-
specified convergence criterion of gradients having magnitude less than
$10^{-5}$) for a collection of $\delta$ ranging from $-0.75\ \text{Bohr}$ to
$0.75\ \text{Bohr}$. The fidelity is close to $1$ for all $\delta$, but the
derivative of the fidelity with respect to $\delta$ is non-zero and decreasing
for large $\delta$.

Figure 1: (A) Fidelity $|\langle\psi(R_{0}+\delta e_{H})|\psi_{0}(R_{0}+\delta
e_{H})\rangle|^{2}$ as a function of $\delta$: the distance the hydrogen atoms
are displaced from equilibrium. Here we consider the fidelity between the
ground state of the BeH2 molecule and an approximation prepared by an
optimized circuit yielded from ADAPT-VQE. (B) The derivative of the fidelity,
computed using a central finite-difference.

Now, we briefly develop some intuition as to how this issue can be resolved by
considering gradients of gates with respect to derivatives of the Hamiltonian
$H(R)$. Consider an adaptive procedure, where gates $\\{G_{a}\\}$ are selected
to build an ansatz $U$ which yields an approximation to the ground state of
$H(R_{0})$, based on the value of gradients of the form

$\frac{\partial}{\partial\theta_{a}}\bra{\psi}G_{a}(\theta_{a})^{\dagger}H(R_{0})G_{a}(\theta_{a})\ket{\psi},$
(7)

where $|\psi\rangle$ is a state constructed during the previous iteration of
the adaptive procedure. If the value of this gradient is large, then the gate
$G_{a}$ is selected for the circuit ansatz which is used to prepare ground
state of $H(R_{0})$. Eventually, after iterating through all gates in
$\\{G_{a}\\}$, we construct the ansatz $U$ which yields
$|\psi(R_{0})\rangle\approx|\psi_{0}(R_{0})\rangle$.

Now, for each $G_{a}$, consider the function

$\Delta_{a}(R)=\frac{\partial}{\partial\theta_{a}}\langle\psi(R_{0})|G_{a}(\theta_{a})^{\dagger}H(R)G_{a}(\theta_{a})|\psi(R_{0})\rangle|_{\theta_{a}=0}.$
(8)

As in the original adaptive procedure, $\Delta_{a}(R)$ can be thought of as
measuring the effect of the gate $G_{a}$ on an optimization procedure which
attempts to minimize the expectation value of $H(R)$, starting from the state
$|\psi(R_{0})\rangle$. In general, we will expect $\Delta_{a}(R)$ to be small
at $R=R_{0}$ and likely also in a neighbourhood around $R$, as adding any gate
that was already selected or not previously selected will likely not have a
large effect on improving the accuracy of the approximate ground state
$|\psi(R_{0})\rangle$. However, these gates could still provide small
corrections to the state. As was indicated above, this is precisely what we
should consider to ensure that first-order state derivatives are accurate.
Thus, we turn our attention to the derivatives of $\Delta_{a}(R)$ evaluated at
$R=R_{0}$,

$\frac{\partial\Delta_{a}(R_{0})}{\partial
R_{j}}=\frac{\partial}{\partial\theta_{a}}\langle\psi(R_{0})|G_{a}(\theta_{a})^{\dagger}\frac{\partial
H(R_{0})}{\partial R_{j}}G_{a}(\theta_{a})|\psi(R_{0})\rangle.$ (9)

If a derivative $\partial\Delta_{a}(R_{0})/\partial R_{j}$ is large, this
indicates that a gate $G_{a}$ contributes to the preparation of
$|\psi(R_{0})\rangle$ and $|\psi(R)\rangle$ for $R$ close to $R_{0}$ with
significantly different magnitudes, relative to the distance from $R_{0}$ to
$R$. If we add all gates of this form to the ansatz, our hope is that this
will lead to more accurate state derivatives, which will in turn lead to more
accurate energy derivatives.

This treatment only takes into account first-order differences, so a natural
subsequent question to ask is whether we have to consider higher-order
derivatives of $\Delta_{a}(R)$ to prepare an ansatz which yields accurate
higher-order energy derivatives. To answer this question in the affirmative,
in the following sections we take a more rigorous approach, via Taylor series
expansions of the Hamiltonian $H(R)$.

###  II.1 Taylor Series Expansions of Hamiltonians

Keeping the importance of Hamiltonian derivatives in mind, we turn our
attention to finding sufficient conditions for a circuit $U(\theta)$ that can
approximate ground states $|\psi_{0}(R)\rangle$ in a neighbourhood around
$R_{0}$ by Taylor-expanding $H(R)$ around $R_{0}$. In this form, we
approximate $H(R)$ in terms of its derivatives at $R_{0}$. Let $H_{n}(R)$ be
the $n$-th order Taylor approximation of $H$ centred at $R_{0}$:

$\displaystyle H_{n}(R)=H(R_{0})+\displaystyle\sum_{j=1}^{N}\frac{\partial
H(R_{0})}{\partial R_{j}}(R_{j}-R_{0j})+\cdots$
$\displaystyle+\frac{1}{n!}\displaystyle\sum_{j_{1},\ldots,j_{n}=1}^{N}\frac{\partial^{n}H(R_{0})}{\partial
R_{j_{1}}\cdots\partial
R_{j_{n}}}(R_{j_{1}}-R_{0j_{1}})\cdots(R_{j_{n}}-R_{0j_{n}}).$ (10)

The Taylor expansion $H_{n}(R)$ provides a good approximation of $H(R)$ when
$R$ is close to $R_{0}$, provided we truncate at sufficiently high order.
Thus, it is reasonable to suspect that for sufficiently large $n$, the
eigenvalues of $H_{n}$ closely approximate those of $H$ for a small
neighbourhood around $R_{0}$. As a result, the energy derivatives computed
with respect to the two Hamiltonians are equal to some order. This is
formalized in the following theorem:

######  Theorem 1.

Let $|\phi_{0}(R)\rangle$ be the ground state of $H_{n-1}(R)$. Suppose that
there exists $\delta>0$ such that a circuit $U(\theta)$ can prepare the state
$|\phi_{0}(R)\rangle$ for all $R$ such that $0<\|R-R_{0}\|<\delta$. Then it
holds that

$\frac{\partial^{p}E(R_{0})}{\partial R_{j_{1}}\cdots\partial
R_{j_{p}}}=\frac{\partial^{p}\tilde{E}(R_{0})}{\partial
R_{j_{1}}\cdots\partial R_{j_{p}}},$ (11)

for all $p$ such that $0\leq p\leq n$, and all choices of coordinates,
$\\{R_{j_{1}},\ ...,\ R_{j_{p}}\\}$.

The proof is given in Appendix A.

The implication of this result is crucial for our purposes: circuits that can
compute accurate energy eigenstates of Taylor approximations of the
Hamiltonian, in some neighbourhood around $R_{0}$, can be used to calculate
accurate energy derivatives. Importantly, if the circuit can prepare accurate
eigenstates for a Taylor expansion of order $n-1$, then it is possible to
calculate accurate energy derivatives up to order $n$. In the case of $n=1$,
our result implies that to compute accurate first-order derivatives, it
suffices that the circuit can prepare accurate ground states of the
Hamiltonian $H(R_{0})$ at the equilibrium geometry.

This supports the observation that adaptive circuits are effective for
computing first-order derivatives Delgado et al. (2021). However, the
sufficient condition of Theorem 1 indicates that the ability to prepare ground
states at equilibrium may not be enough for higher-order derivatives. In such
cases, it is helpful to design circuits that can prepare ground states of
higher-order Taylor approximations in the neighbourhood around equilibrium.
This in turn can be achieved by focusing on the derivatives of the
Hamiltonian, which determine the coefficients in the Taylor expansion. This is
the main insight behind the tailgating algorithm that we describe next.

##  III Tailgating

Figure 2: An illustration of the tailgating procedure for the case of second-
order energy derivatives. (A) We begin with a circuit which can prepare good
approximations of the ground state energy of $H(R)$, along with a pool of
gates and Hamiltonian derivatives. (B) We perform an adaptive selection of
gates by iteratively applying each gate in the pool to the end of the circuit,
and computing derivatives of this particular gate with respect to Hamiltonian
derivatives. We add all gates with non-negligible gradients to the new,
tailgated circuit. (C) This circuit can then be used to compute accurate
second-order energy derivatives. Note that to calculate the desired energy
derivatives, we evaluate derivatives of the added gates at
$\theta_{p}=\theta_{q}=0$.

Similar to the expression in Eq. (7), consider an adaptive procedure where
gates $\\{G_{a}\\}$ are selected based on the value of gradients of the form

$\frac{\partial
C_{n}(\theta)}{\partial\theta_{k}}:=\frac{\partial}{\partial\theta_{k}}\bra{\psi}G_{a}(\theta)^{\dagger}H_{n-1}(R)G_{a}(\theta)\ket{\psi}.$
(12)

If the value of these gradients is large for all $R$ around $R_{0}$, then the
gate $G_{a}$ can be helpful in preparing ground states of the Taylor
approximation $H_{n-1}(R)$ in this region. From Theorem 1, this implies that
it can aid in calculating accurate energy gradients. Rather than performing
this calculation for a large collection of points $R$, we instead employ Eq.
(II.1) to express $H_{n-1}(R)$ as a linear combination of Hamiltonian
derivatives. For example, consider the case $n=2$, where we have

$\displaystyle\frac{\partial C_{2}(\theta)}{\partial\theta_{k}}$
$\displaystyle=\frac{\partial}{\partial\theta_{k}}\bra{\psi}G_{a}(\theta)^{\dagger}H(R_{0})G_{a}(\theta)\ket{\psi}$
$\displaystyle+\sum_{j=1}^{N}(R_{j}-R_{0j})\frac{\partial}{\partial\theta_{k}}\bra{\psi}G_{a}(\theta)^{\dagger}\frac{\partial
H(R_{0})}{\partial R_{j}}G_{a}(\theta)\ket{\psi}.$ (13)

There are contributions due to both $H(R_{0})$ and $\frac{\partial
H(R_{0})}{\partial R_{j}}$, meaning we should select gates depending on their
gradient with respect to both terms. This selection procedure will capture all
of the gates which have non-negligible gradients relative to $H_{2}(R)$. The
same holds true for higher-order $H_{n}(R)$, where we compute expectation
values with respect to higher-order derivatives of $H(R)$ at $R_{0}$ as well.

A naive approach based on the gate-selection strategy outlined above is to
evaluate gradients with respect to all of the necessary Hamiltonian
derivatives at each step of the adaptive procedure, and then optimize the
resulting circuit. This can be unnecessarily costly if there are many
additional gates arising from the Hamiltonian derivatives that would otherwise
be absent when considering only $H(R_{0})$, as in standard adaptive schemes.

An improvement comes from realizing that only the gates that are needed to
prepare approximate ground states need to be optimized. From Theorem 1, it
suffices to add gates such that the circuit is _capable_ of preparing ground
states of $H_{n-1}(R)$, but there is no need to optimize them to do so. These
gates can simply be added at the end of the circuit with their initial
parameters set to zero. They play a role only when computing the expressions
for energy gradients, which may require shifting their parameters away from
zero, as we will have to compute derivatives with respect to circuit
parameters (see Eq. (19) and Eq. (20)). Below we describe the tailgating
procedure in full:

  1. 1.

Given some Hamiltonian $H(R)$, identify a variational circuit $U(\theta)$
which can approximate the ground state of $H(R)$ at $R_{0}$. Let
$U(\theta^{*}(R))$ denote the optimized circuit.

  2. 2.

Let $\mathcal{H}$ be the collection of Hamiltonian derivatives of $H(R)$ up to
order $n-1$, evaluated at $R_{0}$. Similarly, let $\\{G_{a}\\}$ be a pool of
parameterized quantum gates and for each $G_{a}$, define
$V_{a}(\theta)=G_{a}(\theta)U(\theta^{*}(R_{0}))$. Then for each Hamiltonian
derivative $\frac{\partial^{m}H(R_{0})}{\partial j_{1}\cdots\partial
j_{m}}\in\mathcal{H}$, compute the magnitudes of the derivatives

$\displaystyle\frac{\partial}{\partial\theta_{k}}\Big{\langle}\psi(R_{0})\Big{|}V_{a}(\theta)^{\dagger}\frac{\partial^{m}H(R_{0})}{\partial
j_{1}\cdots\partial j_{m}}V_{a}(\theta)\Big{|}\psi(R_{0})\Big{\rangle},$ (14)

evaluated at $\theta=0$. If at least one is larger than some $\epsilon>0$, add
$G_{a}$ to the list $L$.

  3. 3.

Let $Q(\theta)=\prod_{G_{a}\in L}G_{a}(\theta_{a})$. Define the new
variational circuit $U_{\text{tailgated}}$ to be
$U_{\text{tailgated}}(\theta,\theta^{\prime})=Q(\theta)U(\theta^{\prime})$,
and new optimal parameters to be $(0,\theta^{*}(R))$.

Assuming $Q(0)=\mathbbm{1}$, we have

$U_{\text{tailgated}}(0,\theta^{*}(R))|0\rangle=U(\theta^{*}(R))|0\rangle=|\psi(R)\rangle,$
(15)

so the new circuit $U_{\text{tailgated}}$ approximates the ground state just
as the original circuit did. However, using the previous justification, this
circuit is designed to also be capable of preparing ground states of
$H_{n-1}(R)$ in a neighbourhood around $R_{0}$, and hence also accurate $n$-th
order energy derivatives.

We refer to this procedure as tailgating, as we are effectively adding a
sequence of gates to the end of a circuit $U$, to get a new circuit
$U_{\text{tailgated}}$, but are not optimizing any of their corresponding
variational parameters. The entire tailgating procedure is summarized in
Figure 2.

##  IV Numerical examples

| Original Circuit | Tailgated Circuit | GAMESS (FCI)
---|---|---|---
| BeH Frequencies
$\omega_{1}$ | $5088.77\ \text{cm}^{-1}$ | $2568.98\ \text{cm}^{-1}$ | $2569.52\ \text{cm}^{-1}$
$\omega_{2}$ | $2301.41\ \text{cm}^{-1}$ | $2300.68\ \text{cm}^{-1}$ | $2298.31\ \text{cm}^{-1}$
$\omega_{3}$ | $1018.53\ \text{cm}^{-1}$ | $784.72\ \text{cm}^{-1}$ | $780.1\ \text{cm}^{-1}$
$\omega_{4}$ | $1018.53\ \text{cm}^{-1}$ | $784.71\ \text{cm}^{-1}$ | $780.1\ \text{cm}^{-1}$
| H2O Frequencies
$\omega_{1}$ | $9901.80\ \text{cm}^{-1}$ | $3845.51\ \text{cm}^{-1}$ | $3812.60\ \text{cm}^{-1}$
$\omega_{2}$ | $3606.33\ \text{cm}^{-1}$ | $3605.40\ \text{cm}^{-1}$ | $3569.82\ \text{cm}^{-1}$
$\omega_{3}$ | $2047.07\ \text{cm}^{-1}$ | $2043.55\ \text{cm}^{-1}$ | $2036.99\ \text{cm}^{-1}$
Table 1: Calculated normal mode frequencies of BeH2 and H2O using the minimal
STO-3G basis set. The parameter $\omega_{j}$ denotes the $j$-th normal mode
frequency. Note that in the quantum methods, the BeH2 ground state was
prepared with fidelity $>0.9998$ and the H2O ground state was prepared with
fidelity $>0.9999$, but only the tailgated circuit gives accurate values in
both cases.

To highlight the effectiveness of the tailgating procedure, we study a
practical example: computing the normal mode vibrational frequencies of
molecules. Specifically, we consider BeH2 and H2O at their equilibrium
coordinates. Each molecular Hamiltonian is of the form

$H(R)=\displaystyle\sum_{pq}h_{pq}(R)c_{p}^{\dagger}c_{q}+\frac{1}{2}\displaystyle\sum_{pqrs}h_{pqrs}(R)c_{p}^{\dagger}c_{q}^{\dagger}c_{r}c_{s},$
(16)

where $R$ is a set of nuclear coordinates, $c^{\dagger}_{p}$ and $c_{p}$ are
the fermionic creation and annihilation operators, acting on the $p$-th
orbital, and

$h_{pq}(R)=\displaystyle\int dr\
\phi_{p}^{*}(r)\left(-\frac{\nabla^{2}}{2}-\displaystyle\sum_{I}\frac{Z_{I}}{|r-R_{I}|}\right)\phi_{q}(r),$
(17)

$h_{pqrs}(R)=\displaystyle\int dr_{1}dr_{2}\
\frac{\phi_{p}^{*}(r_{1})\phi_{q}^{*}(r_{2})\phi_{r}(r_{2})\phi_{s}(r_{1})}{|r_{1}-r_{2}|},$
(18)

are the one and two-electron integrals yielded from a set of molecular
orbitals $\phi_{p}(r)$, usually obtained by using the Hartree-Fock method
Arrazola et al. (2021b). Note that these molecular orbitals depend implicitly
on $R$. The normal mode frequencies of a molecule are precisely the
eigenvalues of the energy Hessian – the matrix of second-order energy
derivatives with respect to each of the nuclear coordinates,
$\frac{\partial^{2}E(R)}{\partial R_{i}\partial R_{j}}$.

We use the PennyLane library for quantum differentiable programming Bergholm
et al. (2018); Arrazola et al. (2021b) to simulate an adaptive circuit-
building procedure, followed by VQE, to compute approximate ground-state
energies at the equilibrium geometry. More specifically, we run an adaptive
gate selection procedure using a pool of gates composed of all admissible
single and double excitation gates for each particular molecule Arrazola et
al. (2021a), and optimize the resulting circuit by minimizing the energy
expectation value of the Hamiltonian, at geometry $R_{0}$, with gradient
descent. The coordinates $R_{0}$ are the equilibrium geometries of the
molecules obtained at the level of Hartree-Fock with the computational
chemistry package GAMESS Barca et al. (2020). All calculations are performed
using the minimal STO-3G basis set.

This process leads to short-depth circuits with optimized parameters that
prepare approximate ground states. We then employ the tailgating procedure to
adaptively add new gates to the circuit to make it suitable for second-order
energy derivatives, which are calculated with the analytic formula derived in
Ref. Mitarai et al. (2020),

$\frac{\partial^{2}\tilde{E}(R)}{\partial R_{i}\partial
R_{j}}=\displaystyle\sum_{a}\frac{\partial\theta_{a}^{*}(R)}{\partial
R_{i}}\frac{\partial}{\partial\theta_{a}}\Big{\langle}\psi(\theta)\Big{|}\frac{\partial
H(R)}{\partial
R_{j}}\Big{|}\psi(\theta)\Big{\rangle}\biggr{\rvert}_{\theta=\theta^{*}(R)}\\\
+\Big{\langle}\psi(R)\Big{|}\frac{\partial^{2}H(R)}{\partial R_{i}\partial
R_{j}}\Big{|}\psi(R)\Big{\rangle}.$ (19)

The derivatives of the optimal circuit parameters $\theta^{*}(R)$ with respect
to the nuclear coordinates are given by the response equation Mitarai et al.
(2020):

$\displaystyle\sum_{a}\frac{\partial\theta^{*}_{a}(R)}{\partial
R_{i}}\frac{\partial^{2}\tilde{E}(\theta,R)}{\partial\theta_{b}\partial\theta_{a}}\biggr{\rvert}_{\theta=\theta^{*}(R)}=-\frac{\partial}{\partial\theta_{b}}\frac{\partial\tilde{E}(R)}{\partial
R_{i}},$ (20)

where, by the Feynman-Hellman theorem, the first-order energy derivative is
given by

$\frac{\partial\tilde{E}(R)}{\partial
R_{i}}=\left\langle\psi(R)\left|\frac{\partial H(R)}{\partial
R_{i}}\right|\psi(R)\right\rangle.$ (21)

From here, we calculate the normal mode frequencies for the molecules,
comparing the values yielded from utilizing the original variational circuit
to a tailgated version of the circuit. To find the accurate values of the
normal mode frequencies for reference, we use GAMESS to perform full
configuration-interaction (FCI). The results of the calculations are
summarized in Table I. The frequencies yielded from the tailgated circuit are
in much closer agreement with the value given by GAMESS FCI, for both BeH2 and
H2O. The full details of the numerics can be found at
https://github.com/XanaduAI/tailgating.

In addition to the normal mode frequency calculations, we can also return to
plotting fidelity curves (as in Fig. 1) to gain a more intuitive understanding
of the effects of tailgating circuits. In Fig. 3, we provide a comparison of
the curves of fidelity $|\langle\psi(R_{0}+\delta
e_{H}^{\prime})|\psi_{0}(R_{0}+\delta e_{H}^{\prime})\rangle|^{2}$ for both
the cases where $|\psi(R)\rangle$ is prepared by a non-tailgated and a
tailgated circuit. In this example, we consider the molecular Hamiltonian
$H(R)$ corresponding to H${}_{3}^{+}$, where $\delta e_{H}^{\prime}$ is a
displacement vector corresponding to vertical stretching (in opposite
directions) of two of the hydrogen atoms from their equilibrium configurations
at $R_{0}$. As in Fig. 1, we plot the fidelity as a function of $\delta$.

Figure 3: A plot showing the fidelity $|\langle\psi(R_{0}+\delta
e_{H}^{\prime})|\psi_{0}(R_{0}+\delta e_{H}^{\prime})\rangle|^{2}$ between the
true ground state of H${}_{3}^{+}$, $|\psi_{0}(R)\rangle$, and
$|\psi(R)\rangle$ prepared by both non-tailgated and tailgated circuits, at
$R=\delta e_{H}^{\prime}$ for a range of $\delta$, which specifies the
displacement of the hydrogens from equilibrium.

As can be seen in the figure, the tailgating procedure effectively “flattens”
the curve of fidelities between the prepared state and the true state, which
suggests that derivatives of the difference $|\psi(R_{0}+\delta
e_{H}^{\prime})\rangle-|\psi_{0}(R_{0}+\delta e_{H}^{\prime})\rangle$ will be
small. This is exactly the property we desire to compute accurate energy
derivatives, as was explained in Sec. II.

##  V Conclusion

This work outlines possible shortcomings of adaptively-prepared variational
quantum circuits when computing higher-order energy derivatives, and proposes
a solution which relies on the procedure of tailgating. We first discussed the
problem of computing accurate energy derivatives beyond the first-order with
variational quantum circuits, and provide intuition as to why this problem
arises. Specifically, we argued that errors in derivatives appear because the
relative degree of accuracy to which some variational circuits can prepare
correct ground states is large.

Using this intuition, we then propose sufficient conditions for a variational
circuit to yield correct $n$-th order energy derivatives in Theorem 1. The
result we prove about $n$-th order energy derivatives and $(n-1)$-th order
Taylor expansions of $H(R)$ may find more general use in other algorithms
related to high-order energy derivatives. Using Theorem 1, we justify the
tailgating procedure, in which extra gates are added to the end of an already-
optimized variational circuit, and are used when computing derivatives of the
circuit with respect to its parameters.

Finally, we provide numerical examples demonstrating that tailgating yields
better results in certain cases, compared to standard adaptive procedures. We
consider the concrete problem of computing normal mode frequencies of BeH2 and
H2O molecules, which requires knowledge of second-order energy derivatives
with respect to the atomic coordinates of each molecule. Comparing tailgated
and non-tailgated circuits to the true values obtained from classical methods,
tailgating can lead to a substantial increase in the accuracy of the
calculations. Normal mode frequencies are critical for understanding the
vibrational and thermodynamic properties of molecules, as they determine
vibronic spectra, vibrational partition functions, and other observable
quantities. To use variational circuits for chemistry tasks that involve
computing high-order energy derivatives, we argue that a procedure such as
tailgating is valuable to ensure the accuracy of the calculations while
maintaining the lower cost arising from adaptive procedures.

We note that our work does not address other important challenges in quantum
algorithms for quantum chemistry, which still face many obstacles before they
can become competitive with existing classical methods. In particular,
estimating high-order energy derivatives using variational algorithms suffers
from a more pronounced version of the “measurement problem”, which refers to
the often prohibitive number of circuit executions needed to estimate accurate
expectation values. Thus, while we tackle issues regarding the quality and
gate count of variational circuits, other obstacles remain to be overcome.
Ultimately, we hope that tailgating will emerge as a useful member of the
toolbox of techniques that scientists can use to construct variational
circuits for particular classes of problems in chemistry, and beyond.

##  VI Acknowledgements

We thank Nathan Wiebe and John Sipe for valuable conversations relating to
this project. This research was partially funded by a MITACS Accelerate Grant.

## References

  * McArdle et al. (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Reviews of Modern Physics 92, 015003 (2020).
  * Aspuru-Guzik et al. (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
  * Quantum et al. (2020) G. A. Quantum, Collaborators*†, F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, et al., Science 369, 1084 (2020).
  * 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, et al., Nature chemistry 2, 106 (2010).
  * Cao et al. (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, et al., Chemical reviews 119, 10856 (2019).
  * O’Brien et al. (2021) T. E. O’Brien, M. Streif, N. C. Rubin, R. Santagati, Y. Su, W. J. Huggins, J. J. Goings, N. Moll, E. Kyoseva, M. Degroote, et al., arXiv:2111.12437 (2021).
  * Su et al. (2021) Y. Su, D. W. Berry, N. Wiebe, N. Rubin, and R. Babbush, PRX Quantum 2, 040332 (2021).
  * Low and Chuang (2019) G. H. Low and I. L. Chuang, Quantum 3, 163 (2019).
  * 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, et al., Nature Reviews Physics pp. 1–20 (2021).
  * Kandala et al. (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * 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, Nature communications 5, 1 (2014).
  * Higgott et al. (2019) O. Higgott, D. Wang, and S. Brierley, Quantum 3, 156 (2019).
  * Nakanishi et al. (2019) K. M. Nakanishi, K. Mitarai, and K. Fujii, Physical Review Research 1, 033062 (2019).
  * McClean et al. (2017) J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. De Jong, Physical Review A 95, 042308 (2017).
  * Parrish et al. (2019) R. M. Parrish, E. G. Hohenstein, P. L. McMahon, and T. J. Martínez, Physical review letters 122, 230401 (2019).
  * Delgado et al. (2021) A. Delgado, J. M. Arrazola, S. Jahangiri, Z. Niu, J. Izaac, C. Roberts, and N. Killoran, Physical Review A 104, 052402 (2021).
  * Ibe et al. (2020) Y. Ibe, Y. O. Nakagawa, N. Earnest, T. Yamamoto, K. Mitarai, Q. Gao, and T. Kobayashi, arXiv preprint arXiv:2002.11724 (2020).
  * Azad and Singh (2022) U. Azad and H. Singh, Chemical Physics 558, 111506 (2022).
  * Sokolov et al. (2021) I. O. Sokolov, P. K. Barkoutsos, L. Moeller, P. Suchsland, G. Mazzola, and I. Tavernelli, Physical Review Research 3, 013125 (2021).
  * Mitarai et al. (2020) K. Mitarai, Y. O. Nakagawa, and W. Mizukami, Physical Review Research 2, 013129 (2020).
  * Huh et al. (2015) J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean, and A. Aspuru-Guzik, Nature Photonics 9, 615 (2015).
  * Sparrow et al. (2018) C. Sparrow, E. Martín-López, N. Maraviglia, A. Neville, C. Harrold, J. Carolan, Y. N. Joglekar, T. Hashimoto, N. Matsuda, J. L. O’Brien, et al., Nature 557, 660 (2018).
  * Jahangiri et al. (2020) S. Jahangiri, J. M. Arrazola, N. Quesada, and A. Delgado, Physical Chemistry Chemical Physics 22, 25528 (2020).
  * Jahangiri et al. (2021) S. Jahangiri, J. M. Arrazola, and A. Delgado, The Journal of Physical Chemistry Letters 12, 1256 (2021).
  * Stober et al. (2022) S. T. Stober, S. M. Harwood, D. Trenev, P. K. Barkoutsos, T. P. Gujarati, and S. Mostame, Physical Review A 105, 012425 (2022).
  * Grimsley et al. (2019) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nature communications 10, 1 (2019).
  * Tang et al. (2021) H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2, 020310 (2021).
  * Bilkis et al. (2021) M. Bilkis, M. Cerezo, G. Verdon, P. J. Coles, and L. Cincio, arXiv preprint arXiv:2103.06712 (2021).
  * Claudino et al. (2020) D. Claudino, J. Wright, A. J. McCaskey, and T. S. Humble, Frontiers in Chemistry 8, 1152 (2020).
  * Arrazola et al. (2021a) J. M. Arrazola, O. Di Matteo, N. Quesada, S. Jahangiri, A. Delgado, and N. Killoran, arXiv preprint arXiv:2106.13839 (2021a).
  * Arrazola et al. (2021b) J. M. Arrazola, S. Jahangiri, A. Delgado, J. Ceroni, J. Izaac, A. Száva, U. Azad, R. A. Lang, Z. Niu, O. Di Matteo, et al., arXiv:2111.09967 (2021b).
  * Bergholm et al. (2018) V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, M. S. Alam, S. Ahmed, J. M. Arrazola, C. Blank, A. Delgado, S. Jahangiri, et al., arXiv preprint arXiv:1811.04968 (2018).
  * Barca et al. (2020) G. M. J. Barca, C. Bertoni, L. Carrington, D. Datta, N. De Silva, J. E. Deustua, D. G. Fedorov, J. R. Gour, A. O. Gunina, E. Guidez, et al., The Journal of Chemical Physics 152, 154102 (2020), ISSN 0021-9606, 1089-7690, URL http://aip.scitation.org/doi/10.1063/5.0005188.

##  Appendix A Proof of Theorem 1

In this section, we present a proof of Theorem 1, which was introduced in the
main text. We begin by proving a lemma. Note in the proof that follows, for
notational convenience, we define

$\Bigg{|}\frac{\partial^{n}\phi(R_{0})}{\partial R_{j_{1}}\cdots\partial
R_{j_{n}}}\Bigg{\rangle}:=\frac{\partial^{n}|\phi(R_{0})\rangle}{\partial
R_{j_{1}}\cdots\partial R_{j_{n}}}.$ (22)

######  Lemma 1.

Let $H(R)$ be a parameterized Hamiltonian, $R_{0}$ a set of parameters, and
$|\phi(R)\rangle$ a parameterized state vector. Define

$C(R)=\langle\phi(R)|H(R)|\phi(R)\rangle.$ (23)

If $|\phi(R_{0})\rangle$ is an eigenstate of $H(R_{0})$, then any $n$-th order
derivative of $C$ at $R_{0}$ can be written as a sum of inner products
involving $0$-th to $(n-1)$-th order derivatives of $|\phi(R)\rangle$
evaluated at $R_{0}$.

###### Proof.

It can be checked that the only $n$-th order derivative of $|\phi(R)\rangle$
in an $n$-th order derivative of $C$ appears in a term of the form

$2\text{Re}\left\langle\ \frac{\partial^{n}\phi(R_{0})}{\partial
R_{j_{1}}\cdots\partial
R_{j_{n}}}\Bigg{|}H(R_{0})\Bigg{|}\phi(R_{0})\right\rangle.$ (24)

Hence, the problem of writing an $n$-th order derivative in terms of $0$-th to
$(n-1)$-th order derivatives of $|\phi(R)\rangle$ is reduced to writing the
above expression in terms of such derivatives. We know that

$H(R_{0})|\phi(R_{0})\rangle=C(R_{0})|\phi(R_{0})\rangle,$ (25)

as $|\phi(R_{0})\rangle$ is an eigenvector. Hence,

$\displaystyle 2\text{Re}\left\langle\ \frac{\partial^{n}\phi(R_{0})}{\partial
R_{j_{1}}\cdots\partial
R_{j_{n}}}\Bigg{|}H(R_{0})\Bigg{|}\phi(R_{0})\right\rangle=$ $\displaystyle
2\text{Re}\ C(R_{0})\left\langle\ \frac{\partial^{n}\phi(R_{0})}{\partial
R_{j_{1}}\cdots\partial R_{j_{n}}}\Bigg{|}\psi(R_{0})\right\rangle.$ (26)

Finally, since $|\phi(R)\rangle$ is normalized, we have
$\langle\phi(R)|\phi(R)\rangle=1$ for all $R$, which implies that any $k$-th
order derivative, for $k\geq 1$, of this inner product will be equal to $0$.
It follows that

$\displaystyle\frac{\partial^{n}\langle\phi(R_{0})|\phi(R_{0})\rangle}{\partial
R_{j_{1}}\cdots\partial R_{j_{n}}}$ $\displaystyle=2\text{Re}\left\langle\
\frac{\partial^{n}\phi(R_{0})}{\partial R_{j_{1}}\cdots\partial
R_{j_{n}}}\Bigg{|}\phi(R_{0})\right\rangle+S_{n-1}(\phi(R_{0}))=0,$ (27)

where $S_{n-1}(\phi(R_{0}))$ is a sum of inner products involving $0$-th to
$(n-1)$-th order derivatives of $|\phi(R_{0})\rangle$. Rearranging, and using
Eq. (26) and Eq. (27), gives that the expression in Eq. (24) is equal to
$-2\text{Re}[C(R_{0})S_{n-1}(\phi(R_{0}))]$: an expression involving only
$0$-th to $(n-1)$-th order derivatives of $|\phi(R_{0})\rangle$. This
completes the proof. ∎

######  Corollary 1.

If for a circuit $U$, the corresponding $p$-th order derivatives of
$|\psi(R)\rangle$ match those of $|\psi_{0}(R)\rangle$ for $0\leq p\leq n-1$
at $R_{0}$, then $U$ can yield accurate $n$-th order energy derivatives at
$R_{0}$.

###### Proof.

Recall the definitions of $|\psi(R)\rangle$ and $\tilde{E}(R)$ from Sec. II.
Since, when $n\geq 1$, both $|\psi(R)\rangle$ and $|\psi_{0}(R)\rangle$
satisfy the conditions of Lemma 1, it follows that both
$\frac{\partial^{p}\tilde{E}(R_{0})}{\partial R_{j_{1}}\cdots\partial
R_{j_{p}}}$ and $\frac{\partial^{p}E(R_{0})}{\partial R_{j_{1}}\cdots\partial
R_{j_{p}}}$ can be written in terms of inner products of $0$-th to $(p-1)$-th
order derivatives of $|\psi(R)\rangle$ and $|\psi_{0}(R)\rangle$ respectively.
By assumption, we have

$\frac{\partial^{p}|\psi(R_{0})\rangle}{\partial R_{j_{1}}\cdots\partial
R_{j_{p}}}=\frac{\partial^{p}|\psi_{0}(R_{0})\rangle}{\partial
R_{j_{1}}\cdots\partial R_{j_{p}}}$ (28)

for all $p$ from $1$ to $n-1$, and all choice of coordinates, $\\{R_{j_{1}},\
\ldots,\ R_{j_{p}}\\}$. Thus, the representations of the $n$-th order
derivatives of $\tilde{E}$ and $E$ in terms of such derivatives must be equal
as well. ∎

######  Lemma 2.

Let $|\phi(R)\rangle$ be the ground state of $H_{n-1}(R)$. The derivatives of
$|\psi_{0}(R)\rangle$ and $|\phi(R)\rangle$ agree to order $n-1$ when
evaluated at $R_{0}$.

###### Proof.

Note that

$H(R)=H_{n-1}(R)\\\
+\frac{1}{n!}\displaystyle\sum_{j_{1},\ldots,j_{n}=1}^{N}\frac{\partial^{n}H(R_{0})}{\partial
R_{j_{1}}\cdots\partial
R_{j_{n}}}(R_{j_{1}}-R_{0j_{1}})\cdots(R_{j_{n}}-R_{0j_{n}})\\\ +\cdots$ (29)

If $H(R)$ is a non-degenerate Hamiltonian with normalized eigenvectors
$|v_{j}(R)\rangle$ and corresponding eigenvalues $\lambda_{j}(R)$, then

$\frac{\partial|v_{j}(R)\rangle}{\partial R_{k}}=\displaystyle\sum_{i\neq
j}\frac{\left\langle v_{j}(R)\left|\frac{\partial H(R)}{\partial
R_{k}}\right|v_{i}(R)\right\rangle}{\lambda_{j}(R)-\lambda_{i}(R)}|v_{i}(R)\rangle$
(30)

(see Appendix B), and

$\frac{\partial\lambda_{j}(R)}{\partial R_{k}}=\left\langle
v_{j}(R)\left|\frac{\partial H(R)}{\partial
R_{k}}\right|v_{j}(R)\right\rangle,$ (31)

from the Feynman-Hellmann theorem. It follows from Eq. (30) that
$\frac{\partial^{n}|v_{j}(R)\rangle}{\partial R_{k_{1}}\cdots\partial
R_{k_{n}}}$ can be written as a sum of terms involving $(n-1)$-th order
derivatives of the eigenvectors and eigenvalues, and $n$-th order derivatives
of the Hamiltonian. In addition, it follows from Eq. (31) that $(n-1)$-th
order derivatives of the eigenvalues can be written in terms of $(n-1)$-th
order derivatives of the Hamiltonian and $(n-2)$-th order derivatives of the
eigenvectors.

Thus, continuing inductively, it follows that
$\frac{\partial^{n-1}|v_{j}(R)\rangle}{\partial R_{k_{1}}\cdots\partial
R_{k_{n-1}}}$ can be written in terms of the eigenvectors $|v_{i}(R)\rangle$,
eigenvalues $\lambda_{i}(R)$, and $0$-th to $(n-1)$-th order derivatives of
the Hamiltonian.

Since $H(R_{0})=H_{n-1}(R_{0})$, so their eigenvectors and eigenvalues at
$R_{0}$ are the same, and from Eq. (29) the derivatives of the two
Hamiltonians agree to the $(n-1)$-th order at $R_{0}$, it follows that the
derivatives of $|\psi_{0}(R)\rangle$ and $|\phi(R)\rangle$ agree to order
$n-1$ at $R_{0}$. ∎

Now, we are able to prove the desired result.

Theorem 1. Let $|\phi(R)\rangle$ be the ground state of $H_{n-1}(R)$. Suppose
that there exists $\delta>0$ such that a circuit $U(\theta)$ can prepare the
state $|\phi(R)\rangle$ for all $R$ such that $0<\|R-R_{0}\|<\delta$. Then it
holds that

$\frac{\partial^{p}E(R_{0})}{\partial R_{j_{1}}\cdots\partial
R_{j_{p}}}=\frac{\partial^{p}\tilde{E}(R_{0})}{\partial
R_{j_{1}}\cdots\partial R_{j_{p}}}$ (32)

for all $p$ such that $0\leq p\leq n$, and all choices of coordinates,
$\\{R_{j_{1}},\ ...,\ R_{j_{p}}\\}$.

###### Proof.

Recall that $|\psi(R)\rangle$ is the state such that $U(\theta(R))$ closest
approximates $|\psi_{0}(R)\rangle$ (with respect to the vector $2$-norm).
Thus, we have that

$\||\psi(R)\rangle-|\psi_{0}(R)\rangle\|\leq\||\phi(R)\rangle-|\psi_{0}(R)\rangle\|$
(33)

when $0\leq||R-R_{0}||<\delta$, as $U$ can also prepare $|\phi(R)\rangle$. It
follows immediately that the derivatives of $|\psi(R)\rangle$ and
$|\psi_{0}(R)\rangle$ will agree to order $n-1$ at $R_{0}$, as those of
$|\phi(R)\rangle$ and $|\psi_{0}(R)\rangle$ agree. We then apply Corollary 1
to arrive at the result. ∎

##  Appendix B Derivatives of Eigenvectors

In this section, we briefly outline how to compute derivatives of eigenvectors
corresponding to a non-degenerate Hamiltonian. Suppose $H(R)$ is a
Hamiltonian, such that

$H(R)|v_{i}(R)\rangle=\lambda_{i}(R)|v_{i}(R)\rangle,$ (34)

where we assume $\lambda_{i}(R)$ and $|v_{i}(R)\rangle$ are smooth functions
such that the set of eigenvectors $|v_{i}(R)\rangle$ forms an orthonormal
basis, which is guaranteed by the spectral theorem. Taking the derivative of
both sides of Eq. (34) an re-arranging, we get

$\left[\lambda_{i}(R)-H(R)\right]\frac{\partial|v_{i}(R)\rangle}{\partial
R_{k}}=\left[\frac{\partial H(R)}{\partial
R_{k}}-\frac{\partial\lambda_{i}(R)}{\partial R_{k}}\right]|v_{i}(R)\rangle.$
(35)

Taking the inner product of both sides of the above equation with
$|v_{j}(R)\rangle$, it follows that

$(\lambda_{i}(R)-\lambda_{j}(R))\Big{\langle}v_{j}(R)\Big{|}\frac{\partial
v_{i}(R)}{\partial R_{k}}\Big{\rangle}=\\\
\Big{\langle}v_{j}(R)\Big{|}\frac{\partial H(R)}{\partial
R_{k}}\Big{|}v_{i}(R)\Big{\rangle}-\frac{\partial\lambda_{i}(R)}{\partial
R_{k}}\langle v_{j}(R)|v_{i}(R)\rangle.$ (36)

In the case that $i\neq j$, then $\langle v_{j}(R)|v_{i}(R)\rangle=0$, and we
can re-arrange to get

$\Big{\langle}v_{j}(R)\Big{|}\frac{\partial v_{i}(R)}{\partial
R_{k}}\Big{\rangle}=\frac{\Big{\langle}v_{j}(R)\Big{|}\frac{\partial
H(R)}{\partial
R_{k}}\Big{|}v_{i}(R)\Big{\rangle}}{\lambda_{i}(R)-\lambda_{j}(R)}$ (37)

In the case that $i=j$, then since $\langle v_{i}(R)|v_{i}(R)\rangle=1$ for
all $R$, we will get

$\frac{\partial}{\partial R_{k}}\langle
v_{i}(R)|v_{i}(R)\rangle=2\text{Re}\left[\Big{\langle}v_{i}(R)\Big{|}\frac{\partial
v_{i}(R)}{\partial R_{k}}\Big{\rangle}\right]=0$ (38)

In general, we can assume that $\Big{\langle}v_{i}(R)\Big{|}\frac{\partial
v_{i}(R)}{\partial R_{k}}\Big{\rangle}$ is real, as we can always multiply
either state by a non-physical global phase. Hence, this inner product will be
equal to $0$. Therefore, since the $|v_{j}(R)\rangle$ form a basis, we will
have

$\frac{\partial|v_{i}(R)\rangle}{\partial R_{k}}=\displaystyle\sum_{j\neq
i}\frac{\Big{\langle}v_{j}(R)\Big{|}\frac{\partial H(R)}{\partial
R_{k}}\Big{|}v_{i}(R)\Big{\rangle}}{\lambda_{i}(R)-\lambda_{j}(R)}|v_{j}(R)\rangle.$
(39)