Exponentially more precise quantum simulation of fermions in the configuration interaction representation

We present a quantum algorithm for the simulation of molecular systems that is asymptotically more efficient than all previous algorithms in the literature in terms of the main problem parameters. As in Babbush et al (2016 New Journal of Physics 18, 033032), we employ a recently developed technique for simulating Hamiltonian evolution using a truncated Taylor series to obtain logarithmic scaling with the inverse of the desired precision. The algorithm of this paper involves simulation under an oracle for the sparse, first-quantized representation of the molecular Hamiltonian known as the configuration interaction (CI) matrix. We construct and query the CI matrix oracle to allow for on-the-fly computation of molecular integrals in a way that is exponentially more efficient than classical numerical methods. Whereas second-quantized representations of the wavefunction require  ̃ ( N ) qubits, where N is the number of single-particle spin-orbitals, the CI matrix representation requires  ̃ ( η ) qubits, where η ≪ N is the number of electrons in the molecule of interest. We show that the gate count of our algorithm scales at most as  ̃ ( η 2 N 3 t ) .


I. INTRODUCTION
The first quantum algorithm for quantum chemistry was introduced nearly a decade ago [1].That algorithm was based on the Trotter-Suzuki decomposition, which Lloyd and Abrams first applied to quantum simulation in [2,3].The Trotter-Suzuki decomposition has been used in almost all quantum algorithms for quantum chemistry since then [4][5][6][7][8][9][10][11][12][13][14][15][16][17], with the exception of the adiabatic algorithm detailed in [18], the variational quantum eigensolver approach described in [19][20][21], and in our prior papers using the Taylor series technique [22,30].Recently, there has been substantial renewed interest in these algorithms due to the low qubit requirement compared with other algorithms such as factoring, together with the scientific importance of the electronic structure problem.This led to a series of papers establishing formal bounds on the cost of simulating various molecules [13][14][15][16][17].
Whereas qubit requirements for the quantum chemistry problem seem modest, using arbitrarily high-order Trotter formulas, the tightest-known upper bound on the gate count of the second-quantized, Trotter-based quantum simulation of chemistry is O(N 8+o(1) t/ o(1) ) [23,24] 1 , where N is the number of spin-orbitals and is the required accuracy.However, using significantly more practical Trotter decompositions, the best known gate complexity for this quantum algorithm is O(N 9 t 3 / ) [14].Fortunately, recent numerics suggest that the scaling for real molecules is closer to O(N 6 t 3 / ) [15] or O(Z 3 N 4 t 3 / ) [17], where Z is the largest nuclear charge in the molecule.Still, the Trotter-based quantum simulation of many molecular systems remains a costly proposition [25,26].
In Ref. [22], we introduced two novel quantum algorithms for chemistry based on the truncated Taylor series simulation method of [27], which are exponentially more precise than algorithms using the Trotter-Suzuki decomposition.Our first algorithm, referred to as the "database" algorithm, was shown to have gate count scaling as O(N 4 H t). Our second algorithm, referred to as the "on-the-fly" algorithm, was shown to have the lowest scaling of any approach to quantum simulation previously in the literature, O(N 5 t).Both of these algorithms use a second-quantized representation of the Hamiltonian; in this paper we employ a more compressed, first-quantized representation of the Hamiltonian known as the configuration interaction (CI) matrix.We also analyze the on-the-fly integration strategy far more rigorously, by making the assumptions explicit and rigorously deriving error bounds.Our approach combines a number of improvements: • a novel 1-sparse decomposition of the CI matrix (improving over that in [10]), • a self-inverse decomposition of 1-sparse matrices as introduced in [28], • the exponentially more precise simulation techniques of [27], • and the on-the-fly integration strategy of [22].
The paper is outlined as follows.In Section II, we summarize the key results of this paper, and note the improvements presented here over previous approaches.In Section III, we introduce the configuration basis encoding of the wavefunction.In Section IV, we show how to decompose the Hamiltonian into 1-sparse unitary matrices.In Section V, we use the decomposition of Section IV to construct a circuit which provides oracular access to the Hamiltonian matrix entries, assuming access to sample(w) from [22].In Section VI, we review the procedures in [27] and [22] to demonstrate that this oracle circuit can be used to effect a quantum simulation which is exponentially more precise than using a Trotter-Suzuki decomposition approach.In Section VII, we discuss applications of this algorithm and future research directions.

II. SUMMARY OF RESULTS
In our previous work [22], simulation of the molecular Hamiltonian was performed in second quantization using Taylor series simulation methods to give a gate count scaling as O(N 5 t).In this work, we use the configuration interaction representation of the Hamiltonian to provide an improved scaling of O(η 2 N 3 t).This result is summarized by the following Theorem.
Theorem 1 Using atomic units in which , Coulomb's constant, and the charge and mass of the electron are unity, we can write the molecular Hamiltonian as where R i are the nuclear coordinates, r j are the electron coordinates, and Z i are the nuclear atomic numbers.Consider a basis set of N spin-orbitals satisfying the following conditions: 1. each orbital takes significant values up to a distance at most logarithmic in N , 2. beyond that distance the orbital decays exponentially, 3. the maximum value of each orbital, and its first and second derivatives, scale at most logarithmically in N , 4. and the value of each orbital can be evaluated with complexity O(1).
Evolution under the Hamiltonian of Eq. (1) can be simulated in this basis for time t within error > 0 with a gate count scaling as O(η 2 N 3 t), where η is the number of electrons in the molecule.
We note that these conditions will be satisfied for most, but not all, quantum chemistry simulations.To understand the limitations of these conditions, we briefly discuss the concept of a model chemistry (i.e. standard basis set specifications) and how model chemistries are typically selected for electronic structure calculations.There are thousands of papers which study the effectiveness of various basis sets developed for the purpose of representing molecules [29].These model chemistries associate specific orbital basis functions with each atom in a molecule.For example, wherever Nitrogen appears in a molecule a model chemistry would mandate that one add to the system certain basis functions which are centered on Nitrogen and have been pre-optimized for Nitrogen chemistry; different basis functions would be associated with each Phosphorus, and so on.In addition to convenience, the use of standardized model chemistries helps chemists to compare different calculations and reproduce results.
Within a standard model chemistry, orbital basis functions are almost always represented as linear combinations of pre-fitted Gaussians which are centered on each atom.Examples of such model chemistries include Slater Type Orbitals (e.g.STO-3G), Pople Basis Sets (e.g.6-31G*) and correlation consistent basis sets (e.g.cc-DVTZ).We note that all previous studies on quantum algorithms for quantum chemistry in an orbital basis have advocated the use of one of these models.Simulation within any of these model chemistries would satisfy the conditions of our theorem because the basis functions associated with each atom have maximum values, derivatives and distances beyond which each orbital decays exponentially.
Similarly, when molecular instances grow because more atoms are added to the system it is standard practice to perform these progressively larger calculations using the same model chemistry and the conditions of Theorem 1 are satisfied.For instance, in a chemical series such as progressively larger Hydrogen rings or progressively longer alkane chains or protein sequences, these conditions would be satisfied.We note that periodic systems such as conducting metals might require basis sets (e.g.plane waves) violating the conditions of Theorem 1.When systems grow because atoms in the molecule are replaced with heavier atoms, the orbitals do tend to grow in volume and their maximum values might increase (even within a model chemistry).However, there are only a finite number of elements on the periodic table so this is irrelevant for considerations of asymptotic complexity.Finally, we point out that these conditions do not hold if the simulation is performed in the canonical molecular orbital basis, but this is not a problem for our approach since the Hartree-Fock state can easily be prepared in the atomic orbital basis at cost that is quadratic in the number of spin-orbitals.We discuss this procedure further in Section III.
The simulation procedure of Ref. [27] requires a decomposition of the Hamiltonian into a weighted sum of unitary matrices.In [22], we decomposed the molecular Hamiltonian in such a way that all the coefficients were integrals, i.e.
where the H are unitary operators, and the w ( z) are determined by the procedure.We then showed how one could evolve under H while simultaneously computing these integrals.In this paper, we investigate a different representation of the molecular Hamiltonian with the related property that the Hamiltonian matrix elements H αβ can be expressed as integrals, or a sum of a limited number of integrals.We decompose the Hamiltonian into a sum of one-sparse Hamiltonians, each of which has only a single integral in its matrix entries.We then decompose the Hamiltonian by discretizing the integrals and then further decompose the Hamiltonian into a sum of self-inverse operators, H ,ρ .Using this decomposition, we construct a circuit called select(H) which selects and applies the self-inverse operators so that By repeatedly calling select(H), we are able to evolve under H with an exponential improvement in precision over Trotter-based algorithms.
The CI matrix is a compressed representation of the molecular Hamiltonian that requires asymptotically fewer qubits than all second-quantized algorithms for chemistry.Though the CI matrix cannot be expressed as a sum of polynomially many local Hamiltonians, a paper by Toloui and Love [10] demonstrated that the CI matrix can be decomposed into a sum of O(N 4 ) 1-sparse Hermitian operators, where N is the number of spin-orbitals.We provide in this paper a new decomposition of the CI matrix into a sum of O(η 2 N 2 ) 1-sparse Hermitian operators, where η N is the number of electrons in the molecule.This new decomposition enables our improved scaling.Using techniques introduced in [28], we further decompose these 1-sparse operators into unitary operators which are also self-inverse.As a consequence of the self-inverse decomposition, the Hamiltonian is an equally weighted sum of unitaries.select(H) requires the ability to compute the entries of the CI matrix; accordingly, we can use the same strategy for computing integrals on-the-fly that was introduced in [22], but this time our Hamiltonian is of the form in Eq. (3).
Using this approach, the simulation of evolution over time t then requires O(η 2 N 2 t) calls to select(H).To implement select(H), we make calls to the CI matrix oracle as described in Section V, which requires O(N ) gates.This scaling is due to using a database approach to computing the orbitals, where a sequence of N controlled operations is performed.This causes our overall approach to require O(η 2 N 3 t) gates.As in [10], the number of qubits is O(η) rather than O(N ), because the compressed representation stores only the indices of occupied orbitals, rather than occupation numbers of all orbitals.To summarize, our algorithm with improved gate count scaling of O(η 2 N 3 t) proceeds as follows: 1. Represent the molecular Hamiltonian in Eq. (1) in first quantization using the CI matrix formalism.This requires selection of a spin-orbital basis set, chosen such that the conditions in Theorem 1 are satisfied.
2. Decompose the Hamiltonian into sums of self-inverse matrices approximating the required molecular integrals via the method of Section IV.
3. Query the CI matrix oracle to evaluate the above self-inverse matrices, which we describe in Section V.
4. Simulate the evolution of the system over time t using the method of [27], which is summarized in Section VI.

III. THE CI MATRIX ENCODING
The molecular electronic structure Hamiltonian describes electrons interacting in a nuclear potential that is fixed under the Born-Oppenheimer approximation.Except for the proposals in [9][10][11][12]30], all prior quantum algorithms for chemistry use second quantization.While in second quantization antisymmetry is enforced by the fermionic anticommutation relations, in first quantization the wavefunction itself is explicitly antisymmetric.The representation of Eq. (1) in second quantization is where the operators a † i and a j in Eq. ( 5) obey antisymmetry due to the fermionic anti-commutation relations, The one-electron and two-electron integrals in Eq. ( 5) are where (throughout this paper), r j represents the position of the j th electron, and ϕ i ( r j ) represents the i th spin-orbital when occupied by that electron.To ensure that the integrand in Eq. ( 7) is symmetric, we can alternatively write the integral for h ij as The second-quantized Hamiltonian in Eq. ( 5) is straightforward to simulate because one can explicitly represent the fermionic operators as tensor products of Pauli operators, using either the Jordan-Wigner transformation [31,32] or the Bravyi-Kitaev transformation [33][34][35].
With the exception of real-space algorithms described in [9,30], all quantum algorithms for chemistry represent the system in a basis of N single-particle spin-orbital functions, usually obtained as the solution to a classical mean-field treatment such as Hartree-Fock [36].However, the conditions of Theorem 1 only hold when actually performing the simulation in the atomic orbital basis2 (i.e. the basis prescribed by the model chemistry).The canonical Hartree-Fock orbitals are preferred over the atomic orbitals because initial states are easier to represent in the basis of Hartree-Fock orbitals.These orbitals are actually a unitary rotation of the orthogonalized atomic orbitals prescribed by the model chemistry.This unitary basis transformation takes the form φi and κ is anti-Hermitian.For κ and U , the quantities κ ij and U ij respectively correspond to the matrix elements of these operators in the basis of spin orbitals.It is a consequence of the Thouless theorem that this unitary transformation is efficient to apply.The canonical Hartree-Fock orbitals and κ are obtained by performing a self-consistent field procedure to diagonalize a mean-field Hamiltonian for the system which is known as the Fock matrix.Because the Fock matrix describes a system of non-interacting electrons it can be expressed as the following N by N matrix: The integrals which appear in the Fock matrix are defined by Eq. ( 7) and Eq. ( 8).Importantly, the canonical orbitals are defined to be the orbitals which diagonalize the Fock matrix.Thus, the integrals in the definition of the Fock matrix are defined in terms of the eigenvectors of the Fock matrix so Eq. ( 12) is a recursive definition.The canonical orbitals are obtained by repeatedly diagonalizing this matrix until convergence with its own eigenvectors.The Hartree-Fock procedure is important because the Hartree-Fock state (which is a product state in the canonical basis with the lowest η eigenvectors of the Fock matrix occupied and the rest unoccupied) has particularly high overlap with the ground state of H.As stated before, the conditions of Theorem 1 do not apply if we represent the Hamiltonian in the basis of canonical orbitals.But this is not a problem for us because we can still prepare the Hartree-Fock state in the basis of orthgonalized atomic orbitals (which do satisfy the conditions) and then apply the operator U = e −κ to our initial state at cost O(N 2 ).Note that the use of a local basis has other advantages, as pointed out in [16].In particular, in the limit of certain large molecules, use of a local basis allows one to truncate terms from the Hamiltonian so that there are O(N 2 ) terms instead of O(N 4 ) terms.However, Theorem 1 exploits an entirely different property of basis locality which does not require any approximation from truncating terms.
The spatial encoding of Eq. ( 5) requires Θ(N ) qubits, one for each spin-orbital; under the Jordan-Wigner transformation, the state of each qubit indicates the occupation of a corresponding spin-orbital.Many states representable in second quantization are inaccessible to molecular systems due to symmetries in the Hamiltonian.For instance, molecular wavefunctions are eigenstates of the total spin operator so the total angular momentum is a good quantum number, and this insight can be used to find a more efficient spatial encoding [11,12].Similarly, the Hamiltonian in Eq. ( 5) commutes with the number operator, ν, whose expectation value gives the number of electrons, η, Following the procedure in [10], our algorithm makes use of an encoding which reduces the number of qubits required by recognizing η as a good quantum number.Conservation of particle number implies there are only ξ = N η valid configurations of these electrons, but the second-quantized Hilbert space has dimension 2 N , which is exponentially larger than ξ for fixed η.We work in the basis of Slater determinants, which are explicitly antisymmetric functions of both space and spin associated with a particular η-electron configuration.We denote these states as |α = |α 0 , α 1 , • • • , α η−1 , where α i ∈ {1, . . ., N } and α ∈ {1, . . ., N η }.We emphasize that α i is merely an integer which indexes a particular spin-orbital function ϕ αi ( r).While each configuration requires a specification of η occupied spin-orbitals, there is no sense in which α i is associated with "electron i" since fermions are indistinguishable.Specifically, where the bars enclosing the matrix in Eq. ( 14) denote a determinant.Because determinants have the property that they are antisymmetric under exchange of any two rows, this construction ensures that our wavefunction obeys the Pauli exclusion principle.We note that although this determinant can be written equivalently in different orders (e.g. by swapping any two pairs of orbital indices), we avoid this ambiguity by requiring the Slater determinants to only be written in ascending order of spin-orbital indices.
The representation of the wavefunction introduced in [10] uses η distinct registers to encode the occupied set of spin-orbitals, thus requiring Θ(η log N ) = O(η) qubits.However, it would be possible to use a further-compressed representation of the wavefunction based on the direct enumeration of all Slater determinants, requiring only Θ(log ξ) qubits.When using very small basis sets (such as the minimal basis), it will occasionally be the case that the spatial overhead of Θ(N ) for the second-quantized algorithm is actually less than the spatial complexity of our algorithm.However, for a fixed η, the CI matrix encoding requires exponentially fewer qubits.

IV. THE CI MATRIX DECOMPOSITION
The molecular Hamiltonian expressed in the basis of Slater determinants is known to chemists as the CI matrix.Elements of the CI matrix are computed according to the Slater-Condon rules [36], which we will express in terms of the one-electron and two-electron integrals in Eq. (7) and Eq.(8).In order to motivate our 1-sparse decomposition, we state the Slater-Condon rules for computing the matrix element by considering the spin-orbitals which differ between the determinants |α and |β [36]: 1.If |α and |β contain the same spin-orbitals {χ i } η i=1 then we have a diagonal element 2. If |α and |β differ by exactly one spin-orbital such that |α contains spin-orbital k where |β contains spinorbital , but otherwise contain the same spin-orbitals {χ i } η−1 i=1 , then 4. If |α and |β differ by more than two spin-orbitals, These rules assume that α and β have the list of occupied orbitals given in a corresponding order, so all corresponding occupied orbitals are listed in the same positions.In contrast, we will be giving the lists of occupied orbitals in ascending order.In order to use the rules, we therefore need to change the order of the list of occupied orbitals.In changing the order of the occupied orbitals, there is a sign flip on the state for an odd permutation.This sign flip needs to be included when using the above rules.
In general, there is no efficient way to decompose the CI matrix into a polynomial number of tensor products of Pauli operators.It is thus inefficient to directly simulate this Hamiltonian in the same fashion with which we simulate local Hamiltonians.However, the CI matrix is sparse and there exist techniques for simulating arbitrary sparse Hamiltonians.A d-sparse matrix is one which contains at most d nonzero elements in each row and column.As discussed in [10,13], the Slater-Condon rules imply that the sparsity of the CI matrix is Because N is always greater than η, we find that the CI matrix is d-sparse where d ∈ O(η 2 N 2 ).This should be compared with the second-quantized Hamiltonian which is also d-sparse, but where d ∈ O(N 4 ).Our strategy here parallels the second-quantized decomposition, but works with the first-quantized wavefunction.This decomposition is explained in four steps, as follows.
B. Further decompose each of these 1-sparse matrices into 1-sparse matrices with entries proportional to a sum of a constant number of molecular integrals.
C. Decompose those 1-sparse matrices into sums approximating the integrals in Eqs. ( 8) and (9).D. Decompose the integrands from those integrals into sums of self-inverse matrices.
FIG. 1. Example of the 1-sparse coloring, where i is the position of the occupied orbital in α that must be moved.(a) i = 2, p = 4 is sufficient to determine β from α, as well as to determine α from β.(b) i = 2, p = 5 is sufficient to determine β from α, but not the reverse: subtracting p = 5 from β2, β3, or β4 all give different valid values for αi = α2.The spacing condition means that we would need to give the position of the occupied orbital for β instead.

A. Decomposition into 1-sparse matrices
In order to decompose the molecular Hamiltonian into 1-sparse matrices, we require a unique and reversible graph coloring between nodes (Slater determinants).We introduce such a graph coloring here, with the details of its construction and proof of its properties given in Appendix A. The graph coloring can be summarized as follows.
1. Perform the simulation under σ x ⊗ H, where σ x is the Pauli x matrix, in order to create a bipartite Hamiltonian of the same sparsity as H.
2. Label the left nodes α and the right nodes β.We seek a procedure to take α to β, or vice versa, with as little additional information as possible, and without redundancy or ambiguity.
3. Provide an 8-tuple γ = (a 1 , b 1 , i, p, a 2 , b 2 , j, q) which determines the coloring.The coloring must uniquely determine α given β or vice versa.Using the 8-tuples, proceed via either Case 1, 2, 3, or 4 in Appendix A to determine the other set of spin-orbitals, using an intermediate list of orbitals χ.The 4-tuples (a 1 , b 1 , i, p) and (a 2 , b 2 , j, q) each define a differing orbital.For a single difference, we can set p = 0, and for no differences, we can set p = q = 0.
The basic idea is that we give the positions i and j of those orbitals which differ in α, as well as by how much the occupied orbital indices shift, which we denote by p and q.This allows us to determine β from α.However, it does not allow us to unambiguously determine α from β.To explain how to resolve this ambiguity, we consider the case of a single differing orbital.We will denote by i the position of the differing orbital in α, and by k the position of the differing orbital in β.
Consider the example in Figure 1(a): given i which is the position in α, the position k in β can be immediately determined.But given β, multiple potential positions of occupied orbitals would need to be tested to see if they put the occupied orbital in position i = 2 in α.In this case, given β there is only one orbital which can be shifted to position 2 for α so the position in β is unambiguous.Now consider Figure 1 The difference between the two cases is that in Figure 1(a) there is a larger spacing between orbitals for β, whereas in Figure 1(b) there is a larger spacing for α.More specifically, for Figure 1(a) the spacing between α 1 and α 3 is 3, whereas the spacing between β 2 and β 4 is larger at 5. For Figure 1(b) the spacing between α 1 and α 3 is 5, whereas the spacing between β 2 and β 4 is smaller at 2. It is the spacing between the occupied orbitals adjacent to the one that is moved that should be compared.
For the situation in Figure 1(b), rather than specifying the position in α we should specify the position in β to resolve the ambiguity.The bit a determines whether we are specifying the position in α or in β; this is done depending on the relative spacing of the adjacent occupied orbitals in the two.However, this spacing condition does not completely resolve the ambiguity: there are potentially two different choices for the occupied orbital.The choice is made by the bit b.The coloring for the two differing orbitals is done by doing this twice with an intermediate list of occupied orbitals χ.
There are O(η 2 N 2 ) possible colors: there are two possible choices of each of the bits a 1 , a 2 , b 1 , and b 2 , η choices each of i and j, and N choices each of p and q.
B. Decomposition into hij and h ijk Each 1-sparse matrix from Section IV A is associated with some 8-tuple γ = (a 1 , b 1 , i, p, a 2 , b 2 , j, q).However, without further modification, some of these 1-sparse matrices have entries given by a sum over a number of molecular integrals that grows with η, namely the diagonal terms as in Eq. ( 16), and the single-orbital terms as in Eq. ( 17).Here, we further decompose those matrices into a sum of 1-sparse matrices H γ , which have entries proportional to the sum of a constant number of molecular integrals, in order to remove this changing upper bound.
We want to have a new set of 1-sparse matrices, each with entries corresponding to a single term in the sum over molecular integrals.To be more specific, the combinations of γ correspond to terms in Eq. ( 16) to Eq. ( 18) as follows.
1.If p = q = 0, this indicates that we have a diagonal 1-sparse matrix.In Eq. ( 16), the entries on the diagonal would be a sum of O(η 2 ) terms.As we have freedom in how to use i and j, we use these to give terms in the sum.When i = j for p = q = 0, we take the 1-sparse matrix to have diagonal elements given by h ii .If i < j for p = q = 0 we take the 1-sparse matrix to have diagonal entries h χiχj χiχj − h χiχj χj χi .We do not allow tuples γ such that i > j for p = q = 0 (alternatively we could just give zero in this case).The overall result is that the sum over i and j for the 1-sparse matrices for γ with p = q = 0 yields the desired sum in Eq. ( 16).
2. Next, if p = 0 and q = 0, then this indicates that we have a 1-sparse matrix with entries where α and β differ by only one spin-orbital.According to Eq. ( 17), each entry would normally be a sum of O(η) terms.Instead, when p = 0 and q = 0, we use the value of i to index terms in the sum in Eq. ( 17), though we only yield a nonzero result when i is in the Slater determinant.In particular, the 1-sparse matrix has entries h kχi χi − h kχiχi .We allow an additional value of i to indicate a 1-sparse matrix with entries h k .Then the sum over 1-sparse matrices for different values of i gives the desired sum Eq. ( 17).We do not allow γ such that q = 0 but p = 0.
3. Finally, if both p and q are nonzero, then we have a 1-sparse matrix with entries where α and β differ by two orbitals.In this case, there is no sum in Eq. ( 18), so there is no additional decomposition needed.
Combining these three steps we find that the decomposition into 1-sparse matrices H γ can be achieved with the indices (a 1 , b 1 , i, p, a 2 , b 2 , j, q).Thus, there are O(η 2 N 2 ) terms without any redundancies.Note that sorting of the spin-orbital indices requires only O(η) gates, which is less than the number of complexity of evaluating the spinorbitals.In the following sections, we denote the total number of terms given by the above decomposition by Γ, and the sum over H γ yields the complete CI matrix,

C. Discretizing the integrals
Next we consider discretization of the integrals for h ij and h ijk .In [27] it is shown how to simulate Hamiltonian evolution with an exponential improvement in the scaling with 1/ , as compared to methods based on Trotter formulas.In this approach, the time-ordered exponential for the evolution operator is approximated by a Taylor series up to an order K.The time t is broken into r segments, and the integrals are discretized in the following way on each segment: where T is the time-ordering operator.In our case the Hamiltonian does not change in time, so the time-ordering is unimportant.
The Hamiltonian is expanded as a sum of H γ as in Eq. ( 21), and each of those terms has matrix entries that can be given in the form of an integral as In cases where H αβ γ corresponds to h ij , the integral is over a three-dimensional region, and where H αβ γ corresponds to h ijk the integral is over a six-dimensional region, so z represents six parameters.
Ideally, each integral can be truncated to a finite domain D with volume V. Using a set of grid points z ρ , we can approximate the integral by The complexity will then be logarithmic in the number of points in the sum, µ, and linear in the volume times the maximum value of the integrand.
In practice the situation is more complicated than this.That is because the integrals are all different.As well as the dimensionality of the integrals (three for h ij and six for h ijk ), there will be differences in the regions that the integrals will be over, as well as some integrals being in spherical polar coordinates.To account for these differences, it is better to write the discretized integral in the form The Hamiltonian H γ can then be written as the sum As discussed in [22], the discretization is possible because the integrands can be chosen to decay exponentially [36].The required properties of the orbitals are given in Theorem 1.Here we present a more precise formulation of the required properties, and provide specific results on the number of terms needed.We make the following three assumptions about the spin-orbitals ϕ .
1.There exists a positive real number ϕ max such that, for all spin-orbital indices and for all r ∈ IR 3 , 2. For each spin-orbital index , there exists a vector c ∈ IR 3 (called the center of ϕ ) and a positive real number where α is some positive real constant.
3. For each spin-orbital index , ϕ is twice-differentiable and there exist positive real constants γ 1 and γ 2 such that and for all r ∈ IR 3 .
Note that α, γ 1 and γ 2 are dimensionless constants, whereas x max has units of distance, and ϕ max has the same units as ϕ .The conditions of Theorem 1 mean that ϕ max and x max grow at most logarithmically with the number of spin-orbitals.Note that we use x max in a different way than in [22], where it was the size of the cutoff on the region of integrals, satisfying x max = O(log(N t/ )).Here we take x max to be the size scale of the orbitals independent of t or , and the cutoff will be a multiple of x max .We also assume that x max is bounded below by a constant, so the first and second derivatives of the spin-orbitals grow no more than logarithmically as a function of the number of spin-orbitals.We next define notation used for the integrals for h ij and h ijk .These integrals are and for any choices of D 0 , D 1,q ⊆ IR 3 and D 2 ⊆ IR 6 .Thus and Using the assumptions on the properties of the orbitals, we can bound the number of terms needed in a Riemann sum that approximates each integral to within a specified accuracy, δ (which is distinct from the accuracy of the overall simulation, ).These bounds are summarized in the following three lemmas.
Lemma 1 Let δ be any real number that satisfies where Then S (0) ij IR 3 can be approximated to within error δ using a Riemann sum with terms, where the terms in the sum have absolute value no larger than Lemma 2 Let δ be any real number that satisfies where Then S (1,q) ij IR 3 can be approximated to within error δ using a Riemann sum with terms, where the terms in the sum have absolute value no larger than Lemma 3 Let δ be any real number that satisfies where Then S (2) ijk IR 6 can be approximated to within error δ using a Riemann sum with terms, where the terms in the sum have absolute value no larger than The conditions in Eqs. ( 36), ( 40) and ( 44) are just used to ensure that we are considering a reasonable combination of parameters, and not for example a very large allowable error δ or a small value of x max .We prove these Lemmas in Appendix B. Specifically, we prove Lemma 1 in Appendix B 2, Lemma 2 in Appendix B 3 and Lemma 3 in Appendix B 4. In discretizing these integrals it is important that the integrands are Hermitian, because we need H γ,ρ to be Hermitian.The integrands of these integrals are not Hermitian as discretized in the way given in the proofs in Appendix B. This is because the regions of integration are chosen in a non-symmetric way.For example, the region of integration for S (0) ij is chosen centered on the orbital ϕ i , so the integrand is not symmetric.It is simple to symmetrize the integrands, however.For example, for S (0) ij we can add (S (0) ji ) * and divide by two.That ensures that the integrand is symmetric, with just a factor of two overhead in the number of terms in the sum.
As a consequence of these Lemmas, we see that the terms of any Riemann sum approximation to one of the integrals that define the Hamiltonian coefficients h ij and h ijk have absolute values bounded by where µ is the number of terms in the Riemann sum and δ is the desired accuracy of the approximation.Here we have taken Z q to be O(1).

D. Decomposition into self-inverse operators
The truncated Taylor series strategy introduced in [27] requires that we can represent our Hamiltonian as a weighted sum of unitaries.To do so, we follow a procedure in [28] which shows how 1-sparse matrices can be decomposed into a sum of self-inverse matrices with eigenvalues ±1.Specifically, we decompose each ℵ γ,ρ into a sum of M ∈ Θ max γ,ρ ℵ γ,ρ max /ζ 1-sparse unitary matrices of the form where ζ is the desired precision of the decomposition.First, we construct a new matrix ℵ γ,ρ by rounding each entry of ℵ γ,ρ to the nearest multiple of 2 ζ, so that We decompose each C γ,ρ into C γ,ρ max 1-sparse matrices, indexed by m, with entries in {0, −2, 2}, as follows: Finally, we remove zero eigenvalues by further dividing each C γ,ρ,m into two matrices C γ,ρ,m,1 and C γ,ρ,m,2 with entries in {0, −1, +1}.For every all-zero column β in C γ,ρ,m , we choose α so that (α, β) is the location of the nonzero entry in column β in the original matrix H γ,ρ .Then the matrix C γ,ρ,m,1 has +1 in the (α, β) position, and C γ,ρ,m,2 has −1 in the (α, β) position.Thus, we have decomposed each H γ into a sum of 1-sparse, unitary matrices with eigenvalues ±1.
We now use a simplified notation where corresponds to the triples (s, m, γ), and ℵ ,ρ ≡ C γ,ρ,m,s .We denote the number of values of by L, and can write the Hamiltonian as a sum of O(N 4 µM ) unitary, 1-sparse matrices That is, the decomposition is of the form in Eq. ( 2), but in this case W is independent of .To summarize, we decompose the molecular Hamiltonian into a sum of self-inverse matrices in four steps: 1. Decompose the molecular Hamiltonian into a sum of 1-sparse matrices using the bipartite graph coloring given in Appendix A, summarized in Section IV A.
2. Decompose these 1-sparse matrices further, such that each entry corresponds to a single term in the sum over molecular integrals.This does not change the number of terms, but simplifies calculations.
3. Discretize the integrals over a finite region of space, subject to the constraints and bounds given in [22].
4. Decompose into self-inverse operators by the method proposed in [28].
This decomposition gives an overall gate count scaling contribution of O(η 2 N 2 ).

V. THE CI MATRIX ORACLE
In this section, we discuss the construction of the circuit referred to in our introduction as select(H), which applies the self-inverse operators in a controlled way.As discussed in [22], the truncated Taylor series technique of [27] can be used with a selection oracle for an integrand which defines the molecular Hamiltonian.This method will then effect evolution under this Hamiltonian with an exponential increase in precision over Trotter-based methods.For clarity of exposition, we describe the construction of select(H) in terms of two smaller oracle circuits which can be queried to learn information about the 1-sparse unitary integrands.This information is then used to evolve an arbitrary quantum state under a specific 1-sparse unitary.
The first of the oracles described here is denoted as Q col and is used to query information about the sparsity pattern of a particular 1-sparse Hermitian matrix from Eq. ( 21).The second oracle is denoted as Q val and is used to query information about the value of integrands for elements in the CI matrix.We construct Q val by making calls to a circuit constructed in [22] where it is referred to as "sample(w)".The purpose of sample(w) is to sample the integrands of the one-electron and two-electron integrals h ij and h ijk in Eq. ( 8) and Eq. ( 9).The construction of sample(w) in [22] requires O(N ) gates.
The oracle Q col uses information from the index γ.The index γ is associated with the indices (a 1 , b 1 , i, p, a 2 , b 2 , j, q) which describe the sparsity structure of the 1-sparse Hermitian matrix H γ according to the decomposition in Section IV B. Q col acts on a register specifying a color |γ as well a register containing an arbitrary row index |α to reveal a column index |β so that the ordered pair (α, β) indexes the nonzero element in row α of H γ , From the description in Section IV B, implementation of the unitary oracle Q col is straightforward.
To construct select(H) we need a second oracle that returns the value of the matrix elements in the decomposition.This selection oracle is queried with a register | = |s |m |γ which specifies which part of the 1-sparse representation we want, as well as a register |ρ which indexes the grid point ρ and registers |α and |β specifying the two Slater determinants.Specifically, the entries in the tuple identify the color (γ) of the 1-sparse Hermitian matrix from which the 1-sparse unitary matrix originated, which positive integer index (m ≤ M ) it corresponds to in the further decomposition of ℵ γ,ρ into C γ,ρ,m , and which part it corresponds to in the splitting of C γ,ρ,m into C γ,ρ,m,s (where s ∈ {1, 2}).
As a consequence of the Slater-Condon rules shown in Eqs. ( 16), ( 17), (18), and ( 19), Q val can be constructed given access to sample(w), which samples the integrand of the integrals in Eqs. ( 8) and ( 9) [22].Consistent with the decomposition in Section IV B, the i and j indices in the register containing γ = (i, p, j, q) specify the dissimilar spin-orbitals in |α and |β that are needed in the integrands defined by the Slater-Condon rules; therefore, the determination of which spin-orbitals differ between |α and |β can be made in O(log N ) time (only the time needed to read their values from γ).As sample(w) is comprised of O(N ) gates, Q val has time complexity O(N ) and acts as where H αβ ,ρ is the value of the matrix entry at (α, β) in the self-inverse matrix H ,ρ .When either |α or |β represents an invalid Slater determinant (with more than one occupation on any spin-obital), we take H αβ ,ρ = 0 for α = β.This ensures there are no transitions into Slater determinants which violate the Pauli exclusion principle.The choice of H αα ,ρ for invalid α will not affect the result, because the state will have no weight on the invalid Slater determinants.
Having constructed the column and value oracles, we are finally ready to construct select(H).This involves implementing 1-sparse unitary operations.The method we describe is related to the scheme presented in [37] for evolution under 1-sparse Hamiltonians, but is simplified due to the simpler form of the operators.As in Eq. ( 4), select(H) applies the term H ,ρ in the 1-sparse unitary decomposition to the wavefunction |ψ .Writing |ψ = α c α |α , we require that select(H) first call Q col to obtain the columns, β, corresponding to the rows, α, for the nonzero entries of the Hamiltonian: Now that we have the row and column of the matrix element, we apply Q val which causes each Slater determinant to accumulate the phase factor k α = H αβ ,ρ = ±1: Next, we swap the locations of α and β in order to complete multiplication by the 1-sparse unitary, Finally we apply Q col again but this time β is in the first register.Since Q col is self-inverse and always maps |α |b to |α |b ⊕ β and |β |b to |β |b ⊕ α , this allows us to uncompute the ancilla register.
Note that this approach works regardless of whether the entry is on-diagonal or off-diagonal; we do not need separate schemes for the two cases.The circuit for select(H) is depicted in Figure 2.
The circuit implementing select(H), which applies the term H ( zρ) labeled by = (γ, m, s) in the unitary 1-sparse decomposition to the wavefunction |ψ .

VI. SIMULATING HAMILTONIAN EVOLUTION
The simulation technique we now discuss is based on that of Ref. [27].We partition the total time t into r segments of duration t/r.For each segment, we expand the evolution operator e −iHt/r in a Taylor series up to order K, The error due to truncating the series at order K is bounded by In order to ensure the total simulation has error no greater than , each segment should have error bounded by /r.Therefore, provided r ≥ H t, the total simulation will have error no more than for Using our full decomposition of the Hamiltonian from Eq. (51) in the Taylor series formula of Eq. (58), we obtain The sum in Eq. (61) takes the form where U is close to unitary and the V j are unitary.Note that in contrast to [22], β j ≥ 0, consistent with the convention used in [27].Our simulation will make use of an ancillary "selection" register It is convenient to encode k in unary, as |k = |1 k 0 K−k , which requires Θ(K) qubits.Additionally, we encode each | k in binary using Θ(log L) qubits and each |ρ k in binary using Θ(log µ) qubits.We denote the total number of ancilla qubits required as J, which scales as To implement the truncated evolution operator in Eq. (62), we wish to prepare a superposition state over j, then apply a controlled V j operator.Following the notation of Ref. [27], we denote the state preparation operator as B, and it has the effect FIG. 3. The circuit for B as described in Eq. ( 64).An expression for θ k is given in Eq. ( 66).
where λ is a normalization factor.We can implement B by applying Hadamard gates (which we denote as Hd) to every qubit in the | υ and |ρ υ registers, in addition to K (controlled) single-qubit rotations applied to each qubit in the |k register.We need to prepare the superposition state over |k Using the convention that R y (θ k ) := exp[−i θ k σ y /2], this can be prepared by applying R y (θ 1 ) followed by where The Hadamard gates are applied to each of the qubits in the 2K remaining components of the selection register This requires O(K log(µL)) gates; parallelizing the Hadamard transforms leads to circuit depth O(K) for B. A circuit implementing B is shown in Figure 3.
We then wish to implement an operator to apply the V j which is referred to in [22,27] as select(V ), This operation can be applied using O(K) queries to a controlled form of the oracle select(H) defined in Section V.One can apply select(H) K times, using each of the | υ and |ρ υ registers.Thus, given that select(H) requires O(N ) gates, our total gate count for select(V ) is O(KN ).Table I lists relevant parameters along with their bounds, in terms of chemically relevant variables.A strategy for implementing the evolution operator in Eq. ( 62) becomes clear if we consider the combined action of B followed by select(V ) on state |ψ , where |Φ is a state for which the selection register ancilla qubits are orthogonal to |0 ⊗J .Accordingly, we define a projection operator P so that Using these definitions, we follow the procedure for oblivious amplitude amplification of [27] to deterministically execute the intended unitary.To this end, we define the oblivious amplitude amplification operator We use the amplification operator G in conjunction with the projector onto the empty ancilla state P so that The sum of the absolute values of the coefficients in the self-inverse decomposition of the Hamiltonian in Eq. ( 51) is ζLµ.If we choose the number of segments as r = ζLµt/ ln(2), then our choice of K as in Eq. (60) ensures that U − U r max ∈ O( /r), and hence [27] Then the total error due to oblivious amplitude amplification on all segments will be O( ).Therefore, the complexity of the total algorithm is r times the complexity of implementing select(V ) and B. While we implement B with gate count O(K log(µL)), our construction of select(H) has gate count O(N K).
The gate count of our entire algorithm depends crucially on r.Above we have taken r ∈ O(ζLµt) where As a result, we may bound r as As a consequence of Lemma 1, Lemma 2 and Lemma 3, µ max γ,ρ ℵ γ,ρ max can be replaced with where ϕ max is the maximum value taken by the orbitals, and x max is the scaling of the spatial size of the orbitals.
To relate to δ, in Section IV B the Hamiltonian is broken up into a sum of O(η 2 N 2 ) terms, each of which contains one or two of the integrals.Therefore, the error in the Hamiltonian is O(δη 2 N 2 ).The Hamiltonian is simulated for time t, so the resulting error in the simulation will be O(δtη 2 N 2 ).To ensure that the error is no greater than , we should therefore choose δ = Θ( /(tη 2 N 2 )).Since we are considering scaling with large η and N , δ will be small and the conditions Eq. ( 36), Eq. ( 40) and Eq.(44) will be satisfied.In addition, the conditions of Theorem 1 mean that ϕ max and x max are logarithmic in N .Hence one can take, omitting logarithmic factors, The complexity of B does not affect the scaling, because it is lower order in N .Therefore, our overall algorithm has gate count as stated in Theorem 1.This scaling represents an exponential improvement in precision as compared to Trotter-based methods.However, we suspect that the actual scaling of these algorithms is much better for real molecules, just as has been observed for the Trotter-based algorithms [15,17].Furthermore, the approach detailed here requires fewer qubits than any other approach to quantum simulation of chemistry in the literature.

VII. DISCUSSION
We have outlined a method to simulate the quantum chemistry Hamiltonian in a basis of Slater determinants using recent advances from the universal simulation literature.We find an oracular decomposition of the Hamiltonian into 1-sparse matrices based on an edge coloring routine first described in [10].We use that oracle to simulate evolution under the Hamiltonian using the truncated Taylor series technique described in [28].We discretize the integrals which define entries of the CI matrix, and use the sum of unitaries approach to effectively exponentially compress evaluation of these discretized integrals.
Asymptotic scalings suggest that the algorithms described in this paper series will allow for the quantum simulation of much larger molecular systems than would be possible using a Trotter-Suzuki decomposition.Recent work [13][14][15][16][17] has demonstrated markedly more efficient implementations of the original Trotter-Suzuki-based quantum chemistry algorithm [1,38]; similarly, we believe the implementations discussed here can still be improved upon, and that numerical simulations will be crucial to this task.
Finally, we note that the CI matrix simulation strategy discussed here opens up the possibility of an interesting approach to adiabatic state preparation.An adiabatic algorithm for quantum chemistry was suggested in second quantization in [18] and studied further in [39].However, those works did not suggest a compelling adiabatic path to take between an easy-to-prepare initial state (such as the Hartree-Fock state) and the ground state of the exact Hamiltonian.We note that one could start the system in the Hartree-Fock state, and use the CI matrix oracles discussed in this paper to "turn on" a Hamiltonian having support over a number of configuration basis states which increases smoothly with time.
to the division of that graph into Γ sets of disjoint graphs of degree 1, this edge coloring represents a decomposition of the Hamiltonian into Γ 1-sparse matrices.Aharonov and Ta-Shma show a procedure for accomplishing the 1-sparse decomposition of any arbitrary d-sparse matrix using Θ(d2 ) terms by coloring an arbitrary graph of degree d with Θ(d 2 ) colors.This result was tightened from Θ(d 2 ) terms to d 2 terms in [28].Importantly, Aharonov and Ta-Shma also showed how these Hamiltonians can be efficiently simulated using an oracular scheme based on the Trotter-Suzuki decomposition.Toloui and Love used this result to show how the CI matrix can be efficiently simulated under Trotter-Suzuki decomposition with O(N4 ) colors [10].
We provide an improved 1-sparse decomposition into O(η 2 N 2 ) terms.For convenience of notation, we denote the occupied spin-orbitals for |α by α 1 , . . ., α η , and the occupied spin-orbitals for |β by β 1 , . . ., β η .We also drop the bra-ket notation for the lists of orbitals (Slater determinants); that is, we denote the list of occupied orbitals for the left portion of the graph by α, and the list of occupied orbitals for the right portion of the graph by β.We require both these lists of spin-orbitals to be sorted in ascending order.According to the Slater-Condon rules, the matrix element between two Slater determinants is zero unless the determinants differ by two spin-orbitals or less.Thus, two vertices (Slater determinants) in the Hamiltonian graph are connected if and only if they differ by a single occupied orbital or two occupied orbitals.
In order to obtain the decomposition, for each color (corresponding to one of the resulting 1-sparse matrices) we need to be able to obtain β from α, and vice versa.Using the approach in [28], we take the tensor product of the Hamiltonian with a σ x operator.That is, we perform the simulation under the Hamiltonian σ x ⊗ H, which is bipartite and has the same sparsity as H.The σ x operator acts on the ancilla register that determines whether we are in the left (α) or right (β) partition of the graph.We do this without loss of generality as simulation under H can be recovered from simulation under σ x ⊗ H using the fact that e −i(σx⊗H)t |+ |ψ = |+ e −iHt |ψ [28].
In order for the graph coloring to be suitable for the quantum algorithm, for any given color we must have a procedure for obtaining β given α, and another procedure for obtaining α given β.For this to be a valid graph coloring, the procedure must be reversible, and different colors must not give the same β from α or vice versa.
To explain the decomposition, we will first consider how it works for α and β differing by only a single spin-orbital occupation.We are given a 4-tuple (a, b, , p), where a and b are bits, is a number in the sorted list of occupied orbitals, and p is a number that tells us how many orbitals the starting orbital is shifted by.Our notation here differs slightly from that in Section IV, where i and j were used in place of to represent the positions of the two orbitals which differed: here we will use i and j for a different purpose.To simplify the discussion, we do not perform the addition modulo N , and instead achieve the same effect by allowing p to take positive and negative values.If adding p takes us beyond the list of allowable orbitals, then the matrix element returned is zero, and the list of occupied orbitals is unchanged (corresponding to a diagonal element of the Hamiltonian).We will also use the convention that α 0 = β 0 = 0 and α η+1 = β η+1 = N + 1.These values are not explicitly stored, but rather are dummy values to use in the algorithm when goes beyond the range 1, . . ., η.
The register a tells us whether the is for α or β.To simplify the discussion, when a = 0 we take i = , and when a = 1 we take j = .In either case, we require that β j = α i + p, but in the case a = 0 we are given i and need to work out j, whereas in the case a = 1 we are given j and need to work out i.In particular, for a = 0 we just take α i and add p to it.Then j is the new position in the list β, so β j = α i + p.
The general principle is that, if we are given i for α and need to determine j for β, we require that β j+1 − β j−1 ≥ α i+1 − α i−1 , (i.e. the spacing between orbitals is larger in β than in α).Alternatively, if we were given j for β and needed to determine a corresponding i for α, we would require β j+1 − β j−1 < α i+1 − α i−1 (i.e. the spacing between orbitals is larger in α than in β).If the inequality is not consistent with the value of a (i.e.we are proceeding in the wrong direction), then the matrix element for this term in the decomposition is taken to be zero (in the graph there is no line of that color connecting the nodes).This procedure allows for a unique connection between nodes, without double counting.
The reason for requiring these inequalities is that the list of orbitals with a larger spacing will have less ambiguity in the order of occupied orbitals.To reduce the number of terms in the decomposition, we are only given i or j, but not both, so we either need to be able to determine j from i given β, or i from j given α.When the spacing between the occupied orbitals for β is larger, if we are given β and i there is less ambiguity in determining j.In particular, when β j+1 − β j−1 ≥ α i+1 − α i−1 , there can be at most two values of j that could have come from i, and the bit b is then used to distinguish between them.
There are four different cases that we need to consider.
Next we explain the procedure for each of these cases in detail.In the following we use the terminology "invalid" to indicate that we need to return α = β and a matrix element of zero.
1. Given β and need to determine α; a = 0. We are given β, but is the position in the list of occupied orbitals for α.We do not know which is the β j to subtract p from, so we loop through all values as follows to find a list of candidates for α, α(k) .We define this as a procedure so we can use it later.
procedure FindAlphas k = 0 For j = 1, . . ., η: Subtract p from β j and check that this yields a valid list of orbitals, in that β j − p does not yield an orbital number beyond the desired range, or duplicate another orbital.That is: If ((β j − p ∈ {1, . . ., N }) ∧ (∀j ∈ {1, . . ., η} : β j − p = β j )) ∨ (p = 0) then Sort the list of orbitals to obtain α(k) , and denote by i the new position of β j −p in this list of occupied orbitals.Check that the new value of i corresponds to , and that the spacing condition for a = 0 is satisfied, as follows.
After this procedure there is a list of at most two candidates for α, and k will correspond to how many have been found.Depending on the value of k we perform the following: That is, if we have two possibilities for α, then we use b to choose between them.If there is only one, then we only return that one if b = 0 to avoid obtaining two colors that both link α and β.
2. Given α and need to determine β; a = 0. We are given α, and = i is the position of the occupied orbital in α that is changed.We therefore add p to α i and check that it gives a valid list of orbitals.Not only this, we need to check that we would obtain α if we work backwards from the resulting β.If ((α i + p ∈ {1, . . ., N }) ∧ (∀i ∈ {1, . . ., η} : We sort the new list of occupied orbitals to obtain a candidate for β, denoted β.We next check that the spacing condition for a = 0 is satisfied.If ( βj+1 − βj−1 ≥ α i+1 − α i−1 ) then Perform the procedure FindAlphas to find potential candidates for α that could be obtained from β.There can only be 1 or 2 candidates returned from this procedure.
)) then return β = β else return invalid else return invalid else return invalid 3. Given α and need to determine β; a = 1.This case is closely analogous to the case where we need to determine α from β, but a = 0. We are given α, but is the position in the list of occupied orbitals for β.We do not know which is the α i to add p to, so we loop through all values as follows to find a list of candidates for β, β(k) .We define this as a procedure so we can use it later.
procedure FindBetas k = 0 For i = 1, . . ., η: Add p to α i and check that this yields a valid list of orbitals, in that α i + p does not yield an orbital number beyond the desired range, or duplicate another orbital.That is: Hence α i −1 + p is not in the interval (β j−1 , β j+1 ), and therefore cannot give β j .Therefore there can be no third ambiguous value, in the same way as above for a = 0. Hence the single bit b is again sufficient to distinguish between any ambiguous values, and enables us to determine β given α, p, and j.
We now consider the requirement that the procedure is reversible.In particular, Case 1 needs to be the reverse of Case 2, and Case 3 needs to be the reverse of Case 4. Consider starting from a particular β and using the method in Case 1.We have shown that the procedure FindAlphas in Case 1 can yield at most two potential candidates for α, and then one is chosen via the value of b.For the resulting α, adding p to α i will yield the original set of occupied orbitals β.Moreover, the inequality β j+1 − β j−1 ≥ α i+1 − α i−1 must be satisfied (otherwise Case 1 would yield invalid).
If Case 1 yields β from α, then Case 2 should yield β given α.Case 2 simply adds p to α i (where i is given), which we know should yield β. The method in Case 2 also performs some checks, and outputs invalid if those fail.These checks are: 1.It checks that β is a valid list of orbitals, which must be satisfied because we started with a valid β.

It checks that β
which must be satisfied for Case 1 to yield α instead of invalid.
3. It checks that using Case 1 on β would yield α, which must be satisfied here because we considered initially using Case 1 to obtain α from β.
Thus we see that, if Case 1 yields α from β, then Case 2 must yield β from α.
Going the other way, and starting with α and using Case 2 to find β, a result other than invalid will only be provided if Case 1 would yield α from that β.Thus we immediately know that if Case 2 provides β from α, then Case 1 will provide α from β.This means that the methods for Cases 1 and 2 are the inverses of each other, as required.Via exactly the same reasoning, we can see that the methods in Cases 3 and 4 are the inverses of each other as well.
Next, consider the question of uniqueness.The color will be unique if we can determine the color from a pair α, β.Given α and β, we will see that all the occupied orbitals are identical, except one.Then the occupied orbitals for α and β which are different will be i and j, respectively.We can then immediately set p = β j − α i for the color.We can then compare β j+1 − β j−1 and α i+1 − α i−1 .
If β j+1 − β j−1 ≥ α i+1 − α i−1 then for the color a = 0 and = i.We can then find how many ambiguous values of α there would be if we started with β.If α was obtained uniquely from β, then we would set b = 0 for the color.If there were two ambiguous values of α that could be obtained from β, then if the first was correct we would set b = 0, and if the second were correct then we would set b = 1.
If β j+1 − β j−1 < α i+1 − α i−1 then for the color a = 1 and = j.We can then find how many ambiguous values of β there would be if we started with α.If β was obtained uniquely from α, then we would set b = 0 for the color.If there were two ambiguous values of β that could be obtained from α, then if the first was correct we would set b = 0, and if the second were correct then we would set b = 1.In this way we can see that the pair α, β yields a unique color, and therefore we have a valid coloring.
So far we have considered the case where α and β differ by just one orbital for simplicity.For cases where α and β differ by two orbitals, the procedure is similar.We now need to use the above reasoning to go from α to β through some intermediate list of orbitals χ.That is, we have one set of numbers (a 1 , b 1 , 1 , p) that tells us how to find χ from α, then a second set of numbers (a 2 , b 2 , 2 , q) that tells us how to obtain β from χ.
First, it is easily seen that this procedure is reversible, because the steps for going from α to χ to β are reversible.Second, we need to be able to determine the color from α and β.First, we find the two occupied orbitals for α and β that differ.Call the different occupied orbitals for α i 1 and i 2 , and the different orbitals for β j 1 and j 2 (assume in ascending order so the labels are unique).Then there are four different ways that one could go from α to β, through different intermediate states χ.
To resolve this ambiguity we require that the color is obtained by assuming the first alternative that α i1 → β j1 then α i2 → β j2 .Then α and β yield a unique color.This also requires a slight modification of the technique for obtaining α from β and vice versa.First the color is used to obtain the pair α, β, then it is checked whether the orbitals were mapped as in the first alternative above.If they were not, then invalid is returned.
To enable us to include the matrix elements where α and β differ by a single orbital or no orbitals with a coloring by an 8-tuple γ = (a 1 , b 1 , 1 , p, a 2 , b 2 , 2 , q), we can also allow p = 0 (for only one differing) and p = q = 0 (for α = β).The overall number of terms in the decomposition is then O(η 2 N 2 ).sum in order to give the bound on µ in the lemma.We also give a bound on the absolute value of each term in the Riemann sum using the value of x specified in the first stage.
Third stage: In the final stage of each proof, we bound the total error via the triangle inequality as Our choice of x then ensures that the error is bounded by δ.
To be more specific about the approach in the second stage, we partition D into regions T , and the Riemann sum approximates the integral over each T with the value of the integrand multiplied by the volume of T .The error due to this approximation is bounded by observing the following.Suppose f : IR d → IR is once-differentiable and f max is a bound on its first derivative.If r T is any element of T , we will seek to bound the error of the approximation where vol(T ) is the d-dimensional hypervolume of the set T .The error of this approximation is which can be bounded as follows: where We will choose the points r T in the centers of the regions T , so that where diam(T ) := max r1, r2∈T The Riemann sum approximations we define will then take the form and the error of this approximation is which can be bounded via the triangle inequality as We prove the Lemma by deriving exact formulae for Λ µ,x ( c) in the cases c ≤ x and c > x and then deriving bounds on these formulae that have simpler functional forms.
To derive exact formulae for Λ µ,x ( c), we use the Laplace expansion where R ,m and I ,m refer to the regular and irregular solid spherical harmonic functions, respectively, and r ≥ c.That is to say, and where are the spherical harmonics (see §14.30(i) in [40]), P m are the associated Legendre polynomials, and θ and φ are respectively the polar and azimuthal angles of r.Via Eq. ( 8) of §14.30(ii) in [40], we have and where δ a,b denotes the Kronecker delta.If c ≤ x: Eq. (B3) follows from the fact that (1 + z)e −z < 2e −z/2 for all z > 0.
If c > x: where we use the fact that (z 2 + 2z + 2)e −z < e −z/2 for any z > 0 and the fact that which follows from Eq. (B22).This gives us Eq. (B2) for c > x.In the case that c ≤ x, and, since (z 2 + z)e −z < 4e −z/2 for all z > 0, we have 16π Therefore the bound Eq. (B2) holds for c ≤ x as well.

Proof of Lemma 1
Our proof for Lemma 1 roughly follows the three stages presented in Appendix B 1 a.Here we give the proof in summary form and relegate some of the details to the later subsections.

a. First stage for Lemma 1
The first part of the proof corresponds to the first stage discussed in Appendix B 1 a.We choose The condition Eq. ( 36) ensures that x 0 ≥ x max .We show in Appendix B 2 e that the error due to this truncation can be bounded as b. Green's identity for Lemma 1 The next part of the proof is specific to Lemma 1, and is not one of the general stages outlined in Appendix B 1 a.The integral is given in the form with a second derivative of an orbital, which means that to bound the error we would need additional bounds on the third derivatives of the orbitals.We have not assumed such bounds, so we would like to reexpress the integral in terms of first derivatives before approximating it as a Riemann sum.We have already truncated the domain, though, so we will obtain terms from the boundary of the truncated domain.
We reexpress the integral via Green's first identity, which gives where dV and d S are the volume and oriented surface elements, respectively, and ∂D 0 is the boundary of D 0 .The reason why we do not make this change before truncating the domain is that we have not made any assumptions on the rate of decay of the derivatives of the orbitals.We define We show (in Appendix B 2 f) that c. Second stage for Lemma 1 Next we consider the discretization into a Riemann sum for Lemma 1.We define so that The Riemann sum is then where, for every triple of integers k = (k 1 , k 2 , k 3 ) such that 0 ≤ k 1 , k 2 , k 3 < N 0 , we define and Thus we have partitioned D 0 into µ = N 3 0 equal-sized cubes T (0) k that overlap on sets of measure zero.The expression in Eq. (38) of Lemma 1 then follows immediately.
Each term of R 0 satisfies where the second inequality follows from Eq. (29).Using the value of x 0 in Eq. (B28) in Eq. (B38), each term in the sum has the upper bound on its absolute value (corresponding to Eq. ( 39) in Lemma 1) We show (in Appendix B 2 g) that d. Third stage for Lemma 1 In the final part of the proof of Lemma 1 we show that the total error is properly bounded.By the triangle inequality, we have Riemann . (B41) We therefore arrive at a simple bound on the total error: where To ensure that δ total ≤ δ, we should have We can satisfy this inequality with x 0 given by Eq. (B28).This last step completes our proof.The remainder of this subsection gives the details for some of the steps above.Observe first that where the last inequality follows from the fact that B x ( c i ) ⊂ D 0 .Using this fact together with assumptions 2 and 3 from Section IV C, we have We simplify this expression by expressing r in polar coordinates with center c i .After integrating over the solid angles, we have we have Green for Lemma 1 Using Eq. (B30) and Eq.(B32) we have Then using Eq. ( 28) and Eq. ( 29) gives us We further observe that r − c i ≥ x for all r ∈ ∂D 0 , and the cube with side length 2x has surface area 24x 2 , giving where we have noted 12z 2 e −z < 26e −z/2 for all z > 0.
g. Bounding δ (0) Riemann First we bound the derivative of the integrand.We use the chain rule, the triangle inequality, Eq. ( 29) and Eq. ( 30) to find We have for all k.Using Eq. (B15) and Eq.(B33), and noting that 1/ z ≤ 1/z for any z ∈ IR, we have x 0 x max .(B56)

Proof of Lemma 2
For this proof, the discretization into the Riemann sum will be performed differently depending on whether spinorbital i is considered distant from or nearby to nucleus q.If the nucleus is far from the spin-orbital, the singularity in the integrand is not inside our truncated domain of integration and we need not take special care with it.Otherwise, we can remove the singularity by defining spherical polar coordinates centred at the nucleus.In each case, we select different truncated integration domains and therefore different Riemann sums.
We focus on the centre of spin-orbital i for simplicity; in principle, the centre of spin-orbital j could also be taken into account.We again start by truncating the domain of integration.We select The condition in Eq. (40) ensures that x 1 ≥ x max .We regard the spin-orbital as distant from the nucleus if If so, we use the truncated domain Otherwise, we use D 1,q,singular := B 4x1 R q . (B60) We define δ (1,q,non-singular) trunc (1,q,singular) trunc (1,q) trunc := max δ (1,q,non-singular) trunc , δ (1,q,singular) trunc and show in Appendix B 3 e that δ (1,q) trunc < Now we consider in the discretization of the integral for the case that R q − c i ≥ √ 3x 1 + x max , so orbital i can be regarded as distant from the nucleus.We set and define two different Riemann sums containing µ = N 3 1 terms.We also use this expression for N 1 in the case that the spin-orbital is near the nucleus.Using our value of x 1 in Eq. (B57), Then, noting that µ = N 3 1 is the number of terms in either Riemann sum, we obtain the bound on µ in Eq. (42) of Lemma 2.
We approximate S (1,q) ij (D 1,q,non-singular ) with the sum where, for every triple of integers k = (k 1 , k 2 , k 3 ) such that 0 ≤ k 1 , k 2 , k 3 < N 1 , we define and Thus we have partitioned D 1,q,non-singular into N 3 1 equal-sized cubes T (1,q,non-singular) k that overlap on sets of measure zero.Each term of R 1,q,non-singular satisfies where we have used Eq. ( 27) and the fact that R q − r ≥ x max for every r ∈ D 1,q,non-singular .This upper bound is no greater than Now substituting our value of x 1 from Eq. (B57) shows that no term has absolute value greater than (corresponding to Eq. ( 43) in Lemma 2) We show in Appendix B 3 f that δ (1,q,non-singular) Riemann ) c. Second stage for Lemma 2 with spherical polar coordinates Next we consider discretization of the integral for the case where R q − c i < √ 3x 1 + x max , so orbital i is nearby the nucleus.We express where we define s := ( r − R q )/(4x 1 ) and Here we use θ and φ to refer to the polar and azimuthal angles, respectively, of the vector s.Note that the singularity in the nuclear Coulomb potential has been absorbed into the spherical polar volume form s 2 sin θ ds dθ dφ.For every triple of natural numbers k = (k s , k θ , k φ ) such that 0 ≤ k s , k θ , k φ < N 1 , we define Thus our Riemann sum approximation is where we emphasize that vol T (1,q,singular) k = T (1,q,singular) is not the volume of T (1,q,singular) k when considered as a subset of IR 3 under the usual Euclidean metric.The reason for this discrepancy is that we absorbed the Jacobian introduced by switching from Cartesian to spherical polar coordinates into the definition of f 1 .Thus we are integrating f 1 with respect to the volume form ds dθ dφ, not s 2 sin θ ds dθ dφ.The terms of R 1,q,singular are bounded by where the inequality follows from Eq. (27).Again this expression is upper bounded by Eq. (B71), so substituting our value of x 1 from Eq. (B57) gives the upper bound in Eq. ( 43).We show in Appendix B 3 g that δ (1,q,singular) Riemann We again finish the proof by showing that the total error is bounded by δ.From Eq. (B6), we have δ (1,q,non-singular) total (1,q,singular) total We have given a bound that holds simultaneously for both δ (1,q,non-singular) trunc and δ (1,q,singular) trunc , and we have given a bound for δ (1,q,singular) Riemann that is larger (as a function of x) than our bound for δ (1,q,non-singular) Riemann . We are therefore able to assert that the error of our Riemann sum approximation, no matter which we choose, is always bounded above by where We have found two different upper bounds on the magnitudes of the terms in the Riemann sums given in Eqs.(B70) and (B80).Finally, we note that by substituting our value of x 1 from Eq. (B57), this expression is upper bounded by δ.This last step completes our proof of Lemma 2. The remainder of this subsection gives the details for some of the steps above.
e. Bounding δ (1,q) trunc for Lemma 2 Note first that B x1 ( c i ) ⊂ D 1,q,non-singular and B x1 ( c i ) ⊂ D 1,q,singular .To see the latter, note that we only consider D 1,q,singular in the case that R As we have for any D such that B x ( c i ) ⊂ D, we may compute where the function Λ is as defined in Eq. (B1) and the inequality follows from Eq. (28).By Lemma 4, in the case R q − c i > x 1 we have where the second inequality follows from x 1 ≥ x max .In the case R q − c i ≤ x 1 we can use We can add the bounds to find, in general, that ) f. Bounding δ (1,q,non-singular) Riemann for Lemma 2 Following Appendix B 1 a, we note that vol (D 1,q,non-singular for each k.We can bound the derivative of the integrand using the product rule and the triangle inequality as follows: where the last inequaity follows from Eq. ( 27) and Eq. ( 29), as well as the fact that R q − r ≥ x max for any r ∈ D 1,q,non-singular .From Eq. (B15) and Eq.(B65), and noting 1/ z ≤ 1/z, we have ) g. Bounding δ (1,q,singular) Riemann for Lemma 2 Recalling that we are using a non-standard metric to evaluate the volumes and diameters of sets, we find vol (D 1,q,singular ) = 2π 2 (B96) and diam T (1,q,singular) By Eq. (B15), it remains to find a bound on the derivative of f 1 .Throughout this subsection, we write f max for this bound.
To bound this derivative, we consider the gradient in three different ways.First there is ∇, which is the gradient with respect to the unscaled position coordinates.Second there is ∇ s , which is the gradient with respect to the spherical polar coordinates, but just taking the derivatives with respect to each coordinate.That is, We use this because we are treating the coordinates like they were Euclidean for the discretized integral.Third, there is the usual gradient in spherical polar coordinates, Because we consider s ∈ [0, 1], the components of the gradient ∇ s are upper bounded by the components of the gradient ∇ s .Therefore The restriction on the magnitude of the gradient in Eq. ( 29) holds on the usual gradient in spherical polar coordinates.This means that Using these results, we have f max = ∇ s ϕ * i 4x 1 s + R q ϕ j 4x 1 s + R q s sin θ ≤ ϕ j 4x 1 s + R q s sin θ ∇ s ϕ * i 4x 1 s + R q + ϕ * i 4x 1 s + R q s sin θ ∇ s ϕ j 4x 1 s + R q + ϕ * i 4x 1 s + R q ϕ j 4x 1 s + R q ∇ s [s sin θ] Thus we have the bound We now can give a bound for our approximation to S (1,q) ij (D 1,singular ).Using the above definitions of f max , vol (D 1,singular ) and diam T (1,q,singular) k , we have δ (1,q,singular) Riemann ≤ 1 2 16x 2 1 Z q f max diam T k vol (D 1,singular ) ≤ 16π 2 5π 2 + 1 Using Eq. (B65) and noting 1/ z ≤ 1/z, we have δ (1,q,singular) Riemann where we have used x 1 ≥ x max .

Proof of Lemma 3
As in Appendix B 3, we separate our proof into two cases, depending on whether the singularity of the integrand is relevant or not.If the orbitals i and j are distant, then the singularity is unimportant and we can use rectangular coordinates.If these orbitals are nearby, then we use spherical polar coordinates to eliminate the singularity from the integrand.We do not consider the distance between the orbitals k and in order to simplify the analysis.Note that B x2 ( c i ) × B x2 ( c j ) is a subset of both D 2,non-singular and D 2,singular .The former is immediately apparent.To see the latter, observe that c i − c j < 2 √ 3x 2 + x max implies that the maximum possible value of r 1 − r 2 for any r 1 ∈ B x ( c i ) and r 2 ∈ B x ( c j ) is 2 √ 3x 2 + 2x 2 + x max ≤ (2 √ 3 + 3)x 2 = ζx 2 .Therefore, δ We use the fact that (z 2 + 2z + 2)e −z < 4e −z/2 for any z > 0 to find , it only remains to find a bound on the derivative of the integrand.
To bound the derivative of the integrand, we first find bounds on the gradients of the numerator and the denominator separately.The gradient of the numerator can be bounded using the product rule and triangle inequality, as well as Eq. ( 27) and Eq. ( 29): . (B145) The gradient of the denominator can be computed directly: Again by the product rule and the triangle inequality, The last inequality follows from our assumption that c i − c j > 2 √ 3x 2 + x max , which implies that the distance between C x2 ( c i ) and C x2 ( c j ) is greater than x max .Therefore, δ where we have again used the product rule and the triangle inequality and, in the last inequality, Eq. ( 27).We have also used the bounds on the gradient operator ∇ t in the same was as in Appendix B 3 g.We note that ∇ t (t sin θ) ≤ √ 2, ∇ s η jk x 2 s − ζx 2 t + c i ≤ 2xγ 1 ϕ 2 max /x max (as above) and (b): multiple positions in β could lead to position 2 in α.
a. First stage for Lemma 2 ) b. Second stage for Lemma 2 with Cartesian coordinates ) d. Third stage for Lemma 2

) 3 .
If |α and |β differ by exactly two spin-orbitals such that occupied spin-orbital i in |α is replaced with spinorbital k in |β , and occupied spin-orbital j in |α is replaced with spin-orbital in |β , then Table II lists relevant operators and their gate counts.

TABLE I .
Taylor series simulation parameters and bounds