Efficient Two-Electron Ansatz for Benchmarking Quantum Chemistry on a Quantum Computer

Quantum chemistry provides key applications for near-term quantum computing, but these are greatly complicated by the presence of noise. In this work we present an efficient ansatz for the computation of two-electron atoms and molecules within a hybrid quantum-classical algorithm. The ansatz exploits the fundamental structure of the two-electron system, and treating the nonlocal and local degrees of freedom on the quantum and classical computers, respectively. Here the nonlocal degrees of freedom scale linearly with respect to basis-set size, giving a linear ansatz with only $\mathcal{O}(1)$ circuit preparations required for reduced state tomography. We implement this benchmark with error mitigation on two publicly available quantum computers, calculating accurate dissociation curves for 4- and 6- qubit calculations of ${\rm H}_\textrm{2}^{}$ and ${\rm H}_\textrm{3}^+$.


I. INTRODUCTION
Quantum computers possess a natural affinity for quantum simulation and can transform exponentially scaling problems into polynomial ones [1][2][3]. Quantum supremacy, the ability of a quantum computer to surpass its classical counterpart on a designated task with lower asymptotic scaling, is potentially realizable for the simulation of quantum many-electron systems [4,5]. Work over the previous decade has been towards this goal with a focus on calculating the energy of small molecules and exploring strategies to leverage emerging quantum technologies, especially those designed to correct or mitigate quantum errors [6][7][8]. In this paper we introduce an efficient ansatz for a two-electron quantum-mechanical system that can be employed as a benchmark for assessing the capabilities and accuracy of quantum computers. The twin goals of the work are: (i) to present a quantum-computing benchmark based on the correlated but polynomial scaling two-electron problem, solvable on classical computers, that can be used to assess the accuracy of quantum computers, and (ii) to develop an efficient ansatz for solving the two-electron problem on quantum computers, based on an effective partitioning of the computational work between classical and quantum computers that is applicable to more general N -electron molecular systems.
The two-electron density matrix (2-DM) of any twoelectron system can be expressed as a functional of its oneelectron reduced density matrix (1-RDM) and a set of phase factors. This representation of the 2-DM has important connections to natural-orbital functional theories and geminalbased theories in quantum chemistry [9][10][11][12][13][14][15][16]. It offers a natural separation between the nonlocal and local fermionic degrees of freedom in the system [17], scaling linearly and polynomially respectively, and can be leveraged in a variational hybrid quantum-classical algorithm. The entangled, nonlocal degrees are treated on the quantum computer while the local degrees are treated on the classical computer, leading to an efficient simulation of the system.
For a quantum algorithm to exhibit quantum supremacy, obtaining the solution classically will be impractical except * For correspondence, damazz@uchicago.edu for cases that are close to the classical limits of feasibility [4,5]. For some problems such as prime factorization, the solution can be quickly verified, but for many-body quantum systems this is not the case [18][19][20][21][22]. Possessing 'easy', classically solvable problems to implement and verify will be crucial to evaluate the performance of quantum devices and error mitigation schemes [23]. Our proposed quantum-classical hybrid algorithm targets only the necessary entanglement needed on the quantum computer, scales linearly with respect to basis size, and has O(1) circuit preparations, making it an ideal benchmark for molecular simulation. We highlight this by evaluating this ansatz through the computation of H 2 and H + 3 on two generations of publicly available quantum devices.

II. THEORY
Because the representation of the 2-DM in terms of the 1-RDM and phase factors has been well studied elsewhere [11,[24][25][26], we present in section IIA only the aspects of the theory that are relevant to the quantum-classical hybrid algorithm in section IIB. We also discuss the preparation of the linearscaling ansatz in section IIC and practical error mitigation techniques in section IID that are employed in the benchmarks in section III.

A. Structure of the Two-Electron System
For a two-electron system the energy is given as the trace of the Hamiltonian and the density matrix: where H and D are 2r × 2r with 2r being the rank of the one-electron basis set and in which the wavefunction expansion coefficients g ij are elements of the coefficient matrix G. From the antisymmetric nature of fermions, G must be a skew-symmetric matrix, and from a theorem by Zumino [26], G must have a blockdiagonal formG with 2 × 2 matrices G i : The block-diagonal form of G in Zumino's theorem defines an orbital basis set with a natural pairing of the orbitals where we denote the indices of an orbital and its pair by i and i ′ , respectively. The 2-DM in Zumino's basis has only nonzero elements of the form:D The 1-RDM, containing the one-body information, can be obtained from the 2-DM by contraction: 1Di where all 2-DM elements in the contraction vanish except when k = i ′ . Because the 1-RDM is diagonal in Zumino's basis set, we find that Zumino's basis set is identically the natural-orbital basis set and that the occupations n i are the natural orbital occupations. The paired orbitals i and i ′ , we observe, have equal occupations n i and n i ′ . For a system of two electrons with S z = 0, these paired orbitals share the same spatial component with different spin components, denoted by convention as α and β. This decomposition can be viewed as a particular case (N = 2) of a more general result derived by Schmidt [27], and later by Carlson and Keller [28].
The importance of this decomposition as a quantum computing ansatz for two electrons will be manifest below.

B. Variational Hybrid Algorithm
Hybrid quantum-classical algorithms with a variational eigenvalue solver are among the most promising algorithms for near term applications [6,[29][30][31]. For our approach we utilize a variational eigensolver on the quantum device but apply it only to the optimization of the 2-DM in the naturalorbital basis set. Optimization of the natural orbitals by orbital rotations is performed with polynomial scaling on the classical computer. In this manner we are able to partition the nonlocal and local degrees of freedom between the quantum and classical calculations respectively.
Using the structure given in Eqs.
(3)- (7), we see that to evaluate theD matrix, we need only: (1) the natural orbital occupations, measured as the orbital populations on the quantum computer, and (2) the phase corresponding to the natural orbital coefficients, which for a real wavefunction, is simply the parity of the term that can easily be measured on a quantum computer. That is, we need the phase ξ ii ′ wherẽ In general, the phases can be measured through tomography on the quantum computer of certain terms ofD requiring only O(1) additional circuit preparations. The details of the specific ansatz for the tomography are discussed in the next section.
After convergence criteria in the optimization ofD on the quantum computer is satisfied through gradient-free optimization (see Appendix A), we optimize the energy on the classical computer through one-body unitary transformations of the Hamiltonian. Specifically, we optimize the orbitals through a series of Givens rotations. These complementary quantum and classical optimizations are sequentially repeated until the energy and 2-DM converge.

C. Preparation of the Efficient Quantum Ansatz
To create a state of the form in Eq. (3) on the quantum computer, we need to implement double excitations from orbitals ii ′ to kk ′ . If we consider an initial wavefunctionG 0 from a standard Hartree-Fock calculation, the ansatz to generate a genericG is:G whereī = i + 1, and tīī ′ ii ′ is an antihermitian, antisymmetric two-body matrix with nonzero elements corresponding to an excitation between i, i ′ and (i + 1), (i + 1) ′ . The operators acting onG 0 can be easily expressed in second quantization as shown in the Appendix. The ansatz is a subset of the unitary coupled cluster (UCC) [32] or antihermitian contracted Schrödinger equation (ACSE) ansatz [33]. From there we perform a Jordan-Wigner transformation (though others may be utilized) which yields an exponential of Pauli strings that are implementable on a quantum device as strings of CNOT gates [34]. Additionally, the implementation naturally requires only a nearest neighbor connectivity among qubits.
Finally, the tomography of the state involves only the measurement of the orbital occupations for a given qubit in the computational basis as well as the sign. Because there are only r − 1 phase terms since the last phase is equivalent to a global phase, we only require the tomography of a linear number of sequential terms inD of the formD jj ′ ii ′ with j = i + 1. Because terms likeD jj ′ ii ′ andD (j+2)(j+2) ′ (i+2)(i+2) ′ are qubit-wise commuting, we can measure r/2 terms simultaneously, leading to a constant number of circuit preparations. In this work we evaluated H 2 and H + 3 using the Jordan-Wigner transformation in 4-and 6-qubit cases. There is only one phase term in the 4-qubit calculations of H 2 , which we measured by direct tomography on the quantum computer, whereas the two phases of the 6-qubit calculation of H + 3 were computed by optimization on the classical computer to avoid degradation from noise on the quantum computer.

D. Error Mitigation Strategies
Even if we model a two-electron system on a quantum computer and construct the state through the above tomography, we may find that our occupations n i and n i ′ do not match for a given i, which implies a violation of the fermion statistics. Because the two-electron ansatz in Section IIB formally guarantees a two-electron wavefunction, any deviation in pure-state N -representability (up to sampling errors) is due to errors on the quantum computer [11,35,36].
The effect of errors on current quantum computers can easily influence the N -representability of a system [11,37] with the extent somewhat depending on the fermionic mapping. For a compact mapping, S z and N will typically remain constant, but for more general mappings, this is not the case. Other errors can also accumulate, making it difficult to reach certain extrema of the set of density matrices. To address this, we use a projective technique where we map the set of accessible points onto the ideal set of points (see the Appendix). We achieve this by finding an affine transformation A that maps from the accessible but error-prone set S ′ to the ideal set S.
For a general mapping, it is easy for the quantum system to violate N and S z . By utilizing a symmetry verification technique [20,21], along with the structure of our tomography requiring only measurements of diagonal terms, we can filter out results which do not obey the correct N and S z values which effectively projects the resulting state into an eigenstate of the chosen operator. This can be extended to other operators S which commute with the Hamiltonian: As will be seen in the results, the symmetry verification is useful in bringing the results back to the set of all two-electron states, and then the projection restores the equality of the two pairing-related sets of occupations, n i and n i ′ .

III. RESULTS
Using the two-electron ansatz, we first treat the molecular dissociation of H 2 in a minimal Slater-type-orbital-expandedin-three-Gaussians (STO-3G) basis set of two electrons in four orbitals. The quantum algorithm is implemented on both the 5-and 14-qubit devices, denoted as ibm-5 and ibm-14, representing two generations of superconducting quantum devices by IBM. With the Jordan-Wigner transformation [38] the system can be represented with 4-qubits though more compact mappings are certainly possible. Note in this basis only a single excitation is possible. Figure 1 shows the potential energy curve of the H 2 molecule, computed with full error mitigation.
With the help of error mitigation techniques, both quantum devices are able to capture the dissociation of the molecule, achieving mhartree accuracy across the spectrum of states, leaving differences in the devices somewhat unclear. Inspection of the device calibration (Appendix A, double excitation are shown in Figure 2. The ibm-5 device is seen on the top row, and the ibm-14 on the bottom row, and we also show the effect of the symmetry verification in correcting the occupations. While the ibm-5 device maintains continuity with respect to t 22 ′ 11 ′ , it has distinct problems. First, we observe (see top left insert) correlated measurement error in the set of qubit occupation {n i } which cause an inversion in the expected relationship between n 1 and n 2 among the i (α) occupations. Second, we find that symmetry verification is effective in increasing the differentiability of the two states (note for a decohered system, n 1 and n 2 would be identically 0.5), yet it is not able to correct for the reversal in ordering seen on the qubits, and somewhat reinforces this error for n 2 .
On ibm-14, while there is still a contraction of the occupation numbers across the range of t 22 ′ 11 ′ , the obtained curve is continuous, and the two sets of occupations correspond to expected values. With symmetry verification, we obtain nearly the ideal occupations, which otherwise are far from spanning the full spectrum of occupations. The effect of our further correction is to stretch the sinusoidal curves in Fig. 2 so that their maxima and minima are 1.0 and 0.0 respectively, hence it is not shown here.
To quantify the effects of symmetry verification, we calculate the area between the two changing orbital occupations for both {n i } and {n i ′ } (denoted as V i and V i ′ , where V i = n 2 − n 1 , and V i ′ = n 2 ′ − n 1 ′ , ) subject to different symmetries, as well as the uncertainty in measurement after each symmetry is applied, and show these in Table 1. The maximum and minimum values for this metric would be 2 and 0,  representing fully error-free and fully decohered states respectively. In each case, application of multiple symmetries serves to increase the resulting 'reach' of the state, without increasing variance with respect to decreased measurement counts.
Computations with 6 orbitals are performed only on the ibm-14 device since more than 5 qubits are required. Using simplifications seen in Nam et al. [39], we are able to construct a gate with 8 CNOT gates which still requires only neighboring connections (see Appendix C). Due to the longer depth of the circuit, the effects of noise are more pronounced and we find it difficult to reliably measure the phase of the 2-DM terms required in Eq. (8) and note that these are not always continuous. To show the overall effect of errors on the 6-qubit system on the local occupations, we present a scan of possible symmetry verified 1-RDMs over a range of the parametrized entangling gates in Fig. 3. Additional details regarding the computation are provided in the Appendix. While we do not show the occupations in terms of the parameters, the effect of the aggregate errors for this case is again to shrink the portion of the hyperplane accessible to the quantum device.
Expectedly, the obtained results differ greatly depending on the qubits and available connectivity of the quantum device, though in general we still observe a degree of continuity in the local 1-RDM properties. By dealing with the phases classically, we are able to calculate the dissociation curve for triangular H + 3 in Fig. 4. Again, we are able to obtain chemically accurate energies across the dissociation curve, although there was difficulty in sampling the uncorrelated Hartree-Fock state which is a vertex of the polytope. The error mitigation techniques we use also extend the capabilities of noise-limited quantum computers, which otherwise do not span the ideal Nrepresentability of the state [40].

IV. DISCUSSION AND CONCLUSION
In this work we present an ansatz for two-electron quantum systems which can be implemented on near-term and future quantum computers. Applying this on two public-access quantum computers highlights the successes and differences of two generations of quantum computers, as well as the difficulties which must be overcome in approaching more complicated systems. We also show that using error mitigation strategies we are able to simulate both H 2 and H + 3 to high accuracy. The proposed ansatz can be readily applied to twoelectron atoms and molecules in larger basis sets with similar types of error mitigation. The gate sequence proposed here can be applied as a generic ansatz, removing the need for long expansions of the required exponential operators.
The ansatz is efficient mainly in regards to its scalability with respect to other methods of state preparation. Unitary coupled cluster, which for a two-electron system only needs to be expressed with single and double excitations (i.e., UCCSD), and which additionally has only one occupied orbital for the α and β excitations, still will have O(r 2 ) terms in the ansatz. Furthermore, quantum tomography would additionally require the measurement of the αβ block of the 2-DM, which naively has O(r 4 ) terms, and hence, the method would scale as O(r 6 ) where the depth d and the number of measurements contribute factors of r 2 and r 4 , respectively. Other methods based on the propagation of the Hamiltonian could be implemented, but they also have large costs and would not necessarily lead to the same advantages that result from the structure of the 2-DM.
From the results, it is clear one cannot rely solely on the energy or other external molecular properties to investigate the integrity of the quantum device, particularly in comparing the performance of the two-qubit devices on the 4-qubit calculation. While averaged or localized metrics related to qubit depolarization, dephasing, or bit-flip errors are often used as indicators of performance of the quantum device, they may not translate directly to the fidelity of a simulated fermionic system, particularly with multi-qubit or environmental effects. Looking at these metrics in conjunction with the physical properties of benchmark problems like the ground-state energy of two-electron atoms and molecules will yield greater insights into the fidelity of a quantum device and its needs for error mitigation or correction.
The method of symmetry verification for error mitigation is useful in that its different forms are low cost and can easily correct flagrant faults in the output, such as particle count (or more generally, parity) and the projected spin symmetries. Others forms of error mitigation related to the reduced density matrices (RDMs) can be implemented as constraints on the tomography from N -representability, or in a form of postcorrection of the two-electron RDM (2-RDM) where the measured 2-RDM is purified through semidefinite programming, which can be applied to arbitrary N -electron systems [40]. The mapping we use can have difficulties when errors begin to change the ordering of occupations for larger and larger systems.
The electron pair itself plays a key role in such phenomena as superconductivity and bonding, and yet the exact energy of a two-electron system itself cannot be solved exactly with known methods. Such a problem shows the essence of the electron-electron interaction as well as some of the complexities of electron correlation and quantum mechanics. The theory in this work could also be seen as a subset of more complex geminal-based wavefunction methods, which appear in classical electronic structure theory where the electron pair is treated as the fundamental unit to improve the accuracy beyond the mean-field approximation [9][10][11][12][13][14][15][16][41][42][43]. The exploitation of the structure of the 2-DM in the natural-orbital basis set can be extended to more general pairing 2-RDMs for efficient implementations of pairing (geminal) theories or natural orbital functional theories on quantum computers.
Here we show that the properties of the two-electron system lead to an ansatz which is well suited for use in a hybrid quantum-classical approach, where degrees of freedom that would increase exponentially with N are treated on the quantum computer, while non-exponentially increasing degrees of freedom are treated on the classical computer. The treatment of the orbital rotations on the classical computer can be generalized to N -electron molecular systems where the orbital rotations can be used to implement active-space methods where orbital rotations are used to optimize the correlation in a modest subset of total orbitals known as the active orbitals. Generally, such orbital-rotation algorithms including active-space self-consistent-field algorithms [44,45] could assist in achieving quantum supremacy by further lowering resource requirements on the quantum computer as well as tomography and measurement costs. The two-electron system also is useful in its own right as a benchmark for molecular simulation on a quantum computer, where it serves as a simple yet effective way to assess the performance of an arbitrary quantum device. The present work is clearly applicable for the current state of noisy quantum computers, but will continue to be relevant as improving generations of quantum computers are developed. The two-electron ansatz, albeit polynomially scaling even on a classical computer, can serve as a powerful benchmark for quantum computers due to the availability of accurate results from classical computers and the requirements shared by its solution and the solution of exponentially scaling many-electron problems. The electronic structure package PySCF [46] was used to obtain the one-and two-electron integrals and to perform restricted Hartree-Fock and full configuration interaction (FCI) calculations.
For the quantum computation we used the IBM Quantum Experience devices Yorktown (Sparrow, 5-qubits) and Melbourne (Albatross, 14-qubits), available online. The former has triangular-type coupling between qubits, and the latter has square-type coupling between qubits. These cloud accessible quantum devices are fixed-frequency transmon qubits with coplaner waveguide resonators [47,48]. The quantum information software development kit QISKIT was used to interface with the device. We include the calibration data in Table II. For the 6-qubit case, 2 11 measurements were obtained, and we used a simple Nelder-Mead simplex method with the Han TABLE II. Calibration data for the ibm-5 and ibm-14 devices during benchmarking. U2 and U3 represent the errors for single qubit unitaries containing one and two X π/2 pulses and two and three frame changes respectively. RO represents the readout error, and we have the standard T1 depolarization and T2 dephasing times.
[j] specifies the target qubit with control qubit i, and we report the error. initial simplex for the 6-qubit case on the quantum computer. Classically orbital rotations were performed with Givens rotations, with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm being utilized. Convergence criteria were more strict for the 4-qubit case, with convergence between the 2 D and 2 K steps being 1 mH. In the 4-and 6-qubit cases, we chose the best of 2 runs as the optimal results. While this did not make a difference for most points, for those it did, the difference in energies for the 2 runs were usually significant, indicating that noise had led the optimization into some local minima.
where 1 ≤ j ≤ r and H(x) is the Heaviside step function. The vertices for ∧ 2 H 4 are given as v 2 1 = (1, 0), v 2 2 = ( 1 2 , 1 2 ). The vertices for ∧ 2 H 6 are then: v 3 1 . These form a r − 1 dimensional hyperplane in the r dimensional subspace for the two half sets, respectively. The practical effect of symmetry verification here (mostly from the S z application) is to project noisy points onto the plane. The effect of the N-representability application then is to map the measured points to the extreme points of V n ′ . This can be visualized in Fig. (III) by mapping the accessible triangular plane to the black outlined plane.
The unitary coupled cluster term as seen in Fig. 5, but in a simplified form. Note the target qubit for the rotation here is not q3 but q1. Importantly, the circuit require only 8 CNOT operations, and still has the same connectivity requirements. The gate H ′ is defined as S † HS (applied left to right).
For points close to the edges, one can imagine that with significant non-coherent error the projected points might lie outside of the space. We account for this by re-projecting these points into the polytope according to the closest edge. A simple semi-definite program would also suffice, though we take a more geometric approach. For mild errors we find that this method is satisfactory.

Appendix C: Second-Quantization Treatment of Entangling Gates and Phase
The exponential operator in (9) has a readily recognizable form in second quantization. Utilizing creation and annihilation operators in the natural orbital basis the operator: where t ii ′ jj ′ is a scalar, and anti symmetric with respect to swapping lower or upper indices, can be used to construct a two-body unitary operator U : Only one exponential Pauli term is needed to excite the initial double excitation. After that, a simplification was performed akin to that of Nam et al. [39]. For a system of two electrons, because do not consider number or spin changing operations, using the inverse Jordan-Wigner mapping we can simplify the total number of Pauli terms needed in the exponential to two.
To measure the phase of the terms in Eq. (8), we are interested in tomography of the 2-DM elements. Using the Jordan-Wigner transformation we can approximate the total operator as: Where we utilize the fact that any non-number conserving elements will be contribute 0. Due to the limited amount of sign terms we need to measure, and the fact every other excitation in the linear sequence will be completely commuting, we can prepare one circuit with only X i terms, and another with alternating X and Y pairs which give sign terms of either (C3) or (C5), yielding only 2 circuit preparations and a O(1) complexity. Note, obtaining all of the proper terms also scales as O(1), taking no more than 8 additional circuits. This method, however, did not yield significant increases in the accuracy in our computations, and so for our optimizations we used the approximate circuit. One issue with a direct measurement of the sign is that while the diagonal elements we measure can be symmetry verified, the N and S z operators do not necessarily commute with the Pauli terms (despite commuting with the operator as a whole). [N , X 1 X 2 Y 3 Y 4 ] = 0. To partially address this, we attempted to use in-line symmetry measurements, where we cast the Pauli measurement onto an ancilla qubit, allowing by propagating through the circuit to the corresponding entangling gate. This requires explicit connectivity requirements and attention to the layout, as a change in Pauli basis between the applied entangler and the required Pauli term. One issue we found was that while the sign information obtained in this manner was coherent and exhibited the proper behavior and change in sign with respect to t 22 ′ 11 ′ , there was a phase difference on the ancilla which resulted in the sign information being shifted from the magnitude of the occupations. A direct measurement of the Pauli terms, while not allowing for symmetry verification, still contained the correct qualitative information, and so was utilized for the 4-qubit case. For the 6-qubit case, obtaining reliable information throughout the longer optimization requirements was more difficult, and so to ease the demands on the quantum computer, we mapped the sign of the elements to the sign of the ideal function generated by our entanglers (a simple product of sine and cosine functions).