Theory of analytical energy derivatives for the variational quantum eigensolver

The variational quantum eigensolver (VQE) and its variants, which is a method for finding eigenstates and eigenenergies of a given Hamiltonian, are appealing applications of near-term quantum computers. Although the eigenenergies are certainly important quantities which determines properties of a given system, their derivatives with respect to parameters of the system, such as positions of nuclei if we target a quantum chemistry problem, are also crucial to analyze the system. Here, we describe methods to evaluate analytical derivatives of the eigenenergy of a given Hamiltonian, including the excited state energy as well as the ground state energy, with respect to the system parameters in the framework of the VQE. We give explicit, low-depth quantum circuits which can measure essential quantities to evaluate energy derivatives, incorporating with proof-of-principle numerical simulations. This work extends the theory of the variational quantum eigensolver, by enabling it to measure more physical properties of a quantum system than before and to explore chemical reactions.


Introduction
The variational quantum eigensolver (VQE) has attracted much attention as a potential application of near-term quantum computers [1].The VQE is an iterative algorithm to construct a quantum circuit that outputs eigenstates and eigenenergy of a Hamiltonian which describes the system under consideration.Originally, the method was devised for finding the ground state of a system [1].It has subsequently been extended for excited states by several proposals [2][3][4][5].From the generated eigenstates, one can measure its associated physical quantities, such as the particle densities and transition amplitudes between the different eigenstates.
Though eigenenergy and associated particle density are certainly important quantities, the wavefunction of a quantum system has valuable information besides those.Quantum chemistry calculations, which would be one of the most promising applications of the VQE, often utilize such information.Among such, we focus on energy derivatives in this work.Many time-independent physical/chemical properties can be defined as derivatives of the energy [6][7][8][9].For example, the first derivatives of the energy with respect to nuclear coordinates give us the forces acting on atoms, which can be utilized for the task of locating energy extrema on the potential energy surface (i.e., geometry optimization) [10].The second-order derivatives give the force-constant matrix that not only helps to locate and verify transition states but also allows us to compute vibrational frequencies and partition functions within the harmonic approximation [10,11].The energy derivatives with respect to external fields have to do with various spectroscopy: Intensities of infrared and Raman spectroscopy are proportional to the cross derivatives with respect to vibrational normal modes and external electric fields [11]; NMR chemical shifts can be obtained using the cross derivatives with respect to nuclear spin and magnetic fields [12,13].Computing such derivatives is a core part of simulations or analysis of molecular spectra.If quantum chemical methods were unable to extract those properties, presumably they would not be widely used as it is today.
A simple way to compute energy derivatives is to use the finite difference method and calculate them numerically.This approach, however, suffers from high computational costs as well as numerical errors and instabilities [7,9].Say, the number of energy points needed to evaluate the forces increases linearly to the number of atoms.This high computational cost makes the numerical approach impractical in many cases.Moreover, with near-term quantum devices, where noise is inevitable, the numerical difference approach would give us poor results.The other way -analytical approach -is therefore vital.The theory and program codes of analytical energy derivatives indeed support the high practicality (and popularity) of today's molecular electronic structure theory [9].Methods to calculate the derivatives of excited state energy on classical computers has also been widely developed , but still suffers from its high computational cost and relatively low accuracy.The task of computing such derivatives on quantum computers has been addressed in the traditional methods which utilize the quantum phase estimation [35], but not for the VQE.
In this work, we derive the analytical formulae and explicit quantum circuits to address the task of measuring the energy derivatives.More specifically, we describe the methods to obtain the derivatives of the energy with respect to the system parameters up to the third order, from which one can extract the physical properties.We also present a method to extract the derivatives of excited state energy based on the technique presented in Refs.[3,4].Proof-of-principle numerical simulations are also performed to verify the correctness of the derived equation and circuits.The presented methods extend the applicability of the VQE by enabling it to evaluate more physical properties than before.

Variational quantum eigensolver
Here, we briefly review the algorithm of the VQE.In the VQE, we construct a parameterized quantum circuit U (θ) and the corresponding ansatz state |ψ(θ) = U (θ) |0 ⊗n , where |0 ⊗n is an initialized n-qubit state and θ = (θ 1 , . . ., θ N θ ) ∈ R N θ is a vector of parameters implemented on the circuit with N θ being the number of them.A set of parameters θ are variationally optimized so that the expectation value E(θ) = ψ(θ)| H |ψ(θ) of a given Hamiltonian H is minimized.At the optimal point, one naturally expects for all a.We denote such an optimal point by θ * .Let The condition of Eq. ( 1) reduces to Higher derivatives of the wave function will be denoted by to shorten the notation.
As stated in the introduction, many physical properties of a quantum system are calculated from the derivative of the energy with respect to the system parameter in a given Hamiltonian.The parameter can be, for example, the coordinates of atoms or the electric field applied to the system.We denote such parameters by an N x -dimensional vector x ∈ R Nx .In this case, both of the Hamiltonian H and the optimal parameter θ * of the wave function at the specific value of x is also a function of x, which will be denoted by H(x) and θ * (x), respectively.Corresponding to the change of this problem setting, we redefine the energy as, ( Let the optimal ground state energy be E * (x), that is, E * (x) = E(θ * (x), x).In the following sections, we show the analytical forms of the derivatives such as ∂E * ∂x i and ∂ 2 E * ∂x i ∂x j , which are the essential quantities for extracting physical properties of the target system.

Analytical expression of derivatives
The derivation of the formulae presented in this section is in the Appendix for completeness, or you can also refer to Ref. [36].In the following, we use the following notation, and likewise for terms similar to this.

Derivatives of ground state energy
The analytical expressions for the derivatives of ground state energy are the following, where we assumed ∂E(θ * (x),x) ∂θ = 0. Note that, in general, the formulae for the d-th derivative of E * (x) contains θ-detrivatives up to the d-th in the form of ∂ q ∂θ q ∂ d−q E ∂x d−q for q = 1, . . ., d.Also, Wigner's (2n + 1)-rule [36] ensures that x-derivatives of the optimal parameter, θ * (x), up to the n-th are sufficient for calculating (2n + 1)-th derivative of E * (x).In other words, for d-th derivative of E * (x) we only need d/2 -th derivative of θ * (x), where y is the floor function denoting the greatest integer less than or equal to y.The term ∂θ * a (x) ∂x i in the above equation can be calculated by solving the response equation, which we write down in the next subsection.

Derivatives of optimal parameters
The first and second derivatives of the optimal parameter, θ * (x), can be obtained from the following response equation. where, By Wigner's (2n + 1)-rule [36], the above equations are enough to obtain derivatives of the energy up to the 5th-order.
4 Measurement and calculation of derivatives of ground state In this section, we describe the methodology for calculating the derivatives of the ground state energy, whose analytical expressions are shown in the previous section, when given an optimal circuit parameter θ * (x) at some x that gives the local minimum of E(θ, x).

Notations and assumptions
In the VQE, we target a Hamiltonian which acts on an n-qubit system and is decomposed into a sum of Pauli strings, P = {I, X, Y, Z} ⊗n , as follows, where h P (x) ∈ R. h P (x) is assumed to be non-zero only for O(poly(n)) terms.For a quantum chemistry problem, the original Hamiltonian has the following form, where c † i and c i are the fermion creation and annihilation operators acting on the i-th orbital, respectively.Eq. ( 14) is converted to the form of Eq. (13) by, for example, Jordan-Wigner or Braviy-Kitaev transformation [37,38].Note that because we always work in the second quantization formalism, the effect of the change of the molecular orbital corresponding to the change of molecular geometry is totally absorbed in the coefficients h(x).Therefore, the change of the molecular orbital does not explicitly appear in the following discussion.
To calculate the energy derivatives of such Hamiltonian, first of all, we assume that the derivatives of Hamiltonian, ∂H ∂x i , ∂ ∂x i ∂H ∂x j , and so on, can be calculated by the classical computer, e.g., the conventional libraries of quantum chemistry calculation.In other words, we are able to calculate the quantities such as ∂h P (x) ∂x j .For quantum chemisry problems given in terms of Hartree-Fock orbitals, these calculations correspond to solving the coupled perturbed Hartree-Fock equation [30,36].
Notice that under this assumption, we only need to consider how to evaluate quantities which involve the differentiation with respect to θ such as , can be evaluated with the exactly same procedure as the usual VQE, i.e., we can measure the expectation value of each Pauli term which appears in ∂x k .After the measurement of all quantities which appear in Eqs. ( 7), (8), and ( 9), one can compute the energy derivative using a classical computer by summing up the terms.
The parameterized quantum state, |ψ(θ) , is constructed by applying a parameterized unitary matrix, that is, a parameterized quantum circuit, U (θ) to an initialized state; |ψ(θ) = U (θ) |0 ⊗n .Following [39], we assume U (θ) to be a product of unitary matrices each with one parameter, We define each unitary U a (θ a ) to be generated by a generator G a as U a = e iθaGa , which can be decomposed into a sum of Pauli strings; where g a,µ ∈ R and P a,µ ∈ P. We often use this form of parametrized quantum circuits, for example see Ref. [40][41][42][43].

Overview of the algorithm
The presented algorithm in this work for evaluating the d-th derivative is the following: 1. Perform the VQE and obtain the optimal parameter θ * (x).
Let us recall N H and N θ be the number of terms in the target Hamiltonian and the number of the VQE parameters.The cost to evaluate d-th derivatives by the above algorithm is roughly O(N H N d θ ) when we ignore the cost of the task performed on a classical computer.This is because the most time-consuming part of the algorithm is Step 3 and 4 on a quantum computer, which take O(N H N d θ ) as shall be clear in Secs.4.3-4.5 and 5.

Measurement of
The key quantities for evaluating ∂ ∂x i ∂E * ∂x j and higher order derivatives are Note that the first order derivative, ∂E ∂θa can be evaluated by the method presented in Ref. [44].First, we show how to measure |∂ a ψ(θ) can be expressed as, and For convinience, we define, Then, In Fig. 1 (a), we show quantum circuits for evaluation of each of the terms in the above equation, which is a variant of the circuit presented in Ref. [39].From the outputs of the circuits in Fig. 1 (a), each term of Eq. ( 21) can be calculated by, The strategy proposed in Ref. [45] gives low-depth versions of the circuits to measure the same quantities.We show the low-depth quantum circuit in Fig. 1 (b).From the measurement of Q (a,µ,±),(b,ν,±) with the circuit in Fig. 1 (b), we can evaluate each term of Eq. ( 21) by the following formula, Notice that the low-depth version doubles the number of the circuit runs, and therefore, if the device can execute the circuit in Fig. 1 (a) while maintaining sufficient overall fidelity, it is advantageous to utilize Fig. 1 (a) and Eq.(22).Otherwise the low-depth version should be used to obtain meaningful result.
This circuit is a variant of the one presented in Ref. [39].In the figure, Z anc is the Pauli Z operator acting only on the ancilla qubit, and s ∈ {0, 1}.See Eq. ( 22) for detailed procedure.(b) Low-depth version of (a), derived with the strategy presented in Ref. [45].In the figure , R ± a,µ = exp(±iπP a,µ /4).See Eq. ( 23) for detailed procedure.
Each term can be evaluated with the circuit in Fig. 2 (a): The circuits can also be reduced to the low-depth version using the same strategy [45].The low-depth circuit for The circuit can evaluate each term of the summation by the following formula.Notice the pattern in the signs associating Q (a,µ,±),(b,ν,±),(c,ρ,±) , that is, the signs are determined in the parity of ± appearing in the subscript.Higher order derivatives of E with respect to θ can be evaluated by following the strategy described in this and the previous section.
Similar to the previous case, the low-depth version doubles the number of the circuit runs, and therefore, the same argument about the tradeoff between the noise and the number of the circuit runs also applies in this case.9), can be measured using the same strategy.

Computational cost
Here, we give a comparison between the presented algorithm and the numerical differentiation to take d-th energy derivatives.We specifically discuss the cost for each step in the overall algorithm presented in 4.2 for quantum chemistry problem, where a problem Hamiltonian and its derivatives always has O(n 4 ) terms as expressed by Eq. ( 14), as it is the main target of the proposed algorithm.This analysis should be easy to be extended to more general cases.We safely assume that classical computational cost involved in algorithms is much smaller than the quantum ones and thus can be ignored, since usual classical computation is much faster than sampling from the NISQ devices.
First, we ignore the cost of performing the VQE (Step 1), since both of the numerical differentiation and the presented algorithm need this step.Step 2, 5, and 6 is merely classical, and therefore, we also ignore the cost for this part.Step 3 of the presented algorithm requires us to run a quantum computer O(n 4 / 2 ) times to estimate each term up to an additional error of .For Step 4, it takes O(n ) runs.Note that even for the calculation of the terms involving x-and θ-differentiation a the same time, such as ∂ ∂θa ∂E ∂x i , we only need to perform measurement of the θ-derivatives of each term, such as , in the Hamiltonian on the quantum device and then combine them by ijkl where x-differentiated h(x) is computed classically.Therefore, N x does not appear in the number of NISQ runs.Overall, the cost for the quantum part of the computation scales as O(n 4 N d θ / 2 ).Let us compare this cost against the numerical differentiation.More specifically, we consider the case of evaluating the second derivative by the formula, where h > 0. Let Ẽ * (x) be the estimate of E * (x) obtained by measuring the state |ψ(θ | ≤ E , the precision of Eq. ( 27) is [46], For classical computation where the source of E is mainly the round-off error, the second term of the right-hand side is usually negligibly small for a decent h.On the other hand, for the VQE, this term is the leading factor for the precision, since Ẽ * (x) has to be calculated by sampling.To achieve | Ẽ * (x) − E * (x)| ≤ E with high probability, we need to run the quantum computer for O(1/ 2 E ) times.Therefore, if we want to achieve the precision of for , and overall cost is O(n 4 N 2 x /( 2 h 4 )).This scaling with respect to the error and the dimension of x is clearly worse than that of the presented method, whose scaling is O(1/ 2 ).

Derivatives of excited state energy
The generation of excited states can be a powerful application of the VQE, because the classical computation, despite the recent significant improvement in the theory and the computational power, still suffers in the calculation of them [47].Among the several proposals [2][3][4][5] to generate excited states with the VQE, we adopt the one proposed in Refs.[3,4] to compute the derivatives of the excited state energy.
The algorithm [3,4] works as follows.First, we find the ground state of the given Hamiltonian H 0 (x), which we denote by |ψ (0) (θ (0) (x)) , where θ (0) (x) is the optimal parameter for the ground state (θ * (x) in the previous sections).Then, we iteratively define a Hamiltonian where |ψ (r) (θ (r) (x)) is the ground state of H r (x).If β s is sufficiently large, each H r (x) has r-th excited state of the original Hamiltonian, H 0 (x), as its ground state.Therefore, by finding the ground state of each H r (x), one can generate the series of excited states of H 0 (x).We assume |ψ (r) (θ) = U (r) (θ) |0 , where U (r) (θ) has the same structure as U (θ) in previous sections In this algorithm, it is also assumed that we have a device which can measure the overlap between |ψ (r) (θ) and |ψ (s) (ϕ) , that is, we assume we can measure | ψ (r) (θ)|ψ (s) (ϕ) | 2 .Let the expectation value of H r (x) with respect to the state |ψ (r) (θ) be E r (θ, x); E r (θ, x) = ψ (r) (θ)|H r (x)|ψ (r) (θ) .We define the optimal energy by E * r (x) = E r (θ (r) (x), x).The task is to compute the derivatives such as ∂ 2 E * r ∂x i ∂x j .Since E * r is the ground state energy for H r , the formulae in the previous sections, Eqs. ( 7), (8), and ( 9) can be adapted for this task.The only difference from that of the actual ground state is that the derivative of the Hamiltonian, such as ∂ 2 Hr ∂x i ∂x j , cannot be computed classically.For example, the expression of the first derivative of the Hamiltonian is, In the analytical expression for ∂E * r ∂x i (Eq.( 7)), ∂Hr ∂x i (x) appears as the expectation value with respect to |ψ (r) (θ (r) (x)) ; ψ (r) (θ (r) (x))| ∂Hr ∂x i (x)|ψ (r) (θ (r) (x)) .As for ∂ 2 E * r ∂x i ∂x j (Eq.( 8)), it appears as Re ψ (r) (θ (r) (x))| ∂Hr ∂x i (x)|∂ a ψ (r) (θ (r) (x)) .The quantities that cannot be computed classically in Eqs. ( 7) to ( 9) are the terms which involves the inner product between the states, such as Re[ However, notice that if the condition ψ (r) (θ (r) (x))|ψ (s) (θ (s) (x)) = 0 holds for all x, which we naturally expect at the optimal parameter, we obtain by differentiating with respect to x the both hand side of the equation, We can therefore ignore the inner-product terms for the evaluation of derivatives of excited state energy and utilize the same procedure as the ground state energy in this case.Likewise, the inner-product terms that appear in the higher-order derivatives can also be ignored when the orthogonality condition ψ (r) (θ (r) (x))|ψ (s) (θ (s) (x)) = 0 is satisfied.The so-called subspace search VQE method [2], which guarantees the orthogonality condition of the resultant states, can be advantageous to obtain such condition.

Numerical simulation
We provide proof-of-principle numerical simulations using the electronic Hamiltonian of the hydrogen molecule.The Hamiltonians are calculated with the open source library PySCF [48] and OpenFermion [49].The simulation of quantum circuits are performed with Qulacs [50].

Approximation of the potential energy surface
First, we perform a simulation of the VQE and the methods described in Sec. 4 with the Hamiltonian of a hydrogen molecule, calculated with the STO-3G basis set, at the bonding distance r = 0.735 Å.We use the ansatz circuit shown in Fig. 3, which is a variant of so-called hardware efficient circuit [41].The result of the simulation is shown as Fig. 4. The ansatz  used in this simulation could achieve the exact ground state, which is called full configuration interaction (Full-CI) state in the context of chemistry, and therefore, we could draw the harmonic and the third-order approximation of the Full-CI energy as shown in Fig. 4. The harmonic approximation can be used to calculate the vibrational spectra of a molecule, and the third order approximation can be utilized for more accurate description of the vibration.

Continuous determinination of the optimal parameter
The response equation, Eq. ( 10), can be used to determine the optimal paramter θ * (x) from the one at the slightly different system parameter, θ * (x + δx), that is, to the first order, holds up to the additive error of O(δx 2 ).One can iteratively use the equation above, which resembles the Euler method, to obtain the optimal parameter θ Euler (x) ≈ θ * (x) from some θ * (x 0 ) in a range of x around x 0 where the error term is sufficiently small.We demonstrate this parameter update in Fig. 5, by plotting the energy expectation value of |ψ(θ Euler (r)) where θ Euler (r) is determined with the iterative use of Eq. ( 32), starting from the optimal parameter at r = 1.5 Å with the step size δr = 0.02 Å. Eq. ( 32) resembles the Euler method, therefore we refer to the update of θ according to Eq. ( 32) as so in Fig. 5.It can be seen that in the range where r is sufficiently close to the optimized point, r = 1.5 Å, the parameter determined from Eq. ( 32) is almost optimal, but as r goes far from the optimized point, they deviate fast from the optima.This is because, as the error of the each update accumulates, the state becomes a non-variational state, that is, ∂E ∂θ = 0, which results in the breakdown of the response equation, Eq. ( 10), where we assumed ∂E ∂θ = 0.This method should be useful for drawing potential energy surfaces using the VQE, because it can reduce the time for the optimization of the circuit parameter by predicting the optimal parameter from the one determined with the slightly different system.We belive it can also be combined with the parameter interpolation approach proposed in Ref. [51].

Conclusion
We have described a methodology for computing the derivatives of the ground state and the excited state energy.We have shown the straight-forwardly constructed quantum circuit to measure the quantities necessary for the calculation of energy derivatives, and also the low-depth version of the circuit.The low-depth protocol for the measurement makes them suitable for the use in NISQ devices.We also performed numerical simulations of the proposed method, which validates the correctness.This work enables to extract the physical properties of the quantum system under investigation, thus widening the application range of the NISQ devices.
Note added -During the final revision of this article, we have become aware of recent article [52] which also describes a methodology for evaluating energy derivatives.Their work is based on perturbative sum-over-state approach and differs from our methods which provides a way to evaluate analytic derivatives.

A Derivative of the optimal parameter
We derive the expression for the derivatives of the optimal parameter, such as ∂θ * a (x) ∂x i , by Taylor expansion.Assume that at x = α, we know the optimal parameter θ * (α) We perform Taylor expansion of |ψ(θ * (x)) and H(x) around α.For H(x), we obtain, For |ψ(θ * (x)) , we can expand it as follows.
|ψ(θ * (α + x)) When grouped by the order of x, We can derive a similar expression for |∂ a ψ(θ * (α + x)) .Now, we use the condition of Eq. ( 1) to derive the expression of the derivatives of θ * (α).Plugging Eqs. ( 33) and (34) into Eq.( 1) and imposing the coefficient of each order in x to be zero, we get the analytical expression for the derivatives of θ * (α).In the following two subsections, we derive the expression for the first and second derivatives of θ * (α).

A.1 First derivative
For the first order in x, we get Note that the first term of Eq. ( 35) includes the Hessian of E(θ, x), that is, Also, notice that, We define ∇ d θ ∂E ∂x i likewise.Finally, we obtain, which is Eq. ( 10) of the main text.Note that we expect the matrix ∂ ∂θ b ∂E(θ * (α),α) ∂θa to be positive definite because θ * (α) is a local minimum, and thus Eq. ( 38) is solvable to obtain ∂θ * a (α) ∂x i .

A.2 Second derivative
The second derivative, ∂ ∂x i ∂θa ∂x j , is derived from the second order in x of Eq. ( 1).We have, from Eq. ( 1), for all i, j and c, ∂θ b can be used to greatly simplify the above equation, which gives us, This is equivalent to Eq. ( 11) B Derivatives of the ground state energy

B.1 First derivative
The first derivative of the energy is calculated as, ( This expression (Eq.( 43)) can be simplified so as to avoid the explicit calculation of which is Eq.( 9) of the main text.

Figure
Figure 1: Quantum circuit to evaluate ∂
be measured with the same protocol as described in the previous sections; one can substitute h Q with ∂h Q ∂x i .With this substitution, Eqs.(21) and (24) give us the analytical formula for ∂ ∂θa ∂ ∂θ b • • • ∂ ∂θc ∂E ∂x i , where each term in the summation can be evaluated with the same procedure.Also, ∂ ∂θa ∂ ∂x i ∂E ∂x k , which appears in Eq. (

Figure 3 :
Figure 3: Ansatz used in simulations.R x and R y is single-qubit x and y rotation gate, respectively.The parameters θ are implemented as rotation angles of R x and R y .

Figure 4 :
Figure 4: The harmonic and third-order approximation of the energy curve of the hydrogen molecule at the bonding distance, determined by the simulation of the VQE and the proposed method.

Figure 5 :
Figure5: The evolution of the energy expectation value when the circuit parameters are evolved according to Eq. (32), starting from the optimized parameter at interatomic distance of 1.5 Å.
and first term of the above equation vanishes by Eq. (1) of the main text.Thus, we get, Here, we derive the expression for the second derivative.This is Eq.(8) of the main text.