Quantum chemistry calculations on a trapped-ion quantum simulator

Quantum-classical hybrid algorithms are emerging as promising candidates for near-term practical applications of quantum information processors in a wide variety of fields ranging from chemistry to physics and materials science. We report on the experimental implementation of such an algorithm to solve a quantum chemistry problem, using a digital quantum simulator based on trapped ions. Specifically, we implement the variational quantum eigensolver algorithm to calculate the molecular ground state energies of two simple molecules and experimentally demonstrate and compare different encoding methods using up to four qubits. Furthermore, we discuss the impact of measurement noise as well as mitigation strategies and indicate the potential for adaptive implementations focused on reaching chemical accuracy, which may serve as a cross-platform benchmark for multi-qubit quantum simulators.


I. INTRODUCTION
Quantum simulators [1][2][3][4][5] are widely recognized for their promise to help solve challenging problems in a range of fields, including physics, chemistry and materials science. First suggested by Richard Feynman in the early 1980s [6], quantum simulators aim to harness controlled quantum evolutions in order to simulate other quantum systems. Later refined to a formalized notion of a universal quantum computer [7], a quantum simulator's main purpose is to turn the exponential scaling of resources needed to simulate quantum systems on classical computers into a more favorable polynomial overhead on a quantum machine. Experimentally pioneered at the turn of the millennium by nuclear magnetic resonance [8], neutral atom [9] and trapped-ion [10] systems, quantum simulation has attracted much attention over the last decade and might well hold the key to unlock significant advances in our understanding of many-body physics [11].
In the analog approach to quantum simulation, the interactions of a quantum system are engineered to closely match the Hamiltonian of the system of interest. Starting from a well-defined initial state, the simulator evolves under a tailored Hamiltonian to a final state, and measurements are taken that reveal the sought-after information about the system of interest. Prominent examples * These authors contributed equally to this work. † cornelius.hempel@gmail.com are experiments realizing Bose-or Fermi-Hubbard models with ultracold atoms [12], or long-range spin models using either up to 20 fully controlled [13][14][15][16][17][18] or over 100 entangled ions [19][20][21].
The digital approach to quantum simulation employs computational gates to approximate the time evolution of arbitrary local Hamiltonians. This can be done in a freely programmable way [7] that is amenable to quantum error correction [22,23]. With a versatility similar to classical computers, digital quantum simulators offer the prospect to shed light on properties of systems that are fundamentally different from the simulator. Lasercooled ensembles of trapped ions that simulate highenergy physics processes are an example [24]. The first multi-qubit demonstrations of the technique on scalable platforms were implemented on trapped ions [25] and a superconducting system [26].
To enable quantum simulations with a level of complexity that is beyond the reach of classical computers, multiple avenues are being pursued to scale up universal quantum processors realized on platforms such as trapped ions or superconducting qubits [27,28]. As experimental progress paves the way towards the realization of larger scale devices, it is important to develop algorithms that make efficient use of these continually improving quantum resources. This is where quantumclassical hybrid algorithms have recently emerged, which leverage the unique capabilities of quantum devices by incorporating them in classical numerical calculations.
In this article, we use a digital quantum simulator based on trapped ions to experimentally investigate one such algorithm, the variational quantum eigen-solver [29][30][31] for the calculation of molecular ground state energies. We begin with a review of the underlying methods of quantum chemistry as well as the overall algorithmic approach and briefly summarize previous experimental implementations. We then introduce our experimental realization and present results for two example molecules, molecular hydrogen and lithium hydride. Finally, we discuss the results and sources of error with respect to a threshold known as chemical accuracy and provide suggestions for further improvements of the algorithm.

II. SIMULATING QUANTUM CHEMISTRY
A. Approach and specific steps While quantum computers have been shown to offer exponential speedups for a variety of problems, a particularly compelling application is the quantum computation of molecular energies [32][33][34]. Efficient quantum simulations of classically intractable instances of the associated electronic structure problem promise breakthroughs in our understanding of basic chemistry and could revolutionize research into new materials, pharmaceuticals, and industrial catalysts. In the last few years, many efforts have sought to develop new algorithms [35][36][37][38][39][40][41] and better implementation strategies [42][43][44][45][46][47][48][49] for these simulations. Among the explored approaches, the variational quantum eigensolver (VQE) algorithm [29][30][31] has been shown experimentally to be resilient to some systematic errors [50] from imperfect control pulses, making it a promising candidate for quantum simulations in the near-term [51].
The central problem of theoretical chemistry is to compute the lowest energy eigenvalue of the molecular electronic structure Hamiltonian. The eigenstates of this Hamiltonian determine almost all of the properties of interest in a molecule or material, and as the energy gap between the ground and first excited state is often much larger than 25.7 meV (k B T at room temperature), the ground state is of particular interest. To arrive at the standard form of this Hamiltonian used in quantum computation, one begins from a collection of nuclear charges Z i and a number of electrons in the system for which the corresponding Hamiltonian is written as in atomic units ( = 1). Here the positions, masses, and charges of the nuclei are R i , M i , Z i , and the positions of the electrons are r i . This form of the Hamiltonian and its real-space discretization is often referred to as the first-quantized formulation of quantum chemistry and generally enforces the fermionic nature of the electron through an anti-symmetric wavefunction. Several approaches have been developed for treating this form of the problem on a quantum computer [35,[66][67][68][69][70]; however, the focus of this work is on the second-quantized formulation.
To reach the second quantized formulation, one typically first approximates the nuclei as fixed classical point charges under the Born-Oppenheimer approximation and chooses a basis φ i in which to represent the electronic wavefunction. Often, one chooses a basis of N molecular orbitals, constructed as a linear combination of atomic orbitals (LCAO), which are computed using a meanfield procedure known in chemistry as the Hartree-Fock method [71]. The atomic orbital basis functions are derived from variations of hydrogen-like atomic orbitals for different values of Z. They are numerically optimized to match desired physical properties across a range of systems and to be compatible with systematic improvement [72]. Usually the basis functions are expressed as sums of Gaussian functions rather than the original Slater type orbitals to enhance the efficiency of integral evaluation. This choice is convenient for small problems, but not mandatory and much work has been done recently to improve the Hamiltonian representation for electronic structure problems [40,41].
The second-quantized formulation of quantum chemistry leads to the (electronic) Hamiltonian in which the wavefunction's anti-symmetry is enforced through the anti-commutation relations of the fermion creation and annihilation operators a † i and a j . The goal now is to compute the energy of electrons interacting in the fixed external potential of the atomic nuclei. Encoding the electron's spatial and spin coordinate as σ i = (r i , s i ), the two sets of scalar coefficients in Eq. (2) can be calculated via In order to implement the Hamiltonian on a qubitbased quantum simulator, the specific fermionic Hamiltonians need to be transformed to spin Hamiltonians. The most common schemes for this transformation are the Jordan-Wigner transformation [73,74] and the Bravyi-Kitaev transformation [68,[74][75][76][77][78][79]. Characteristically, the Jordan-Wigner (JW) transformation leads to N -local Hamiltonians and the Bravyi-Kitaev (BK) transformation leads to log(N )-local Hamiltonians. In this work, we explore both approaches to arrive at a spin Hamiltonian H, which can be implemented with qubits on a quantum simulator.
VQE algorithms for quantum chemistry typically seek to prepare the ground state of the target system for a particular geometric configuration specified by R i . The energy landscape created by this geometry for the electrons, derives from the integrals in Eqs. (3) above and will be captured by scalar values c throughout the remainder of the article. The calculation then proceeds along the following four steps, visually summarized in Fig. 1.

A digital quantum simulator is initialized in a sim-
ple state |ϕ(0) which represents a good classical approximation to the ground state of Hamiltonian H. When using a molecular orbital basis, this initialization is particularly straightforward as the classical mean-field state (most often a Hartree-Fock solution [71]) is a product state.

2.
A quantum circuit implementing the unitary operation U ( θ 0 ) is applied to |ϕ(0) mapping the initial state to a parameterized "ansatz" state |ϕ( θ 0 ) = U ( θ 0 ) |ϕ(0) . The operation U ( θ 0 ) and its parameter vector θ 0 are chosen based on known structure in the target system, which is usually obtained from classical approximation methods.
3. One measures the expectation value of the energy H = ϕ( θ 0 )| H |ϕ( θ 0 ) of the prepared ansatz state. This makes use of the form of the Hamiltonian H = c H , where H are tensor products of Pauli matrices and c the above mentioned scalars that were pre-calculated for a given internuclear configuration R i . Repeated rounds of state preparation and measurement of the individual terms H allow one to provide an estimate for the expectation value H = c H . In this context, this step is often referred to as Hamiltonian averaging [30,31]  Steps of a VQE-based quantum chemistry calculation. A classical preprocessing stage (A) translates the chemical problem to a classical description of a quantum circuit for unitary operator U ( θ) that can be directly implemented on a quantum processor. (B) As outlined in the main text above, the four steps of the iterative variational optimization are performed in a loop combining quantum resources with a numerical search running on a classical computer. Finally, in (C), the full potential energy surface (PES) of the molecule in question, including both nuclear and electronic contributions, is reconstructed from the data.
It is worth noting that, when applied to classical (diagonal) Hamiltonians, the above strategy is the basis of the quantum approximate optimization algorithm [80].
In order to obtain the total molecular energy of the potential energy surface, the mutual Coulomb energy of the nuclei, separated out in the Born-Oppenheimer approximation of the classical preprocessing step before, has to be added to the VQE result at every configuration R i .
A good ansatz unitary U ( θ) is typically chosen to balance the strengths of available experimental resources with insights about the structure of the fermionic system under study. For instance, it has been argued that a parameterized Trotter approximation [7,81], to the adiabatic state preparation algorithm may serve as an effective ansatz [40,41,51]. In this article, we parameterize our variational ansatz based on insights from a cornerstone method of modern electronic structure theory known as coupled cluster [82]. While widely regarded the "gold standard" of accuracy for the classical study of strongly correlated systems, a drawback of the coupled cluster ansatz is that the method is non-unitary and there is no way to ensure that solutions are prop-erly normalized. This can lead to significant errors when the ground state has what is referred to as multireference character, i.e. when it has significant support on more than one classical configuration within the chosen basis. To address this problem, a unitary variant of the coupled cluster method was proposed in [83]. However, while theoretically immune to some of the problems of the traditional method, unitary coupled cluster (UCC) calculations cannot be performed efficiently on classical computers. Thus, the classical study of UCC involves additional approximations which are also problematic in some cases [84]. However, as first discussed in [30], one can perform finite order UCC calculations efficiently on a quantum computer without any approximation [31,57].
UCC essentially suggests a series of fermionic operators that one should evolve under in order to prepare the ground state. For instance, one can think of the unitary coupled cluster ansatz in a Trotter approximation as where T γ − T † γ are anti-Hermitian fermionic excitation operators. The approximation associated with using only a single Trotter step for the exponential is generally acceptable since the coupled cluster amplitudes θ γ are usually quite small. As discussed in [57], one can scalably deploy classical coupled cluster calculations in order to get a reasonable initial guess of the amplitudes θ 0 . Furthermore, due to selection rules involving the symmetry point groups of molecular orbitals and suitable quantum numbers of the Hamiltonian, most of the coupled cluster amplitudes associated with operators that one might naively anticipate will need to be simulated, actually are zero. A vast majority of the other amplitudes are typically extremely small [57], such that one may choose a threshold value for the θ 0 provided by the classical coupled cluster method to further reduce the number of relevant e θγ (Tγ −T † γ ) terms.

B. Previous demonstrations
Experimental demonstrations of quantum chemistry algorithms have been performed on architectures ranging from quantum photonics [30,85,86], NV centers [61], ion traps [62,87] and NMR platforms [88,89] to superconducting qubits [50,63,64]. The first realizations date back to 2010, when a photonic experiment [85] and a closely related NMR experiment [88] simulated the hydrogen molecule in the simplest basis, using a method that required considerable classical precomputation. A similar experiment computing the dissociation curve of the HeH + cation was performed using a nitrogen vacancy based system in 2015 [61] and the first VQE experiment investigated the same molecule in a minimal basis in 2013 [30] using a photonic implementation.
The first scalable quantum chemistry simulation was carried out on a superconducting platform in 2016 [50].
Here scalable means that the required classical precomputation resources increase in the same way that they would for an arbitrarily large problem. The work in [50] used a one-parameter variational ansatz and performed VQE in post-processing to calculate the ground state energy of H 2 , reaching chemical accuracy for the binding energy (deviations ≤ 1.6 × 10 −3 Hartree). The authors also implemented the phase estimation algorithm for quantum chemistry [32], which did not reach chemical accuracy due to the degradation of gate fidelity in the absence of error correction on their device. As both experiments were performed on the same chip using gates of the same fidelity, these results provide a first validation that the VQE approaches we focus on can deliver superior results on current hardware.
The most advanced implementation of the VQE algorithm in terms of molecular size and resource scaling was published 2017 in the work reported in [63], which simulated three molecules H 2 , LiH and BeH 2 on a superconducting qubit platform. The authors also performed a theoretical investigation of the resource requirements in terms of circuit depth for each molecule and found a significantly more favorable scaling for devices with allto-all qubit connectivity, validating the potential utility of ion trap hardware with its native all-to-all connectivity for these types of simulations. Most recently a superconducting implementation of a modified VQE algorithm was used to also determine excited states of the H 2 molecule [64].
Ion-trap implementations of quantum simulation applied to quantum chemistry have so far been limited to working with a single qubit. The work in reported in [62] used a single 171 Yb + ion and mapped the Hamiltonian of the molecule to four of its internal energy levels, which is not a scalable procedure and provides no advantage over classical computation [90]. However, the other steps of the VQE algorithm were performed scalably, i.e. as they would be for a different state mapping. Recently, the same system was used to explore vibronic molecular spectra [87] employing the bosonic motional degree of freedom of a single trapped ion.
VQE is not the only method that reduces the experimental overhead of phase estimation. Recently, Bayesian approaches to phase estimation have been proposed [91]. This method has been implemented in a photonic simulation of minimal basis H 2 system -reaching chemical accuracy [86], however with scalability properties comparable to [85].

III. EXPERIMENTAL IMPLEMENTATION
The trapped ion system used to implement our digital quantum simulation consists of a linear Paul trap in which a variable number of 40 Ca + ions are electrically confined for many days at a time. Each ion stores a quantum bit that is encoded in a pair of Zeeman states chosen from their 4S 1/2 electronic ground and 3D 5/2 metastable states. In the following, these electronic states will be labeled as qubit states |1 and |0 , respectively.
The qubits are manipulated via a set of global and tightly-focussed, addressed laser beams. Single qubit gate operations are implemented via a three-pulsesequence that combines two global qubit rotations with an intermediate addressed laser pulse manipulating only the targeted qubit [92]. Multi-qubit entangling operations are realized through laser-driven interactions that are mediated by the collective motional modes of the ions within their common trapping potential.
After mapping the UCC excitation operators in Eq. (4) to tensor products of Pauli operators, we implement each operator using appropriate quantum circuits, specifically illustrated for each case below. The circuits' general structure combines two multiqubit entangling gates, each realizing the unitary , with a single qubit operation that encodes a single element of the parameter vector θ in Fig. 1 and a particular θ γ in Eq. (4), respectively. The subscript MS stands for Mølmer and Sørensen, the originators [93,94] of this type of entangling gate and φ = π/2 for the creation of an entangled state, when starting in any state of the computational basis (z). As a unit, this circuit building block [95,96] can effectively realize arbitrary many-body interactions as required by the Pauli operators H resulting from the transformed UCC operators.

The estimation of the expectation values of individual
Pauli operators H is carried out through a projective measurement in the logical z-basis. This measurement is implemented using laser-induced, state-dependent fluorescence which we detect on a CCD camera, allowing for simultaneous readout of all qubits (cf. Appendix A). To determine H with terms in bases other than z, we employ suitable global and single qubit operations to rotate these into the measurement basis prior to the projective measurement. We now discuss the individual steps and results for the two example molecules investigated experimentally.

IV. MOLECULAR HYDROGEN
The simplest possible neutral molecule is formed by two hydrogen atoms. We aim to calculate the potential energy surface of its ground state as function of the internuclear distance R. This will be done in two ways, (1) by scanning the entire parameter space associated with the UCC operator of the spin-transformed molecular Hamiltonian H and (2) by running the VQE algorithm for different values of R separately, leading to only a sparse sampling of this parameter space. Due to its simplicity, H 2 has become a standard example for experimental implementations of quantum chemistry algorithms.
A. Encoding the problem A common molecular basis set ( Figure 1, step A), is that of Slater-type orbitals (STO), represented by linear combinations of, e.g., three Gaussian functions (STO-3G). In this minimal basis set, each hydrogen atom contributes a 1s atomic orbital. Molecular orbitals are then formed by adding the corresponding two atomic wavefunctions. In-phase addition results in a lower energy σ 1s bonding orbital with increased electron density between the nuclei and out-of-phase addition in a σ * 1s anti-bonding orbital associated with a depletion in electron density between them (Figure 2.a).
Using this minimal basis set, a classical Hartree-Fock calculation allows us to pre-compute the relevant molecular orbitals for each geometric configuration specified by the respective internuclear separation R. This numerical step determines the scalars c (R) by solving the integrals in Eq. (3) capturing the spatial and spin coordinates. It also yields a product state solution for the molecular wavefunction, which is used as an initial guess |ϕ(0) in Figure 1's step B. Following a series of theoretical considerations detailed in Appendix B, we determine the only relevant unitary coupled cluster operator for single and double excitations (UCCSD), which can be expressed as, Here, θ is the parameter we aim to optimize, the indices correspond to the minimal set of orbitals and a (a † ) to fermionic annihilation (creation) operators. Using the Jordan-Wigner transformation the fermionic operators of Eq. (5) can be mapped to Pauli operators acting on four qubits as with the initial Hartree-Fock state |ϕ(0) JW = |0011 . Alternatively, after employing the Bravyi-Kitaev transformation, one can make use of the fact that after the mapping only qubits 0 and 2 are affected by operators other than σ z and I, allowing the number of required qubits to be reduced to two [50]. The resulting unitary acts on the Hartree-Fock ansatz state |ϕ(0) BK = |01 . Each transformed UCC operator is implemented using the corresponding circuit shown in Figure 2.b. They only differ in the number of required qubits and the projective measurements prescribed by the Pauli operators H in the respectively transformed Hamiltonians. Each circuit contains two single shot entangling gates "MS" that make use of the all-to-all connectivity of our device, making it particularly resource efficient for the larger register size of the four qubit implementation. Following steps outlined in Appendix B, the effective Hamiltonian under the BK transformation becomes where the coefficients c were all derived in this classical pre-processing step and c 0 specifically captures the spatially fixed nuclei's Coulomb potential. It correspondingly mandates three projective measurement settings to obtain the associated set of expectation values In the JWtransformed case, we are to obtain 14 expectation values

B. Results
With only a single parameter θ in both circuits, it is possible to efficiently scan the complete parameter space. Figure 2.c compares experimental results of such a parameter scan for the two qubit case under the BK transformation with theoretical predictions from a noise-free simulation of the circuit. In this way, measurements of a limited number of values for θ can be extrapolated to estimate the entire parameter range via Gaussian process regression [50], or by fitting sinusoidal functions to either the expectation values or the resulting energy landscape, as was done here. Once the expectation values H (θ) are known over the entire parameter space spanned by θ, one can calculate the molecular energies over the given range of internuclear distances R yielding A ground state potential energy curve is then obtained by performing the energy minimization procedure in post-processing, effectively generating an experimental reference for the average performance of our simulator. The resulting energies are usually expressed in Hartree in chemistry, where 1 Ha = 2 /mea0 ≈ 27.2 eV, with Planck's constant , the mass of an electron m e and the Bohr radius a 0 .
Under the influence of the MS gate, both the Hartree-Fock reference states |01 and |0011 are transformed through a decoherence-free subspace [97][98][99] that is protected against correlated dephasing -a leading source of error in trapped ion qubit implementations. To experimentally investigate the impact of this decoherence channel on the results, we also implement a simulation that operates outside of the protected subspace. This is achieved by performing a basis rotation at the level of the Hamiltonian, which results in sign changes of the c (R) coefficients and a different initial state. In the Bravyi-Kitaev case, this corresponds to simply changing the ansatz function to |ϕ BK (0) = |11 , while still maintaining the same UCCSD operator and circuit implementation. More details on the construction of this particular decoherence free subspace (DFS) on our architecture can be found in reference [100].
In total, we implement four cases in our demonstration, whose results are depicted in Figure 3.a and b on the next page. They include a DFS-protected and an unprotected implementation of the Bravyi-Kitaev mapping with a two-qubit MS gate fidelity of 99(3)% as well as a Jordan-Wigner implementation using a four-qubit MS gate with either 97(4)% or 93(3)% fidelity. In each case, the MS gate entanglement fidelity was estimated from population averages and the parity contrast relating to the coherence of the respective Bell or GHZ state generated from |11 or |1111 at the conclusion of the operation [101][102][103]. We observe that the absolute energy values calculated from the parameter scans and illustrated in Figure 3.a are shifted to larger values both under increased correlated dephasing error as well as reduced gate fidelities. In particular, we observe that a two-qubit ansatz state that is not protected against correlated dephasing yields results similar to those of a DFS-protected fourqubit ansatz state. The latter is significantly more sensitive to correlated dephasing due to the larger number of qubits involved, highlighting the benefit of employing decoherence-free subspaces in algorithmic implementations. As measurements in chemistry generally refer to energy differences as opposed to absolute values, it is common to translate the potential energy curves to their nominal reference value at large separation R as illustrated in Figure 3.b. This depiction more clearly reveals the respective upshift in energy with respect to the calculated binding energy (or well depth) and the simultaneously occurring shift in the position of the energy minimum towards larger internuclear distances.
We proceed to implement the VQE algorithm in full at five different internuclear separations R, yielding the results shown in Figure 3.c. For each configuration R, we start at a random initial value of θ 0 , prepare |ϕ BK (θ 0 ) , measure the expectation values H (θ 0 ) corresponding to the terms in the Hamiltonian and pass the measurement results to a Nelder-Mead simplex algorithm running on a classical computer. Here, the energy is calculated according to Eq. (7) before a new value for θ is suggested for the next iteration. In parallel, we execute a noise-free simulation of the circuit at each iteration in order to monitor the convergence towards the theoretically expected energy.
Generally, the algorithm converges in simulation and experiment. Residual energy fluctuations seen in the experimental results (cf. Figure 9.a) are related to (1) measurement errors and (2) noisy gate operations. The first source of error stems from quantum projection noise (QPN) [104] that scales proportional to 1 /r for r repetitions of the circuit (here, r ≤ 1000 for the VQE points). The second source of error is related to the experimental environment and manifests itself, e.g., in laser intensity fluctuations and transient electrical noise coupling to the motion of the ions. Both introduce a loss of fidelity in digital quantum simulations as shown in our earlier work [25] and can be mitigated through technical improvements.
As a result, we did not fix the number of VQE iterations to a specific value, but instead implemented a sinusoidal fit to the 1D parameter space of the energy explored throughout all iterations (cf. Figure 9.b), with each point weighted according to the QPN contributions from its constituent expectation value measurements. Figure 3.c shows each run's result superimposed on the previously discussed parameter scan. Error bars for the VQE points are derived from the above fitting procedure. In some cases, e.g. R = 0.6 shown in the Figure 3.c inset, the simplex algorithm appears to get stuck, which likely is the result of a premature reduction in the step-width of θ caused by the noise sources discussed above. We return to this effect in the next section.

V. LITHIUM HYDRIDE
We now increase the complexity by turning to a heteronuclear molecule with four electrons and aim to sim-z(π) z(α) ulate the ground state energy of lithium hydride (LiH). This requires the introduction of additional variational parameters and thereby increases the circuit depth. LiH is also a natural small-molecule example and was previously simulated using four superconducting qubits in [63]. We implement its simulation using three ion qubits.

A. Encoding the problem
We again begin with the classical preprocessing step, choosing a minimal basis set of Slater-type orbitals represented by linear combinations of 6 Gaussian functions (STO-6G). A Hartree-Fock calculation then leads us to determine an energy-ordered molecular orbital basis into which the molecules' 4 electrons are filled sequentially (Figure 4.a). The corresponding complete UCC ansatz, truncated to single and double excitations (UCCSD), yields 32 single and 168 double excitation operators. Under a Bravyi-Kitaev transformation, a direct implementation of this ansatz would require 12 qubits. Using the same approximation that was employed in the four qubit implementation of [63], we reduce the number of required qubits by identifying an active space of two electrons in three spatial orbitals, making an arbitrary choice in the degenerate subspace. Such active spaces average out (freeze) the core electrons not thought to be involved strongly in the correlations responsible for bonding [105].
Lastly, we take an additional step in order to identify the dominant contributions to the active space at an efficient level of classical theory known as configuration interaction singles and doubles (CISD). In this case, two singlet excitations are found to dominate and their corresponding unitary coupled cluster formulation is approximated as where α, β are two components of the vector θ and correspond to the parameters of the variational optimization. After mapping the above fermionic representation to Pauli operators using the Bravyi-Kitaev (BK) transformation, only three qubits are acted on non-trivially, i.e. in a way that changes state populations. This is related to the previously observed fact that BK computational basis states naturally reflect certain spin symmetries, allowing for a more compact representation in some cases than the corresponding Jordan-Wigner mapping [50,68]. The problem may hence be reduced to its action on these three qubits and simulated by selecting appropriate subterms, acting on the initial state |ϕ(0) = |111 . The final operator we implement has the form Further details of this derivation can be found in Appendix C 1. The quantum circuit corresponding to Eq. (8), which acts on subsets of the qubit register, is shown in Fig. 4.b. We choose to implement it in the experiment using the circuit shown in Fig. 4.c, which is based on a refocusing technique [95]. Here, an addressed π-phase shift between two half-entangling MS gates effectively decouples the addressed qubit from the two remaining qubits: these become entangled, while the phase-shifted qubit does not obtain any correlations with the rest of the register. Other implementation strategies range from other algorithmic solutions such as spectro- The theoretical LiH potential energy surface calculated for the minimal basis set (black) is shown in comparison with experimentally obtained results, offset to overlap at maximum distance R to better illustrate the well-depth differences in the grid scan reconstructions. In absolute values, all measured data is above the FCI calculation. The data points result from sampling the energy landscape H(R) α,β using the VQE algorithm or the parameter scan shown in a) and fitting the explored space with a Gaussian process regression (GPR) based machine learning algorithm (red dashed line) or a 2D quadratic fit (blue solid line). The error bars are obtained from the fits with the underlying data weighted by quantum projection noise. The slight "kink" close to R = 3.5Å is due to the interplay of rounding errors introduced in the fitting routine with small deviations originating in our active space approximation.

B. Results
We first perform a parameter scan to establish a baseline for the performance of our system before implementing the VQE algorithm for select points. Moving sequentially over a section of the two-dimensional parameter space bound by α = [1.5, 6], β = [2,5] in a gridlike pattern, we perform three rounds of projective measurements at each setting to determine expectation values for each term in the BK transformed Hamiltonian: Fig. 12 for data). The results are combined via Eq. (7) in order to calculate the energy landscape for each internuclear separation R. An example for R = 1.6Å is shown in Figure 5.a) with the experimentally measured parameter space superimposed on a theoretical calculation of the full range.
In order to reconstruct the full potential energy curve of the electronic ground state from this data, we investigate two approaches: (1) a two-dimensional quadratic fit to the energy minimum and (2) a Gaussian process regression (GPR) fit. The fit minima from all energy surfaces of the different internuclear separations R finally yield the potential energy curves in Fig. 5.b, which we again use as experiment-based reference.
We now implement an iterative VQE procedure similar to the one described in the case of H 2 above, but with an important modification. Numerical simulations, detailed in Appendix C 3, and experiments reveal that the previously observed convergence failure of the bare Nelder-Mead search algorithm in the presence of noise has a significant impact on the accuracy of the energies obtained in each VQE run. To combat this effect, we switch the optimization to a hybrid algorithm [108] that also incorporates an element of simulated annealing by introducing random perturbations, sampled from a distribution D, that are added to the cost function in Eq. (7). In this way, the VQE algorithm is forced to continuously sample the surroundings of the minimum as shown in Fig. 5.a without converging any further. The larger number of samples effectively allows for a more precise estimation of the minimum's location through a fit to the data. We heuristically choose D to be on the order of the energy error caused by quantum projection noise such that the perturbations only become dominant in the vicinity of the minimum. For example, at R = 1.6Å the energy error from quantum projection noise after 500 repetitions of the experiment is be- we proceed for another 10 − 20 iterations before stopping the outer VQE loop. We then evaluate the data using (1) a two-dimensional quadratic fit to a subset of the VQE iterations (four standard deviations from the median; illustrated in Fig. 14) or (2) a Gaussian process regression fit. Both sets of results are shown in Fig. 5.b in comparison with the ideal theoretical result for our basis set. While resulting in a smooth potential energy surface, the GPR-based fit appears to systematically underestimate the binding energy, which highlights the impact the chosen data evaluation method has in the final step of a VQE algorithm.

VI. DISCUSSION
To put the results presented above into perspective, it is useful to establish a benchmark against which they can be compared. For quantum chemistry calculations, the widely used concept of "chemical accuracy" constitutes such a point of reference.
Chemists are often particularly interested in free energy landscapes which provide mechanistic insight into chemical events of significant practical importance such as drug binding, catalysis and material properties. Free energies are obtained from partition functions which can be sampled using Monte Carlo or Molecular Dynamics simulations, which assume the ability to compute solutions to the electronic structure problem. The resulting energy landscapes must be extremely accurate as chemical rates are exponentially sensitive to changes in free energy and thus changes in potential energies. This sensitivity can be seen from the Erying equation [109] for chemical rates which is proportional to e −β∆G ‡ /β, where ∆G ‡ is the difference in free energy between reactants and the transition state, and β is the inverse temperature in atomic units. At room temperature and atmospheric pressure, an energy error in ∆G ‡ of 43.3 × 10 −3 eV (1.6 × 10 −3 Hartree) translates to a chemical rate that is wrong by a factor of ten [82]. With respect to this threshold, it follows that the expectation values of most operators H in the chemical Hamiltonian need to be determined at a similar level of precision. In practical terms, this translates to a minimum number of measurement averages required to overcome the intrinsic quantum projection noise in the estimation of each Pauli string H . Figure 6 illustrates this requirement for the case of the H 2 molecule, assuming no other sources of noise and knowledge of the final value of each constituent expectation value at every point.
The graph shows that one would need to repeat each of the measurements at least 15000 times in order to reduce the intrinsic measurement noise below chemical accuracy. Given that every repetition incurs the overhead of state initialization, state preparation and measurement, the number of repetitions has a strong impact on run-time. For example, on our trapped ion system the time needed for one repetition is on the order of 20 ms, requiring at least 5 minutes of averaging to ensure a measurement limit at chemical accuracy in each VQE iteration. In the current setup this limitation stems from both state initialization through laser cooling and the duration of quantum gate operations. However, ways to overcome these constraints have been demonstrated in a variety of experiments already [110][111][112], promising order-of-magnitude advances in speed and thereby much shorter runtimes. In addition to technical improvements on the hardware side, an adaptive measurement strategy might also alleviate the resource needs with respect to the required number of averages. This could entail varying the number of repetitions throughout the VQE iterations based on the observed energy changes in each step, gradually increasing the measurement precision as the algorithm converges towards the energy minimum.
In the proof-of-principle experiments above we have focused on probing the effects of system noise from both a limited number of measurements and decoherence effects rather than achieving chemical accuracy. The parameter scans for H 2 (LiH) were taken by averaging 100 (500) repetitions of the corresponding circuit, which for the case of LiH already corresponded to a significant time overhead due to the two-dimensional nature of the parameter space. As a consequence, there is a large uncertainty associated with the energy error of the extracted potential energy curves. The effect of errors resulting from decoherence, seen in Figure 3.a, was modeled for the case of molecular hydrogen using the Bravyi-Kitaev transformation with help of the quantum chemistry package OpenFermion [113]. The result is shown in Figure 7 and largely matches the upshift in energy observed in the experiment.
The VQE method appears to be robust to calibration errors in some of the gates as evidenced by the shift in the parameter value θ in Fig. 2.c, also seen in [50]. However, similar to standard quantum tomography, the method still relies on accurate determination of the measurement operators H , which in turn are bound to local gate operations that rotate the desired basis for each qubit into the measurement basis. Errors at this point in the circuit will be folded into the results and cannot directly be recovered from. Possible ways around this issue might be found both on the quantum and classical side. Quantum control techniques tailored to experimental errors allow for the suppression of certain types of gate errors [114]. An algorithmic modification might be the addition of a second variational optimization loop parametrizing the pre-measurement rotations, which is entered once the nominal VQE outer loop has converged to a desired level. In this way, the measurement basis can be corrected for offsets in the same way as the calibration errors inside the quantum circuit. Lastly, the quantum circuit might be extended to include ancilla qubits, which allow for error detection based on stabilizers that indicate whether the computational subspace was still maintained at circuit completion [115]. In this way, post selection might be used to raise the effective fidelity of non-error corrected quantum devices to a level that makes them useful for quantum chemistry calculations in the near-term future.

VII. CONCLUSION
In summary, we have performed the first multi-ion quantum simulation of quantum chemistry. In our simulation of molecular hydrogen, we have employed both the Bravyi-Kitaev as well as the Jordan-Wigner transformation in the mapping from fermions to qubits and varied the respective Hartree-Fock reference state to reveal the impact of decoherence. Unstable behavior of a classical Nelder-Mead optimization routine was circumvented in our simulation of lithium hydride by amending the variational optimization with simulated annealing and subsequently employing a quadratic fit to estimate the location of the minimum in the two-dimensional energy landscape. When scaling up to larger parameter spaces, this step of the classical data evaluation routine will likely become an even more important issue.
Despite a significant increase in experiment runtime, it was still possible to sequentially scan the two-dimensional parameter space of LiH. However, the addition of additional excitation operators for more accurate energy determinations and the scaling up to more complex molecules reveals the true power of the VQE approach: sparse and adaptive sampling of a potentially very high dimensional energy landscape. While the number of different measurements needed to estimate the energy of an arbitrary molecule can be large, recent work shows how using structure fundamental to the fermionic nature of particles can reduce this number by an order of magnitude or more for minimal additional complexity [59].
To deliver on the promise of near-term usefulness, further work is needed on both the quantum and classical aspects of the VQE algorithm. The mitigation of errors in short depth quantum circuits has increasingly come into focus [55,59,64,[116][117][118] and the use of dynamic error suppression techniques at the physical gate level promises improved fidelities of the practically implemented operations. A recent numerical study [119] specifically investigates the impact of errors in the simulation of the ground state energy of H 2 and LiH and illustrates the expected impact on reaching chemical accuracy.
Furthermore, there are a number of recent proposals in the literature which appear to significantly reduce the resources required to scale up these simulations. For instance, in [40] authors introduce a new class of basis functions for the simulation of electronic structure. They quadratically reduce the number of terms in the Hamiltonian from O(N 4 ) to O(N 2 ), which also reduces naive bounds on the number of measurements from O(N 8 ) to O(N 4 ) [51]. This speeds up data collection by making the number of circuit repetitions substantially more practical. The number of measurements required can also be improved further using recent techniques for enforcing n-representability conditions (known constraints on the geometry of fermionic states) from [59].
The most significant challenge in realizing variational algorithms in the near-term remains the circuit depth required. A naive application of unitary coupled cluster requires a number of gates that scales as O(N 4 ) assuming arbitrary connectivity between qubits. This strongly suggests that more scalable variational circuits will be needed if we are to approach classically intractable calculations without error-correction. One alternative to uni-tary coupled cluster is to use an unstructured variational circuit, as demonstrated recently in [63]. Recent work in [120], however, has shown that such strategies become exponentially expensive as a consequence of concentration of measure in random quantum circuits. Nevertheless, there has been significant recent progress in realizing low depth structured variational ansatze; e.g., in [41] the authors introduce variational circuits based on Trotterized adiabatic state preparation which can be implemented with linear gate depth even on a device with extremely limited (e.g. next-neighbor) qubit connectivity. In [121], a similar ansatz is introduced, also with linear gate depth which is equivalent to the fermionic swap network from [41] with a different variational parameter initialization.
It is far from a forgone conclusion that we will one day solve classically intractable problems in quantum chemistry without error-correction. Yet, by using these ever improving algorithms combined with error-mitigation strategies and simpler Hamiltonian representations, it is certainly the case that we should be able to push forward and obtain reasonably accurate quantum simulations of molecules with tens of spin-orbitals (qubits) in the near future. In comparison to other quantum simulations with less stringent requirements, the notion of chemical accuracy in this context provides a clear benchmark to determine the point at which support by quantum processors transitions into a useful regime. In this way, metrics like the deviation in well-depth or the worst case deviation along the entire potential energy surface, captured by the so-called non-parallel error, might in fact provide a convenient way to measure multi-qubit performance in different architectures, which will help assess technological progress. The qubits used for this experiment are encoded in Zeeman sublevels of a metastable D and an S ground state of 40 Ca + ions confined in a macroscopic linear Paul trap described in detail in [122] and shown in Figure 8.a. In particular, we associate |0 = |↑ = |3d 2 D 5/2 (m = +3/2) and |1 = |↓ = |4s 2 S 1/2 (m = +1/2) .
Both electronic levels are connected via an optical quadrupole transition at 729 nm, which is used to drive the qubits via an ultra-stable laser with ∼ 1 Hz linewidth. The same laser mediates multi-qubit interactions via the joined motional modes of the ions in their trapping potential. An auxiliary transition at 397 nm between 4s 2 S 1/2 and 4p 2 P 1/2 is used for laser cooling, state preparation via optical pumping and state detection. Each experimental repetition consists of laser cooling, ground state initialization (sideband cooling to the motional ground state of the ion strings center of mass mode and optical pumping to the electronic state |1 ), execution of the desired gate sequence and finally parallel state detection for each qubit via an electron multiplying CCD camera (Figure 8.b). Every sequence is repeated at least 100 times to calculate an expectation value from the observed probabilities. The molecular hamiltonian for H 2 in the STO-3G minimal basis can be mapped from its fermionic form (2) to qubits yielding where the coefficients f i depend on the internuclear separation (R) between the hydrogen atoms and are derived from the integrals in (3). The reference state for the calculation corresponds to the Hartree-Fock solution |ϕ HF = |0001 . We observe that the terms in the Hamiltonian only act with the identity and σ z operations on qubits 1 and 3. This fact allows us to rewrite our reference state as |ϕ HF = |0 1 |0 3 ⊗ |0 2 |1 0 . As qubits 1 and 3 will not experience population changes under the Hamiltonian, we can reduce Eq. (B1) to an effective Hamiltonian acting on two qubits, with reference state |ϕ HF = |01 : Here we have relabeled qubits 0 and 2 as 0 and 1. The coefficients c i are now given by: This reduction of the problem for the hydrogen molecule was first noted in [50], developed into a general method in [68] and used in superconducting implementations of several problems in [63].

Jordan-Wigner (JW) transformation
The molecular Hamiltonian under the Jordan-Wigner transformation gets mapped to with coefficients c i again derived from the integrals (3). Under the Jordan-Wigner transformation all the qubits are used to store occupation numbers while in in the Bravyi-Kitaev transformation even qubits store occupations and odd qubits keep track of the parity of all the qubits with smaller indices. Hence, four qubits are needed to encode the ansatz state |ϕ HF = |0011 .

Application of Unitary Coupled Cluster to H2
For H 2 in the minimal basis, the second quantized formulation of the Unitary Coupled Cluster operator for single and double excitations corresponds to where θ 23 01 is the coupled cluster amplitude that is variationally optimized. Note that in this case the singleexcitation operators are effectively incorporated in the basis we are using (i.e. the single excitations rotate the basis and do not need to be applied explicitly in the circuit). Using the BK mapping this operator is expressed as follows: As the terms in U (θ 23 01 ) only act with the identity and σ z on qubits 1 and 3, the operator can be reduced to 2 )I + i sin( θ 23 This allows us to define the ansatz via the U (θ 23 01 ) operator simply as Note that this form is only valid when the operator acts on the reference state |01 .

Implementation of unitary coupled cluster operator using MS gates
In order to implement the UCC operator (B3) above using Mølmer-Sørensen gates, we employ a technique first demonstrated in Müller et al. [95], formulae 10-12. using the fact that A and B do not commute and therefore must anticommute we obtain: Specifically, for α = π/4 and the case above 5. Data for molecular hydrogen theoretical simulation and appears to converge after 15 iterations. In order to maintain the same analysis that was used to accommodate cases in which the optimization routine got stuck, the corresponding energy minimum is extracted by fitting a sinusoidal function to the parameter space explored by all iterations (see Fig 9.b). The vertical displacement between theory simulation and experiment fits are likely due to noise. The horizontal displacement is likely due to calibration drifts, leading to a different α value having the smallest energy.

Decoherence simulation
To understand the effect of various decoherence channels, we performed a simulation of the entire circuit of the two-qubit H 2 experiment using the open source framework OpenFermion [113].
We assume that throughout the execution of the entire state preparation circuit, the qubits experience dephasing, e.g. due to magnetic field fluctuations induced by the environment. We model the dephasing via an i.i.d. channel of the form is a Kraus map and p d is the probability for a single phase flip. In our simulations, we applied the dephasing channel to all the qubits after the application of each gate in the circuit of Figure 2.a, with probability p d = 1.0 − exp(−T g /T 2 ), where T g is the physical time of the gate and T 2 is the dephasing time. We employed T 2 = 40 ms, as determined from Ramsey experiments on a single ion. In addition to dephasing, we model the effect of errors in the MS gates using a two-qubit depolarizing channel. Here, we consider all single-and two-qubit errors with the same probability. For a 2-qubit MS gate, the noise is described by the quantum operation where p M S is the probability of a MS depolarizing error, and Λ a and Λ α correspond to the set of indexes for the active ions and Pauli matrices, respectively. For the 2-qubit MS gate, there are 15 possible Pauli errors (6 single-qubit and 9 two-qubit), resulting in the prefactor 1/15. The probability p M S is related to the fidelity of the gate as F = 1 − 14 15 p M S [123]. The experimental fidelity estimated for the 2-qubit MS gate is 0.99. In the simulation, the two-qubit depolarizing channel is applied with probability p M S after the MS gate. Figures 7 and 10 display the results of the simulation for H 2 under the BK mapping. Our simulation appears to account for the observed experimental errors along a significant portion of the energy curve. We observe an uneven upshift of the energy values and effectively reduces the estimated well depth, respectively binding energy. While our simulation are close in magnitude to the observed results, we note that other factors such as faulty measurement operators (related to basis rotations and detection fidelity) could also contribute to the discrepancies. Internuclear distance R (Å) the natural orbital basis (NMO). NMOs are obtained by diagonalizing the exact one-electron reduced density matrix (1-RDM) of the system and are ordered by natural orbital occupation numbers (NOONs). It has been shown that NMOs with small NOONs or NOONs close to fulloccupancy have a negligible effect in the electron correlation and therefore can be discarded [124]. Consequently, approximate NMOs and NOONs, obtained from perturbation theory or truncated configuration interaction (CI) calculations, are usually employed to reduce the computational cost of more involved correlated calculations and for the selection of active spaces [125].
After the BK transformation, the Hamiltonian for LiH comprises 193 terms with amplitudes larger than 10 −10 Hartree. A VQE simulation using the UCC ansatz truncated to single and double excitations requires 12 qubits and involves 32 single excitation operators and 168 double excitation operators (without imposing spin constrains). To reduce the number of excitation operators in the calculation, we employed the NOONs derived from a configuration interaction calculation with single and double excitations (CISD) in order to select an appropriate active space. Fig. 11 shows the NOONs for the six molecular orbitals of LiH, calculated using CISD for four different internuclear separations. Based on the variations in the NOONs, we establish orbitals 1 to 4 as an appropriate active space.
A reasonable choice of excitation operators in the selected active space would be the singlet double excita- tions from orbital 1 to orbitals 2, 3 and 4, respectively. Based on the NOONs, we expect the amplitude of the excitation operator from orbital 1 to orbital 2 to be largest. Similarly, we expect excitation operators from orbital 1 to orbitals 3 and 4 to have the same or similar amplitudes. Due to constrains in circuit depth, we consider only excitations from orbital 1 to orbitals 2 and 3, explicitly the operators: a † 5 a † 4 a 3 a 2 −a † 2 a † 3 a 4 a 5 and a † 7 a † 6 a 3 a 2 −a † 2 a † 3 a 6 a 7 , where a † i (a i ) denote the creation (annihilation) operator in the i-th spin-orbital. In our notation, spin-orbitals with odd (even) indices correspond to spin-up (spindown) electrons, with indices starting at 0. Using the BK mapping, these operators can be expressed as The initial state of the simulation, corresponding to the Hartree-Fock wavefunction, is the state |000000000101 , which simplifies to |000001 in the active space required for the selected excitation operators. A full-simulation of LiH with the two double excitation operators listed above would require at most 32 MS gates for a single Trotter step if all the subterms were going to be implemented.
To make the simulation affordable in the current device we approximated each of the operators using only the first subterm, corresponding to σ X 2 σ Y 4 and σ X 2 σ Y 6 , respectively. The absolute error in the total energy introduced by this approximation is smaller than chemical accuracy within the basis set used, when compared to the full configuration interaction (FCI) reference solution. The two selected subterms can be implemented using MS gates and single qubit rotations as follows where R z represents a rotation around the Z-axis and U S M S is an MS gate acting on the set of qubits S. As the entangling operations involve only qubits 2, 4 and 6, we can efficiently construct an effective Hamiltonian involv-ing only operations on these qubits. The corresponding 3-qubit Hamiltonian has the form H =c 0 I + c 1 Z 0 + c 2 Z 1 + c 3 Z 2 + c 4 Z 1 Z 0 + c 5 Z 2 Z 0 + c 6 Z 2 Z 1 + c 7 X 1 X 0 + c 8 Y 1 Y 0 + c 9 X 2 X 0 + c 10 Y 2 Y 0 + c 11 X 2 X 1 + c 12 Y 2 Y 1 .
The applied pulse sequence to realize the Hamiltonian in Eq. (C1) is shown in Figure 4.b acting on initial state |111 .

Parameter scans for LiH
The 2D LiH energy landscape H(R) α,β , as shown in the inset of Fig. 5.a), is obtained in the following way: We apply the pulse sequence, shown in Fig. 4.c) to the initial state |ϕ(0) = |111 . The MS gate fidelity for three ions reaches 97(3)% and the experiment is performed with 100 repetitions. The energy depends on two parameters α, β, both corresponding to a rotation angle on σ z , applied on one of three ions. These angles are scanned in the range of α ∈ [1.5, 6] , β ∈ [2, 5] with resolutions of 0.1 and 0.15 respectively. The measured expectation values are shown in Fig. 12. The 2D energy landscapes H(R) α,β , as shown in Fig. 5.a), are finally calculated by combining the measured expectation values according to Eq. (C1) and adding the contribution from the Coulomb interaction of the nuclei H nuc (R) (Eq. 7).
The final molecular energy curve H(R) for LiH, as shown in Fig. 5.b, is derived by fitting two-dimensional functions on the 2D energy landscapes and extracting the minimum value.

VQE runs for LiH
In our first implementation of the VQE we employed the Nelder-Mead optimization routine in the same way as for molecular hydrogen. In the experimental results shown in Figure 13.a it appears that the optimization algorithm becomes trapped in a local minimum, preventing it from converging any further.
To investigate this phenomenon, we simulate the experiment, including quantum projection noise, and run the a full VQE simulation 10 times at a fixed separation R using the same values for initial guess of parameters α and β. The simulation results are shown in Figure 13.b. For better visibility, we plot only the last 20 points for each run, after which the algorithm converges to a "stable" position. The black dot represents the optimal combination of parameters, predicted by the theory. One finds, that each of the 10 simulated runs is trapped in a different local minimum, leading to a deviation from the ideal energy values. We conclude that the basic Nelder-Mead algorithm is not suitable for this optimization problem. We therefore instead employ a hybrid of the Nelder-Mead direct search algorithm and simulated annealing theory [108] to obtain a better estimate of the energy H(R) α,β in the noisy environment. We now fit a quadratic function to the 2D landscape and extract the minimum value as H(R) . Fig. 14 shows an example of fitting the 2D quadratic function E(α, β) = m + (c · α − a) 2 + (d · β − b) 2 , with the minimum energy value m, to a subset of the VQE iterations, taken to include those within 4 standard deviations from the median. We perform VQE runs for two different internuclear separations R = 1.6, R = 2.75, with 3 ion MS gate fidelities of 98(5)%. The results are shown in Fig. 5.b, together with the parameter scan discussed in the previous section.