Calculating energy derivatives for quantum chemistry on a quantum computer

TL;DR

This paper introduces two new methods for calculating energy derivatives in quantum chemistry on quantum computers, enabling geometry optimization and property estimation with experimental validation.

Contribution

It presents novel quantum algorithms for energy derivatives, including quantum phase estimation and response approximation, with experimental implementation on a quantum processor.

Findings

Achieved the first quantum-processor-based geometry optimization of H2.

Estimated H2 polarizability with 2% relative error.

Provided error bounds and computational scalings for the methods.

Abstract

Modeling chemical reactions and complicated molecular systems has been proposed as the `killer application' of a future quantum computer. Accurate calculations of derivatives of molecular eigenenergies are essential towards this end, allowing for geometry optimization, transition state searches, predictions of the response to an applied electric or magnetic field, and molecular dynamics simulations. In this work, we survey methods to calculate energy derivatives, and present two new methods: one based on quantum phase estimation, the other on a low-order response approximation. We calculate asymptotic error bounds and approximate computational scalings for the methods presented. Implementing these methods, we perform the world's first geometry optimization on an experimental quantum processor, estimating the equilibrium bond length of the dihydrogen molecule to within 0.014 Angstrom of the full configuration interaction value. Within the same experiment, we estimate the polarizability of the H2 molecule, finding agreement at the equilibrium bond length to within 0.06 a.u. (2% relative error).

Figures

Tables

Table 1. Table 1: Calculated performance of the energy derivative estimation methods suggested in this work. The computation time scaling as a function of the estimation error ϵ italic-ϵ \epsilon is given, alongside an approximate range of scalings with respect to the system size N sys subscript 𝑁 sys N_{\mathrm{sys}} . Calculation details, full approximation details, and intermediate steps are given in Sec.  IV.2 . Highlighted scalings correspond to the methods used in experiment in Sec.  II.3 . † † {\dagger} denotes performance estimates when using a VQE for state preparation, ∗ * denotes performance when using QPE for state preparation, ∗ ⁣ ∗ ** denotes performance if the amplitude amplification technique of  Knill  et al. ( 2007 ) is used (which requires a VQE for state preparation), and ▽ ▽ \triangledown denotes performing differentiation on a quantum computer with the methods of  Kassal and Aspuru-Guzik ( 2009 ); Jordan ( 2005 ) (which requires the ability to call the Hamiltonian as a quantum oracle as a function of the system parameters). Details of errors not accounted for are given in Sec.  A (basis set error), Sec.  IV.3.4 (state error), Sec.  IV.5.3 (resolution error), and Sec.  IV.6 (truncation error). We further note whether methods require phase estimation (which requires long coherence times).

	Hellmann–Feynman		PPE		ETA		Direct
	(first order)
Time scaling	${ϵ}^{-2}$	$†$	${ϵ}^{-2}$	*	${ϵ}^{-2}$	$†$	${ϵ}^{-d+1}$	$†$
with $ϵ$	${ϵ}^{-2}$	*			${ϵ}^{-2}$	*	${ϵ}^{-\frac{d+2}{2}}$	*
(fixed $N_{\mathrm{sys}}$)	${ϵ}^{-1}$	**			${ϵ}^{-1}$	**	${ϵ}^{-1}$	$▽$
Time scaling	$N_{\mathrm{sys}}^4-N_{\mathrm{sys}}^{13}$	$†$	$N_{\mathrm{sys}}^4-N_{\mathrm{sys}}^{17.5}$	*	$N_{\mathrm{sys}}^8-N_{\mathrm{sys}}^{21}$	$†$	$N_{\mathrm{sys}}^4-N_{\mathrm{sys}}^{13}$	$†$
with $N_{\mathrm{sys}}$	$N_{\mathrm{sys}}^4-N_{\mathrm{sys}}^{15}$	*			$N_{\mathrm{sys}}^8-N_{\mathrm{sys}}^{23}$	*	$N_{\mathrm{sys}}^2-N_{\mathrm{sys}}^7$	*
(approx, fixed $ϵ$)	$N_{\mathrm{sys}}^{3.5}-N_{\mathrm{sys}}^{11}$	**			$N_{\mathrm{sys}}^{6.5}-N_{\mathrm{sys}}^{17}$	**	$N_{\mathrm{sys}}^2-N_{\mathrm{sys}}^7$	$▽$
Error sources	Basis set error		Basis set error		Basis set error		Basis set error
unaccounted for	State error	($†$,**)	Resolution error		Truncation error		State error	($†$,**)
Req. $\hat{H}$ sim.	$\ast$, $\ast \ast$		All		$\ast$, $\ast \ast$		$\ast$,$▽$

Table 2. Table 2: Intermediate steps for the approximate calculation of the computational complexity as a function of the system size give in Tab.  1 . Symbols are defined in Sec.  IV.3.1 ( N U subscript 𝑁 U N_{\mathrm{U}} , N T subscript 𝑁 T N_{\mathrm{T}} , n m subscript 𝑛 m n_{\mathrm{m}} ), Sec.  IV.3.2 ( T U subscript 𝑇 U T_{\mathrm{U}} ), Sec.  IV.4 ( K 𝐾 K , A 0 subscript 𝐴 0 A_{0} , t 𝑡 t ), Sec.  IV.5.3 ( Δ Δ \Delta ) and Sec.  IV.3.3 ( T P subscript 𝑇 P T_{\mathrm{P}} ).

	Hellmann–Feynman		PPE		ETA		Direct
	(first order)
Error scaling ($ϵ$)	$N_T^{\frac{1}{2}}{n}_m^{-\frac{1}{2}}$	$†$	$N_U^{\frac{d}{2}}\times$	*	$N_E{N}_T^{\frac{1}{2}}{n}_m^{-\frac{1}{2}}$	$†$	${\left(N_T{n}_m^{-1}\right)}^{\frac{1}{d+1}}$	$†$
	$N_U^{\frac{1}{2}}{N}_m^{-\frac{1}{2}}{A}_0^{-\frac{1}{2}}$	*	$(d^{\frac{d}{2}}{\delta }^{d-2}{\Delta }^{\frac{d}{2}}{K}^{-1}{t}^{-1}{A}_0^{-1}$		$N_E{N}_U^{\frac{1}{2}}{N}_m^{-\frac{1}{2}}{A}_0^{-\frac{1}{2}}$	*	${\left(K{A}_{0}t\right)}^{\frac{-2}{d+2}}$	*
	$N_U^{\frac{1}{2}}{K}^{-1}$	**	$+{\delta }^{d-1}{\Delta }^{\frac{d}{4}}{N}_m^{-\frac{1}{2}}{A}_0^{-\frac{1}{2}})$		$N_E{N}_U^{\frac{1}{2}}{K}^{-1}$	**	${\left(K{A}_{0}t\right)}^{-1}$	$▽$
Time scaling	$N_T{n}_m{T}_P$	$†$	$d{N}_U^{d}K{T}_U$	*	$N_T{N}_E^2{n}_m{T}_P$	$†$	$d{N}_T{n}_m{T}_P$	$†$
	$N_{U}K{T}_U$	*			$N_U{N}_E^{2}K{T}_U$	*	$dK{T}_U$	*
	$N_{U}K{T}_P$	**			$N_U{N}_E^{2}K{T}_P$	**	$2^{d}K{T}_U$	$▽$
Time scaling	$N_T^2{T}_P$	$†$	$\times T_U$	*	$N_T^2{N}_E^4{T}_P$	$†$	$d{N}_T^2{T}_P$	$†$
at fixed $ϵ$	$A_0^{-2}{N}_U^2{T}_U$	*	$(d^{\frac{3}{2}}{\delta }^{d-2}{t}^{-1}{\Delta }^{\frac{d}{2}}{A}_0^{-1}{N}_U^{\frac{3d}{2}-1}$		$N_U^2{N}_E^4{A}_0^{-2}{T}_U$	*	$d{t}^{-1}{A}_0^{-1}{T}_U$	*
	$N_U^{\frac{3}{2}}{T}_P$	**	$+d{\delta }^{2d-2}{\Delta }^{\frac{d}{2}-1}{A}_0^{-\frac{1}{2}}{N}_U^{2d-1})$		$N_U^{3/2}{N}_E^3{T}_P$	**	$2^d{t}^{-1}{A}_0^{-1}{T}_U$	$▽$

Table 3. Table 3: Summary of approximations used to derive the scaling laws with N sys subscript 𝑁 sys N_{\mathrm{sys}} in Tab.  1 .

	Lower bound	Upper bound
$T_P$	$N_{\mathrm{sys}}^2$	$N_{\mathrm{sys}}^5$
$T_U$	$1$	$N_{\mathrm{sys}}^5$
$N_T$	$N_{\mathrm{sys}}$	$N_{\mathrm{sys}}^4$
$N_U$	$N_{\mathrm{sys}}$	$N_{\mathrm{sys}}^4$
$N_E$	$N_{\mathrm{sys}}$	$N_{\mathrm{sys}}^2$
$t^{-1}$	$N_{\mathrm{sys}}$	$N_{\mathrm{sys}}$
$A_0^{-1}$	$N_{\mathrm{sys}}$	$N_{\mathrm{sys}}$
${\delta }^{-1}$	1	1
${\Delta }^{-1}$	1	1

Equations

$$\hat{H}|{\Psi }_j⟩=E_j|{\Psi }_j⟩,$$

$$D_{{\lambda }_1,{\lambda }_2,\dots }^{d_1,d_2,\dots }=\frac{{\partial }^d{E}_0({\lambda }_1,{\lambda }_2,\dots )}{{\partial }^{d_1}{\lambda }_1,{\partial }^{d_2}{\lambda }_2,\dots },$$

$$\hat{H}=\sum _{pq}{h}_{pq}{\hat{E}}_{pq}+\frac{1}{2}\sum _{pqrs}{g}_{pqrs}\left({\hat{E}}_{pq}{\hat{E}}_{rs}-{\delta }_{q,r}{\hat{E}}_{ps}\right),$$

$$\frac{\partial E_0}{\partial \lambda }=⟨{\Psi }_0|\frac{\partial \hat{H}}{\partial \lambda }|{\Psi }_0⟩.$$

$$\frac{{\partial }^2{E}_0}{\partial {\lambda }_1\partial {\lambda }_2}=⟨{\Psi }_0|\frac{{\partial }^2\hat{H}}{\partial {\lambda }_1\partial {\lambda }_2}|{\Psi }_0⟩$$

$$+\sum _{j\ne 0}2\mathrm{Re}\left[⟨{\Psi }_0|\frac{\partial \hat{H}}{\partial {\lambda }_1}|{\Psi }_j⟩⟨{\Psi }_j|\frac{\partial \hat{H}}{\partial {\lambda }_2}|{\Psi }_0⟩\right]\frac{1}{E_0-E_j},$$

$${\mathcal{A}}_1(j)$$

$${\mathcal{A}}_2$$

$$f_1(E_0;E_j)=\frac{2}{E_0-E_j},f_2=1.$$

$$D=$$

$$\times f_{\mathcal{A}}(E_0;E_{j_1},\dots ,E_{j_{X_{\mathcal{A}}-1}}),$$

$$\mathcal{A}(j_1$$

$$=⟨{\Psi }_0|{\hat{P}}_{X_{\mathcal{A}}}\prod _{x=1}^{X_{\mathcal{A}}-1}\left(|{\Psi }_{j_x}⟩⟨{\Psi }_{j_x}|{\hat{P}}_x\right)|{\Psi }_0⟩.$$

$$|{\chi }_j⟩={\hat{E}}_j|{\Psi }_0⟩.$$

$$S_{j,k}^{(\mathrm{QSE})}=⟨{\chi }_j|{\chi }_k⟩,$$

$$H_{j,k}^{(\mathrm{QSE})}=⟨{\chi }_j|\hat{H}|{\chi }_k⟩,$$

$${\hat{H}}^{(\mathrm{QSE})}{\vec{v}}^{(j)}={\tilde{E}}_j{\hat{S}}^{(\mathrm{QSE})}{\vec{v}}^{(j)}\to |{\tilde{\Psi }}_j⟩=\sum _l{\vec{v}}_l^{(j)}|{\chi }_l⟩.$$

$$\frac{{\partial }^{2}E}{\partial {\lambda }^2}$$

$$=\frac{1}{\delta \lambda }\left(\frac{\partial E}{\partial \lambda }(\lambda +\delta \lambda /2)-\frac{\partial E}{\partial \lambda }(\lambda -\delta \lambda /2)\right)+O(\delta {\lambda }^2).$$

$$h_{pq}=\int \mathrm{d}\mathbf{r}{\varphi }_p^{\ast }(\mathbf{r})\left(-\frac{1}{2}{\nabla }_{\mathbf{r}}+V(\mathbf{r})\right){\varphi }_q(\mathbf{r}),$$

$$g_{pqrs}$$

$$\frac{\partial \hat{H}}{\partial \lambda }\approx \frac{\hat{H}(\lambda +\delta \lambda /2)-\hat{H}(\lambda -\delta \lambda /2)}{\delta \lambda },$$

$$\frac{{\partial }^2\hat{H}}{\partial {\lambda }^2}\approx \frac{\hat{H}(\lambda +\delta \lambda )+\hat{H}(\lambda -\delta \lambda )-2\hat{H}(\lambda )}{\delta {\lambda }^2},$$

$$\rho=\left(\begin{array}[]{cc}p_{0}&p_{+}e^{i\phi}\\ p_{+}e^{-i\phi}&p_{1}\end{array}\right),$$

$$2p_{+}{e}^{i\varphi }$$

$$+$$

$$|{\tilde{\Psi }}_0⟩=\sum _j{a}_j|{\Psi }_j⟩.$$

$$2\Vert \hat{H}{\Vert }_{\mathrm{S}}(1-|a_0{|}^2)\ge |E_0-{\tilde{E}}_0|$$

$$\ge \delta (1-|a_0{|}^2)\ge 0,$$

$$|{\tilde{D}}_0-D_0|=|\sum _{i,j>0}{a}_i^{\ast }{a}_j⟨{\Psi }_i|\hat{D}|{\Psi }_j⟩$$

$$+2\mathrm{Re}\sum _j{a}_j^{\ast }{a}_0⟨{\Psi }_j|\hat{D}|{\Psi }_0⟩+(|a_0{|}^2-1)D_0|.$$

$$\sum _{i,j}{a}_i^{\ast }{a}_j\le (\sum _i|a_i{|}^2{)}^{1/2}(\sum _j|a_j{|}^2{)}^{1/2}=(1-|a_0{|}^2),$$

$$|{\tilde{D}}_0-D_0|\le 2\Vert \hat{D}{\Vert }_{\mathrm{S}}\left[(1-|a_0{|}^2)+|a_0|\sqrt{1-|a_0{|}^2}\right].$$

$$|{\tilde{D}}_0-D_0|\le 2\Vert \hat{D}{\Vert }_{\mathrm{S}}\left(\frac{|{\tilde{E}}_0-E_0{|}^{1/2}}{{\delta }^{1/2}}+\frac{|{\tilde{E}}_0-E_0|}{\delta }\right).$$

Full text

I Introduction

Quantum computers are at the verge of providing solutions for certain classes of problems that are intractable on a classical computer Preskill (2018). As this threshold nears, an important next step is to investigate how these new possibilities can be translated into useful algorithms for specific scientific domains. Quantum chemistry has been identified as a key area where quantum computers can stop being science and start doing science Lloyd (1996); Reiher et al. (2017); Cao et al. (2018); McArdle et al. (2018a). This observation has lead to an intense scientific effort towards developing and improving quantum algorithms for simulating time evolution Abrams and Lloyd (1997); Zalka (1998) and calculating ground state energies Aspuru-Guzik et al. (2005); Kitaev (1995); Peruzzo et al. (2014); McClean et al. (2016) of molecular systems. Small prototypes of these algorithms have been implemented experimentally with much success Peruzzo et al. (2014); O’Malley et al. (2016); Kandala et al. (2017); Santagati et al. (2018); Colless et al. (2018). However, advances over the last century in classical computational chemistry methods, such as density functional theory (DFT) Dreizler and Gross (1990), coupled cluster (CC) theory Shavitt and Bartlett (2009), and quantum Monte-Carlo methods Booth et al. (2009), set a high bar for quantum computers to make impact in the field.

The ground and/or excited state energy is only one of the targets for quantum chemistry calculations. For many applications one also needs to be able to calculate the derivatives of the molecular electronic energy with respect to a change in the Hamiltonian Jensen (2007); Norman et al. (2018). For example, the energy gradient (or first-order derivative) for nuclear displacements is used to search for minima, transition states, and reaction paths Schlegel (2011) that characterize a molecular potential energy surface (PES). They also form the basis for molecular dynamics (MD) simulations to dynamically explore the phase space of the system in its electronic ground state Marx and Hutter (2009) or, after a photochemical transition, in its electronically excited state Tully (1990). While classical MD usually relies on force-fields which are parameterized on experimental data, there is a growing need to obtain these parameters on the basis of accurate quantum chemical calculations. One can easily foresee a powerful combination of highly accurate forces generated on a quantum computer with machine learning algorithms for the generation of reliable and broadly applicable force-fields Behler (2016). This route might be particularly important in exploring excited state PES and non-adiabatic coupling terms, which are relevant in describing light-induced chemical reactions Li and Liu (2014); Curchod et al. (2013); Faraji et al. (2018). Apart from these perturbations arising from changing the nuclear positions, it is also of interest to consider the effect that small external electric and/or magnetic fields have on the molecular energy. These determine well-known molecular properties, such as the (hyper)polarizability, magnetizability, A- and g-tensors, nuclear magnetic shieldings, among others.

Although quantum algorithms have been suggested to calculate derivatives of a function represented on a quantum register Jordan (2005); Gilyén et al. (2017); Dallaire-Demers et al. (2018); Schuld et al. (2018); Harrow and Napp (2019), or of derivatives of a variational quantum eigensolver (VQE) for optimization purposes Romero et al. (2018); Guerreschi and Smelyanskiy (2017), the extraction of molecular properties from quantum simulation has received relatively little focus. To the best of our knowledge only three investigations; in geometry optimization and molecular energy derivatives Kassal and Aspuru-Guzik (2009), molecular vibrations McArdle et al. (2018b), and the linear response function Roggero and Carlson (2018); have been performed to date.

In this work, we survey methods for the calculation of molecular energy derivatives on a quantum computer. We calculate estimation errors and asymptotic convergence rates of these methods, and detail the classical pre-and post-processing required to convert quantum computing output to the desired quantities. As part of this, we detail two new methods for such derivative calculations. The first involves simultaneous quantum phase and transition amplitude (or propagator) estimation, which we name ’propagator and phase estimation’ (PPE). The second is based on truncating the Hilbert space to an approximate (relevant) set of eigenstates, which we name the ’eigenstate truncation approximation’ (ETA). We use these methods to perform geometry optimization of the H2 molecule on a superconducting quantum processor, as well as its response to a small electric field (polarizability), and find excellent agreement with the full configuration interaction (FCI) solution.

II Main

Let $\hat{H}$ be a Hamiltonian on a $2^{N_{\mathrm{sys}}}$-dimensional Hilbert space (e.g. the Fock space of an $N_{\mathrm{sys}}$-spin orbital system), which has eigenstates

[TABLE]

ordered by the energies $E_{j}$. In this definition, the Hamiltonian is parametrized by the specific basis set that is used and has additional coefficients $\lambda_{1},\lambda_{2},\ldots$, which reflect fixed external influences on the electronic energy (e.g. change in the structure of the molecule, or an applied magnetic or electric field). An $d$th-order derivative of the ground state energy with respect to the parameters $\lambda_{i}$ is then defined as:

[TABLE]

where $d=\sum_{i}d_{i}$. As quantum computers promise exponential advantages in calculating the ground state $E_{0}$ itself, it is a natural question to ask how to efficiently calculate such derivatives on a quantum computer.

II.1 The quantum chemical Hamiltonian

A major subfield of computational chemistry concerns solving the electronic structure problem. Here, the system takes a second-quantized ab initio Hamiltonian, written in a basis of molecular spinors $\{\phi_{p}(\mathbf{r})\}$ as follows:

[TABLE]

where $\hat{E}_{pq}=\hat{c}_{p}^{\dagger}\hat{c}_{q}$ and $\hat{c}^{{\dagger}}_{p}$ ($\hat{c}_{p}$) creates (annihilates) an electron in the molecular spinor $\phi_{p}$. With equation (3) relativistic and non-relativistic realizations of the method only differ in the definition of the matrix elements $h_{pq}$ and $g_{pqrs}$ Visscher (2002). A common technique is to assume pure spin-orbitals and integrate over the spin variable. As we want to develop a formalism that is also valid for relativistic calculations, we will remain working with spinors in this work. Adaptation to a spinfree formalism is straightforward, and will not affect computational scaling and error estimates.

The electronic Hamiltonian defined above depends parametrically on the nuclear positions, both explicitly via the nuclear potential and implicitly via the molecular orbitals that change when the nuclei are displaced.

II.2 Asymptotic convergence of energy derivative estimation methods

In this section, we present and compare various methods for calculating energy derivatives on a quantum computer. In Tab. 1, we estimate the computational complexity of all studied methods in terms of the system size $N_{\mathrm{sys}}$ and the estimation error $\epsilon$. We also indicate which methods require quantum phase estimation, as these require longer coherence times than variational methods. Many methods benefit from the amplitude estimation algorithm of Knill et al. (2007), which we have included costings for. We approximate the scaling in $N_{\mathrm{sys}}$ between a best-case scenario (a lattice model with a low-weight energy derivative and aggressive truncation of any approximations), and a worst-case scenario (the electronic structure problem with a high-weight energy derivative and less aggressive truncation). The lower bounds obtained here are competitive with classical codes, suggesting that these methods will be tractable for use in large-scale quantum computing. However, the upper bounds will need reduction in future work to be practical, e.g. by implementing the strategies suggested in McClean et al. (2016); Romero et al. (2018); McClean et al. (2017a).

For wavefunctions in which all parameters are variationally optimized, the Hellmann–Feynman theorem allows for ready calculation of energy gradients as the expectation value of the perturbing operator Kassal and Aspuru-Guzik (2009); Wecker et al. (2015):

[TABLE]

This expectation value may be estimated by repeated measurement of a prepared ground state on a quantum computer (Sec. IV.3.1), and classical calculation of the coefficients of the Hermitian operator $\partial\hat{H}/\partial\lambda$ (Sec. A). If state preparation is performed using a VQE, estimates of the expectation values in Eq. 4 will often have already been obtained during the variational optimization routine. If one is preparing a state via QPE, one does not get these expectation values for free, and must repeatedly measure the quantum register on top of the phase estimation routine. Such measurement is possible even with single-ancilla QPE methods which do not project the system register into the ground state (see Sec. IV.5.1). Regardless of the state preparation method, the estimation error may be calculated by summing the sampling noise of all measured terms (assuming the basis set error and ground state approximation errors are minimal).

The Hellmann–Feynman theorem cannot be so simply extended to higher-order energy derivatives. We now study three possible methods for such calculations. The propagator and phase estimation (PPE) method uses repeated rounds of quantum phase estimation to measure the frequency-domain Green’s function, building on previous work on Green’s function techniques Dallaire-Demers and Wilhelm (2016); Bauer et al. (2016); Roggero and Carlson (2018). We may write an energy derivative via perturbation theory as a sum of products of path amplitudes $A$ and energy coefficients $f_{A}$. For example, a second order energy derivative may be written as

[TABLE]

allowing us to identify two amplitudes

[TABLE]

and two corresponding energy coefficients

[TABLE]

The generic form of a $d$-th order energy derivative may be written as

[TABLE]

where $X_{\mathcal{A}}$ counts the number of excitations in the path. As this is different from the number of responses of the wavefunction, $X_{\mathcal{A}}$ does not follow the $2n+1$ rule; rather, $X_{\mathcal{A}}\leq d$. The amplitudes $\mathcal{A}$ take the form 111Higher order ($d\geq 4$) amplitudes will eventually contain terms corresponding to disconnected excitations, which then are products of multiple terms of the form of Eq. 10. Our procedure may be extended to include these contributions; we have excluded them here for the sake of readibility.

[TABLE]

These may be estimated simultaneously with the corresponding energies $E_{j_{x}}$ by applying rounds of QPE in between excitations by operators $\hat{P}_{x}$ (Sec. IV.5.2). One may then classically calculate the energy coefficients $f_{A}$, and evaluate Eq. 9. Performing such calculation over an exponentially large number of eigenstates $|\Psi_{j_{x}}\rangle$ would be prohibitive. However, the quantum computer naturally bins small amplitudes of nearby energy with a resolution $\Delta$ controllable by the user. We expect the resolution error to be smaller than the error in estimating the amplitudes $\mathcal{A}(j_{1},\ldots,j_{X_{\mathcal{A}}-1})$ (Sec. IV.5.4); we use the latter for the results in Tab. 1.

In lieu of the ability to perform the long circuits required for phase estimation, one may approximate the sum over (exponentially many) eigenstates $|\Psi_{j}\rangle$ in Eq. 9 by taking a truncated set of (polynomially many) approximate eigenstates $|\tilde{\Psi}_{j}\rangle$. We call such an approximation the eigenstate truncation approximation, or ETA for short. However, on a quantum computer, we expect both to better approximate the true ground state $|\Psi_{0}\rangle$, and to have a wider range of approximate excited states Parrish et al. (2018); McClean et al. (2017a); Santagati et al. (2018); Higgott et al. (2018); Endo et al. (2018). In this work, we focus on the quantum subspace expansion (QSE) method of McClean et al. (2017a). This method proceeds by generating a set of $N_{\mathrm{E}}$ vectors $|\chi_{j}\rangle$ connected to the ground state $|\Psi_{0}\rangle$ by excitation operators $\hat{E}_{j}$,

[TABLE]

This is similar to truncating the Hilbert space using a linear excitation operator in the (classical) equation of motion coupled cluster (EOMCC) approach Bartlett (2011). The $|\chi_{j}\rangle$ states are not guaranteed to be orthonormal; the overlap matrix

[TABLE]

is not necessarily the identity. To generate the set $|\tilde{\Psi}_{j}\rangle$ of orthonormal approximate eigenstates, one can calculate the projected Hamiltonian matrix

[TABLE]

and solve the generalized eigenvalue problem:

[TABLE]

Regardless of the method used to generate the eigenstates $|\tilde{\Psi}_{j}\rangle$, the dominant computational cost of the ETA is the need to estimate $N_{\mathrm{E}}^{2}$ matrix elements. Furthermore, to combine all matrix elements with constant error requires the variance of each estimation to scale as $N_{\mathrm{E}}^{-2}$ (assuming the error in each term is independent). This implies that, in the absence of amplitude amplification, the computational complexity scales as $N_{\mathrm{E}}^{4}$. Taking all single-particle excitations sets $N_{\mathrm{E}}\propto N_{\mathrm{sys}}^{2}$. However, in a lattice model one might consider taking only local excitations, setting $N_{\mathrm{E}}\propto N_{\mathrm{sys}}$. Further reductions to $N_{\mathrm{E}}$ will increase the systematic error from Hilbert space truncation (Sec. IV.6), although this may be circumvented somewhat by extrapolation.

For the sake of completeness, we also consider here the cost of numerically estimating an energy derivative by estimating the energy at multiple points;

[TABLE]

Here, the latter formula is preferable if one has direct access to the derivative in a VQE via the Hellmann–Feynman theorem, whilst the former is preferable when one may estimate the energy directly via QPE. In either case, the sampling noise (Sec. IV.3.1 and Sec. IV.4) is amplified by the division of $\delta\lambda$. This error then competes with the $O(\delta\lambda^{2})$ finite difference error, the balancing of which leads to the scaling laws in Tab. 1. This competition can be negated by coherently copying the energies at different $\lambda$ to a quantum register of $L$ ancilla qubits and performing the finite difference calculation there Kassal et al. (2008); Jordan (2005). Efficient circuits (and lower bounds) for the complexity of such an algorithm have not been determined, and proposed methods involve coherent calculation of the Hamiltonian coefficients on a quantum register. This would present a significant overhead on a near-term device, but with additive and better asymptotic scaling than the QPE step itself (which we use for the results in Tab. 1).

II.3 Geometry optimization on a superconducting quantum device

To demonstrate the use of energy derivatives directly calculated from a quantum computing experiment, we first perform geometry optimization of the diatomic H2 molecule, using two qubits of a superconducting transmon device. (Details of the experiment are given in Sec. IV.7.) Geometry optimization aims to find the ground state molecular geometry by minimizing the ground state energy $E_{0}(\mathbf{R})$ as a function of the atomic co-ordinates $R_{i}$. In this small system, rotational and translational symmetries reduce this to a minimization as a function of the bond distance $R_{\rm H-H}$ In Fig. 1, we illustrate this process by sketching the path taken by Newton’s minimization algorithm from a very distant initial bond distance ($R_{\rm H-H}=1.5$Å). At each step of the minimization we show the gradient estimated via the Hellman–Feynman theorem. Newton’s method additionally requires access to the Hessian, which we calculated via the ETA (details given in Sec. IV.7). The optimization routine takes $5$ steps to converge to a minimum bond length of $0.749$Å, within $0.014$Å of the target FCI equilibrium bond length (given the chosen STO-3G basis set). To demonstrate the optimization stability, we performed $100$ simulations of the geometry optimization experiment on the quantumsim density-matrix simulator O’Brien et al. (2017), with realistic sampling noise and coherence time fluctuations (details given in Sec. IV.8). We plot all simulated optimization trajectories on Fig. 1, and highlight the median $(R_{\rm H-H},E(R_{\rm H-H}))$ of the first $7$ steps. Despite the rather dramatic variations between different gradient descent simulations, we observe all converging to within similar error bars, showing that our methods are indeed stable.

To study the advantage in geometry optimization from direct estimation of derivatives on a quantum computer, we compare in Fig. 2 our performance with gradient-free (Nelder-Mead) and Hessian-free (conjugate gradient, or CG) optimization routines. We also compare the performance of Newton’s method with an approximate Hessian from Hartree-Fock (HF) theory. All methods converge to near-identical minima, but both Newton methods converge about twice as fast as the raw CG method, which in turn converges about twice as fast as Nelder-Mead. The density-matrix simulations predict that the ETA method Hessians provide less stable convergence than the HF Hessians; we attribute this to the fact that the HF Hessian at a fixed bond distance does not fluctuate between iterations. The density-matrix simulations also predict the CG method performance to be on average much closer to the Newton’s method performance. However, we expect the separation between gradient and Hessian-free optimization routines to become more stark at larger system sizes, as is observed typically in numerical optimization Nocedal and Wright (2006).

To separate the performance of the energy derivative estimation from the optimization routine, we study the error in the energy $E$, the Jacobian $J$ and Hessian $K$ given as $\epsilon_{A}={|A_{\mathrm{FCI}}-A_{\mathrm{expt}}|}$, $(A=E,J,K).$ In Fig. 3, we plot these errors for different bond distances. For comparison we additionally plot the error in the HF Hessian approximation. We observe that the ETA Hessian is significantly closer than the HF-approximated Hessian to the true value, despite the similar performance in geometry optimization. The accuracy of the ETA improves at large bond distance, where the HF approximation begins to fail, giving hope that the ETA Hessian will remain appropriate in strongly correlated systems where this occurs as well.

II.4 Polarizability estimation

A key property to model in quantum chemistry is the polarizability, which describes the tendency of an atom or molecule to acquire an induced dipole moment due to a change in an external electric field $\vec{F}$. The polarizability tensor may be calculated as $\alpha_{i,j}=\left.\frac{\partial E(\vec{F})}{\partial F_{i}\partial F_{j}}\right|_{\vec{F}=0}$ 222The first-order derivative $\partial E/\partial F_{i}$ gives the dipole moment, which is also of interest, but is zero for the hydrogen molecule.. In Fig. 4, we calculate the $z$-component of the polarizability tensor of H2 in the ETA, and compare it to FCI and HF polarizability calculations on a classical computer. We observe good agreement to the target FCI result at low $R_{\rm H-H}$, finding a $0.060$ a.u. ($2.1\%$) error at the equilibrium bond distance (including the inaccuracy in estimating this distance). However our predictions deviate from the exact result significantly at large bond distance ($R_{\rm H-H}\gtrsim 1.2$ Å). We attribute this deviation to the transformation used to reduce the description of H2 to a two-qubit device (see Sec. IV.7), which is no longer valid when adding the dipole moment operator to the Hamiltonian. To confirm this, we classically compute the FCI polarizability following the same transformation (which corresponds to projecting the larger operator onto a $2$-qubit Hilbert space). We find excellent agreement between this and the result from the quantum device across the entire bond dissociation curve. This implies that simulations of H2 on a $4$-qubit device should match the FCI result within experimental error.

III Conclusion

In this work, we have surveyed possible methods for estimating energy gradients on a quantum computer, including two new techniques of our own design. We have estimated the computational complexity of these methods, both in terms of the accuracy required for the result and the size of the studied system. We have demonstrated the use of these methods on a small-scale quantum computing experiment, obtaining the equilibrium bond length of the H2 molecule to $0.014$Å ($2\%$) of the target Full-CI value, and estimating the polarizability at this bond length to within $0.060$ a.u. ($2.1\%$) of the same target.

Our methods do not particularly target the ground state over any other eigenstate of the system, and so can be used out-of-the-box for gradient estimation for excited state chemistry. They hold further potential for calculating frequency-domain Green’s functions in strongly correlated physics systems (as PPE estimates the gradient through a Green’s function calculation). However, throughout this work we made the assumption that the gap $\delta$ between ground and excited state energies was sufficiently large to not be of concern (namely that $\delta\propto N_{\mathrm{sys}}^{-1}$). Many systems of interest (such as high-temperature superconductors) are characterized by gap closings in the continuum limit. How this affects PPE is an interesting question for future work. Further investigation is also required to improve some of the results drawn upon for this work, in particular reducing the number of measurements required during a VQE and improving amplitude estimation during single-round QPE.

IV Methods

IV.1 Classical computation

The one- and two-electron integrals defining the fermionic Hamiltonian in Eq. 3 are obtained from a preliminary HF calculation that is assumed to be easily feasible on a classical computer. In non-relativistic theory the one-electron integrals are given by

[TABLE]

where $V(\mathbf{r})$ is the electron-nuclear attraction potential from fixed nuclei at positions $\mathbf{R}_{i}$. The two-electron integrals are given by,

[TABLE]

For simplicity we used a finite difference technique to compute the matrix representations of perturbations corresponding to a change in nuclear coordinates and an external electric field

[TABLE]

and

[TABLE]

where $\delta\lambda=0.001$ corresponds to a small change in $\lambda$. The above (perturbed) quantum chemical Hamiltonians have been determined within the Dirac program DIR and transformed into qubit Hamiltonians using the OpenFermion McClean et al. (2017b) package. This uses the newly-developed, freely-available Senjean OpenFermion-Dirac interface, allowing for the simulation of relativistic quantum chemistry calculations on a quantum computer. While a finite difference technique was sufficient for the present purpose, such schemes are sensitive to numerical noise and have a high computational cost when applied to larger molecular systems. A consideration of the analytical calculation of energy derivatives can be found in the Supplementary Materials.

IV.2 Approximate bound calculation details

In this section we detail our method for calculating the approximate bounds in Table. 1. We first estimate the error $\epsilon$ (Tab. 2, first row; details of the non-trivial calculations for the PPE and ETA methods given in Sec. IV.5 and Sec. IV.6) respectively. Separately, we may calculate the time cost by multiplying the number of circuits, the number of repetitions of said circuits ($n_{\mathrm{m}}$, $N_{\mathrm{m}}$, and $K$ depending on the method), and the time cost of each circuit (Tab. 2, second row). (This assumes access to only a single quantum processor, and can in some situations be improved by simultaneous measurement of commuting terms, as discussed in Sec. IV.3.1.) We then choose the scaling of the number of circuit repetitions as a function of the other metaparameters to fix $\epsilon$ constant (Tab. 2, third row). We finally substitute the lower and upper bounds for these metaparameters in terms of the system size as stated throughout the remaining sections. For reference, we summarize these bounds in Tab. 3.

IV.3 Quantum simulation of the electronic structure problem

IV.3.1 Preliminaries

To represent the electronic structure problem on a quantum computer, we need to rewrite the fermionic creation and annihilation operators $\hat{c}^{{\dagger}}_{i}$, $\hat{c}_{i}$ in terms of qubit operators (e.g. elements of the Pauli basis $\mathbb{P}^{N}=\{I,X,Y,Z\}^{\otimes N}$). This is necessarily a non-local mapping, as local fermionic operators anti-commute, while qubit operators commute. A variety of such transformations are known, including the Jordan-Wigner Jordan and Wigner (1928); Ortiz et al. (2001), and Bravyi-Kitaev Bravyi and Kitaev (2002) transformations, and more recent developments Verstraete and Cirac (2005); Seeley et al. (2012); Whitfield et al. (2016); Steudtner and Wehner (2018); Setia et al. (2018); Setia and Whitfield (2018); Tranter et al. (2018); O’Brien et al. (2018).

After a suitable qubit representation has been found, we need to design quantum circuits to implement unitary transformations. Such circuits must be constructed from an appropriate set of primitive units, known as a (universal) gate-set. For example, one might choose the set of all single-qubit operators, and a two-qubit entangling gate such as the controlled-NOT, C-Phase, or iSWAP gates Nielsen and Chuang (2002). One can then build the unitary operators $e^{i\theta\hat{P}}$ for $\hat{P}\in\mathbb{P}^{N}$ exactly (in the absence of experimental noise) with a number of units and time linear in the size of $\hat{P}$ Whitfield et al. (2011). (Here, size refers to the number of non-identity tensor factors of $\hat{P}$.) Optimizing the scheduling and size of these circuits is an open area of research, but many improvements are already known Hastings et al. (2014).

Transformations of a quantum state must be unitary, which is an issue if one wishes to prepare e.g. $\partial\hat{H}/\partial\lambda|\Psi_{0}\rangle$ on a quantum register ($\partial\hat{H}/\partial\lambda$ is almost always not unitary). To circumvent this, one must decompose $\partial\hat{H}/\partial\lambda$ as a sum of $N_{\mathrm{U}}$ unitary operators, perform a separate circuit for each unitary operator, and then combine the resulting measurements as appropriate. Such a decomposition may always be performed using the Pauli group (although better choices may exist). Each such decomposition incurs a multiplicative cost of $N_{\mathrm{U}}$ to the computation time, and further increases the error in any final result by at worst a factor of $N_{\mathrm{U}}^{1/2}$. This makes the computational complexities reported in Tab. 2 highly dependent on $N_{\mathrm{U}}$. The scaling of $N_{\mathrm{U}}$ with the system size is highly dependent on the operator to be decomposed and the choice of decomposition. When approximating this in Tab. 3 we use a range between $O(N_{\mathrm{sys}})$ (which would suffice for a local potential in a lattice model) to $O(N_{\mathrm{sys}}^{4})$ (for a two-body interaction).

To interface with the outside world, a quantum register needs to be appropriately measured. Similarly to unitary transformations, one builds these measurements from primitive operations, typically the measurement of a single qubit in the $Z$ basis. This may be performed by prior unitary rotation, or by decomposing an operator $\hat{O}$ into $N_{\mathrm{T}}$ Hermitian terms $\hat{O}_{i}$ (which may be measured separately). $N_{\mathrm{T}}$ differs from $N_{\mathrm{U}}$ defined above, as the first is for a Hermitian decomposition of a derivative operator and the second is for a unitary decomposition. Without a priori knowledge that the system is near an eigenstate of a operator $\hat{O}$ to be measured, one must perform $n_{\mathrm{m}}$ repeated measurements of each $\hat{O}_{i}$ to estimate $\langle\hat{O}\rangle$ to an accuracy $\propto n_{\mathrm{m}}^{-1/2}N_{\mathrm{T}}^{1/2}$. As such measurement is destructive, this process requires $n_{\mathrm{m}}N_{\mathrm{T}}$ preparations and pre-rotations on top of the measurement time. This makes the computational costs reported in Tab. 2 highly dependent on $N_{\mathrm{T}}$. The scaling of $N_{\mathrm{T}}$ with the system size $N_{\mathrm{sys}}$ is highly dependent on the operator $\hat{O}$ to be measured and the choice of measurements to be made McClean et al. (2016); Romero et al. (2018). In Tab. 3, we assume a range between $O(N_{\mathrm{sys}})$ and $O(N_{\mathrm{sys}}^{4})$ to calculate the approximate computation cost. This is a slight upper bound, as terms can be measured simultaneously if they commute, and error bounds may be tightened by accounting for the covariance between non-commuting terms McClean et al. (2016). The details on how this would improve the asymptotic scaling are still lacking in the literature however, and so we leave this as an obvious target for future work.

Throughout this text we require the ability to measure a phase $e^{i\phi}$ between the $|0\rangle$ and $|1\rangle$ states of a single qubit. (This information is destroyed by a measurement in the $Z$ basis, which may only obtain the amplitude on either state.) Let us generalize this to a mixed state on a single qubit, which has the density matrix Nielsen and Chuang (2002)

[TABLE]

where $p_{0}+p_{1}=1$, and $0\leq p_{+}\leq\sqrt{p_{0}p_{1}}\leq 0.5$. If one repeatedly performs the two circuits in Fig. 5 (which differ by a gate $R=I$ or $R=R_{Z}=e^{i\frac{\pi}{4}Z}$), and estimates the probability of a final measurement $m=0,1$, one may calculate

[TABLE]

We define the circuit element $M_{T}$ throughout this work as the combination of the two circuits to extract a phase using this equation. As above, the error in estimating the real and imaginary parts of $2p_{+}e^{i\phi}$ scales as $n_{\mathrm{m}}^{-1/2}$.

IV.3.2 Hamiltonian Simulation

Optimal decompositions for complicated unitary operators are not in general known. For the electronic structure problem, one often wants to perform time evolution by a Hamiltonian $\hat{H}$, requiring a circuit for the unitary operator $U=e^{i\hat{H}t}$. For a local (fermionic or qubit) Hamiltonian, the length $T_{\mathrm{U}}$ of the required circuit is polynomial in the system size $N_{\mathrm{sys}}$. However, the coefficient of this polynomial is often quite large; this depends on the chosen Hamiltonian, its basis set representation, the filling factor $\eta$ (i.e. number of particles), and whether additional ancilla qubits are used McArdle et al. (2018a); Cao et al. (2018). Moreover, such circuits usually approximate the target unitary $U$ with some $\tilde{U}$ with some bounds on the error $\epsilon_{\mathrm{H}}=\|U-\tilde{U}\|_{\mathrm{S}}$. This bound $\epsilon_{\mathrm{H}}$ is proportional to the evolution time $t$, providing a ‘speed limit’ for such simulation Berry et al. (2007). For the electronic structure problem, current methods achieve scaling between $O(N_{\mathrm{sys}}^{2})$ Low and Wiebe (2018) and $O(N_{\mathrm{sys}}^{6})$ Hastings et al. (2014); Motzoi et al. (2017) for the circuit length $T_{\mathrm{U}}$, assuming $\eta\propto N_{\mathrm{sys}}$ (and fixed $t$, $\epsilon$). (When $\eta$ is sublinear in $N_{\mathrm{sys}}$, better results exist Babbush et al. (2018a).) The proven $O(N_{\mathrm{sys}}^{6})$ scaling is an upper bound, and most likely reduced by recent work Motta et al. (2018); Campbell (2018). For simpler models, such as the Hubbard model, scalings between $O(1)$ and $O(N_{\mathrm{sys}})$ are available O’Brien et al. (2018); Kivlichan et al. (2019). As we require $t\propto N_{\mathrm{sys}}^{-1}$ for the purposes of phase estimation (described in Sec. IV.4), this scaling is reduced by an additional factor throughout this work (though this cannot reduce the scaling below $O(1)$). For Tab. 3, we use a range of $T_{\mathrm{U}}=O(1)$ and $T_{\mathrm{U}}=O(N_{\mathrm{sys}}^{5})$ when approximating the scaling of our methods with the system size.

IV.3.3 Ground state preparation and measurement

A key requirement for our derivative estimation methods is the ability to prepare the ground state $|\Psi_{0}\rangle$ or an approximation to it on the system register. Various methods exist for such preparation, including QPE (see Sec. IV.4), adiabatic state preparation Wu et al. (2002), VQE Peruzzo et al. (2014); McClean et al. (2016), and more recent developments Motta et al. (2019); Kyriienko (2019). Some of these preparation methods (in particular adiabatic and variational methods) are unitary, whilst others (phase estimation) are projective. Given a unitary preparation method, one may determine whether the system remains in the ground state by inverting the unitary and measuring in the computational basis (Sec. IV.3.1). By contrast, such determination for QPE requires another phase estimation round, either via multiple ancilla qubits or by extending the methods in Sec. IV.4. Unitary preparation methods have a slight advantage in estimating expectation values of unitary operators $\hat{U}$; the amplitude amplification algorithm Knill et al. (2007) improves convergence of estimating $\langle\hat{U}\rangle$ from $\epsilon\propto T^{-1/2}$ to $\epsilon\propto T^{-1}$ (in a total computation time $T$). However, this algorithm requires repeated application of the unitary preparation whilst maintaining coherence, which is probably not achievable in the near-term. We list the computation time in Tabs. 1 and 2 both when amplitude amplification is (marked with $**$) and is not available.

Regardless of the method used, state preparation has a time cost that scales with the total system size. For quantum phase estimation, this is the time cost $KT_{\mathrm{U}}$ of applying the required estimation circuits, where $K$ is the total number of applications of $e^{i\hat{H}t}$ Higgins et al. (2009). The scaling of a VQE is dependent on the variational ansatz chosen McClean et al. (2016); Romero et al. (2018). The popular UCCSD ansatz for the electronic structure problem has a $O(N_{\mathrm{sys}}^{5})$ computational cost if implemented naively. However, recent work suggests aggressive truncation of the number of variational terms can reduce this as far as $O(N_{\mathrm{sys}}^{2})$ Romero et al. (2018). We take this as the range of scalings for our approximations in Tab. 1.

IV.3.4 Systematic error from ground state approximations (state error)

A variational quantum eigensolver is not guaranteed to prepare the true ground state $|\Psi_{0}\rangle$, but instead some approximate ground state

[TABLE]

In general we expect $|a_{0}|^{2}$ to be relatively large, although this may not be the case for systems with a small gap $\delta$ to nearby excited states. One may place very loose bounds on the error this induces in the energy:

[TABLE]

where here $\|\hat{H}\|_{\mathrm{S}}$ is the spectral norm of the Hamiltonian (its largest eigenvalue). (Note that while in general $\|\hat{H}\|_{\mathrm{S}}$ is difficult to calculate, reasonable approximations are usually obtainable.) As $|\tilde{\Psi}_{0}\rangle$ is chosen to minimize the approximate energy $\tilde{E}_{0}$, one expects to be much closer to the smaller bound than the larger. For an operator $\hat{D}$ (such as a derivative operator $\partial\hat{H}/\partial\lambda$) other than the Hamiltonian, cross-terms will contribute to an additional error in the expectation value $D_{0}=\langle|\hat{D}|\rangle$:

[TABLE]

One can bound this above in turn using the fact that

[TABLE]

which leads to

[TABLE]

Combining this with the error in the energy gives the bound

[TABLE]

This ties the error in our derivative to the energy in our error, but with a square root factor that unfortunately slows down the convergence when the error is small. (This factor comes about precisely because we do not expect to be in an eigenstate of the derivative operator.) Unlike the above energy error, we cannot expect this bound to be loose without a good reason to believe that the orthogonal component $\sum_{j>0}a_{j}|\Psi_{j}\rangle$ has a similar energy gradient to the ground state. This will often not be the case; the low-energy excited state manifold is usually strongly coupled to the ground state by a physically-relevant excitation, causing the energies to move in opposite directions. Finding methods to circumvent this issue are obvious targets for future research. For example, one could optimize a VQE on a cost function other than the energy. One could also calculate the gradient in a reduced Hilbert space (see Sec. IV.6) using eigenstates of $\hat{H}^{(\rm QSE)}+\epsilon\hat{D}^{(\rm QSE)}$ with small $\epsilon$ to ensure the coupling is respected.

IV.4 Quantum Phase Estimation

Non-trivial measurement of a quantum computer is of similar difficulty to non-trivial unitary evolution. Beyond learning the expectation value of a given Hamiltonian $\hat{H}$, one often wishes to know specific eigenvalues $E_{i}$ (in particular for the electronic structure problem, the ground and low-excited state energies). This may be achieved by QPE Kitaev (1995). QPE entails repeated hamiltonian simulation (as described above), conditional on an ancilla qubit prepared in the $|+\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ state. (The resource cost in making the evolution conditional is constant in the system size.) Such evolution causes phase kickback on the ancilla qubit; if the system register was prepared in the state $\sum_{j}a_{j}|\Psi_{j}\rangle$, the combined (system plus ancilla) state evolves to

[TABLE]

Repeated tomography at various $k$ allows for the eigenphases $E_{j}$ to be inferred, up to a phase $E_{j}t+2\pi\equiv E_{j}t$. This inference can be performed directly with the use of multiple ancilla qubits Kitaev (1995), or indirectly through classical post-processing of a single ancilla tomographed via the $M_{T}$ gate of Fig. 5 Knill et al. (2007); Higgins et al. (2009); Svore et al. (2013); Kimmel et al. (2015); Wiebe and Granade (2016); O’Brien et al. (2019).

The error in phase estimation comes from two sources; the error in Hamiltonian simulation and the error in the phase estimation itself. The error in Hamiltonian simulation may be bounded by $\epsilon_{\mathrm{H}}$ (as calculated in the previous section), which in turn sets the time for a single unitary $T_{U}$. Assuming a sufficiently large gap to nearby eigenvalues, the optimal scaling of the error in estimating $E_{j}$ is $A_{j}^{-1}t^{-1}K^{-1}$ (where $A_{j}=|a_{j}|^{2}$ is the amplitude of the $j$th eigenstate in the prepared state). Note that the phase equivalence $E_{j}t+2\pi=E_{j}t$ sets an upper bound on $t$; in general we require $t\propto\|\hat{H}\|_{S}$, which typically scales with $N_{\mathrm{sys}}$. (This scaling was incorporated into the estimates of $T_{\mathrm{U}}$ in the previous section.) The scaling of the ground state amplitude $A_{0}$ with the system size is relatively unknown, although numeric bounds suggest that it scales approximately as $1-\alpha N_{\mathrm{sys}}$ Reiher et al. (2017) for reasonable $N_{\mathrm{sys}}$, with $\alpha$ a small constant. Approximating this as $N_{\mathrm{sys}}^{-1}$ implies that $K\propto N_{\mathrm{sys}}^{2}$ is required to obtain a constant error in estimating the ground state energy.

The error in estimating an amplitude $A_{j}$ during single-ancilla QPE has not been thoroughly investigated. A naive least-squares fit to dense estimation leads to a scaling of $n_{\mathrm{m}}^{-1/2}k_{\max}^{-1/3}$, where $n_{\mathrm{m}}$ is the number of experiments performed at each point $k=1,\ldots,k_{\max}$. One requires to perform controlled time evolution for up to $k_{\max}\approx\max(t^{-1},A_{j}^{-1})$ coherent steps in order to guarantee separation of $\phi_{j}$ from other phases. To obtain a constant error, one must then perform $n_{\mathrm{m}}\propto k_{\max}^{-\frac{2}{3}}\propto\max(A_{j}^{\frac{2}{3}},t^{\frac{2}{3}})$ measurements at each $k$, implying $K=n_{\mathrm{m}}k_{\max}\propto\max(A_{j}^{-\frac{1}{3}},t^{-\frac{1}{3}})$ applications of $e^{i\hat{H}t}$ are required. For the ground state ($A_{j}=A_{0}$), this gives scaling $K\propto N_{\mathrm{sys}}^{\frac{1}{3}}$. By contrast, multiple-ancilla QPE requires $N_{\mathrm{m}}$ repetitions of $e^{i\hat{H}t}$ with $k_{\max}=\max(A_{0}^{-1},t^{-1})$ to estimate $A_{0}$ with an error of $(A_{0}(1-A_{0})N_{\mathrm{m}}^{-1})^{1/2}$. This implies that $N_{\mathrm{m}}\propto A_{0}$ measurements are sufficient, implying $K\propto N_{\mathrm{m}}A_{0}^{-1}$ may be fixed constant as a function of the system size for constant error in estimation of $A_{0}$. Though this has not yet been demonstrated for single-round QPE, we expect it to be achievable and assume this scaling in this work.

IV.5 The propagator and phase estimation method

In this section, we outline the circuits required and calculate the estimation error for our newly developed PPE method for derivative estimation.

IV.5.1 Estimating expectation values with single-ancilla QPE

Though single-ancilla QPE only weakly measures the system register, and does not project it into an eigenstate $|\Psi_{j}\rangle$ of the chosen Hamiltonian, it can still be used to learn properties of the eigenstates beyond their eigenvalues $E_{j}$. In particular, if one uses the same ancilla qubit to control a further unitary operation $\hat{U}$ on the system register, the combined (system plus ancilla) state evolves from Eq. 29 to

[TABLE]

The phase accumulated on the ancilla qubit may then be calculated to be

[TABLE]

Note that the gauge degree of freedom is not present in Eq. 31; if one re-defines $|\Psi_{j}\rangle\rightarrow e^{i\theta}|\Psi_{j}\rangle$, one must send $a_{j}\rightarrow e^{-i\phi_{j}}a_{j}$, and the phase cancels out. One may obtain $g(k)$ at multiple points $k$ via tomography of the ancilla qubit (Fig. 5). From here, either Prony’s method or Bayesian techniques may be used to extract phases $\omega_{j}\approx E_{j}t$ and corresponding amplitudes $\alpha_{j}\approx\sum_{j^{\prime}}a_{j}a_{j^{\prime}}^{*}\langle\Psi_{j^{\prime}}|\hat{U}|\Psi_{j}\rangle$ O’Brien et al. (2019). The amplitudes $\alpha_{j}$ are often not terribly informative, but this changes if one extends this process over a family of operators $U$. For instance, if one chooses $U=e^{ik^{\prime}\hat{H}t}\hat{V}e^{ik\hat{H}t}$ (with $\hat{V}$ unitary), an application of Prony’s method on $k$ returns amplitudes of the form

[TABLE]

from which a second application of Prony’s method obtains phases $\omega_{j^{\prime}}=E_{j^{\prime}}t$ with corresponding amplitudes

[TABLE]

Each subsequent application of QPE requires taking data with $U^{k}$ fixed from $k=1,\ldots,K$ to resolve $K$ individual frequencies (and corresponding eigenvalues). However, if one is simply interested in expectation values $\langle\Psi_{j}|\hat{V}|\Psi_{j}\rangle$ (i.e. when $j=j^{\prime}$), one may fix $k=k^{\prime}$ and perform a single application of Prony’s method, reducing the number of circuits that need to be applied from $O(K^{2})$ to $O(K)$ (see Fig. 6). The error in the estimator $\alpha_{j,j}$ (Eq. 33) may be bounded above by the error in the estimator $\alpha_{j}$ (Eq. 32). However, to estimate $\langle\Psi_{j}|\hat{V}|\Psi_{j}\rangle$ from Eq. 33, one needs to divide by $A_{j}$. This propagates directly to the error, which then scales as $A_{j}^{-1/2}N_{\mathrm{m}}^{-1/2}$. Thus constant error in estimating $\langle\Psi_{j}|\hat{V}|\Psi_{j}\rangle$ is achieved if $K\propto N_{\mathrm{sys}}$.

IV.5.2 PPE circuits

As presented, the operator $\hat{V}$ in Fig. 6 must be unitary. However if one applies additional phase estimation within $\hat{V}$ itself, one can extend this calculation to non-unitary operators, such as those given in Eq. 9. This is similar in nature to calculating the time-domain Green’s function for small $t$ on a quantum computer (which has been studied before in Refs. Wecker et al. (2015); Dallaire-Demers and Wilhelm (2016); Bauer et al. (2016)), but performing the transformation to frequency space with Prony’s method instead of a Fourier transform to obtain better spectral resolution. It can also be considered a generalization of Ref. Roggero and Carlson (2018) beyond the linear response regime. To calculate a $X$th order amplitude (Eq. 10), one may set

[TABLE]

which is unitary if the $\hat{P}_{x}$ are chosen to be a unitary decomposition of $\partial\hat{H}/\partial\lambda_{x}$. In Fig. 7, we show two circuits for the estimation of a second order derivative with $\hat{P}=\partial\hat{H}/\partial\lambda_{1}$, $\hat{Q}=\partial\hat{H}/\partial\lambda_{2}$ (or some piece thereof). The circuits differ by whether QPE or a VQE is used for state preparation. If QPE is used for state preparation, the total phase accumulated by the ancilla qubit over the circuit is

[TABLE]

Applying Prony’s method at fixed $k_{1}$ will obtain a signal at phase $2tE_{0}$ with amplitude

[TABLE]

A second application of Prony’s method in $k_{1}$ allows us to obtain the required amplitudes

[TABLE]

and the eigenvalues $\omega_{j}\approx E_{j}t$, allowing classical calculation of both the amplitudes and energy coefficients required to evaluate Eq. 9. If a VQE is used for state preparation instead, one must post-select on the system register being returned to $|\vec{0}\rangle$. Following this, the ancilla qubit will be in the state

[TABLE]

with an accumulated phase $g(k)=\alpha_{0,0}(k)$ (where $\alpha_{0,0}$ is as defined above). Here, $|\Psi_{0}^{(\rm VQE)}\rangle$ is the state prepared by the VQE unitary (which may not be the true ground state of the system). Both methods may be extended immediately to estimate higher-order amplitudes by inserting additional excitations and rounds of QPE, resulting in amplitudes of the form $\alpha_{0,0}(k_{1},\ldots,k_{X})$. To explore this in more detail, in App. C we apply this method to a simple toy system.

We note that the VQE post-selection does not constitute ‘throwing away data’; if the probability of post-selecting $|\Psi_{0}\rangle$ is $p$, we have

[TABLE]

and as the variance in any term $\alpha_{i,j}(k_{1},\ldots,k_{X})$ scales as $|\alpha_{i,j}(k_{1},\ldots,k_{X})(1-\alpha_{i,j}(k_{1},\ldots,k_{X}))|$, the absolute error in estimating a derivative scales as $p^{1/2}$ (note the lack of minus sign). (Note here that the relative error scales as $p^{-1/2}$.)

IV.5.3 Energy discretization (resolution error)

The maximum number of frequencies estimatable from a signal $g(0),\ldots,g(k)$ is $(k+1)/2$. (This can be seen by parameter counting; it differs from the bound of $k$ for QPE O’Brien et al. (2019) as the amplitudes are not real.) As the time required to obtain $g(k)$ scales at best linearly with $k$ (Sec. IV.3.2), we cannot expect fine enough resolution of all $2^{N_{\mathrm{sys}}}$ eigenvalues present in a $N_{\mathrm{sys}}$-qubit system. Instead, a small amplitude $\mathcal{A}(j_{1},\ldots,j_{X})$ ($|\mathcal{A}(j_{1},\ldots,j_{X})|\leq\Delta$) will be binned with paths $\mathcal{A}^{\prime}(l_{1},\ldots,l_{X})$ of similar energy ($\delta=\max_{x}|E_{j_{x}}-E_{l_{x}}|\ll\Delta$), and labeled with the same energy $E_{B_{x}}\approx E_{j_{x}}\approx E_{k_{x}}$ O’Brien et al. (2019). Here, $\Delta$ is controlled by the length of the signal $g(k)$, i.e. $\Delta\propto\max(k)^{-1}$. This grouping does not affect the amplitudes; as evolution by $e^{ik\hat{H}t}$ does not mix eigenstates (regardless of energy), terms of the form $|\Psi_{j_{x}}\rangle\langle\Psi_{l_{x}}|$ do not appear. (This additional amplitude error would occur if one attempted to calculate single amplitudes of the form $\langle\Psi_{j}|\hat{P}|\Psi_{k}\rangle$ on a quantum device, e.g. using the method in Sec. IV.4 or that of Ref. Roggero and Carlson (2018), and multiply them classically to obtain a $d>1$-th order derivative.) The PPE method then approximates Eq. 9 as

[TABLE]

Classical post-processing then need only sum over the (polynomially-many in $\Delta$) bins $B_{x}$ instead of the (exponentially-many in $N_{\mathrm{sys}}$) eigenstates $|\Psi_{j_{x}}\rangle$, which is then tractable.

To bound the resolution error in the approximation $f_{\mathcal{A}}(E_{0};E_{j_{1}},\ldots,E_{j_{X_{\mathcal{A}}-1}})\rightarrow f_{\mathcal{A}}(E_{0};E_{B_{1}},\ldots,E_{B_{X_{\mathcal{A}}-1}})$, we consider the error if $E_{j}$ were drawn randomly from bins of width $\|\hat{H}\|_{\mathrm{S}}\Delta$ (where $\|\hat{H}\|_{\mathrm{S}}$ is the spectral norm). The energy functions $f$ take the form of $X_{\mathcal{A}}-1$ products of $\frac{1}{E_{j_{x}}-E_{0}}$. If each term is independent, these may be bounded as

[TABLE]

where $\delta$ is the gap between the ground and excited states. Then, as the excitations $\hat{P}$ are unitary, for each amplitude $\mathcal{A}$ one may bound

[TABLE]

Propagating variances then obtains

[TABLE]

Where $N_{\mathcal{A}}$ is the number of amplitudes in the estimation of $D$. As we must decompose each operator into unitaries to implement in a circuit, $N_{\mathcal{A}}\propto N_{\mathrm{U}}^{d}$.

This bound is quite loose; in general we expect excitations $\partial\hat{H}/\partial\lambda$ to couple to low-level excited states, which lie in a constant energy window (rather than one of width $\|\hat{H}\|_{S}$), and that contributions from different terms should be correlated (implying that $N_{\mathcal{A}}$ should be treated as constant here). This implies that one may take $\Delta$ roughly constant in the system size, which we assume in this work.

IV.5.4 Sampling noise error

We now consider the error in calculating Eq. 38 from a finite number of experiments (which is separate to the resolution error above). Following Sec. IV.4 we have that, if QPE is used for state preparation

[TABLE]

If one were to use a VQE for state preparation, the factors of $A_{0}$ would be replaced by the state error of Sec. IV.3.4. We have not included this calculation in Tab. 2 for the sake of simplicity. Then, assuming each term in Eq. 38 is independently estimated, we obtain

[TABLE]

Substituting the individual scaling laws one obtains

[TABLE]

where again $N_{\mathcal{A}}\propto N_{U}^{d}$. This result is reported in Tab. 2.

IV.6 Eigenstate truncation approximation details

In this section, we outline the classical post-processing required to evaluate Eq. 9 in the ETA, using QSE to generate approximate eigenstates. We then calculate the complexity cost of such estimation, and discuss the systematic error in an arbitrary response approximation from Hilbert space truncation.

The chosen set of approximate excited states $|\tilde{\Psi}_{j}\rangle$ defines a subspace $\mathcal{H}^{(\mathrm{QSE})}$ of the larger FCI Hilbert space $\mathcal{H}^{(\mathrm{FCI})}$. To calculate expectation values within this subspace, we project the operators $\hat{O}$ of interest (such as derivatives like $\partial\hat{H}/\partial\lambda$) onto $\mathcal{H}^{(\mathrm{QSE})}$, giving a set of reduced operators $\hat{O}^{(\mathrm{QSE})}$ ($O^{(\mathrm{QSE})}_{i,j}=\langle\tilde{\Psi}_{i}|\hat{O}|\tilde{\Psi}_{j}\rangle$). These are $N_{\mathrm{E}}\times N_{\mathrm{E}}$-dimensional classical matrices, which may be stored and operated on in polynomial time. Methods to obtain the matrix elements $O^{(\mathrm{QSE})}_{i,j}$ are dependent on the form of the $|\tilde{\Psi}_{j}\rangle$ chosen. Within the QSE, one can obtain these by directly measuring McClean et al. (2017a)

[TABLE]

using the techniques outlined in Sec. IV.3.1, and rotating the basis from $\{|\chi_{j}\rangle\}$ to $\{|\tilde{\Psi}_{j}\rangle\}$ (using Eq. 14).

The computational complexity for a derivative calculation within the QSE is roughly independent of the choice of $|\tilde{\Psi}_{j}\rangle$. The error $\epsilon$ may be bounded above by the error in each term of the $N_{\mathrm{E}}\times N_{\mathrm{E}}$ projected matrices, which scales as either $N_{\mathrm{T}}^{1/2}n_{\mathrm{m}}^{-1/2}$ (when directly estimating), $A_{j}^{-1/2}N_{\mathrm{m}}^{-1/2}$ (when estimating via QPE), or $N_{\mathrm{U}}^{1/2}K^{-1}n_{\mathrm{m}}^{-1/2}$ (using the amplitude estimation algorithm). We assume that the $N_{\mathrm{E}}^{2}$ terms are independently estimated, in which case $\epsilon$ scales with $N_{\mathrm{E}}$. In general this will not be the case, and $\epsilon$ could scale as badly as $N_{\mathrm{E}}^{2}$, but we do not expect this to be typical. Indeed, one can potentially use the covariance between different matrix elements to improve the error bound McClean et al. (2017a). As we do not know the precise improvement this will provide, we leave any potential reduction in the computational complexity stated in Tab. 2 to future work. The calculation requires $n_{\mathrm{m}}$ repetitions of $N_{\mathrm{T}}$ circuits for each pair of $N_{\mathrm{E}}$ excitations, leading to a total number of $n_{\mathrm{m}}N_{\mathrm{T}}N_{\mathrm{E}}^{2}$ preparations (each of which has a time cost $T_{\mathrm{P}}$), as stated in Tab 2. (With the amplitude amplification algorithm, the dominant time cost comes from running $O(N_{\mathrm{T}}N_{\mathrm{E}}^{2})$ circuits of length $KT_{\mathrm{P}}$.)

Regardless of the method of generating eigenstates, the ETA incurs a systematic truncation error from approximating an exponentially large number of true eigenstates $|\Psi_{j}\rangle$ by a polynomial number of approximate eigenstates $|\tilde{\Psi}_{j}\rangle=\sum_{l}\tilde{A}_{j,l}|\Psi_{l}\rangle$. This truncation error can be split into three pieces. Firstly, an excitation $\hat{P}|\Psi_{0}\rangle$ may not lie within the response subspace $\mathcal{H}^{(\mathrm{QSE})}$, in which case the pieces lying outside the space will be truncated away. Secondly, the term $\hat{P}|\tilde{\Psi}_{j}\rangle\langle\tilde{\Psi}_{j}|\hat{Q}$ may contain terms of the form $\hat{P}|\Psi_{j}\rangle\langle\Psi_{l}|\hat{Q}$, which do not appear in the original resolution of the identity. Thirdly, the approximate energies $\tilde{E}_{j}$ may not be close to the true energies $E_{j}$ (especially when $|\tilde{\Psi}_{j}\rangle$ is a sum of true eigenstates $|\Psi_{l}\rangle$ with large energy separation $E_{j}-E_{l}$). If one chooses excitation operators $\hat{E}_{j}$ in the QSE so that $\hat{P}=\sum_{j}p_{j}\hat{E}_{j}$, one completely avoids the first error source. By contrast, if one chooses a truncated set of true eigenstates $|\tilde{\Psi}_{j}\rangle=|\Psi_{j}\rangle$, one avoids the second and third error sources exactly. In App. D we expand on this point, and place some loose bounds on these error sources.

IV.7 Experimental methods

The experimental implementation of the geometry optimization algorithm was performed using two of three transmon qubits in a circuit QED quantum processor. This is the same device used in Ref. Sagastizabal et al. (2019) (raw data is the same as in Fig.1(e) of this paper, but with heavy subsequent processing). The two qubits have individual microwave lines for single-qubit gating and flux-bias lines for frequency control, and dedicated readout resonators with a common feedline. Individual qubits are addressed in readout via frequency multiplexing. The two qubits are connected via a common bus resonator that is used to achieve an exchange gate,

[TABLE]

via a flux pulse on the high-frequency qubit, with an uncontrolled additional single-qubit phase that was cancelled out in post-processing. The exchange angle $\theta$ may be fixed to a $\pi/6000$ resolution by using the pulse duration (with a $1~{}\mathrm{ns}$ duration) as a rough knob and fine-tuning with the pulse amplitude. Repeat preparation and measurement of the state generated by exciting to $|01\rangle$ and exchanging through one of $41$ different choices of $\theta$ resulted in the estimation of $41$ two-qubit density matrices $\rho_{i}$ via linear inversion tomography of $10^{4}$ single-shot measurements per pre-rotation Saira et al. (2014). All circuits were executed in eQASM Fu et al. (2019) code compiled with the QuTech OpenQL compiler, with measurements performed using the qCoDeS Johnson et al. (2016) and PycQED Rol et al. (2016) packages.

To use the experimental data to perform geometry optimization for H2, the ground state was estimated via a VQE Peruzzo et al. (2014); McClean et al. (2016). The Hamiltonian at a given H-H bond distance $R_{\rm H-H}$ was calculated in the STO-3G basis using the Dirac package DIR , and converted to a qubit representation using the Bravyi-Kitaev transformation, and reduced to two qubits via exact block-diagonalization O’Malley et al. (2016) using the Openfermion package McClean et al. (2017b) and the Openfermion-Dirac interface Senjean . With the Hamiltonian $\hat{H}(R_{\rm H-H})$ fixed, the ground state was chosen variationally: $\rho(R_{\rm H-H})=\min_{\rho_{i}}\mathrm{Trace}[\hat{H}(R_{\rm H-H})\rho_{i}]$. The gradient and Hessian were then calculated from $\rho(R_{\rm H-H})$ using the Hellmann–Feynman theorem (Sec. II.2) and ETA (Sec. IV.6) respectively. For the ETA, we generated eigenstates using the QSE, with the Pauli operator $XY$ as a single excitation. This acts within the number conserving subspace of the two-qubit Hilbert space, and, being imaginary, will not have the real-valued H2 ground state as an eigenstate. (This in turn guarantees the generated excited state is linearly independent of the ground state.) For future work, one would want to optimize the choice of $\theta$ at each distance $R_{\rm H-H}$, however this was not performed due to time constraints. We have also not implemented the error mitigation strategies studied in Ref. Sagastizabal et al. (2019) for the sake of simplicity.

IV.8 Simulation methods

Classical simulations of the quantum device were performed in the full-density-matrix simulator *(*quantumsim) O’Brien et al. (2017). A realistic error model of the device was built using experimentally calibrated parameters to account for qubit decay ($T_{1}$), pure dephasing ($T_{2}^{*}$), residual excitations of both qubits, and additional dephasing of qubits fluxed away from the sweet spot (which reduces $T_{2}^{*}$ to $T_{2}^{*,red}$ for the duration of the flux pulse). This error model further accounted for differences in the observed noise model on the individual qubits, as well as fluctuations in coherence times and residual excitation numbers. Further details of the error model may be found in Ref. Sagastizabal et al. (2019) (with device parameters in Tab.S1 of this reference).

With the error model given, $100$ simulated experiments were performed at each of the $41$ experimental angles given. Each experiment used unique coherence time and residual excitation values (drawn from a distribution of the observed experimental fluctuations), and had appropriate levels of sampling noise added. These density matrices were then resampled $100$ times for each simulation.

V Author contributions

New quantum algorithms were designed by TEO. Errors were estimated by TEO, BS, and AD. Experiment was performed by RS and LDC. Density matrix simulations were performed by XBM. Computational chemistry calculations performed by BS and LV. Gradient descent and polarizability calculations from experimental and simulated data were performed by TEO and BS. Relevance to quantum chemistry and suggested applications provided by FB and LV. Paper written by TEO, BS, FB, and LV (with input from all other authors).

VI Competing interests

The authors declare no conflicts of interest.

Acknowledgements.

We would like to thank B.M. Terhal, I. Kassal, V.P. Ostroukh and C.H. Price for useful discussions, M. Singh, M.A. Rol, C.C. Bultink, X. Fu, N. Muthusubramanian, A. Bruno, M. Beekman, N. Haider, F. Luthi, B. Tarasinski and C. Dickel for experimental assistance, and C.W.J. Beenakker and D. Hohl for support during this project. This research was funded by the Netherlands Organization for Scientific Research (NWO/OCW), an ERC Synergy Grant, Shell Global Solutions BV, and IARPA (U.S. Army Research Office grant W911NF-16-1-0071).

Appendix A Analytical derivative of the one- and two-electron integrals

Molecular orbitals are usually obtained as a linear combination of atomic orbitals (LCAO) via HF or another self-consistent field calculation:

[TABLE]

where $\{c_{\mu p}\}$ are the MO coefficients and $\{\chi_{\mu}(\mathbf{r})\}$ the atomic orbitals (AO). The atomic orbitals are usually chosen to be Gaussian functions centered at the nuclear positions $\mathbf{R}_{i}$ and will move with the nuclei if the geometry is modified. An alternative is to express $\{\phi_{p}(\mathbf{r})\}$ in terms of plane-wave type or Gausslet-type basis functions Babbush et al. (2018b); White (2017), which reduce the position dependence or reduce the number of non-zero $g_{pqrs}$ terms, respectively.

With truncated atom-centered basis sets, one needs to consider the so-called ‘wavefunction forces’ (see Ref. Pulay (1997) and references therein) resulting from the basis-set dependence on the nuclear geometry. However, with an appropriately basis-dependent second-quantized Hamiltonian, all changes in the basis can be incorporated into those of the Hamiltonian Surján et al. (1988); Surján (1989); Helgaker et al. (2014)).

Similarly to classical computing, a quantum computer will be limited by the systematic error in the truncated basis set chosen for the problem. Reducing this by directly increasing the number of basis functions is costly (as this defines the system size $N_{\mathrm{sys}}$). Indeed, the errors in correlation energies decay slowly as $\epsilon=O(N_{\mathrm{sys}}^{-1})$ due to the presence of the Coulomb cusp in the wavefunction Helgaker and Klopper (2000). To bypass this slow convergence, the interelectron distance $r_{12}$ can be incorporated explicitly in the wavefunction (such as R12/F12 wavefunctions), thus leading to the so-called explicitly correlated methods Klopper et al. (2006); Kong et al. (2011). In particular for quantum computers, the basis set error must be traded against the number of terms $N_{\mathrm{H}}$ in the Hamiltonian, and the length of the circuits required to prepare ground states or perform phase estimation; optimizing this trade-off is a topic of much debate in the field Kivlichan et al. (2019); Babbush et al. (2018b); Poulin et al. (2017); Berry et al. (2018).

Let us consider the general case in which both the MO coefficients and the Hamiltonian matrix elements depend on the perturbation, and consider first the analytical calculation of these derivatives with respect to a general perturbation $\lambda$. (When $\lambda=\mathbf{R}_{i}$ this yields the Hessian for geometry optimization, and when $\lambda=F_{i}$ this describes an applied electric field for a polarizability calculation.) According to Eq. (49), the one- and two-electron integrals in the MO basis read:

[TABLE]

and

[TABLE]

respectively. Differentiating the above expressions by $\lambda$ leads to

[TABLE]

and

[TABLE]

where

[TABLE]

The matrix $\mathbf{U}^{(1)}(\lambda)$ parametrizes the first-order changes in the MO coefficients and can be obtained by solving the coupled perturbed Hartree–Fock (CPHF) equations Gerratt and Mills (1968); Caves and Karplus (1969); Osamura et al. (1986). When the number of perturbations that need to be evaluated is large, the use of an explicitly $\lambda$-dependent $\mathbf{U}$ matrix in the evaluation of the energy derivative can be avoided with the Z-vector technique of Handy and Schaeffer Handy and Schaefer III (1984). The last terms in Eqs. (A) and (A) depend on the derivatives of the one- and two-electron integrals in the AO basis. Those are simply given by

[TABLE]

and

[TABLE]

respectively. Once all derivatives in Eqs. (A) and (A) are determined, the Hamiltonian

[TABLE]

can be mapped onto a qubit Hamiltonian via whichever mapping was chosen for said Hamiltonian (see Methods section of main text), as these encodings are Hamiltonian- (and thus perturbation-) independent. In some cases (for instance, the application of an electric field) the derivatives of the basis functions $\{\chi_{\mu}(\mathbf{r})\}$ are equal to zero because the basis functions do not depend on the perturbation.

Appendix B Numerical optimization and approximate Hessian calculations

Numerous numerical methods for geometry optimization exist, some gradient-free, some requiring only gradient calculations, and some making use of both gradients and Hessian data Nocedal and Wright (2006); Schlegel (2011). As sampling noise from a quantum computer is typically far larger than the fixed point error on a classical computer, optimization techniques are required to be stable in the presence of this noise. In particular, common implementations of algorithms that numerically estimate gradients tend to construct approximate derivatives by difference approximations, which (as we investigated above) dramatically enhance sampling noise unless care is taken. The Nelder–Mead gradient-free algorithm Nelder and Mead (1965) is a common choice for optimization in quantum algorithms for this reason; as it does not rely on such an approximation, and implementations in scipy Jones et al. (01) prove relatively stable. Gradient- and Hessian-requiring algorithms do not tend to suffer from such instability as gradient-free methods.

In this work, our geometry optimization was reduced to a one-dimensional problem, removing some of the complexity of the task. With more atoms, one need to choose both the direction and the distance to step towards the minima of the energy landscape. Both the CG and Newton’s methods are adjustments to the steepest descent algorithm (which aims to go solely in the direction of the derivative) to account for local curvature. In the absence of any higher order derivatives to assist adjustment, the non-linear CG algorithm weights each direction against traveling in previously-explored directions, and then performs a line-search in this direction (absent additional information that allows an estimation of how far to initially travel). Newton’s method, by comparison, benefits from access to the Hessian, allowing us to choose

[TABLE]

for the direction. One must compensate here for the fact that we wish to minimize, and not maximise, the energy. For our one-dimensional problem this is achieved by taking the absolute value of $\partial^{2}E_{0}/\partial R_{\mathrm{H-H}}^{2}$; for a higher-dimensional problem this is slightly more involved Nocedal and Wright (2006). Regardless, such modified Newton’s methods tend to provide a far more optimal method for estimating higher dimensional functions than Hessian-free methods Nocedal and Wright (2006); Schlegel (2011). We are further able to bound the minimum bond length in our geometry optimization (in particular to $R_{\mathrm{H-H}}>0.3$ Å), which can be of importance for stability as classical methods tend to fail when atoms are unrealistically close together.

For large systems when low accuracy is needed (e.g. at the start of a geometry optimization calculation), one may consider calculating the Hessian via the HF Hamiltonian for the same geometry as a low-cost alternative to explicit calculation on the quantum computer. This is a standard technique for geometry optimization in computational chemistry Fischer and Almlöf (1992). As the Hessian is not used to determine convergence (which depends instead on the size of the gradient), the approximation only affects the convergence rate and stability, rather than the final result. This is even more so for quasi-Newton methods, as the Hessian is updated during the geometry optimization by the estimated gradients, which are more accurate. Calculating the HF Hessian is a standard procedure in most computational programs; for further mathematical information, we refer the reader to Ref. Yamaguchi et al. (1994). This information is additionally obtained ‘for free’ if the $\mathbf{U}^{(1)}$ matrix, required for first-order derivative calculations of the electronic integrals [Eqs. (A) and (A)], has been obtained (and vice-versa).

Appendix C Application of PPE method to toy system

In this section, we detail an example of PPE for a toy Hamiltonian,

[TABLE]

One may immediately calculate eigenvalues and derivatives of this operator to second order

[TABLE]

where the energy of the ground (excited) state $|E_{-}\rangle$ ($|E_{+}\rangle$) is $E_{-}$ ($E_{+}$). We wish to recover these second-order excitations via repeated QPE (using also QPE for state preparation). For simplicity, let us consider calculating the derivatives when $\lambda_{Z}=0$, which makes $|0\rangle$ an equal superposition of the two eigenstates $|E_{\pm}\rangle=\frac{1}{\sqrt{2}}(|0\rangle\pm|1\rangle)$. The required circuit for second-order gradient calculation, using $|0\rangle$ as a starting state, is given in Fig. 8. Here, $\hat{P}$ and $\hat{Q}$ should be set to $X$ or $Z$ depending on which derivative is to be estimated. Let us first consider estimating $\frac{\partial^{2}E_{0}}{\partial\lambda_{Z}^{2}}$, which requires that we set $\hat{P}=\hat{Q}=Z$. One can calculate the system state prior to measurement of the phase of the ancilla qubit to be

[TABLE]

and so the reduced density matrix on the ancilla qubit may be found to be

[TABLE]

The phase measurement $M_{T}$ then obtains a function

[TABLE]

A Fourier transform 333In practice, one will obtain faster convergence here with Prony’s method than a Fourier transform, but the latter is probably more familiar to the reader. in $k_{0}$ at fixed $k_{1}$ will return peaks at $2t\omega_{\pm}=\pm 2t\lambda_{X}$, with complex amplitude (following the notation of Eq. 32 of the main text)

[TABLE]

We further Fourier transform $\alpha_{-,-}(k_{1})$ (in terms of $k_{1}$), as we are interested in the derivative of the ground state. This obtains a peak at $t\omega_{+}=t\lambda_{X}$, with amplitude

[TABLE]

To finish the computation of the gradient, we note that $|\alpha_{-,-}(k_{1}=0)|=a_{-}^{*}a_{-}=\frac{1}{2}$, allowing us to compute $\langle E_{-}|Z|E_{+}\rangle\langle E_{+}|Z|E_{-}\rangle=1$. We then substitute into Eq. 5 of the main text.

[TABLE]

as required. We observe that the above procedure just as easily obtains the gradient of the excited state $|E_{+}\rangle$, as our starting state was a superposition of both ground and excited states.

Let us now repeat the calculation for the other second order derivatives. Running the circuit of Fig. 8 with $\hat{P}=X$, $\hat{Q}=Z$, the system state evolves to

[TABLE]

and the total phase accumulated on the ancilla qubit may be calculated to be

[TABLE]

A Fourier transform of this function in $k_{0}$ obtains spurious peaks at

[TABLE]

but none at $\omega_{-}$ or $\omega_{+}$, implying all amplitudes are zero, and the total derivative $\left.\frac{\partial^{2}E_{-}}{\partial\lambda_{X}\partial\lambda_{Z}}\right|_{\lambda_{Z}=0}$ is [math] as well.

Finally, running the circuit of Fig. 8 with $\hat{P}=\hat{Q}=X$ gives a pre-measurement state of

[TABLE]

and an accumulated phase of

[TABLE]

This time, a Fourier transform at fixed $k_{1}$ obtains a peak at $\omega_{-}=-t\lambda_{X}$ as required, but with an amplitude

[TABLE]

The second Fourier transform in $k_{1}$ will then return a peak at $\omega_{-}$ again (and no peak at $\omega_{+}$). One must then note that the sum over eigenstates in Eq. 5 of the main text is over excited states only, implying that this peak should be excluded, and the final derivative is

[TABLE]

as required.

The considerations in calculating the last two above derivatives demonstrate facets of the PPE resolution error. In practice, one must take sufficient data to resolve any spurious peaks from the true ground state (so as to not confuse these contributions). In a larger system, spurious peaks may appear halfway between any pair of eigenstates, corresponding to scattering processes that do not return the system to the initial state. Interestingly, these peaks will be absent if a VQE is used for state preparation, which suggests the possibility that they might not be fundamental to the estimation technique. One must also take sufficient data to determine whether peaks from the second Fourier transform correspond to excitations from the ground state to the ground state, so that they may be excluded. Both such methods then require a resolution on the order of the gap between ground and excited states.

Appendix D Bounding the systematic error in the ETA

In this section we attempt to expand on the sources of systematic error in the ETA. This approximation revolves around the approximation of the exponentially-many true eigenstates $|\Psi_{j}\rangle$ with a set of polynomially-many approximate eigenstates $|\tilde{\Psi}_{j}\rangle$. Let us write this transformation in the eigenbasis of the original Hamiltonian $\hat{H}$, as the product of a projection $\Pi^{(\mathrm{ETA})}$ into the Hilbert space $\mathcal{H}^{(\mathrm{ETA})}$, followed by a rotation $V^{(\mathrm{ETA})}$ into the eigenbasis of the projected Hamiltonian

[TABLE]

where $\Sigma^{(\mathrm{ETA})}$ contains the eigenvalues of the $|\tilde{\Psi}_{j}\rangle$. One has a freedom in choosing the unitary $V^{(\mathrm{ETA})}$ (as one can always re-label indices on the $|\Psi_{j}\rangle$); for our purposes, we want to choose $V^{(\mathrm{ETA})}$ to maximise the overlap $|\langle\Psi_{j}|\tilde{\Psi}_{j}\rangle|^{2}=|V_{j,j}|^{2}$. Let us now for simplicity focus on the second-order derivative $\partial^{2}E/\partial\lambda_{1}\partial\lambda_{2}$, and fix $\hat{P}=\partial\hat{H}/\partial\lambda_{1}$, $\hat{Q}=\partial\hat{H}/\partial\lambda_{2}$. One may rewrite the term $\sum_{j}\mathcal{A}_{1}(j)f(E_{j})$ (Eq. 6 of the main text) as

[TABLE]

where $F(j,l)=\delta_{j,l}f(E_{j})$, $p_{j}=\langle\Psi_{j}|\hat{P}|\Psi_{0}\rangle$, and $q_{l}=\langle\Psi_{j}|\hat{Q}|\Psi_{0}\rangle$. In the response method, this is replaced by

[TABLE]

The error in the approximation can then be split into:

The error in the projection $\Pi$: in general $\|\Pi\vec{p}\|\leq\|\vec{p}\|$, with equality when $\Pi\vec{p}=\vec{p}$ (in which case this error source does not exist). 2. 2.

The error in the unitary rotation $V$: ideally $V=I$, or at least $V\Pi\vec{p}=\Pi\vec{p}$ (in which case this error source does not exist). 3. 3.

The error in the energy function $\tilde{F}$: ideally $\tilde{F}=F$, or at least $\tilde{F}V\Pi\vec{p}=FV\Pi\vec{p}$ (in which case this error source does not exist).

Generic bounds on these three error sources or measurements thereof are difficult to come by. The first may be measured directly when $\hat{P}$ or $\hat{Q}$ is a unitary operator, as then $\|\vec{p}\|=1$ and $\|\Pi\vec{p}\|$ may be measured. Moreover, for the QSE approximation, if $\hat{P}=\sum_{i}p_{i}\hat{E}_{i}$ (where $\hat{E}_{i}$ are the excitation operators), then $\hat{P}|\Psi_{0}\rangle$ lies within $\mathcal{H}^{(\mathrm{ETA})}$ and this error source is cancelled. The second may be estimated partially by determining the uncertainty in the energies of the approximate eigenstates

[TABLE]

However, this underestimates the uncertainty from eigenstates that are of similar energy (but may be not connected by $\hat{P}\hat{Q}$). One may further estimate this by measuring $\langle\Psi_{0}|\hat{P}\hat{Q}|\Psi_{0}\rangle$ and comparing to the approximation $\langle\Psi_{0}|\hat{P}^{(\mathrm{ETA})}\hat{Q}^{(\mathrm{ETA})}|\Psi_{0}\rangle$ (although this in turn may miss favourable cancellations of different amplitudes). However, if one estimates $|\tilde{\Psi}_{j}\rangle$ via direct approximation methods, this source of error should be similarly negligible. In this case the final source of error is also negligible (if one can both approximate $|\tilde{\Psi}_{j}\rangle=|\Psi_{j}\rangle$ and estimate the energy $E_{j}$). For the QSE, the final source of error is also bounded above by $\sigma^{2}_{|\tilde{\Psi}_{j}\rangle,\hat{H}}$, which is in turn bounded by the maximum energy of the excitation operators $\|[\hat{H},E_{j}]\|$. Making these bounds more precise and tighter for various excited state methods is a clear direction for future research.

Appendix E Amplitude estimation error for single-round QPE

In this section, we extend the analysis of O’Brien et al. (2019)(App.B) to bound the error in estimating an amplitude $A_{j}=|a_{j}|^{2}$ via single-ancilla QPE. As mentioned in the main text, this is a loose bound, and we expect single-ancilla QPE to achieve similar scaling to QPE performed with multiple ancillas. To calculate the bound, let us consider the situation of attempting to estimate a single phase $\phi$ and a corresponding (complex-valued) amplitude $a$ from terms of the function

[TABLE]

for $k=0,\ldots,k_{\max}$. (Note that $k_{\max}$ corresponds to $K$ in O’Brien et al. (2019)) We assume the real and imaginary part of $g(k)$ are estimated independently as $g_{k}=g_{k}^{0}+ig_{k}^{1}$ from $n_{\mathrm{m}}$ single-shot measurements. Following the procedure of O’Brien et al. (2019), one first estimates the phase $\phi$ via least squares as

[TABLE]

From here, one may construct the vector $f_{k}=e^{ik\tilde{\phi}}=f_{k}^{0}+if_{k}^{1}$ and estimate the amplitude $A$ via least squares as

[TABLE]

(We may evaluate the denominator immediately as, regardless of any error in $\phi$, $f_{k}$ will satisfy $|f_{k}|^{2}=1$.) We do not necessarily wish to include all $k\in\{1,\ldots,k_{\max}\}$ in this estimation, as the error in $f_{k}$ at large $k$ is smaller than the error at small $k$. To make things worse, errors in $f_{k}^{a}$ are correlated. It was derived in O’Brien et al. (2019) that this correlation depends only on $g_{k_{\max}}$; adjusting for the fact that $|A|\leq 1$, we obtain

[TABLE]

Propagating variances then obtains

[TABLE]

and a similar scaling for $\mathrm{Var}[\tilde{A}_{1}]$ We then choose $L$ to minimize the variance

[TABLE]

resulting in a final variance scaling as

[TABLE]