Eﬃcient and Accurate Electronic Structure Simulation Demonstrated on a Trapped-Ion Quantum Computer

Quantum computers have the potential to perform accurate and eﬃcient electronic structure calculations, enabling the simulation of properties of materials. However, today’s noisy, intermediate-scale quantum (NISQ) devices have a limited number of qubits and gate operations due to the presence of errors. Here, we propose a systematically improvable end-to-end pipeline to alleviate these limitations. Our proposed resource-eﬃcient pipeline combines problem decomposition techniques for compact molecular representations, circuit optimization methods for compilation, solving the eigenvalue problem on advanced quantum hardware, and error-mitigation techniques in post-processing the results. Using the density matrix embedding theory for compact representation, and an ion-trap quantum computer, we simulate a ring of 10 hydrogen atoms taking into account all electrons equally and explicitly in the electronic structure calculation. In our experiment, we simulated the largest molecular system on a quantum computer within chemical accuracy with respect to total molecular energy calculated by the full CI method. Our methods reduce the number of qubits required for high-accuracy quantum simulations by an order of magnitude in the present work, enabling the simulation of larger, more industrially relevant molecules using NISQ devices. They are further systematically improvable as devices’ computational capacity continues to grow.

(Date: February 23, 2021) Quantum computers have the potential to perform accurate and efficient electronic structure calculations, enabling the simulation of properties of materials. However, today's noisy, intermediatescale quantum (NISQ) devices have a limited number of qubits and gate operations due to the presence of errors. Here, we propose a systematically improvable end-to-end pipeline to alleviate these limitations. Our proposed resource-efficient pipeline combines problem decomposition techniques for compact molecular representations, circuit optimization methods for compilation, solving the eigenvalue problem on advanced quantum hardware, and error-mitigation techniques in postprocessing the results. Using the density matrix embedding theory for compact representation, and an ion-trap quantum computer, we simulate a ring of 10 hydrogen atoms taking into account all electrons equally and explicitly in the electronic structure calculation. In our experiment, we simulated the largest molecular system on a quantum computer within chemical accuracy with respect to total molecular energy calculated by the full CI method. Our methods reduce the number of qubits required for high-accuracy quantum simulations by an order of magnitude in the present work, enabling the simulation of larger, more industrially relevant molecules using NISQ devices. They are further systematically improvable as devices' computational capacity continues to grow.
Electronic structure simulation is an essential tool for understanding chemical properties of molecules. It is a basis for contemporary materials design and drug discovery. Performing accurate electronic structure simulations on classical computers requires a great amount of computational resources. In particular, they grow exponentially with the system size when employing the full configuration interaction (full CI) method, which calculates the exact solution of the electronic Schrödinger equation in a given basis set. Quantum computing, a computing paradigm that leverages the laws of quantum physics, promises to deliver scalable and accurate electronic structure calculations [1,2] beyond the reach of classical computers.
Quantum computing technologies are rapidly advancing, and there has been major progress in simulating molecular systems in the last two decades [3][4][5][6][7][8][9][10][11][12][13]. However, simulating the electronic structure of industrially relevant molecular systems on today's noisy, intermediate-scale quantum (NISQ) devices [14] will require systematically scalable, robust methods that allow a given problem to be represented by a small number of qubits and shallow quantum circuits. Without these methods, the stringent constraints of NISQ devices will inhibit their ability to perform high-accuracy simulations of larger molecular systems, particularly those systems where the electron correlation is strong.
The treatment of electron correlation is necessary for applying electronic structure calculations to make accu-rate predictions about the process of a chemical reaction. Calculating the molecular energy remains a challenge, even for small systems, without limiting the number of configurations or the number of electrons active in performing the electronic structure calculations on NISQ devices.
The attempt to perform the largest calculation of the total energy including electron correlation to date is the simulation of BeH 2 using six qubits for the six-electron problem [7]. Note further that some of the authors of the present manuscript simulate a water molecule using a trapped-ion quantum computer. Using a small number of electron configurations that are known to contribute significantly to the total energy of the molecule, the energy estimates obtained were within the widely used measure for chemical accuracy (1.5936 × 10 −3 hartrees) when compared to classically simulated results [9]. Other recent research [13] reports the simulation of a 12-electron problem of a chain of 12 hydrogen atoms using 12 qubits, setting the record in terms of the largest number of qubits used for a chemistry simulation. However, in this last study, the energy of the molecular systems considered is calculated using the Hartree-Fock method, which does not take electron correlation into account.
One promising direction for alleviating the limitations encountered by researchers in earlier studies is to leverage problem decomposition (PD) techniques. These techniques admit a more compact representation of a molecule, and, therefore, enable the explicit inclusion of more electrons in calculating correlation energies. Problem decomposition techniques decompose a given molecular system into small subsystems, without sacrificing the accuracy of the electronic structure simulation for a wide class of chemical systems [15][16][17][18]. They have the potential to substantially reduce the qubit count requirements for performing electronic structure simulations [19][20][21][22][23][24][25][26][27].
In the present work, we develop a PD pipeline based on the density matrix embedding theory (DMET) [28,29], which can be applied to a wide range of molecules. A quantum algorithm can be used to calculate the electronic structure of each fragment resulting from the decomposition. Although our pipeline is agnostic with respect to the choice of quantum algorithm, we choose the variational quantum eigensolver (VQE) [4] to achieve shallower quantum circuits amenable to NISQ devices. Figure 1(A) shows an illustration of the pipeline. Note that DMET has been applied to quantum simulations of the Hubbard model [19] on quantum computers, but not for molecular systems. We include a trapped-ion quantum computer in our pipeline, as it is a platform that has performed quantum simulation experiments well [6,8,9] and allows the efficient implementation of quantum algorithms via complete qubit connectivity [30][31][32]. We further employ an efficient density matrix purification algorithm to post-process the experimentally determined results and mitigate any residual error, by leveraging the high-quality, two-qubit gates of the trapped-ion hardware.
We use a ring of ten hydrogen atoms to demonstrate the viability of our pipeline. Taking electron correlation into account, we demonstrate the viability considering all electrons equally and explicitly in the electronic structure calculation. We can thus further widen the applicability of electron structure calculations on NISQ devices for larger, industrially relevant molecular systems.
We begin by describing our DMET-based methodology for simulating the electronic structure of a molecular system (see Appendix A for further details). Consider a system described by a second-quantized Hamiltonian H = H 1-e + H 2-e , where H 1-e is the one-electron interaction and H 2-e is the two-electron interaction of the entire molecule. In contrast to simulating the entire system usingĤ, in DMET, the system to be simulated is divided into small fragments. Each fragment is treated as an open quantum system, entangled with its surrounding environment, or bath. Here we use the mean-field, Hartree-Fock (HF) solution of the entire molecular system as a pre-processing step to find local orbitals that we use to fragment the molecule. Then, the following iterative process, which we call the DMET cycle, initiates. Once the cycle terminates, the electronic structure calculation of the entire molecule is complete. The DMET algorithm is shown in Fig. 1(B).
The DMET cycle begins by constructing the bath orbitals for each fragment. The bath orbitals describe the environment that is active for the electronic structure calculation of the fragment, by virtue of Schmidt decom-position [33]. Note that, if the bath is large, the description of the bath can be greatly simplified. With the simple description of the bath, the Hamiltonian H A for each fragment A, along with its specific bath B, is constructed according to the equation where H 1-e,AB denotes the one-electron interaction within and across the fragment and the bath, H 2-e,A denotes the two-electron interaction within the fragment, µ is the chemical potential, and N A is the number of electrons in the fragment. See Appendix B for details. As will be discussed shortly, we use a quantum computer to our advantage in the DMET calculation, efficiently evaluating the minimal expectation value of H A , as well as the number of electrons N A in fragment A.
Once both the energy expectation value and the number of electrons for each fragment are computed, we combine them to compute the total system energy, and check for self-consistency. In particular, we choose to compute the sum of the number of electrons in each fragment and check whether the sum is equal to the total number of electrons in the system. If the sum is within a pre-specified range with respect to the total number of electrons, the DMET cycle terminates. If the sum is not within the specified range, we run the DMET cycle again, with the chemical potential µ updated as the difference between the sum and the total number of electrons.
As an explicit example of calculating the electronic structure of a molecule using DMET, we simulate a ring of 10 hydrogen atoms, H 10 . The molecule is divided into 10 fragments, with one atom per fragment. We allocate two spin-orbitals each for the fragment and the bath. This may be compared to a total of 20 spin-orbitals in the simulation of an entire molecule, which shows a large reduction in the problem size: a 20-qubit problem is reduced to a two-qubit problem. Based on the symmetry exhibited by the system, we simulate only a single twoqubit fragment and use the results of the simulation to infer that the results are identical for every fragment in the DMET cycle described above.
To estimate the expectation value of the fragment energy and the number of particles per fragment, we use VQE [4] with the qubit coupled-cluster (QCC) ansatz [34]. For all calculations we perform, we use symmetry-conserving Bravyi-Kitaev transformation [35,36] to transform from a fermion to a qubit basis. The QCC ansatz operatorÛ (τ ) is specified according to the equationÛ where τ k is a variational parameter, n g is the number of   multi-qubit Pauli operatorsP k , defined aŝ where n q is the number of qubits and X, Y, Z, and I are the Pauli matrices and a single-qubit identity operator, respectively. While the depth of the circuit rapidly increases as the size of the molecule increases, the QCC ansatz admits a low-depth quantum circuit compared to a widely used unitary coupled-cluster single and double ansatz. The details of the QCC ansatz can be found in Appendix C.
We first consider an ansatz state that is a product state of n q arbitrary single-qubit states, following the method described in the work of Ryabinkin et al. [37]. We evaluate the expectation value of the mean-field Hamiltonian with respect to the ansatz state and use VQE, simulated on a classical computer, to minimize the value. The variational parameters that result from the optimization correspond to the optimal wavefunction for the mean field. We next consider a QCC ansatz operatorÛ (τ ) applied to the previously determined, mean-field optimized state. Note that we aim to minimize the expectation value of the fragment Hamiltonian with respect to the ansatz operator parameters τ . We find whichP k to include in our ansatz by computing the derivative of the expectation value with respect to τ k , which can be computed in a straightforward way when τ = 0. We removeP k terms that have small derivative values. Applied to our example molecule H 10 , we find thatP k of XY and Y X have large, identical derivative values. Note that we carry out the computation of the derivatives on a classical computer. We now investigate the performance of DMET with the QCC ansatz, which we denote as DMET-QCC, on a classical computer. Specifically, we consider the potential energy curve of the symmetric expansion (i.e., increasing the bond length R while maintaining it for each pair of neighbouring atoms) of H 10 . We compare the total energies resulting from the two-qubit DMET-QCC ansatz with those obtained from other known methods, such as HF and full CI, calculated classically, to assess the performance of DMET-QCC. We choose five points R = 0.7, 1.0, 1.1, 1.3, and 1.6 A along the potential energy curve for comparison. The total energies per atom of the H 10 molecule are listed in Table I. All DMET-QCC results nearly coincide with the full CI results.
We next describe the simulation of our DMET method on a trapped-ion quantum computer. Instead of variationally optimizing the energies, then running through the DMET cycles, here we focus on the evaluation of the total energy from the quantum simulation for the classically pre-computed optimal parameters. Note that for the DMET cycles, we require XX , Y Y , ZZ , XZ , ZX , Z 0 , Z 1 , X 0 , and X 1 to be simulated, where the subscripts denote the qubit index. See Appendix D for details on the Hamiltonian and the DMET energy expressions used in the experiments.
We implement the quantum simulation circuits on IonQ's 11-qubit trapped-ion quantum computer, which is described in detail elsewhere [38]. In the quantum computer, 15 171 Yb + ions, aligned to form a linear crystal with spacing of about 5 µm, are suspended in a chip trap with a radial pseudopotential frequency of ∼3.1 MHz. We cool the crystal to its motional ground state and use the 11 ions in the middle as qubits. Counter-propagating laser beams capable of illuminating individual ions are used to implement quantum gates, leveraging the ionion coupling mediated by the collective radial motional modes. We read out the quantum state by fluorescing the ions using a detection laser.
We extract the expectation values of the Pauli terms using the three compiled, optimized circuits shown in Fig. 1(C)(2)(a). We use the XX and ZZ circuits for XX and the XZ circuit for XZ and ZX , as the two preoptimization circuits for XZ and ZX reduce to the same circuit upon optimizing compilation. The singlequbit Pauli terms X 0 , X 1 , Z 0 , and Z 1 , where B i denotes the i-th qubit's expectation value in the B ba-sis, are computed using all of the statistics available from the three circuit executions. For instance, for X 0 , we use the results from both XX and XZ . Likewise approaches are used for X 1 , Z 0 , and Z 1 .
The total energies obtained from the quantum computer (blue plotted points), along with the classically computed potential energies (lines), using the reference R values are shown in Fig. 2. The experimentally determined energies nearly coincide with the full CI energies (black line). The full CI values lie within the error bars for every point. With the exception of R = 1.1 A, the experimental energy values are in agreement with the ideal DMET-QCC results to within chemical accuracy (1.5936 × 10 −3 hartrees). See Table I for the precise values of the experimental energies and their error values. The error bars are evaluated from the empirical bootstrapping method detailed in Appendix G. This shows that the DMET-QCC ansatz, employed in a trappedion quantum computer, is capable of describing bondbreaking and bond-forming processes in molecular systems. We now perform classical post-processing, purifying the two-particle, reduced density matrix using McWeeney's density purification technique [39] to see whether we can further improve the total energy evaluation. The two-particle, reduced density matrices obtained from the experiment are iteratively purified as P new pqrs = 3(P old pqrs ) 2 − 2(P old pqrs ) 3 until the convergence criterion of T r(P 2 pqrs − P pqrs ) < 1.0 × 10 −2 is met. See Appendix E for details. This requires us to run an additional circuit, YY, on the trapped-ion quantum computer with three single-qubit and one two-qubit native operations, as shown in Fig. 1(C)(2)(b). See Appendix H for details on the optimizing compilation. In order to capture the nature of the entanglement in our simulation, TABLE I. The calculated total energies per atom (in hartrees), using the HF, full CI, and DMET-QCC methods, for the H10 molecule along the potential energy curve R = 0.7, 1.0, 1.1, 1.3, and 1.6 A of the symmetric expansion. For DMET-QCC, we report the theoretical (T), experimental (E), and post-processed (P) results. For the experimental and the post-processed results, we report the standard deviation of the energies, calculated using bootstrapping. See Appendix G for details. Note that we intentionally report four digits in the experimentally determined and post-processed energies for the purpose of simple comparison.
which is the characteristic phenomenon of quantum simulation, we also investigate two types of quantum entanglement. The Hilbert space of the two-qubit system is split into two subsystems: the first and second qubits, and the fragment and the bath orbitals. In both cases, the experimental results are very close to theoretical predictions and clearly show that the states we observe in the present experiment are entangled. See Appendix F for more details. Figure 2 shows the post-processed total energies (red plotted points). The post-processed energies are within chemical accuracy with the full CI energies for all points calculated. See Table I for the precise values of both the post-processed energies and their respective errors. The high quality of the two-qubit gates in the hardware enables accurate post-processing. See Appendix E for further details. For R = 1.1 A, the post-processed total energy is E = −0.539 ± 0.007 hartrees, which is within chemical accuracy. The two-particle, reduced density matrices with large error values are greatly improved by the density matrix purification, as also seen in a previous work [10]. Note that not all total energies for all points in our work are improved (R = 0.7, 1.0 A) compared to the ideal DMET-QCC results, as the errors in some of the determined energies prior to post-processing are very small.
Using the post-processed result where R = 1.1 A, we achieve chemical accuracy for the total energies with ideal DMET-QCC results for all points. In addition, the postprocessed energies are within chemical accuracy with respect to the full CI energies for all points. This is realized by treating electron correlation for a 10-electron system without using frozen-core approximation (i.e., taking into account all electrons equally and explicitly) in our simulation on a trapped-ion quantum computer.
In this work, we have described a systematic end-toend pipeline developed using a PD technique to reduce the size of an electronic structure simulation to realize the simulation on a trapped-ion NISQ device. Our methodology involves the following. We combine DMET and VQE with a QCC ansatz to compress the quantum simulation circuit for the electronic structure calculation and apply circuit optimization techniques that target trapped-ion quantum computers. We further apply density matrix purification algorithms to mitigate residual errors. We construct the potential energy curve of the ring of 10 hydrogen atoms as a proof of concept to investigate the viability of our pipeline. Our results achieve chemical accuracy with respect to both the ideal simulated results and the full CI results. Our pipeline allows us to reduce the number of qubits required for our simulation, hence enabling us to equally and explicitly consider all electrons in the computation of electronic correlation. Further, our results demonstrate that the approach employed herein can describe bond-breaking and bond-forming in molecular systems. To our knowledge, the current molecular system is the largest chemistry problem simulated on a quantum hardware device to obtain the total molecular energies, including electron correlation.
Computer-aided molecular design in both materials science and the life sciences is key to our attaining a sustainable future through the accurate prediction of chemical reactions, such as the catalytic reaction of organometallic compounds in advanced materials innovation, enzymatic reactions in the life sciences, and the electrochemical reactions implemented in next-generation batteries. In a previous work by some authors of the present paper [26], it is shown that a quantum simulation of a typical, industrially relevant organometallic compound would require over 2000 qubits. Leveraging the pipeline presented herein and selecting the appropriate PD technique could result in the reduction of qubit requirements significantly. A reduction by a factor of five for industrially relevant systems is shown in another work [22]. This sort of reduction could place these simulations within the reach of near-future NISQ devices that have hundreds of qubits. The precise reduction in the amount of required resources can be both algorithm dependent and application dependent, but our current work shows a potential for resource reduction by up to a factor of 10, while preserving the accuracy of simulations. Our pipeline has the capacity to leverage the capabilities of quantum hardware, paving the road toward a molecular design platform enabled by a quantum computer.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The code is available from the corresponding author upon reasonable request.

A. THE SYSTEM HAMILTONIAN
The second-quantized electronic Hamiltonian can be written aŝ where p, q, r, and s are distributed over all spin-orbitals, and a † p and a q are the corresponding creation and annihilation operators. The terms h pq and pq|rs are the one-and two-electron integrals (in chemists' notation), respectively. The evaluation of these integrals, as well as the Hartree-Fock (HF) and full CI calculations, are carried out using PySCF [40]. The minimal basis set MINAO [41] is used in our calculation.
The electronic Hamiltonian is then transformed into the qubitized form by using symmetry-conserving Bravyi-Kitaev transformation (scBK) [35,36]. The qubitized Hamiltonian can be written aŝ where p, q, r,. . . are distributed over all qubits, and σ α p ∈ {σ x , σ y , σ z , I} acts on qubit p. The transformation of the electronic Hamiltonian into a qubit basis is performed using OpenFermion [42].

B. THE DENSITY MATRIX EMBEDDING THEORY
Let us consider a molecular system divided into fragments, in which a small fragment A with N A states is surrounded by a large bath B with N B states. If the wavefunction of the entire system |Ψ is known, Schmidt decomposition [33] can be applied, and we may write where |α i and |β i denote particular many-body bases, and |χ i is the orthogonal set of |χ i . The states |χ i are the states of bath B; however, the number of states coincides with that of fragment A. This shows that, regardless of the size of bath B, only N A states can be entangled with the fragment. This will reduce the size of the problem drastically for large-sized systems. The exact wavefunction of the entire system |Ψ is the eigenfunction of the Hamiltonian for the entire system H. The Hamiltonian H of fragment A embedded in bath B can be defined using a projection operator P , that is, where It is now evident that the electronic structure of the entire system can be described exactly by that of the fragments and their surrounding baths. The electronic structure calculation of the entire system can be solved using this smaller problem. However, the exact wavefunction of the molecular system is usually not known a priori; thus, introduction of an approximation is necessary. The wavefunction of the entire system obtained by low-level, mean-field theory, such as the HF calculation, would be a straightforward approximation of the exact wavefunction of the entire system. Using a low-level wavefunction to construct a bath and solve the reduced problem employing a high-level theory is the principal idea behind the density matrix embedding theory (DMET).
The orbitals of the entire system are transformed into unentangled occupied orbitals, unentangled virtual orbitals, local fragment orbitals, and bath orbitals. The orbital space of the entire system is then greatly reduced, as only the local fragment and bath orbitals are employed for each high-level DMET fragment calculation.
Practically, we need to optimize the embedding of a bath. In DMET, a high-level calculation for each frag-ment is carried out individually until self-consistency has been attained according to a certain criterion: the sum of the one-particle reduced density matrix (1-RDM) of all of the fragments agrees with that of the low-level one for the entire system. The DMET energy is calculated using the 1-RDM and two-particle reduced density matrix (2-RDM). The DMET algorithm used in this work, the single-shot algorithm [29], can be described as follows: 1. Calculate the wavefunction of the entire molecular system using a low-level method and then localize the orbitals to fragment the molecule.
2. Construct the bath orbitals so as to include the surrounding environmental effect.
3. Construct the Hamiltonian of a fragment (including environmental effect) and calculate the wavefunction using a method based on a high-level theory.
4. Calculate the fragment energy and the number of electrons for each fragment with the 1-and 2-RDMs from the wavefunction obtained in Step 3.

5.
Repeat steps 2-4 for each fragment and obtain the total energy and the high-level 1-RDM of the entire molecular system. 6. Repeat steps 2-5 until the sum of the number of electrons in the fragments agrees with the number of electrons for the entire system.
The initial step of a DMET calculation is to perform a mean-field HF calculation for the entire molecule. The localized orbitals are obtained by localizing the canonical orbitals from the HF calculation to determine how to fragment the molecules. The Meta-Löwdin localization scheme [43] is used in this work.
After the molecule is divided, the DMET cycle (steps 2-6) is initiated. The first step in the cycle is to obtain the bath orbitalsâ † r (the active orbitals in the electronic structure calculation used for fragments to describe environmental effects), and the environment density matrix D env,A is calculated as follows, where C represents the molecular orbital coefficients obtained from the mean-field calculation of the entire molecule. Now, the embedding Hamiltonian H emb,A can be constructed (step 3). It can be defined as where the h pq are the one-electron integrals, the pq|rs are the two-electron integrals in chemists' notation, L A is the number of orbitals in the fragment, L B is the number of bath orbitals, L is the number of orbitals in the entire molecule, and p, q, r, and s are general orbital indices. We introduce the chemical potential δµ, which is optimized in the DMET cycle. Once the Hamiltonian H emb,A has been obtained, the electronic structure calculation is performed for fragment A. We employ VQE with the QCC ansatz in this work. Following these calculations, the algorithm constructs the one-and two-particle density matrices, D A and P A , respectively, from the QCC wavefunction Ψ A . The fragment energy E A is calculated (step 4) from the reduced density matrices as Note that only the elements with fragment orbital indices (p ∈ A) are used for calculating the fragment energy.
Here, the 1-and 2-RDMs are defined as and respectively. The number of electrons N A in fragment A are calculated (step 4) as The DMET energy is calculated by summing the fragment energy for each fragment (step 5), which is obtained according to the equation where E nuc is the nuclear repulsion energy. The DMET cycle (steps 2-6) iterates until the number of electrons in the DMET calculation given by converges to the total number of electrons in the molecule N Total . Convergence is achieved by updating the chemical potential δµ according to the equation where a is positive number. For the DMET calculation for the H 10 ring, the high symmetry of the molecular structure allows us to reduce the calculation further. We calculate only a single fragment with the assumption that all of the fragments have both the same energy and number of electrons. The DMET total energy with only one fragment multiplied by 10 (the number of hydrogen atoms) coincides with the energy using all fragments; thus, we treated only one fragment in the DMET calculation and simply multiplied the energy and the number of electrons by 10.

C. THE QUBIT COUPLED-CLUSTER METHOD
An accurate and affordable description of the correlated wave function |Ψ required to evaluate the energy of fragment A according to Eq. (B7) is achieved using the qubit coupled-cluster (QCC) method [34]. Within the QCC approach, a mean-field wave function |Ω is determined and subsequently utilized in a heuristic [37] to construct a unitary operator ansatzÛ (τ ). This operator recovers the missing electron correlation for the mean-field state |Ω and results in the QCC wave function according to the equation A parameterized mean-field wave function is defined as a direct product of n q single-qubit states and can be expressed as is the set of 2n q mean-field parameters. Each single-qubit state is then represented in the qubit computational basis as and is characterized by the Bloch sphere polar and azimuthal angles θ j and φ j , respectively [34]. The meanfield energy functional is given by Minimization of Eq. (C4) with respect to Γ gives the ground-state, mean-field energy E MF and the corresponding optimal set of parameters Γ opt . The QCC unitary operator ansatzÛ (τ ) takes the form whereP k is a multi-qubit Pauli operator defined aŝ τ k is a variational parameter, and n g is the number of Pauli operatorsP k included in the ansatz. The QCC energy functional is defined as is the set of n g variational parameters. The heuristic screening procedure [37] is utilized to construct the set of operators {P k } ng k=1 appearing in Eq. (C5). This heuristic approach [34] relies on the gradient of Eq. (C7) with respect to τ k evaluated using the optimal mean-field wave function (i.e., Γ = Γ opt and τ = 0), the form of which is is quantified for a representativeP k from each group of electron correlation generators contained in the direct interaction set (DIS) [37]. The n g representative generators from the DIS with the largest energy gradient magnitudes are utilized in Eq. (C5). Minimization of the energy functional given by Eq. (C7) with respect to the 2n q +n g parameter set T = Γ ∪ τ gives the ground state energy E QCC and the QCC correlated wave function |Ψ (T opt ) .
The DMET-QCC framework is applied to compute the energy of each point along the H 10 potential energy curve (see Table I and Fig. 2). The pre-optimized quantum circuits utilized for the trapped ion hardware experiments are constructed with T opt = Γ opt ∪ τ opt obtained from classical DMET-QCC simulations of the embedding Hamiltonian for fragment A (see Eq. (B6)). At each point R along the potential energy curve,Ĥ emb,A (R) is a twoqubit operator (n q = 2) after scBK encoding and takes the form The inter-atom, spacing-dependent expansion coefficients ofĤ emb,A (R) in Eq. (C9) are provided in Table II for the different R values considered in this manuscript. Minimization of the mean-field and QCC energy functionals given by Eqs. (C4) and (C7) is performed utilizing the L-BFGS-B solver [44,45] as implemented in the SciPy library [46] with convergence criteria ftol = 10 −8 a.u. and gtol = 10 −6 a.u. For the embedding Hamiltonians specified by Eq. (C9) and Table II, the DIS is characterized by a single flip index [37] F = {0, 1} that gives rise to a set of two electron correlation generators equivalent energy gradient magnitudes. For each DMET-QCC ansatz circuit considered in the present work, a single electron correlation generator is selected (n g = 1) and utilized in Eq. (C5):P 0 = X 0 Y 1 . Figure 3 shows a generic pre-optimized DMET-QCC ansatz circuit and illustrates its explicit dependence on the parameter set T = {θ 0 , θ 1 , φ 0 , φ 1 , τ 0 }. Table III provides the optimal set of parameters T opt obtained from classical DMET-QCC simulations along with the energy gradient magnitudes forP 0 = X 0 Y 1 and the converged QCC energies for each point on the H 10 potential energy curve.
ProjectQ [47] is employed to perform simulation on classical hardware. Details for further optimization of the DMET-QCC ansatz circuits are described in Appendix H.

D. HAMILTONIAN AND DMET TOTAL ENERGY EXPRESSION EMPLOYED IN THE EXPERIMENT
The pre-optimized circuit used for the experiment is obtained from the final iteration of DMET-QCC classical simulation. The VQE in conjunction with the QCC ansatz (VQE-QCC) is employed to minimize the energy of the subsystem in DMET. The Hamiltonians (see also Eq. (B6)) used to minimize the energy in the VQE-QCC calculation for the final iteration of DMET in the qubitized form are as follows: The above Hamiltonians are for the points R = 0.7, 1.0, 1.1, 1.3, and 1.6 A, respectively. These Hamiltonians are used to obtain the pre-optimized circuit, as described in Appendix C. Note that the Hamiltonians shown above are not exactly the same as the Hamiltonian shown in Eq. (B6). The orbitals are transformed into molecular orbitals using the coefficients C optimized in Hartree-Fock calculations of the subsystems.
The DMET energy expression is not the same as that of the Hamiltonians shown above. The total energy of the DMET calculation, shown in Eq. (B7), is calculated from the one-and two-particle reduced density matrices constructed from the wavefunction. The total energy per atom for the five points in the final DMET iteration can be expressed in an analytical form as follows: FIG. 3. Generic DMET-QCC ansatz circuit with explicit dependence on the mean-field and QCC parameters. Both single-qubit states are initialized as |0 . The DMET energy expressions are for the points R = 0.7, 1.0, 1.1, 1.3, and 1.6 A, respectively. The total energy per atom can be calculated using these expressions. Note that when density matrix purification is carried out, the two-particle reduced density matrix is constructed and purified; thus, these equations cannot be used to compute the total energies. The total energies are evaluated using Eq. (B7) with the purified reduced density matrices. Note that the expectation value Y Y is not needed to compute the total energies. When we compute the DMET-QCC energy, we first measure the Pauli operators XX, YY, ZZ, XZ, and ZX in the experiment. An additional experiment is needed that measures Y Y for post-processing.

E. DENSITY MATRIX PURIFICATION
We perform density matrix purification based on McWeeny's purification scheme [39]. This iterative method purifies the two-particle reduced density matrix P pqrs according to The iteration is continued until the convergence criterion is met. In this work, = 1.0×10 −2 is used. Note that the change in the total energy is smaller than a millihartree when we changed the criterion to 1.0 × 10 −7 . The tensor multiplication of P is defined as where the Einstein summation is implied. Note that this purification technique can only be applied to two-electron systems [10,13]. Although we here consider a 10-electron system, the reduced density matrix in our current work can be purified, as the fragment calculation for DMET involves a two-electron system.
In our preliminary investigation, two-particle reduced density matrices with large errors are largely improved by the density matrix purification. However, total energies of all points in this work are not improved, as the errors in them are very small. To understand the performance of the density matrix purification method in this regime, we study the dependency of the purified energies on the accuracy of the individual expectation values. In particular, we explore XZ , ZX , X 0 , and X 1 , whose ideal simulated values are extremely small. We find that the finite values measured via experiment have a significant effect on both the resulting energy and the success of the density matrix purification. Figure 4 shows a simulated analysis of the purified energy estimates and their dependence on XZ + ZX and X 0 + X 1 . In the ideal case, each term along the axes is equal, and contributes equally to the total energy. The experimental values are indicated by the blue dot, and the colour bar represents the error value in the purified energy, with chemical accuracy indicated by the yellow-green-blue region. The area within the dashed red lines shows the region where the purification is effective in improving the experimental results. Here we show two cases. In Fig. 4(a), for the specific case of R = 1.1 A, the purification improves the experimental results (the blue dot is within the dashed red lines), achieving chemical accuracy. Figure 4(b) is an example at the point R = 1.0 A, where the purification process does not improve the energy estimate, and the resulting difference is within the error range of the experimental data. The rest of the experimental points follow this trend. The range where the purification will improve the experimental result (the width between the dashed red lines) varies between (a) and (b), as it is dependent on the accuracy of the rest of the expectation values. This is corroborated by the series of experiments described in this manuscript, where the experimental values which are already within chemical accuracy do not improve after density matrix purification.
We investigate how the accuracy of Y Y , for which a two-qubit circuit is used, affects both the experimental energy and the post-processed energy by manually changing the Y Y expectation value while keeping the other expectation values fixed. We use the point R = 1.1 A, which allows us to succeed in improving the total energy value by employing density matrix purification. The analysis of the energy dependence on Y Y is depicted in Fig. 5. The values without purification (experimentally determined) are indicated in blue, the ideal values are shown using a black line, and the post-processed (purified) energy values are shown plotted in red. The region of energy values within chemical accuracy is shown using a shaded region. The experimental values that are not purified do not depend on Y Y . However, we observe that the value of Y Y affects how the energy is post-processed: the post-processed energy falls out of the chemical accuracy range when Y Y is smaller than 0.07 or larger than 0.23. The present analysis shows that density matrix purification with a large error in Y Y expectation value will not succeed. In this work, however, the trapped-ion hardware provides us with the Y Y expectation value of 0.176, which is accurate enough to bring the post-processed total energy within chemical accuracy.

F. ENTANGLEMENT ENTROPY
In this section, we describe the computation of two types of entanglement entropies. The first is the entanglement between different molecular orbitals (MO) and the second is the entanglement between fragment and bath. The embedding Hamiltonian given by Eq. (B6) is obtained by projecting the total system onto a fragment described by using the basis of fragment and bath orbitals. From this expression, the orbitals are transformed into the MO basis by using the MO coefficients C optimized in the HF calculation of the subsystem. The Hamiltonians used in the experiment are in this MO basis. Therefore, the entanglement between different MOs corresponds to the entanglement between different qubits. One of the motivations for looking at this quantity is to check whether the ground state of the fragment Hamiltonian is an entangled state, despite the fact that the optimized circuits (see Fig. 1(C)) for evaluating the embedding Hamiltonian contain only single-qubit gates. The other type of entanglement entropy, between the fragment and the bath, will capture the structure of the Schmidt decomposition given by Eq. (B1) and, therefore, the properties of DMET more directly.
We consider the entanglement entropies in the fermion picture. The elements of the density matrices we obtain via experimentation are a † p a † r a s a q and a † p a q , where the indices p, q, r, and s take values in the spin-orbitals {1α, 1β, 2α, 2β}. In order to split the Hilbert space based on orbitals, we describe the density matrix ρ in the basis of |1α, 1β, 2α, 2β 1α , 1β , 2α , 2β |. The RDM in our analysis is obtained by taking the trace over the degrees of freedom for orbital 2: It is straightforward to express the RDM in terms of fermionic operators: Note that off-diagonal elements in Eq. (F2) are zero due to the conservation of spins and electron numbers. These symmetries are guaranteed once the scBK transformation has been performed. Therefore, experimental noise will not change the structure of the RDM. Note that while there is a systematic way to expand the two-qubit RDM into the 4-qubit RDM in the sense that there is no additional experimental information required, the values of the entanglement entropies between these two expressions are different.
The entanglement entropy is defined by the equation where p i are the eigenvalues of the RDM.
First, we choose a basis such that 1α and 1β correspond to the first qubit and 2α and 2β correspond to the second qubit, both after scBK transformation. The results show that the entanglement entropies between the first qubit and the second qubit increase as the bond lengths increase.
We perform a transformation such that 1α and 1β R = 0.7 A R = 1.  represent spin-up and spin-down states of the fragment orbitals, respectively, while 2α and 2β represent those of the bath orbitals. To compute the entropies in the fragment-bath basis, we transform the basis into the tensor product of the fragment and the bath Hilbert spaces. In this case, the entanglement entropies decrease as the bond lengths increase. In the MO basis picture, the ground state in the absence of interactions is a product state (i.e., an HF state). Therefore, the entanglement entropies should be zero. Two-body interactions generate entanglement between different orbitals. Intuitively speaking, the Hilbert space is divided based on the energy levels (i.e., the Slater determinants) in this picture. The fact that the entanglement entropies in the MO picture are small reflects the fact that the ground states are close to the HF states.
On the other hand, in the fragment-bath orbital picture, the the Hilbert space is split depending on the orbitals' position in real space. Therefore, the ground state in the fragment-bath basis is highly entangled even when there is no interaction. More explicitly, the MO coefficients C (the transformation matrix between the fragment-bath basis and the MO basis) for R = 0.7 A is C = 0.70710679 0.70710677 −0.70710677 0.70710679 The last term in Eq. (F4) represents the discrete Fourier transformation. The entanglement entropy is close to the maximal value, which is log 4 = 2.

G. THE BOOTSTRAP METHOD
The total energies and their statistical errors are calculated using an empirical bootstrapping method. We follow the procedure described in a previous work [9]. We start from the state preparation and measurement (SPAM)-corrected histograms for each Pauli operator (XX, YY, ZZ, XZ, and ZX), and construct a distribution of total energies. The mean and the standard deviation (σ) are computed from the distribution of energies. The procedure of the bootstrapping method is as follows.
1. Draw a random bootstrap sample of the same size as the original dataset with replacement of the data.
2. Construct a new histogram based on step 1.
3. Compute the expectation value of the Pauli terms using the new histogram.
4. Repeat steps 1-3 for each Pauli term and obtain all expectation values needed to construct reduced density matrices.
5. Construct one-and two-particle reduced density matrices.
6. Calculate the total energy. 7. Repeat steps 1-6 500 times and obtain a distribution of total energies. 8. Calculate the mean and the standard deviation from the distribution of energies constructed in step 7.
We follow this procedure to construct a histogram of possible measurements consistent with the empirical data. The calculated value of σ is represented using error bars. The mean of this distribution is the measured energy, and the 1σ error estimate is the its standard deviation. The difference from the previous work [9] is the construction of the reduced density matrices (step 5), which is required for DMET energy calculation. Note that when we perform density matrix purification, we purify the two-particle reduced density matrix between steps 5 and 6.
We combine the circuit compilation and optimization techniques reported elsewhere [9,48,49]     Parameter specification of the gates that appear in the post-optimization circuits of Y Y in Fig. 1(C)(2)(b). All parameters are in radians.