Robust determination of molecular spectra on a quantum processor

Harnessing the full power of nascent quantum processors requires the efficient management of a limited number of quantum bits with finite lifetime. Hybrid algorithms leveraging classical resources have demonstrated promising initial results in the efficient calculation of Hamiltonian ground states--an important eigenvalue problem in the physical sciences that is often classically intractable. In these protocols, a Hamiltonian is parsed and evaluated term-wise with a shallow quantum circuit, and the resulting energy minimized using classical resources. This reduces the number of consecutive logical operations that must be performed on the quantum hardware before the onset of decoherence. We demonstrate a complete implementation of the Variational Quantum Eigensolver (VQE), augmented with a novel Quantum Subspace Expansion, to calculate the complete energy spectrum of the H2 molecule with near chemical accuracy. The QSE also enables the mitigation of incoherent errors, potentially allowing the implementation of larger-scale algorithms without complex quantum error correction techniques.

Harnessing the full power of nascent quantum processors requires the efficient management of a limited number of quantum bits with finite lifetime. Hybrid algorithms leveraging classical resources have demonstrated promising initial results in the efficient calculation of Hamiltonian ground states -an important eigenvalue problem in the physical sciences that is often classically intractable. In these protocols, a Hamiltonian is parsed and evaluated term-wise with a shallow quantum circuit, and the resulting energy minimized using classical resources. This reduces the number of consecutive logical operations that must be performed on the quantum hardware before the onset of decoherence. We demonstrate a complete implementation of the Variational Quantum Eigensolver (VQE), augmented with a novel Quantum Subspace Expansion, to calculate the complete energy spectrum of the H2 molecule with near chemical accuracy. The QSE also enables the mitigation of incoherent errors, potentially allowing the implementation of larger-scale algorithms without complex quantum error correction techniques.
Quantum computing, the field of physics dedicated to harnessing quantum phenomena to process information, is rapidly progressing along the path from theoretical curiosity to practical technology. Recent experimental progress has been swift, with successful demonstrations of proof-of-concept algorithms on a range of fledgling quantum processors comprised of natural or engineered quantum spins [1-5]. However, even in leading architectures such as superconducting circuits and trapped ions, state-of-the-art systems comprise only few to tens of qubits-the quantum analog of classical bits-and are difficult to control with high precision. For gate-based quantum processors to be competitive with, or outperform their classical counterparts, both qubit number and gate fidelity must increase significantly [6,7]. Indeed, much of the field is currently focused on the design of a multi-qubit architecture capable of demonstrating an unambiguous quantum speedup over classical computers.
Recent theoretical advances suggest that a hybrid approach-judiciously dividing a computation between quantum and classical resources-will likely find utility in specific applications prior to the emergence of universal quantum computation [8][9][10][11]. One such example is calculating the ground-state energy of complex chemical systems, such as is often required in photovoltaics, biological reactions, and catalyst design. Based on the quantum variational principle-that the ground-state wavefunction of any Hamiltonian minimizes the expected energy [12]-an iterative protocol, known as the Variational Quantum Eigensolver (VQE) was developed. This approach uses a classical optimization routine to minimize the expected energy of candidate wavefunctions, using the quantum hardware to evaluate the expected energy. Essentially the VQE leverages the unique capac-ity of even shallow quantum circuits to prepare entangled states from which efficient classical sampling is not known to be possible.
Essential ingredients of the VQE algorithm have recently been demonstrated on a variety of experimental platforms [13][14][15][16][17][18][19]. These initial experiments indicate a robustness to systematic control errors (so-called coherent errors) which would preclude fully quantum calculations, as well as a manageable scaling of quantum circuit depth with Hamiltonian complexity [18,19]. However, work to date has focused primarily on calculating molecular ground-state energies; while extending the VQE approach to find excited states has been demonstrated in the optical domain, it required additional qubits, complicated multi-qubit control, and additional variational searches [18]. Furthermore, while a theoretical scheme for mitigating stochastic, incoherent errors has been proposed, this has yet to be verified experimentally [20]. In this work we demonstrate this extension of the VQE algorithm using a quantum processor comprising two superconducting transmon qubits, with real-time classical optimization augmented by a novel quantum subspace expansion (QSE) protocol. We diagonalize the electronic structure Hamiltonian of the hydrogen molecule over a wide range of nuclear separations and demonstrate the ability of this extended algorithm to calculate excitedstate energies and partially mitigate stochastic error channels, attaining near chemical accuracy (1.6 × 10 −3 H) in the calculated energy spectrum.

General Approach
The electronic structure Hamiltonian, an operator on the space of electronic wavefunctions, is first cast into a form suitable for evaluation on a quantum processor. Specifically, the Hamiltonian is first projected onto a discrete At k = 1 one has the linear response (LR) subspace while at k = N one spans the entire subspace corresponding to the particle number of the reference state, adapted from [20]. set of molecular orbitals-we use the conventional STO-3G basis set [21], which constitutes a so-called minimal set in that it represents the minimum number of orbitals required to represent a given atomic shell. The resulting fermionic Hamiltonian H F is finally mapped onto a two-qubit Hamiltonian H Q (SI Mapping of the H 2 Hamiltonian to qubits). For the hydrogen molecule, H Q takes the form where σ i k is the k th Pauli operator on the i th qubit, and the coefficients {g m }, and thus the total Hamiltonian, depend parametrically on R, the separation between the two hydrogen nuclei. For a given two-qubit state |ψ , prepared on the quantum processor, the expectation H Q is evaluated through repeated measurements of Pauli correlators.
An outline of the VQE algorithm is depicted in Fig. 1A and consists of parameterizing a quantum circuit U ( θ) to prepare an ansatz wavefunction |ψ( θ) , evaluating the expectation ψ( θ)|H Q |ψ( θ) term-wise using a quantum processor, and then using a classical minimization algorithm to update parameters until a minimum, θ min is found. The quantum state |ψ( θ min ) then constitutes an approximation to the ground state of H Q , with an estimated energy of ψ( θ min )|H Q |ψ( θ min ) .
Once the VQE algorithm has converged on an estimate of the ground state wavefunction, the quantum subspace expansion can be applied by measuring additional Pauli correlators that form an approximate matrix representation of H Q within an expanded subspace. This matrix can then be diagonalized classically to yield both lowlying excited-state energies and a refined ground-state energy estimate (Fig. 1C). If the expansion is chosen such that its dimension scales polynomially with system size, this classical matrix calculation is efficient [20]. The effectiveness of the QSE thus requires the existence of such a subspace which captures a significant amount of the excited state support.
We expect that molecular excited energy levels differ from the ground state primarily by excitations which promote a single electron from an occupied to an unoccupied orbital. Therefore to a good approximation, these states are linear combinations of {S 1 : a † i a j |ψ GS }, where a j (a † i ) are fermionic annihilation (creation) operators for the electronic orbitals. While S 1 could serve as a subspace, a more natural choice when working with qubits involves the set of single Pauli flips {P 1 : σ k α |ψ GS | α ∈ {x, y, z}, k ∈ {1, 2}} (Fig. 1D), which we refer to as a linear response expansion. To calculate the matrix elements H ij in the P 1 basis, we use the quantum processor to evaluate the inner products where |ψ GS is taken to be the initial approximate ground state |ψ( θ min ) , found via the VQE routine.
Beyond providing a means of calculating molecular excited state energies, it was conjectured in ref.
[20] that the inclusion of specific measurement operators expanding the number of states under consideration, the QSE could improve the accuracy of the initial VQE ground state estimate. While the VQE can in principle correct for the presence of coherent gate errors, the QSE was thought to additionally correct for incoherent errors, such as dephasing or amplitude damping. As discussed in the results section, we find experimental support for this conjecture.

Experimental Methods
Quantum. The quantum processor we use to evaluate expectation values consists of two superconducting qubits of the transmon variety [22, 23], one of the leading types of superconducting qubits in terms of design simplicity and coherence time. The qubits are initialized in the joint ground state |00 via a heralding measurement [24]. A generating circuit U ( θ) is then used to prepare the desired trial wavefunction (with θ specified by the classical hardware-see next section).
Single-qubit pulses on qubit A and B last 50 and 70 ns respectively, and achieve fidelities of ∼99%. The twoqubit pulse takes up to 310 ns and approaches a fidelity of ∼96%. After state preparation via U ( θ), tomographic reconstruction is used to evaluate H = ij h ij σ i σ j .
A near-quantum-limited traveling wave parametric amplifier [26, 27] enables high-fidelity single-shot measurement of the joint qubit state (SI Experimental Details). The entire sequence, including both state preparation and measurement, comprises less than ∼1.5 µs, below the coherence times of the qubits: 16 µs T 1A , 13.5 µs T * 2A , 12 µs T 1B , 3.5 µs T * 2B . Classical. With the two-qubit processor providing a means to efficiently evaluate H ( θ), the classical computer uses a particle-swarm optimizer (PSO) to find parameter values θ min which minimize this objective function, as shown in Fig. 2A. The PSO approach has two properties useful for this work: it is likely to avoid getting trapped in local minima and it is more robust to noisy objective-function calls [28]. The optimization treats a single point in parameter space as a particle, which has a velocity and position. A swarm of n such particles {|ψ( θ s,i ) , i ∈ [1, n]} (with s the swarm iteration number) is first randomly initialized and then at each iteration, the particles' positions are updated based on both their own energy evaluation and those of others in the swarm (SI Experimental Details). Figure 2B shows how iterating through this loop allows the particles to converge on a set of control parameters that prepares the FIG. 3: H 2 energy spectrum as a function of internuclear distance. Swarm particle energies for each bond length are histogrammed after application of a linear response expansion and Gaussian filter. Energy estimates obtained by a peak finding routine are indicated by dots with theoretically predicted energy levels shown as solid lines. An unphysical spurious state emerges at internuclear distances greater than ∼ 1.2Å due to uncorrected incoherent errors. Inset shows errors in the estimated ground and excited state energies as compared to chemical accuracy (1.6 × 10 −3 Ha).
best approximation of our system's ground state ψ( θ g ) and its associated energy.

Results
We apply our algorithm to the H 2 molecule for 45 internuclear distances between 0.05Å and 3.85Å. As shown in Fig. 2A for a internuclear distance of 1.55Å and a random initialization of 20 swarm particles over θ, we observe good convergence of the control parameters within 12 swarm iterations. Each function evaluation consists of 10,000 acquisitions and lasts on the order of a minute, resulting in a total algorithm run time of approximately four hours. Experimentally optimized parameters show deviation from those that would be expected in the case of idealized gates (SI VQE and Coherent Errors). In particular, while the experimental single-qubit gate amplitudes and two-qubit BSWAP length agree with numerical simulations, the phase of the BSWAP differs significantly, most likely due to an unaccounted for Stark shift during application of the gate. The successful convergence of the algorithm despite this mis-calibration thus provides additional proof of the protocols intrinsic ability to correct for coherent errors.
Plotting the median energy of the swarm as a function of iteration number, we observe a large initial energy error due to the random nature of the particle initialization, followed by an almost monotonic decrease towards the ex-act theoretical value. When calculating an estimate for a new internuclear distance, we exploit the smoothness of the parameter landscape and re-initialize the swarm particles around the minimum found in the preceding run, allowing them to vary by only 5% from their previous optimum values. This results in subsequent runs requiring fewer resources-20 particles and 6 swarm iterations-in order to reach convergence. Once each internuclear separation of interest has been processed, we have an initial approximation for the ground state energy function of the H 2 molecule.
To derive excited states from this approximate ground state, we apply the linear-response QSE to each individually-reconstructed density matrix recovered during the minimization process. The results of applying this expansion are plotted in Fig. 3 where data are binned with 1.5 mHa resolution before convolution with a Gaussian filter (standard deviation of 7.5 mHa). Peak-finding routines are then used to estimate the mean energies for both the corrected ground and excited states. This shows improved robustness for small numbers of swarm iterations as it is less affected by outlying particles in the swarm yet to reach the global minimum. level required to make realistic chemical predictions, is achieved for the ground and highest excited state across a wide range of internuclear distances. Estimates of the second and third excited state energies are generally within an order of magnitude of this level. It is interesting to note that although the ground electronic state wavefunction near equilibrium requires little entanglement to accurately represent, the same is not true of the excited states. The QSE is able to approximate these states with only additional local measurements and efficient classical computation, without an increase in required entanglement on the quantum state of the qubits. Figure 4 shows the deviations from the theoretically expected values for the corrected ground-state energies when using different underlying measurement operators for the applied QSE. Those involving just a single Pauli operator offer sporadic improvement over the uncorrected case, with the σ z operator achieving best results at smaller internuclear separations while the σ x operator is most useful at larger ones. The complete linear-response expansion is able to mitigate incoherent errors for which that the bare VQE algorithm is unable to compensate and produces a reduction in the energy estimate error of almost two orders of magnitude over the entire range.
Note that ideally, the total number of extracted energy levels should be upper-bounded by the dimension of the Hamiltonian. However, if the extant error channels cause the prepared VQE ground state to be sufficiently mixed (for a given set of QSE operators), it is possible to extract additional "spurious" energy levels. Such a spurious state is observed as indicated in Fig. 3 for internuclear distances between ∼ 1.2Å and ∼ 1.7Å. In some cases, these states may be discarded on the basis of continuity of the energy as a function of internuclear distance. Alternatively, these states can be removed by increasing the span of the QSE operators (at the cost of an increased tomographic measurement overhead). The exact conditions for the presence of a spurious state are currently being investigated.

Conclusion
We present a novel extension of the variational quantum eigensolver that only uses a polynomial number of additional tomographic measurements to extract molecular excited states and mitigate incoherent errors on the ground state estimate. With the hydrogen molecule as a test case, we additionally confirm the intrinsic ability of the algorithm to correct for coherent gate errors when pulse properties are optimized directly. Used with classical particle swarm minimization routines well suited to high-dimensional noisy environments, these techniques yield ground-and excited-state energy estimates with near-chemical accuracy. Our results highlight the potential of QSE to significantly reduce the need for more advanced error correction techniques, thereby facilitating practical applications of near-term quantum hardware.  Supplemental Materials: Robust determination of molecular spectra on a quantum processor

SI: Experimental Details
Here we provide details of the quantum processor used in the experiment. The device consists of two superconducting transmon qubits [22] on a single silicon chip, mounted in and coupled to a three-dimensional copper cavity [23]. Each transmon consists of a capacitor shunted by a nonlinear inductance; in our device one of the qubits uses a single Josephson junction as the nonlinear inductance, with a fixed frequency of 3.788 GHz, while the other uses a SQUID loop, allowing for tuning the frequency (via an external magnetic field) from a zero-flux value of roughly 5 GHz to the working frequency of 4.111 GHz. The copper cavity exhibits a resonant frequency of 7.122 GHz with a loaded linewidth κ ext ≈ 8 MHz, set primarily by the coupling (in a reflection geometry) to a 50-ohm environment (see Fig.  S1). The cavity is mounted at the 10 mK stage of a dilution refrigerator. We detect the state of the qubits by using a heterodyne measurement (at heterodyne frequency 11 MHz) of the resonant cavity frequency, exploiting the dispersive shift between the qubits and cavity. Because the two-qubits are coupled to a single cavity, the dispersive shift is roughly equal in magnitude for both qubits, and thus distinguishing the states |01 and |10 with single-shot fidelity is impossible. However, our measurement is able to distinguish the joint two-qubit ground state |00 from all other computational states, which is sufficient for reconstructing the Pauli correlators necessary for evaluating the expectation value H . To evaluate H we first reconstruct the twoqubit density matrix using a set of 32 tomographic measurements (see [29] for details), then calculate the necessary correlators given the density matrix. In future implementations of VQE on larger quantum systems, full tomography will be impossible (due to an exponential scaling of the number of required measurements). Instead, only the necessary correlators will be directly measured. For this reason, our reconstruction of the two-qubit density matrix from the tomographic measurements did not use any method such as maximum-likelihood estimation which enforces physicality (positivity and trace-normalization) on the result.
For the classical optimization routine, we used the Pyswarm package, a Python implementation of particle-swarm optimization.

SI: Mapping of the H 2 Hamiltonian to Qubits
The H 2 Hamiltonian was determined by calculation of the electronic integrals in the standard Gaussian STO-3G basis [21]. For a hydrogen atom, this basis set consists of the single 1s orbital. After determination of the one-and two-electron integrals, the Hamiltonian was projected into a particle conserving manifold defined by a determinantal configuration interaction expansion that does not flip spins. That is, within a molecular orbital basis we define spatial orbitals 1 and 2 with possible spins α and β such that the 4 spin-orbitals in the system can be populated by the second quantized operators a † 1α , a † 1β , a † 2α , a † 2β . Beginning from the reference state defined by a † 1α a † 1β |vac , where |vac is the Fermi vacuum state. This generates the following four basis states that we map to computational basis states explicitly as follows We note that such a reduction to 2-qubits or fewer can also be achieved either through the Bravyi-Kitaev transfor- for each nuclear configuration R, where σ i α is a Pauli operator acting on qubit i from σ i α ∈ {I i , σ i x , σ i y , σ i z }. Due to additional spatial, spin, and time-reversal symmetry in the molecular Hamiltonian, many of the coefficients are 0 for all nuclear configurations R and all are real-valued. As a result, the Hamiltonian may be more compactly expressed as H Q (R) = g 0 (R) + g 1 (R)σ 1 z + g 2 (R)σ 2 z + g 3 (R)σ 1 z σ 2 z + g 4 (R)σ 1 y σ 2 y + g 5 (R)σ 1 FIG. S1: Schematic of the measurement setup used in the experiment. At the 10 mK stage of the dilution refrigerator, the sample is connected in a reflection geometry to the 50 ohm environment from which it receives qubit/cavity pulses. These pulses are generated at room temperature by the electronics shown on the right of the figure, and pass through several stages of attenuation on the way to the sample. To enable high-fidelity measurement of the qubit state, a near-quantum-limited Traveling Wave Parametric Amplifier (TWPA) [26] amplifies the signal after reflection off the cavity. Further amplification is provided by a HEMT amplifier at the 4 K state of the dilution refrigerator, after which the signal is amplified at room temperature, downconverted to a heterodyne frequency of 11 MHz, and digitized by an AlazarTech 9373 ADC. From this data, the qubit state is determined in software.
The exact coefficients used in this work are given in Table I of this SI.

SI: QSE and Choice of Expansion Operators
The choice of operators which act on the ground-state density matrix to form the expanded subspace influences which excited can be extracted. This we show in figure S2. Using only the identity and single Pauli operators (on each qubit) results in only a partial resolution of the low-lying excited states, while the full linear-response is able to resolve the entire spectrum.

SI: QSE with Errors
The quantum subspace expansion (QSE), works by resolving the action of an operator H within a subspace defined by a set of operator {O i }, such as the single fermion excitation operators S 1 or k th order Pauli operators P k defined in the text. This is done by measuring the matrix elements coupling the states generated by these operators H ij = ψ 0 |O † i HO j |ψ 0 as well as the corresponding identity operator in this space, also known as the overlap matrix, The action within this subspace is then used to provide increasingly accurate approximations (as a function of the subspace size) by solving the generalized eigenvalue problem HC = SCE with the matrices H and S for the eigenvectors C and diagonal matrix of eigenvalues E.
Defining the density matrix for the pure state |ψ 0 as ρ 0 = |ψ 0 ψ 0 |, it is easy to see these matrix elements are equivalent are equivalent to . This formulation naturally generalizes to mixed states ρ with rank > 1, which are the case in essentially any real system with incoherent errors, and gives a clear prescription for the measurement of the matrix elements. However, in this notation it is less clear how the measurement correspond to action within a subspace and what this means in the case for mixed states ρ with rank > 1. To clarify these situations, we may alternatively use the vectorization of the density matrix to re-express these matrix elements.
We denote the row-major vectorization [31] of a matrix ρ as |ρ . In this notation, we have that This construction clarifies a number of the mathematical properties, including the hermiticity of the matrices and their dimensionality. In the case of a pure state, the operator ρ T has a single non-zero eigenvalue, and the maximum non-trivial dimension of the space, which is determined by the trace of the identity operator S using normalized operators {O i }, is that of the Hamiltonian. It is important to consider in more detail when the rank of ρ is > 1. In these cases, the dimension of the space is potentially greater than the dimension of the original Hamiltonian. The easiest way to see this is to consider the case of the maximally mixed state ρ = 1 d I, where d is the dimension. In this case, the dimension determined by the identity is the square of that of the original Hamiltonian, which this construction makes clear is the maximal dimension of this problem. Moreover, by properties of the standard tensor product, it is easy to verify that the eigenvalues of H ⊗ I are the eigenvalues of H, but d−fold degenerate. We note that the factor of 1/d is treated by its appearance in the metric matrix S in the generalized eigenvalue problem.
Thus, if one measures a linearly independent, complete set of operators {O i } on the totally mixed state, the resulting eigenvalues will be the spectrum of H with d−fold degeneracies. These additional states represent the different possible expansions from components |ψ i ψ i | with ρ = i |ψ i ψ i | that allowed one to prepare the eigenstates using {O i }.
Consider, however, if the resolution of the operator ρ is incomplete with the given measurements {O i }, then one may have difficulty distinguishing the eigenstate approximations generated from different pure states. This leads to the so-called "spurious states" observed experimentally in this work, which are extra predicted eigenvalues that do not coincide with the eigenvalues of H.
To clarify these effects, consider the following Hamiltonian which has a ground state given by ρ 0 = α|00 + β|11 . In the case of ρ 0 , it is clear that the maximum dimension of this problem is 4 with any set of measurement. Now considering the mixed state generated by a Pauli-X channel that occurs with probability p = 0, 1, ρ = (1 − p)ρ 0 + pσ 1 x ρ 0 σ 1 x . In this case, the operator ρ has as non-trivial eigenvalues p and (1 − p). As an example we choose p = 1 2 such that it has a degenerate non-trivial spectrum of 1 2 and 1 2 . In the case one measures a complete set of operators {O i }, one then finds the eigenvalues of H with a degeneracy of 2 in each case. If we consider only the error in the estimate of the ground state energy, one finds that the operator set {I, σ 1 z σ 2 z } is sufficient to correct it exactly (if applied to the state resulting via acting with the error channel on the ideal ground state, i.e. without minimization on this value). The exact condition for the set of operators that correct errors in the ground state for a given H and error channel and their relation to traditional theories of quantum error correction is an open problem, currently the subject of ongoing research. We conjecture here based on numerical observations that conditions are related to the ability to construct operators within Span({O i }) that both commute with H but not with the error channel E.
To study the case of spurious states, suppose one measures a set of operators with dimension greater than the dimension of the Hamiltonian but not sufficient to resolve ρ.
In this case, one such set is If one measures the Hamiltonian and overlap on state ρ with this basis, one sees examples of the observed behaviors.
First, the non-trivial dimension of the problem is 7, which would be an experimental signature that the measured state is not a pure state but also not the totally mixed state. Second, the eigenspectrum contains the exact spectrum of H, but is not degenerate. Rather it contains 3 erroneous eigenvalues that correspond to the spurious states we define above. Thus the total spectrum is formed from a combination of an exact expansion from one state and a poor expansion from another. If one continues to add operators, the spurious values disappear, replaced by degeneracies in the spectrum on the exact values. It is interesting to note that if one chooses operators capable of correcting these errors, a smaller set such as {I, σ 1 x , σ 1 z , σ 2 z , σ 1 x σ 2 z , σ 1 z σ 2 z } produces the exact spectrum with degeneracies only on the 2nd and 3rd eigenvalues with no spurious states.

SI: VQE and Coherent Errors
The VQE is expected to have an intrinsic ability to correct for coherent gate errors (such as under or over rotations) due to the direct parameterization of the microwave pulse amplitudes/lengths and phases. As a signature of this ability, we plot in Fig. S3 the optimal parameters found by the VQE algorithm (for an internuclear distance of 1.55Å) and compare them to the parameters expected from our initial simulations. The amplitudes of the singlequbit rotations converge to nearly zero, as expected from simulations. The length of the bSWAP gate, which is finite so as to create entanglement between the two qubits, also agrees with simulation. However, the phase of the bSWAP drive differs significantly from the expected value, namely zero. To say this another way, at this internuclear separation, the theoretically-expected ground state wavefunction is a superposition of the states |00 and |11 with equal phases, yet the experimentally prepared ground state clearly exhibits a phase difference in the amplitudes of the |00 and |11 states. This discrepancy is likely due to an uncalibrated Stark shift during the bSWAP gate which is corrected for automatically by the classical minimization routine, which has no knowledge of the true bSWAP unitary transformation. Phases of the single qubit drives are not included on the figure, as the single qubit amplitudes have converged to zero, which renders the phase meaningless.

SI: QSE beyond Linear Response
In this two-qubit example of the QSE, it is straightforward to go beyond the linear-response subspace and include additional measurement operators (see Fig. S4), as a demonstration that further error mitigation is possible. The  FIG. S4: Starting from a large initial internuclear distance, the purely linear response expansion (blue) achieves chemical accuracy in the calculation of the ground state energy until a restart of the data collection run due to technical reasons (vertical black line). From this point onwards the accuracy in linear response estimate is degraded, most likely due to calibration errors in our tomographic reconstructions. Overcoming this calibration drift can be achieved by including additional two-qubit correlators such as σ x σ x in the measurement span (gray). dataset shown in the figure was taken in two separate runs: the first, for internuclear separations greater than 2.6 Angstroms, and the second for separations lower than that value. Technical reasons necessitated a restart of the data collection at that value. In this dataset, the bare VQE ground state error is more than an order of magnitude above the threshold for chemical accuracy. In the initial data run (for internuclear separations greater than 2.6 Angstroms), the linear-response correction is able to bring this error down below the chemical-accuracy threshold. However, after the restart, the linear-response is no longer able to get below chemical accuracy. The likely reason for this is a drift in the gates which effected the tomographic reconstructions. But even though the linear response fails to fall below chemical accuracy, such accuracy can still be achieved over a large range of the separations by including additional operators in the QSE. In the figure, specifically, we show how the addition of the operator σ 1 x σ 2 x dramatically improves the accuracy of the ground-state estimate. as a function of the internuclear distance R. The electronic integrals were calculated in the STO-3G basis and columns with entries that were 0 by symmetry have been omitted.