Exponentially more precise quantum simulation of fermions in second quantization

We introduce novel algorithms for the quantum simulation of fermionic systems which are dramatically more efficient than those based on the Lie–Trotter–Suzuki decomposition. We present the first application of a general technique for simulating Hamiltonian evolution using a truncated Taylor series to obtain logarithmic scaling with the inverse of the desired precision. The key difficulty in applying algorithms for general sparse Hamiltonian simulation to fermionic simulation is that a query, corresponding to computation of an entry of the Hamiltonian, is costly to compute. This means that the gate complexity would be much higher than quantified by the query complexity. We solve this problem with a novel quantum algorithm for on-the-fly computation of integrals that is exponentially faster than classical sampling. While the approaches presented here are readily applicable to a wide class of fermionic models, we focus on quantum chemistry simulation in second quantization, perhaps the most studied application of Hamiltonian simulation. Our central result is an algorithm for simulating an N spin–orbital system that requires  ˜ ( N 5 t ) ?> gates. This approach is exponentially faster in the inverse precision and at least cubically faster in N than all previous approaches to chemistry simulation in the literature.


I. INTRODUCTION
As small, fault-tolerant quantum computers come increasingly close to viability [1][2][3][4] there has been substantial renewed interest in quantum simulating chemistry due to the low qubit requirements and industrial importance of the electronic structure problem.A recent series of papers tried to estimate the resources required to quantum simulate a small but classically intractable molecule [5][6][7][8][9].Although qubit requirements seem modest, initial predictions of the time required were daunting.Using arbitrarily high-order Trotter formulas, the tightest known bound on the gate count of the second quantized, Trotter-based quantum simulation of chemistry is O(N 8 t/ o (1) ) [10,11] 1 , where is the precision required and N is the number of spin-orbitals.However, using significantly more practical Trotter decompositions, the best known gate complexity for this quantum algorithm is O(N 9 t 3 / ) [6].
Fortunately, numerics indicated that the average circuit depth for real molecules may be closer to O(N 6 t 3 / ) [7], or O(Z 3 N 4 t 3 / ) when only trying to simulate ground states, where Z is the largest nuclear charge for the molecule [9].While this improved scaling restores hope that fault-tolerant devices will have an impact on some classically intractable chemistry problems, the Trotter-based quantum simulation of large (e.g.N > 500) molecules still seems prohibitively costly [9,12,13].This limitation would preclude simulations of many important molecular systems, such as those involved in biological nitrogen fixation and high-T c superconductivity [12,13].
The canonical quantum algorithm for quantum chemistry, based on the Trotter-Suzuki decomposition which was first applied for universal quantum simulation in [14,15], was introduced nearly one decade ago [16].This approach was later refined for implementation with a set of universal quantum gates in [17].With the exception of the adiabatic algorithm described in [18] and a classical variational optimization strategy making use of a quantum wavefunction ansatz described in [19], all prior quantum algorithms for chemistry have been based on Trotterization [20][21][22][23][24][25][26][27].
Trotter-Suzuki approaches were also applied to simulation of evolution under sparse Hamiltonians with the entries given by an oracle [28,29].A related problem is the simulation of continuous query algorithms; in 2009, Cleve et al. showed how to achieve such simulation with exponentially fewer discrete queries than Trotterization in terms of 1/ [30].The algorithm of [30] still required a number of ancilla qubits that scaled polynomially in 1/ , but this limitation was overcome in [31] which demonstrated that the ancilla register in [30] could be compressed into exponentially fewer qubits.In [32,33], Berry et al. combined the results of [28][29][30][31] to show exponentially more precise sparse Hamiltonian simulation techniques.A major contribution of [32] was to use oblivi-ous amplitude amplification to make the algorithm from [30,31] deterministic, whereas prior versions had relied on probabilistic measurement of ancilla qubits.An improvement introduced in [33] was to show how to simulate arbitrary Hamiltonians using queries that are not self-inverse (a requirement of the procedure in [32]).We focus on the methodology of [33] which is relatively selfcontained.
The algorithm of [33] approximates the propagator using a Taylor series expansion rather than the Trotter-Suzuki decomposition.By dividing the desired evolution into a number of simulation segments proportional to the Hamiltonian norm, one can truncate the Taylor series at an order which scales logarithmically in the inverse of the desired precision [33].The truncated Taylor series must be expressed as a weighted sum of unitary operators.To simulate the action of this operator, one first initializes the system along with an ancilla register that indexes all terms in the Taylor series sum.The ancilla register is then put in a superposition state with amplitudes proportional to the coefficients of terms in the Taylor series sum.Next, an operator is applied to the system which coherently executes a single term in the Taylor series sum that is selected according to the ancilla register in superposition.Finally, by applying the transpose of the procedure which prepares the ancilla register, one probabilistically simulates evolution under the propagator.The algorithm is made deterministic using an oblivious amplitude amplification procedure from [32].This is the first paper of a two-paper series which applies the algorithm of [33] to quantum chemistry simulation.The algorithms discussed in this paper employ a second quantized encoding of the Hamiltonian, where we dynamically perform the Jordan-Wigner transformation [34,35] on the quantum computer.In the second paper of this series, we use a compressed, first quantized encoding of the wavefunction which requires a number of qubits that scales almost linearly with the number of electrons [36].
In the present paper we develop two new algorithms for the application of the Hamiltonians terms, which we refer to as the "database" algorithm and the "on-thefly" algorithm.In the database algorithm, the ancilla register's superposition state is prepared with amplitudes from a precomputed classical database.In the on-the-fly algorithm, those amplitudes are computed and prepared on-the-fly, in a way that is exponentially more precise than classically possible.

II. OVERVIEW OF RESULTS
The simulation procedure described in [33] assumes the ability to represent the Hamiltonian as a weighted sum of unitaries which can be individually applied to a quantum state.Specifically, we must be able to express the simulation Hamiltonian as where the W γ are complex-valued scalars 2 , the H γ are unitary operators and a mechanism is available for selectively applying the H γ .Using the Jordan-Wigner transformation [34,35] or the Bravyi-Kitaev transformation [37][38][39], the second quantized molecular Hamiltonian can be mapped to a sum of Γ ∈ O(N 4 ) local Hamiltonians.Since these local Hamiltonians are each a tensor product of Pauli operators multiplied by some coefficient, they automatically satisfy the form of Eq. ( 1).We will need a circuit referred to in [33] as select(H) which is queried within the algorithm such that One could construct select(H) by storing all the Pauli strings in a database.However, accessing this data would have time complexity of at least Ω(Γ).Instead, we compute and apply the Pauli strings using O(N ) gates (which can be parallelized to O(1) circuit depth) by dynamically performing the Jordan-Wigner transformation on the quantum hardware.The algorithm of [33] also requires an operator that we refer to as prepare(W ) which applies the mapping where is a normalization factor that will turn out to have significant ramifications for the algorithm complexity.In the first of two algorithms discussed in this paper, we implement prepare(W ) using a database via a sequence of totally controlled rotations at cost O(Γ).Because our first approach uses a database to store classically precomputed values of W γ in order to implement prepare(W ), we refer to the first algorithm as the "database" algorithm.
While we suggest a different strategy in Section III, a database could also be used to construct select(H). 2 The convention of [33] requires that the Wγ are real, non-negative scalars.This treatment remains general as arbitrary phases can be factored into the Hγ .However, we break with that convention and allow the Wγ to take arbitrary complex values.This is done for pedagogical purposes: so that we may separately describe computation of the Hγ and the Wγ for the chemistry Hamiltonian.Consequentially, our Eq.( 41) differs from the analogous equation in [33] by a complex conjugate operator.
That is, a controlled operation is performed which applies H 1 if γ = 1, followed by a controlled operation which performs H 2 if γ = 2, and so forth.This would result in a slightly higher gate count than prepare(W ), because each of the Γ controlled operations must act on O(log N ) qubits even if the Bravyi-Kitaev transformation is used.Nevertheless, this might represent a simpler solution than our construction of select(H) for early experimental implementations.
Our second algorithm involves modifications to the algorithm of [33] which allows us to avoid some of this overhead.We exploit the fact that the chemistry Hamiltonian is easy to express as a special case of Eq. ( 1) in which the coefficients are defined by integrals such as ( Because our approach involves computing integrals onthe-fly, we refer to the second algorithm as the "on-thefly" algorithm.We begin by numerically approximating the integrals as finite Riemann sums such as where z ρ is a point in the integration domain at grid point ρ.Equation ( 6) represents a discretization of the integral in Eq. ( 5) using µ grid points where the domain of the integral, denoted as Z, has been truncated to have total volume V.This truncation is possible because the functions w γ ( z) can be chosen to decay exponentially for the molecular systems usually studied in chemistry.Note that this might not be true for other systems, such as conducting metals.
Our algorithm is effectively able to numerically compute this integral with complexity logarithmic in the number of grid points.It might be thought that this is impossible, because methods of evaluating numerical integrals on quantum computers normally only give a square-root speedup over classical Monte-Carlo algorithms [40].The difference here is that we do not output the value of the integral.The value of the integral is only used to control the weight of a term in the Hamiltonian under which the state evolves.
We construct a circuit which computes the values of w γ ( z ρ ) for the quantum chemistry Hamiltonian with O(N ) gates.We call this circuit sample(w) and define it by its action, where w γ ( z ρ ) is the binary representation of w γ ( z ρ ) using log M qubits.By expanding the W γ in Eq. ( 1) in terms of the easily computed w γ ( z) as in Eq. ( 6), we are able to compute analogous amplitudes to those in Eq. ( 3) in an efficient fashion.Thus, we no longer need the database that characterizes that algorithm.State preparation where the state coefficients can be computed on the quantum computer is more efficient than when they are stored on, and accessed from, a database [41].The worst-case complexity is the square root of the dimension (here it would be O( √ Γµ)), whereas the database state preparation has complexity linear in the dimension (which is O(Γ) for W γ ).Here this would not be an improvement, as we have increased the dimension in the discretization of the integral.
However, the worst-case complexity is only if the amplitudes can take arbitrary values (as this would enable a search algorithm, where the square root of the dimension is optimal [42]).If the amplitudes differ only by phases, then complexity of the state preparation is logarithmic in the dimension.We therefore decompose each w γ ( z) into a sum of terms which differ only by a sign.Then, although the dimension is increased, the complexity of the state preparation is reduced.The decomposition is of the form where In turn, we can express the Hamiltonian as a sum of unitaries weighted by identical amplitudes which differ only by an easily computed sign, As discussed above, the state preparation needed can be performed much more efficiently because the amplitudes are now identical up to a phase.By making a single query to sample(w) and then performing phase-kickback we can implement the operator prepare(w) whose action is is a normalization factor that will turn out to have significant ramifications for the algorithm complexity.Later, we will show that λ ∈ O(N 4 ) and that prepare(w) can be implemented with O(N ) gate count, the cost of a single query to sample(w).
The database algorithm performs evolution under H for time t by making O(Λt) queries to both select(H) and prepare(W ).Because prepare(W ) requires O(Γ) = O(N 4 ) gates, the overall gate count of this approach scales as O(N 4 Λt).To avoid the overhead from prepare(W ), our on-the-fly algorithm exploits a modified version of the truncated Taylor series algorithm which allows for the same evolution by making O(λt) queries to select(H) and prepare(w).As prepare(w) requires O(N ) gates, the gate count for our on-the-fly algorithm scales as O(N λt).
The paper is outlined as follows.In Section III we introduce the second quantized encoding of the wavefunction and construct select(H).In Section IV we review the procedure in [33] to demonstrate our database algorithm which uses select(H) and prepare(W ) to perform a quantum simulation which is exponentially more precise than Trotterization.In Section V we show that one can modify the procedure in [33] to allow for essentially the same result while simultaneously computing the integrals on-the-fly, and show how to implement prepare(w) so as to compute the integrals on-the-fly.In Section VI we bound the errors on the integrals by analyzing the integrands.In Section VII we discuss applications of these results and future research directions.

III. THE HAMILTONIAN ORACLE
The molecular electronic structure Hamiltonian describes electrons interacting in a fixed nuclear potential.Using atomic units in which the electron mass, electron charge, Coulomb's constant and are unity we may write the electronic Hamiltonian as where R i are the nuclei coordinates, Z i are the nuclear charges, and r i are the electron coordinates [43].We 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 [43].Throughout this paper, ϕ i ( r j ) denotes the i th spinorbital occupied by the j th electron which is parameterized in terms of spatial degrees of freedom r j .
In second quantization, antisymmetry is enforced by the operators whereas in first quantization antisymmetry is explicitly in the wavefunction.The second quantized representation of Eq. ( 13) is where the one-electron and two-electron integrals are The operators a † i and a j in Eq. ( 14) obey the fermionic anti-commutation relations, In general, the Hamiltonian in Eq. ( 14) contains O(N 4 ) terms, except in certain limits of very large molecules where use of a local basis and truncation of terms lead to scaling on the order of O(N 2 ) [8].The spatial encoding of Eq. ( 14) requires Θ(N ) qubits, one for each spin-orbital.While fermions are antisymmetric, indistinguishable particles, qubits are distinguishable and have no special symmetries.Accordingly, in order to construct the operator select(H), which applies terms in the second quantized Hamiltonian to qubits as in Eq. ( 2), we will need a mechanism for mapping the fermionic raising and lowering operators in Eq. ( 14) to operators which act on qubits.Operators which raise or lower the state of a qubit are trivial to represent using Pauli matrices, Throughout this paper, σ x j , σ y j and σ z j denote Pauli matrices acting on the j th tensor factor.However, these qubit raising and lowering operators do not satisfy the fermionic anti-commutation relations in Eq. ( 17).To enforce this requirement we can apply either the Jordan-Wigner transformation [34,35] or the Bravyi-Kitaev transformation [37][38][39].
The action of a † j or a j must also introduce a phase to the wavefunction which depends on the parity (i.e.sum modulo 2) of the occupancies of all orbitals with index less than j [38].If f j ∈ {0, 1} denotes the occupancy of orbital j then In general, two pieces of information are needed in order to make sure the qubit encoding of the fermionic state picks up the correct phase: the occupancy of the state and the parity of the occupancy numbers up to j.The Jordan-Wigner transformation maps the occupancy of spin-orbital j directly into the state of qubit j.Thus, in the Jordan-Wigner transformation, occupancy information is stored locally.However, in order to measure the parity of the state in this representation, one needs to measure the occupancies of all orbitals less than j.
Because of this, the Jordan-Wigner transformed operators are N -local, which means that some of the Jordan-Wigner transformed operators are tensor products of up to N Pauli operators.The Jordan-Wigner transformed operators are It would be convenient if we could construct select(H) by applying the Jordan-Wigner transform and acting on the quantum state, one spin-orbital index at a time.For instance, select(H) might control the application of a fermionic operator as follows However, the operators a † j and a j are not unitary because the operators σ + and σ − are not unitary.To correct this problem, we add four qubits to the selection register where each of the four qubits indicates whether the σ x or the ±i σ y part of the σ + and σ − operators should be applied for each of the four fermionic operators in a string such as a † i a † j a k a .For ease of exposition, we define new fermionic operators which are unitary, a † j,q and a j,q where q ∈ {0, 1}, We use these definitions to rewrite the Hamiltonian in Eq. ( 14) so that it is explicitly a weighted sum of unitary Pauli products of the form we require in Eq. ( 1), Inspection reveals that applying the transformations in Eq. ( 27) and Eq. ( 28) to Eq. ( 29) gives the same expression as applying the transformations in Eq. ( 24) and Eq. ( 25) to Eq. (14).By removing factors of 1/2 from both transformation operators and instead placing them in Eq. ( 29), we obtain transformation operators that are always unitary tensor products of Pauli operators.Accordingly, we can implement select(H) in the spirit of Eq. ( 26) by using four additional qubits and the transformation operators in Eq. ( 27) and Eq. ( 28) so that |ijk |q 1 q 2 q 3 q 4 |ψ → |ijk |q 1 q 2 q 3 q 4 a ,q4 |ψ → |ijk |q 1 q 2 q 3 q 4 a k,q3 a ,q4 |ψ → |ijk |q 1 q 2 q 3 q 4 a † j,q2 a k,q3 a ,q4 |ψ → |ijk |q 1 q 2 q 3 q 4 a † i,q1 a † j,q2 a k,q3 a ,q4 |ψ .( A circuit which implements these operators controlled on the selection register is straightforward to construct.Furthermore, the transformation of the terms can be accomplished in O(1) time.Because the Jordan-Wigner transformation is N -local, the number of gates required to actually apply the unitaries in select(H) is O(N ).However, the terms in Eq. ( 27) and Eq. ( 28) are trivial to apply in parallel so that each query takes O(1) time.
Whereas the Jordan-Wigner transformation stores occupancy information locally and parity information Nlocally, the Bravyi-Kitaev transformation stores both parity and occupancy information in a number of qubits that scales as O(log N ) [37][38][39].For this reason, the operators obtained using the Bravyi-Kitaev basis act on at most O(log N ) qubits.It might be possible to apply the Bravyi-Kitaev transformation with O(log N ) gates, which would allow for an implementation of select(H) with O(log N ) instead of O(N ) gates.However, the Bravyi-Kitaev transformation is much more complicated and this would not change the asymptotic scaling of our complete algorithm.The reason for this is because the total cost will depend on the sum of the gate count of select(H) and the gate count of prepare(W ) or prepare(w), and the latter procedures always require at least O(N ) gates.

IV. SIMULATING HAMILTONIAN EVOLUTION
Using the method of [33], Hamiltonian evolution can be simulated with an exponential improvement in precision over Trotter-based methods by approximating the truncated Taylor series of the time evolution operator U = e −iHt .We begin by partitioning the total simulation time t into r segments of t/r.For each of these r segments we perform a Taylor expansion of the propagator and truncate the series at order K, i.e.
where in the second line we have expanded H as in Eq. (1).Notice that if we truncate the series at order  2) applies specified strings of terms, Eq. ( 36) prepare (W ) prepares a superposition of states weighted by coefficients, Eq. ( 3) prepare (β) prepares a superposition of states weighted by coefficients, Eq. ( 37) O (KΓ) W probabilistically performs simulation under H for time t/r, Eq. ( 41) O (KΓ) P projects system onto |0 ⊗J state of selection register, Eq. ( 44) Θ (K log Γ) G amplification operator to implement sum of unitaries, Eq. ( 45) If we wish for the total simulation to have error less than , each segment must have error less than /r.Accordingly, if we set r ≥ H t then our total simulation will have error at most if We now discuss how one can actually implement the truncated evolution operator in Eq. (31).First note that the sum in Eq. ( 31) takes the form where the V j are unitary and U is close to unitary.Our simulation uses an ancillary "selection" register |j = |k |γ 1 • • • |γ K where 0 ≤ k ≤ K and 1 ≤ γ υ ≤ Γ for all υ.We will encode k in unary, which requires Θ(K) qubits, so that |k = |1 k 0 K−k .Additionally, we encode each |γ υ in binary using Θ(log Γ) qubits.While we need K of the |γ υ registers, we note that only k will actually be in use for a given value of |k .The total number of ancilla qubits required for the selection register |j , denoted as J, scales as By making O(K) queries to select(H) from Section IV, we can implement an operator to apply the V j which is referred to in [33] as select(V ), where the V j are defined as in Eq. ( 34).This is equivalent to k applications of select(H), using each of the |γ υ registers, together with k multiplications by −i.In order to obtain k applications of select(H), we may perform a controlled form of select(H) K times, with each successive qubit in the unary representation of k as the control.Given that the gate count for select(H) scales as O(N ), we can implement select(V ) with O(N K) gates.Applying the Pauli strings in parallel leads to circuit depths of O(1) and O(K), respectively.Table I lists relevant parameters along with their bounds in our database algorithm.Table II lists relevant operators and their gate counts in our database algorithm.We will also need an operator that we refer to as prepare(β), which initializes a state, where s is a normalization factor.
To implement prepare(β) we first prepare the state Using the convention that R y (θ) ≡ exp[−i θ σ y /2], we apply R y (θ 1 ) to the first qubit of the unary encoding for FIG. 1.The circuit for prepare(β) as described in Eq. ( 37).An expression for θ k is given in Eq. (39).prepare(W ) is implemented using a precomputed classical database.
k followed by R y (θ k ) to the kth qubit controlled on qubit k − 1 for all k ∈ [2, K] sequentially, where (39) To each of the K remaining components of the selection register |γ 1 • • • |γ K , we apply prepare(W ) once.In principle, we only need to perform prepare(W ) k times, because the registers past k are not used.However, it is simpler to perform prepare(W ) K times, because it does not require control on k.
Using results from [44], prepare(W ) can be implemented with O(Γ) gates by using a classically precomputed database of the Γ molecular integrals.The gate count for prepare(β) thus scales as O(KΓ).However, this construction is naturally parallelized to depth O(K) + O(Γ).A circuit implementing prepare(β) is shown in Figure 1.
The general strategy for implementing the truncated evolution operator in Eq. ( 34) becomes clear if we consider what happens to state |ψ when we apply prepare(β) followed by the operator select(V ), The similarity of this state to the state U |ψ motivates the operator where |Φ is a state with the ancilla qubits orthogonal to |0 ⊗J .Note that in [33], the authors use the convention that all W γ are positive and phases are incorporated into the operators H γ .Since we depart from that convention for reasons described in Section II, the second application of prepare(β) in Eq. ( 41) is the transpose of the first application, in contrast to [33] where the conjugate transpose is used instead.At this point, we choose the number of segments to be Since Λ ≥ H , our choice of K in Eq. ( 33) remains valid.The additional factor of 1/ ln(2) is included to satisfy a requirement of oblivious amplitude amplification as described in [32] so that We now define a projection operator P onto the target state, which has support on the empty ancilla register, We also define the amplification operator, where R = 1 1 − 2P is the reflection operator.With these definitions, we follow the procedure in [33] which uses the oblivious amplitude amplification procedure of [32] to deterministically execute the intended unitary.We use G in conjunction with P to amplify the target state, Recalling the definition of U r in Eq. ( 31), our choice of K in Eq. ( 33) and r = Λt/ ln 2 imply that so that the total error from applying oblivious amplitude amplification to all the segments will again be order .The gate count of the entire algorithm is thus r times the cost of implementing select(V ) plus the cost of implementing prepare(β).Though we implement select(V ) with O(N K) gates, our brute-force construction of prepare(W ) led to a gate count for prepare(β) which scales as O(KΓ).Thus, the total gate count of our database algorithm scales as While this bound suggests an exponentially more precise algorithm than those based on Trotterization, in the remainder of our paper we discuss an even more efficient algorithm with improved dependence on N .

V. EVOLUTION UNDER INTEGRAL HAMILTONIANS
In Section IV we analyzed the database algorithm for quantum simulating chemistry Hamiltonians in a manner that is exponentially more precise than Trotterization.The most costly part of that procedure is the implementation of prepare(W ) as in Eq. ( 3), which prepares a superposition state with amplitudes that are given by integrals over spin-orbitals as in Eq. ( 15) and Eq. ( 16).Instead of classically precomputing these integrals and implementing prepare(W ) with a database, the strategy we introduce is to numerically sample the integrals on-the-fly using the quantum computer.Because of this, we call this the "on-the-fly" algorithm.To accomplish this, we discretize the integrals as sums and design a circuit which returns the integrand of these integrals at particular domain points.The motivation for approximating integrals as sums comes from a direct analogy between the discretization of time in the Taylor series approach for simulating time-dependent Hamiltonians [33] and the discretization of space in Riemann integration.
In [33], the time-ordered exponential is approximated by a Taylor series up to order K, and the integrals are then discretized as follows on each segment, Now let us suppose that our Hamiltonian does not change in time, but instead that the Hamiltonian itself is given as a definite integral over the spatial region Z so that The second quantized Hamiltonian given by Eq. ( 29) is similar to this, except it includes both terms h ij , which are integrals over one spatial coordinate, and h ijk , which are integrals over two spatial coordinates.While those integrands are also defined over all space, the integrands decay exponentially so we can approximate them as definite integrals over the finite region Z, having volume V.
Then we can approximate the integral by where z ρ is a point in the domain Z at the ρ th grid point.
As in Section IV, we begin by dividing t into r seg-ments.We turn our attention to a single segment, where in the second line we have performed a Taylor expansion of the propagator and truncated at order K.In Eq. ( 52), the bolded symbol z indicates a vector of vectors.Like before, if r ≥ H t then the relationship between K and is given by Eq. ( 33).To approximate the integral, we divide it into µ regions of volume V/µ.We now have For the second quantized Hamiltonian, the W γ in Eq. ( 1) are integrals over scalar functions w γ ( z) as in Eq. ( 5).Using this property it is clear that and the segment U r can be expressed as The second quantized Hamiltonian given in Eq. ( 29) is a sum of terms which are integrals over one spatial coordinate and terms which are integrals over two spatial coordinates.This case is easily accounted for by taking z to be a vector including both spatial coordinates, and V to be the product of the volumes for the two coordinates.One can take the terms with the integral over the single spatial coordinate to also be integrated over the second spatial coordinate, and divided by the volume of integration of that second coordinate to give the correct normalization.We may now proceed with the truncated Taylor series simulation as in Section IV.Whereas our database algorithm required prepare(W ) to create a superposition of states weighted by the W γ , as in Eq. ( 3), our on-the-fly algorithm will need to create a superposition of states weighted by the scalar integrands w γ ( z ρ ).
We now show a method for dynamically decomposing each w γ ( z) into a sum of terms which differ only by a sign as in Eq. ( 8) so that each where ζ is the precision of the decomposition.First, we round each entry of w γ ( z) to the nearest multiple of 2 ζ and decompose it in such a fashion that To construct prepare(w), as in Eq. ( 11), we perform logic on the output of sample(w) (introduced in Eq. ( 7), and detailed in Section VI) which determines whether the w γ,m ( z) for a given value of m should be 1 or −1.Since the superposition in Eq. ( 11) must be weighted by the square root of this coefficient, we need to prepare states that either do or do not have a phase factor i so where | = |γ |m |ρ and | w γ ( z ρ ) was obtained from sample(w).The phase factor can be obtained using phase-kickback in the usual way.Then we apply sample(w) a second time to erase the register | w ( z ρ ) .
A single query to this circuit allows for the construction of prepare(w) with the same complexity as sample(w).
Before explaining the integrand circuit we briefly comment on the additional resources required for the Taylor series simulation under a discretized, position-dependent, integrand Hamiltonian.As in the constant Hamiltonian case, we need one register with Θ(K) qubits to encode |k and K registers of Θ(log Γ) qubits to encode |γ 1 • • • |γ k .However, as in the explanation below Eq. ( 11) we also need extra ancilla qubits to store the value of m, the grid point registers, as well as the value registers which are used by the integrand oracle sample(w) in Eq. ( 7).This represents an additional ancilla overhead of Θ(K log(M µ)).
The sources of simulation error are also similar to the constant Hamiltonian case.As we show in Section VI, we can approximate the integrals with discrete sums to precision at a cost that is logarithmic in 1/ .The error due to the discrete sum is controlled by the choice of µ, which we need to select so that the resulting error per segment is less than /r.The most costly integrals (due to the size of their domain) will be the two-electron integrals in Eq. ( 16) which have integrands of the form where x and y represent the three spatial degrees of freedom of two separate electrons.In Section VI, we bound the cost to the quantum algorithm of estimating the corresponding integrals.

VI. THE INTEGRAND ORACLE
In Section V, we showed how one can implement the truncated Taylor series simulation technique by replacing a superposition state having amplitudes given by integrals with a superposition state having amplitudes given by their integrands, as well as a way of decomposing those integrands.We begin this section by constructing a circuit which allows us to sample from the integrands as in Eq. (7).First, we will need a circuit which computes values of the N spin-orbital basis functions ϕ 1 ( z ρ ) to ϕ N ( z ρ ) at z ρ , a real-valued position vector at grid point ρ.The action of each these oracles is where ϕ j ( z ρ ) represents the binary expansion of ϕ j ( z ρ ) using log M qubits.We will need N different circuits of this form, one for each basis function ϕ 1 ( z) to ϕ N ( z).
Usually, the molecular spin-orbital basis functions are represented as sums of Gaussians multiplied by polynomials [43].In that case, the circuit Q ϕj would be a reversible classical circuit that evaluates and sums the Gaussians associated with ϕ j ( z).For example, in the STO-nG basis set, each orbital is a linear combination of n Gaussian basis functions [43].In general, Gaussian functions may be evaluated classically with complexity that is at most polylogarithmic in 1/ [45].The use of Gaussians is a historical precedent used because those functions are simple to integrate on a classical computer.However, the use of a Gaussian basis is not necessarily an optimal strategy for a quantum implementation.We leave open the problem of determining an optimal representation of molecular orbital basis functions for evaluation on a quantum computer and develop a strategy based on the model chemistries used in classical treatments of electronic structure.Next, we combine N different Q ϕj circuits, one for each ϕ j ( z), to construct a circuit Q which allows us to apply any of the N basis functions.This circuit will have depth O(N polylog(N t/ )) and may be constructed as the block diagonal operator Thus, Q is a sequence of Q ϕj circuits with the spinorbital selection completely controlled on a register encoding the basis function index j.There will a factor of log(N t/ ) because the controlled operations need to access O(log N ) qubits for j, as well as O(log(N t/ )) qubits storing the position z.In addition, the circuit needs to perform analytic operations (e.g.calculating exponentials for STO-nG), which will contribute an additional factor polynomial in log(N t/ ).We now discuss how one can use Q to compute the twoelectron integrands in Eq. (58).To avoid singularities that would occur when two electrons occupy the same point in space, we change variables in Eq. ( 58) so that ξ = x − y.With this substitution, the integral becomes Expressing ξ in spherical polar coordinates, with ξ ≡ | ξ|, we have We define the maximum value of any spin-orbital function as ϕ max and the maximum value of its derivative in any direction as ϕ max .In addition, we truncate the integral at a finite distance x max .Now assume that we discretize x in intervals of size δx along each degree of freedom.We can take the maximum value of ξ to be x max , and choose δξ = δx, δθ = δφ = δx/x max .The primary contribution to the complexity is in terms of the number of segments.The maximum value in the integrand of Eq. ( 62) is upper bounded by x max ϕ 4 max .When discretizing the integral, each term in the sum is no larger than O(x max ϕ 4 max δx 4 (δx/x max ) 2 ) = O(ϕ 4 max δx 6 /x max ) and there are O((x max /δx) 6 ) terms.Multiplying these together gives us the contribution of the integral to the scaling of our on-the-fly algorithm, which corresponds to the factor of V max z,γ |w γ ( z)| in Eq. (12).But how do ϕ max and x max scale with N ?The maximum values ϕ max are predetermined by the model chemistry, and hence are independent of N .Determining the appropriate value of x max is a little more complicated.
Because the Hamiltonian is a sum of O(N 4 ) of the integrals, each integral should be approximated within error O( /(N 4 t)) to ensure that the final error is bounded by .Since the functions ϕ j ( z) can be chosen to decay exponentially, x max can be chosen logarithmically in the allowable error .The quantum chemistry problem is always defined within a particular basis, specified as part of a model chemistry [43].The model chemistry prescribes how many spin-orbitals, how many basis functions, and what type of basis functions should be associated with each atom in a molecule.This includes a specification for parameters of the basis functions which impose a particular maximum value ϕ max , as well as a cutoff distance beyond which each ϕ j ( z) is negligibly small.However, the space between basis functions on different atoms must grow as the cube root of N , because the molecular volume will grow as O(N ).This would imply that the value of x max needed scales as Nevertheless, each individual orbital ϕ j ( z) is nonnegligible on a region that grows only as O(log N ) for a given model chemistry.It is therefore advantageous to modify the grid used for the integral so it only includes points where one of the associated orbitals is nonnegligible.This can be performed at unit cost if the center of each spin-orbital function is provided in an additional register when querying the circuit Q.As above, the region where the orbital can be regarded as nonnegligible can be chosen logarithmically in , to ensure that the overall error in the simulation is within error .
To be more specific, the coordinates for x can be chosen to be centered around the center of orbital ϕ i , with the components of x only going up to a maximum value scaling as For ξ, we only wish to take values such that ϕ j ( x − ξ) are non-negligible.Here it should be noted that the spherical polar coordinates are only advantageous if we are in a region where ξ is near zero, where the Cartesian coordinates would have a divergence.In regions where ξ is large, the extra factor of ξ for the integral in spherical polar coordinates increases the complexity.Therefore, if the minimum value of | ξ| such that ϕ j ( x − ξ) is non-negligible is O (log (N t/ )), then the maximum value of | ξ| such that ϕ j ( x− ξ) is non-negligible will also be O (log (N t/ )).Therefore we can use spherical polar coordinates, and obtain scaling as in Eq. ( 63) with x max as in Eq. ( 65).On the other hand, if the minimum value of | ξ| such that ϕ j ( x − ξ) is non-negligible is Ω (log (N t/ )), then we can use Cartesian coordinates, and the division by | ξ| can only lower the complexity.We obtain a contribution to the complexity scaling as O ϕ 4 max x 3 max with x max as in Eq. (65).Here the power of x max is 3 rather than 5, because we divide instead of multiplying by | ξ| as we did spherical polar coordinates.
Next we consider the grid size needed to appropriately bound the error in the discretized integration.The analysis in the case where Cartesian coordinates are used is relatively straightforward.Considering a single block in 6 dimensions with sides of length δx, the value of the integrand can only vary by the maximum derivative of the integrand times δx (up to a constant factor).The error for the approximation of the integral on this cube is therefore that maximum derivative times δx 7 .Then the number of these blocks in the integral is O((x max /δx) 6 ), giving an overall error scaling as x 6 max δx times the maximum derivative of the integrand.
The maximum derivative of the integrand can be bounded in the following way.For the derivative with respect to any component of x, we obtain the derivative of the integrand scaling as where we have used the fact that we are only using Cartesian coordinates for | ξ| = Ω(x max ).For the derivative of the integrand with respect to any component of ξ in the numerator of the integrand, the same scaling is obtained.
For derivatives with respect to components of ξ in the denominator, the scaling is Overall, we therefore bound the error when discretizing in Cartesian coordinates as The analysis for spherical polar coordinates is a little more subtle, but it is largely equivalent if we scale the angular variables.It is convenient to define scaled angular variables θ ≡ x max θ, φ ≡ x max φ. (69) Then the discretization lengths for all variables are δx.
The volume of each block in the discretization is again δx 6 , and there are O((x max /δx) 6 ) blocks.The total error is again therefore the maximum derivative of the integrand multiplied by x 6 max δx.The derivative of the integrand with respect to any component of x is again given by Eq. (66).Multiplication by ξ gives a factor O(x max ), but the change of variables to θ and φ gives division by a factor of x 2 max .The derivative of the integrand with respect to ξ, θ or φ in any of the spin orbitals again gives a factor scaling as in Eq. (66).The derivative of the integrand with respect to ξ or θ in ξ sin(θ /x max ) scales as in Eq. (67).
As a result, regardless of whether Cartesian coordinates are used or spherical polar coordinates, the error due to discretization is bounded as in (68).Thus, to achieve error in the integral no larger than O( /(N 4 t)), we require that The total number of terms in the sum then scales as This is quite large, but because we only need to use a number of qubits that is the logarithm of this, it only contributes a logarithmic factor to the complexity.Because the logarithm scales as O(log(N t/ )), it contributes this factor to the complexity of sample(w).
Given Q, computing the integrand in Eq. ( 62) is straightforward.We need to call Q four times on registers that contain x and ξ.We then multiply those outputs together with the value of ξ sin θ.Next we consider how to construct a circuit for the one-electron integrals in Eq. (15).First, one constructs N additional circuits similar to the ones in Eq. (59) that return ∇ 2 ϕ j ( z) as opposed to ϕ j ( z).These oracles are incorporated into a one-electron version of Q which is called along with a routine to compute the nuclear Coulomb interactions.The one-electron integrals have singularities at the positions of the nuclei.Similar to the two-electron integrals, these singularities can be avoided by using spherical polar coordinates.Each term in the sum over the nuclei should use spherical polar coordinates centered at that nucleus.Selection between the one-electron and two-electron routines is specified by |γ .Putting this together, we can implement sample(w) as in Eq. (7) with O(N log N ) gates, and, as discussed in Section V, prepare(w) has the same complexity.
While the O(N ) gate count of prepare(w) is much less than the O(N 4 ) gate count of prepare(W ), our on-thefly algorithm requires more segments than the database algorithm.Whereas our database algorithm requires r = Λt/ ln(2) segments where Λ is the normalization in Eq. ( 3), our on-the-fly algorithm requires r = λt/ ln(2) segments where λ ∈ Θ(Γϕ 4 max x 5 max ) is the normalization in Eq. (11), which is accounted for in Eq. ( 12) and Eq. ( 63).Thus, by performing the algorithm in Section IV using prepare(w) instead of prepare(W ) and taking r = λt/ ln(2), we see that our on-the-fly algorithm scales as Recall that the O notation indicates that logarithmic factors have been omitted.The full scaling includes a power of the logarithm of 1/ .

VII. DISCUSSION
We have introduced two novel algorithms for the simulation of molecular systems based primarily on the results of [33].Our database algorithm involves using a database to store the molecular integrals; its gate count scales as O(N 8 t).Our on-the-fly algorithm involves computing those integrals on-the-fly; its gate count scales as O(N 5 t).Both represent an exponential improvement in precision over Trotter-based methods which scale as O(N 9 t 3 / ) when using reasonably low-order decompositions, and over all other approaches to date.Specifically, our database algorithm scales like O(N 4 Λt) where we have used the bound Λ ∈ O(N 4 ).
However, we believe this bound is very loose.As discussed in [8,43], the use of local basis sets leads to a number of two-electron integrals that scales as O(N 2 ) in certain limits of large molecules.Accordingly, the true scaling of the database algorithm is likely to be closer to O(N 6 t).It also seems possible that our integration scheme is suboptimal; it is possible that it can be improved by taking account of smaller values of h ijk .
Our asymptotic analysis suggests that these algorithms will allow for the quantum simulation of molecular systems larger than would be possible using Trotter-based methods.However, numerical simulations will be crucial in order to further optimize these algorithms and better understand their scaling properties.Just as recent work showed significantly more efficient implementations of the original Trotterized quantum chemistry algorithm [5][6][7][8][9], we believe the implementations discussed here are far from optimal.Furthermore, just as was observed for Trotterized quantum chemistry [7,9], we believe our simulations might scale much better when only trying to simulate ground states of real molecules.In light of this, numerical simulations may indicate that the scaling for real molecules is much less than our bounds predict.
Using the scaling in Eq. (65), we can bound λ asλ ∈ O Γϕ 4 max x 5 max ∈ O N 4 [log (N t/ )] 5 .(73)so that the overall gate count of the on-the-fly algorithm scales asO N 5 Kt = O N 5 t .

TABLE I .
Database algorithm parameters and bounds