Exploiting fermion number in factorized decompositions of the electronic structure Hamiltonian

Achieving an accurate description of fermionic systems typically requires considerably many more orbitals than fermions. Previous resource analyses of quantum chemistry simulation often failed to exploit this low fermionic number information in the implementation of Trotter-based approaches and overestimated the quantum-computer runtime as a result. They also depended on numerical procedures that are computationally too expensive to scale up to large systems of practical interest. Here we propose techniques that solve both problems by using various factorized decompositions of the electronic structure Hamiltonian. We showcase our techniques for the uniform electron gas, finding substantial (over 100x) improvements in Trotter error for low-filling fraction and pushing to much higher numbers of orbitals than is possible with existing methods. Finally, we calculate the T-count to perform phase-estimation on Jellium. In the low-filling regime, we observe improvements in gate complexity of over 10x compared to the best Trotter-based approach reported to date. We also report gate counts competitive with qubitization-based approaches for Wigner-Seitz values of physical interest.


I. Introduction
There is considerable interest in whether quantum computers -both those available at present, and those under development -can be used to solve problems of scientific and commercial importance. This is particularly evident in the field of quantum simulation of chemical systems -for recent reviews of progress in this area, we direct the reader to Refs. [1][2][3]. Several algorithms have been developed to obtain the eigenstates of chemical systems. These include variational quantum algorithms [4,5] that aim to maximise the limited coherence times of currently available hardware. However, this comes at the cost of introducing heuristic aspects, making it difficult to obtain rigorous performance guarantees. In contrast, approaches based on quantum phase estimation [6,7] provide a route to calculate eigenstates to within a specifiable error, assuming only that we can efficiently prepare approximate eigenstates with sufficiently high overlap with the true eigenstates.
The resources we allocate to a fault tolerant quantum computation will depend on our ability to bound errors in the algorithm; the tighter our error estimates, the fewer resources we will require. Several previous works have estimated the resources required for phase estimation based on product-formula decompositions (also known as Trotterization) [8][9][10][11][12][13][14]. It was recently shown by Su, Huang and Campbell [15] that knowledge about the number of fermions present in a chemical system can be exploited to improve the asymptotic performance of Trotterization. That work introduced an error metric, termed the fermionic seminorm, to bound the Trotter error. This approach uses knowledge of the number of fermions in the system to offset the dependence of the error on the number of orbitals. This effect may be particularly important for applications to chemical systems in realistically sized basis sets, which will need to be large in order to accurately resolve dynamic correlation in the wavefunction. The Su-Huang-Campbell (SHC) bound aimed to find an analytic bound with the best asymptotic complexity. Here we present complementary work that also uses the fermionic seminorm with the goal of developing techniques for numerically obtaining bounds with best performance in practice.
In this work, we introduce three factorized decompositions of the electronic structure Hamiltonian in a plane wave dual basis, and use these in conjunction with the fermionic seminorm to obtain tighter Trotter error bounds in practice. Our approach is inspired by prior work using low-rank decompositions to reduce the number of terms in a Hamiltonian and thereby reduce the gate complexity of quantum algorithms [13,[16][17][18]. However, our use of factorized decompositions is purely computational and optimised for tightest error bounds, with no corresponding change in the execution of the quantum algorithm. A high-level overview of our approach can be found in Section III.
Each of our three factorized decompositions exhibits its own advantage. The spectral decomposition is generally applicable and extends beyond the plane wave dual basis. The cosine decomposition best exploits fermion number information and so performs the most effectively in the low-filling fraction regime. The Cholesky decomposition has the smallest constant factor overhead and so performs best in the medium and half-filling regimes. We discuss these decompositions in detail and compare the resulting Trotter error bounds in Section IV.
These performance observations are supported by numerical results in Section VI, obtained by applying our approach to the uniform electron gas (Jellium) introduced in Section V. In these numerics, we also benchmark against three prior art bounds: the analytic SHC bound described earlier [15]; the fermionic commutator approach used by Kivlichan et al [14]; and a similar Pauli commutator approach where there is anecdotal evidence of good performance (see App. A of Ref. [19]). We report a substantial classical runtime advantage for the calculation of our bounds. The fermionic and Pauli commutator approaches became intractable to calculate at larger spin-orbital number N , so could not be computed beyond N ∼ 200, without access to > 100 GB of RAM. In contrast, it took fewer than 6 hours (using a 3.6GHz c5.2xlarge EC2 instance on AWS) to calculate our new bounds on a 512 spin-orbital instance, using < 16 GB of RAM. One target problem for Trotter methods has been for phase estimation of the ground state energy of the uniform electron gas [14]. Using our improved Trotter error bounds for Jellium, we calculate the T -count for this problem and demonstrate the expected improvements in runtime. We also compare our gate counts to those obtained using qubitization [20], and find comparable results in some parameter regimes of interest.
We present mathematical preliminaries in Section II that are necessary to understand our factorized decompositions and their numerical implementations. We conclude the paper in Section VII with a brief summary of our contributions and a collection of avenues for future work.

A. Fermionic systems and seminorm
The electronic structure Hamiltonian is a widely used model for molecular and material systems where the positions of the nuclei are considered fixed. In an arbitrary basis of N electronic spin-orbitals, the Hamiltonian can be written as where a s is the fermionic annihilation operator on spin-orbital s and the coefficients h pq and h pqrs are defined by integrals over the basis functions [21]. Using the plane wave dual basis given by [22], the number of terms is reduced from O(N 4 ) to O(N 2 ) with the simple form which is split into the electron kinetic, electron-nuclei, and electron-electron terms, respectively. The coefficients T pq , U p , and V pq are defined by integrals over the basis functions, as discussed in Section IV. When simulating time evolution under a Hamiltonian (such as those given above), the error is typically quantified using the spectral-norm distance between the time evolution operator, and the quantum circuit used to approximate it. However, it is possible to use knowledge about the initial state to improve the error bound. In Ref. [15] the fermionic seminorm of an operator X was defined as the maximum transition amplitude of the operator between two states in the η-electron subspace We say an operator X is number preserving if X acting on an η-electron state yields some other η-electron state. It was shown in Ref. [15] that the fermionic seminorm has similar properties to well-known existing norms. For number preserving operators X, Y , we will make use of the following properties: • ||X · Y || η ≤ ||X|| η · ||Y || η (Hölder inequality) • ||λX|| η = |λ| · ||X|| η (for λ ∈ C) • ||X † || η = ||X|| η • ||U XW || η = ||X|| η (for U, X, W number preserving, U, W unitary) We remark that it is a seminorm rather than a norm because it can evaluate to zero for some non-zero operators. For example, for a system with a single fermion, we have n p n q η=1 = 0 (for p = q), but n p n q is a nonzero operator.
It has been shown [14,23] that this approximation has an error given by where W 1 and W 2 are defined as where ||...|| denotes the operator norm (also known as the spectral norm -i.e. the largest singular value of the operator).
In practice, it can be difficult to get a tight value for W 1/2 because of the complexity in evaluating the operator norm of a high-dimensional operator such as [H c , [H b , H a ]]. As such, in aid of numerical expediency, a further relaxation is often made. Each nested commutator is expanded in terms of operators P j with known operator norm ||P j || = 1 so that c>a b>a and then one can bound c>a b>a Common choices include choosing P j as tensor products of Pauli operators, or as fermionic excitation operators (e.g.
. Throughout, we refer to bounds using these relaxations as the Pauli commutator bound and Fermionic commutator bound, respectively. For example, Ref. [14] used the Fermionic commutator bound to estimate the resources for phase estimation in the plane wave dual basis. In our numerical examples, we will benchmark against these prior art bounds.
Ref. [15] showed the commutator bounds can be tightened in the special case where H is a fermionic Hamiltonian and every H j in the Trotter decomposition (Eq. (5)) is number-preserving, so that for second-order Trotter where the operator norm has been replaced by the tighter fermionic semi-norm. Ref. [15] further considered Hamiltonians in the plane wave dual basis (recall Eq. (2)) and a Trotterization where the Hamiltonian is considered as containing two terms; H t = pq T pq a † p a q and H v = pqV pq n p n q . Ref. [15] derived bounds for arbitrary order product formulae, with the second-order result where || . . . || max is the max-norm that represents the largest matrix element in absolute value. A key observation is that the bound depends on η and so captures the expected dependence on the fermion number. The big-O of this result hides the constant factors that are needed for numerical comparisons. For the case of the plane wave dual basis, Eq. (12) reduces to a sum of two terms. In Appendix A we have evaluated these terms, which are given by [ We refer to this as the 'SHC bound' throughout.

III. Improved fermionic seminorm bounds
In this work, we make particular use of the properties of free-fermionic Hamiltonians H(A) := i,j A ij a † i a j . We refer to A as the coefficient matrix of the free-fermionic Hamiltonian. A free-fermionic Hamiltonian can be efficiently diagonalised by diagonalising its coefficient matrix. We can then calculate the fermionic seminorm of a free-fermionic Hamiltonian as where V is a unitary matrix that diagonalises the free-fermionic Hamiltonian, and λ k are the eigenvalues of the coefficient matrix A. This expression can be evaluated using Eq. (3) to give Here we have defined another seminorm · η which takes a coefficient matrix A as its argument. We call this the reduced fermionic semi-norm as the argument is a smaller N -by-N matrix A, rather than the large operator H(A) that is represented by a 2 N -by-2 N matrix. The result of Eq. (17) tells us that for free-fermionic operators H(A) the problem of evaluating the fermionic semi-norm simplifies to the easier problem of evaluating the reduced fermionic seminorm. Evaluating the reduced fermionic semi-norm takes the set λ(A) of eigenvalues of A and finds the subset S ⊂ λ(A) with η elements and largest sum in absolute value. If A is Hermitian this is further simplified, as we can consider the sum of the η largest eigenvalues, and the sum of the η most-negative eigenvalues, and choose the larger absolute value. Therefore, Eq. (17) can be efficiently computed for Hermitian A.
Another useful property involves the commutator of two free-fermionic Hamiltonians itself a free-fermionic Hamiltonian. This has previously been noted and made use of in the context of quantum simulation in Refs. [10,24]. Motivated by these properties of free-fermionic Hamiltonians, we consider decomposing the Hamiltonian as where A, X l , Y l are N × N coefficient matrices, and N is the number of spin-orbitals considered. Decompositions of this form have been considered in the context of quantum computing in Refs. [10,13], where they were obtained by Approach Memory Runtime Exploits η N/2 η = N/2 fermion # rank rank Fermionic commutator [14] O  Figure 1, for a uniform electron gas system with 200 spin-orbitals, and varying electron number. The memory requirement of the Pauli commutator approach was too severe to carry out this second-order Trotter calculation, though we will later present first-order Trotter results for this approach.
eigen/Cholesky decompositions of the tensor h pqrs . This yields the Hamiltonian in a 'single factorised' form [13,16,19]. For example, if we consider the electronic structure Hamiltonian in a Gaussian orbital basis set (described by Eq. (1)) we can apply a spectral decomposition of the tensor h pqrs to write the Hamiltonian in the form (see Appendix B) where L denotes the number of terms in the spectral decomposition, and λ are the corresponding eigenvalues. We consider Trotter decompositions with each term H j in Eq. (4) or Eq. (5) corresponding to some subset of terms from Eq. (20). We show in Appendix B that we can bound the first-order Trotter error with the commutator bound To obtain this form, we made use of commutator identities such as A similar approach was attempted in Ref. [10], however, that work did not explicitly make use of information about the number of electrons in the system, and therefore the result is not tight in the low-filling regime. Depending on the form of h pqrs , other decompositions may be possible. In the following section, we present three decompositions of the plane wave dual basis Hamiltonian, motivated by its simple form, and the analytic expressions available for the Hamiltonian coefficients in this basis. We summarise the main features of these decompositions in Table I. Each of these three decompositions has a particular benefit; the spectral decomposition is the extension of the approach discussed above (and in Appendix B) to the plane wave dual basis, and so is generally applicable to any orbital basis. The Cholesky decomposition performs best in the half-filling regime (η = N/2), while the cosine decomposition performs best in the low-filling regime (η N/2). All of these bounds are more efficient to compute than the existing fermionic and Pauli commutator bounds.

IV. Plane wave dual basis decompositions and Trotter error bounds
The plane wave dual basis electronic structure Hamiltonian given by Eq. (2) describes a system with η electrons in a simulation box of size Ω ∝ L d , where d is the dimensionality of the system, and L is the number of grid points along each side of the box. The spin-orbitals are obtained from a discrete Fourier transform of plane waves. These plane waves are defined by where N is the number of spin-orbital basis functions used, and ν enumerates the N/2 possible distinct momentum vectors of the system. Note that if L is odd, the interval of ν is closed, rather than half-open. The plane wave dual basis resembles a smooth approximation to a grid of delta functions. The coefficients in Eq. (2) are given by [22] T pq = δ σp,σq Here, r p is the position of the orbital centroid corresponding to spatial-orbital p σ p is the spin of the pth spin-orbital (here, we have mapped the vector index p to an integer value by defining an ordering for the spin-orbital basis functions), and R j and ζ j are the position and charge of the jth nucleus in the system. In the plane wave dual basis, the Hamiltonian terms can be partitioned into kinetic and potential terms, respectively We can approximate the time evolution operator by applying the potential terms (which all commute with each other, and so induce no Trotter error), implementing a basis change to plane waves, such that the kinetic term becomes diagonal and can be implemented without Trotter error, and then changing back to the plane wave dual basis (or the equivalent, but starting in the plane wave basis). The second-order Trotter error for these approaches are given by The kinetic and electron-nuclei interaction terms are free-fermionic Hamiltonians. This section presents three ways to decompose the electron-electron interaction term into a sum of products of free-fermionic Hamiltonians such that we can write H v = H(U ) + l H(X l )H(Y l ). We use these decompositions, the aforementioned commutator identities, and the fermionic seminorm properties of free-fermion Hamiltonians to derive expressions for first-and second-order commutator bounds. We calculate the first-order bound here, and refer the reader to Appendix C for calculations of the second-order bounds. The first-order commutator is given by We can simplify the second term using Using the triangle and Hölder inequalities, the fermionic seminorm of the first-order commutator is then upper bounded by A. Chemical potentials When working in a fixed particle number manifold, we can shift the chemical potential of the problem to try and reduce the resulting Trotter error bound. This technique has previously been found to be beneficial in simulations of the Fermi-Hubbard model [24]. We can transform the Hamiltonian to where we have used that n 2 p = n p . This transformation adds a constant C to the diagonal of V pq .

B. Spectral decomposition
In the plane wave dual basis, the electron-electron Coulomb interaction matrix V pq is real symmetric, and therefore admits a spectral decomposition While Eq. (2) corresponds to defining V with V pp := 0, we can also use the chemical potential shift outlined above to set V pp := C. We factorise the Hamiltonian as Here, v i are diagonal N × N coefficient matrices. The first-order bound is given by The second-order bounds are given in Appendix D. This decomposition can be regarded as an instance of the general approach of spectral decomposing tensors h pqrs and we discuss this further in Appendix B.

C. Cholesky decomposition
We can also consider a Cholesky decomposition of the matrix V . The Cholesky decomposition factorises a positive (semi)-definite Hermitian matrix into the product of a lower triangular matrix and its Hermitian conjugate, V = LL † . For the real symmetric matrix V pq , we first shift the chemical potential to make V positive definite. The Cholesky decomposition is then given by We can then factorise the Hamiltonian as where L i are diagonal coefficient matrices such that [L i ] pq = δ pq L qi . The first-order bound is given by The second-order bounds are given in Appendix D. As the Cholesky matrix L is lower triangular, the free-fermionic Hamiltonians H(L i ) become increasingly low rank at higher values of i, suggesting that this decomposition may not fully exploit fermion number.

D. Cosine decomposition
We consider the following decomposition that depends explicitly on the structure of the terms in the matrix V pq . We introduce the shorthand ω p ν := k ν · r p . Applying the double angle formula to Eq. (23) yields We can use this to write In the fixed electron-number manifold, the final term will only contribute a global phase during Hamiltonian simulation, and so can be dropped. We can rewrite H v as where C ν and S ν are diagonal N × N coefficient matrices defined by The first-order bound is given by The second-order bounds are given in Appendix D. We remark that it is possible to further simplify Eq. (38) to The more compact form of this decomposition suggests that it may offer a tighter bound. However, the resulting free-fermionic Hamiltonians p e ±iω p ν n p are non-Hermitian, and so yield operators that are neither Hermitian nor anti-Hermitian when commuted with the Hermitian kinetic operator. The resulting matrices may not be diagonalisable, making it unclear how to efficiently evaluate the fermionic seminorm of the operator. In Table I, we reported that the cosine decomposition is the top-ranked approach in the low-filling fraction regime.

E. Outlook
In the following sections, we will apply these bounds to the 2D uniform electron gas in a plane wave dual basis set. The Hamiltonian for this system is given by Eq. (2), but with U p = 0 ∀p. As a result, all of the commutators containing U can be dropped from the above expressions when considering this system. In the following section, we provide further background on the uniform electron gas. We then present numerical results comparing the Trotter error bounds derived above for the uniform electron gas, which form the basis of the rankings assigned in Table I.

V. Uniform electron gas
The uniform electron gas consists of η electrons in a box of size Ω. We are interested in the properties of this system as it scales to the thermodynamic limit -where η, Ω → ∞, but the electron density ρ = η/Ω stays constant. At zero temperature, the physics of the system depends only on ρ. It is conventional to define a quantity referred to as the Wigner-Seitz radius r s , that represents the average distance between electrons in the simulation cell. For a 3D simulation cell, the Wigner-Seitz radius is given by r s = (3/4πρ) 1/3 (for a 2D simulation cell, r s = 1/πρ). In order to make the system charge neutral, the electrons are immersed in a uniformly distributed sea of positive charge. Consequently, the system is often also referred to as 'Jellium'. The Hamiltonian of the Jellium is given by [25] where the first term represents the kinetic energy of the electrons, the second term describes the Coulomb repulsion of the electrons, the third term is interaction of the electrons with the uniform charge density of the positive background, and the final term is the self-interaction of the background charge. The long range nature of the Coulomb interaction causes divergences in the final two terms as the system scales to the thermodynamic limit. These divergences can be cancelled with a divergence of the opposite sign that arises in the electron-electron interaction term. The length scales are typically rescaled to be measured in Bohr radii ( 2 /me 2 ). When performing calculations on Jellium, we can either consider the real-space formulation of the problem discussed above, or project the Hamiltonian onto a basis set.
In addition to acting as a simple model of interacting electrons, the energy density of Jellium is used to parameterize some of the functionals used in density functional theory [26][27][28]. Although the behaviour of Jellium is well understood in the low [29] and high [25,30,31] density limits, small energy differences in the intermediate regime lead to difficulty in resolving competing phases. This has led to unresolved questions about the existence of a superconducting phase in 2D Jellium [32][33][34], as well as disagreements on the order of 0.7 mHartree per electron between different density functional parametrizations at electron densities of interest [35]. While existing computational techniques, such as quantum Monte Carlo methods, are able to obtain accurate energies of relatively large system sizes, these methods typically introduce an uncontrolled bias. It is conventional to perform calculations on a succession of system sizes, which enables extrapolation to the thermodynamic limit. Extrapolation and correction for finite size effects [35][36][37][38] often accounts for a large amount of the uncertainty present in the values estimated [25].
Quantum Monte Carlo (QMC) methods, in particular, variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC), are the leading techniques for calculating the ground state energy of Jellium. Following the pioneering calculations of Ceperley and Alder [39,40], there have been a number of VMC/DMC calculations on both 3D Jellium [35,38,41,42] and 2D Jellium [43][44][45][46][47] (see Ref. [48] for a review of QMC calculations). Both VMC and DMC are typically performed in real-space, and have been applied to systems with on the order of 10 3 electrons [25].
However, these methods are particularly susceptible to the fermion sign problem. This is typically mitigated by fixing the nodal points of the wavefunction to those of the trial wavefunction. Although this fixed node approximation is believed to work well for the uniform electron gas [45], it introduces an uncontrolled bias that is not systematically improvable. While techniques can be used to mitigate this error, DMC energies for high density (r s ≤ 5) electron gases are thought to possess an error of around 1 mHartree per electron (the fixed node error is believed to be smaller at larger r s values) [35,49]. State-of-the-art DMC calculations require on the order of 10 2 CPU core hours [50].
Calculations have also been performed using full configuration interaction quantum Monte Carlo (FCIQMC) [51], which evolves a population of random walkers using update rules that effectively propagate the wavefunction in imaginary time. FCIQMC is applied to systems that have been projected onto a basis set (typically plane waves for Jellium calculations). While this projection appears to mitigate the fermionic sign problem, it introduces a basis set error that must be eliminated by extrapolation to the continuum limit [52]. The basis set error decays as 1/N , although this may be improved using explicitly correlated methods [53]. FCIQMC formally scales exponentially with the system size, but can in practice achieve bias-free results for small, weakly correlated Jellium systems (e.g. 19 electrons at r s = 1 [35]). The approach is also practical for larger system sizes at high densities; producing more accurate results than DMC in 54 electron systems with r s ≤ 1 [50]. Modern FCIQMC methods require around 10 3 − 10 5 CPU core hours (depending on the value of r s investigated) [49,50]. As r s increases, the correlation present in the system becomes large, which makes FCIQMC methods too costly to converge [49,50].
Calculations can be made more challenging by considering the system at non-zero temperature, which acts as a model for the interiors of stars and planets, or for laser-ignited plasma used in fusion experiments [54,55]. Alternatively, we can consider additional interactions, such as spin-orbit coupling [56].
The uniform electron gas has previously been identified as a candidate system for quantum phase estimation [22] due to the desire to seek accurate, bias-free ground state energies. Existing resource estimates for applying phase estimation to Jellium [14,20] project the Jellium Hamiltonian onto the plane wave dual basis, and so can be directly compared with FCIQMC methods. As discussed above, these calculations must first be extrapolated to the basis set limit, before extrapolation to the thermodynamic limit is performed. Previous estimates [14,20] have only considered the quantum resources required for phase estimation at half-filling (η = N/2). However, the most challenging calculation performed in a realistic study will be that with the largest computationally feasible η value, subject to the constraint that N η. Without this constraint, it will not be possible to perform an accurate extrapolation to the continuum limit. In this work, we explicitly consider this regime of interest, and make use of the fermionic seminorm bounds presented in Section IV to reduce estimates of the Trotter error bound, compared to the state-of-the-art [14].

A. Trotter error comparison
We have numerically evaluated the Trotter error bounds derived in Section IV for 2D uniform electron gas systems with up to 49 electrons in 512 plane wave dual spin-orbitals. These calculations were performed as outlined in Appendix E, with the help of subroutines present in OpenFermion [57], an electronic structure package for quantum computational chemistry 1 .
In Figure 1 we plot the first (W 1 ) and second (W 2 ) order commutator bounds for the Hamiltonian decompositions discussed in this work. We consider a simulation cell resolved with 200 spin-orbitals, and vary the number of electrons in the cell. The Wigner-Seitz radius is set to r s = 5. We fix the electron density, such that the volume of the simulation cell increases proportionally with the number of electrons considered. As the cell volume increases, the Hamiltonian coefficients decrease in magnitude. This effect will contribute to a reduction of the commutator bound. However, increasing the number of electrons in the system also increases the number of eigenvalues considered when taking the fermionic seminorm of the relevant free-fermionic coefficient matrices. This effect increases the commutator bound. The competition between these effects can lead to non-trivial behaviour as the number of electrons is varied -this is particularly evident for the spectral and cosine decomposition bounds. These are the decompositions that maximally exploit the fermionic seminorm, leading to their improved behaviour in the low-filling regime. In contrast, the Pauli and Fermionic commutator bounds receive no benefit from decreasing fermion number η. However, we see that these bounds, as well as our Cholesky bound (which only makes partial use of the fermionic seminorm) perform  well close to half-filling, due to their sensitive dependence on the Hamiltonian coefficients. Although it is masked by the log-scale used in the plots, we observe that for a fixed number of spin-orbitals ||T || ∝ 1/η and ||V || max ∈ O(1), so the first-order SHC bound is proportional to η, and the second-order SHC bound is proportional to η 2 . While the SHC bound exploits the fermionic seminorm, it does not fully exploit the reduction in Hamiltonian coefficient magnitudes at high-filling fractions.
In Figure 2 we again plot the first and second-order commutator bounds, but here keep the fermion number fixed at η = 49 (as well as fixing r s = 5, and the cell volume) and instead vary the number of spin-orbitals used. The secondorder Pauli and Fermionic commutator bounds were only calculated up to 128 and 288 spin-orbitals, respectively, as the memory required for these calculations was prohibitive beyond this point. We have extrapolated the performance of the second-order Pauli bounds to larger N values, as described in Appendix F. Performing simulations with a fixed number of electrons, while increasing the number of spin-orbitals, would enable us to perform extrapolation to the basis-set limit. We observe that close to half-filling, the Cholesky and Pauli bounds outperform all others considered. However, the cosine decomposition performs best in the low-filling fraction regime, due to its increased exploitation of the fermionic seminorm.

B. Phase estimation resource estimates
In this section, we discuss the resources required for performing Trotter-based phase estimation on the uniform electron gas systems discussed in the previous section. Our cost estimates focus on the number of logical qubits and T & Toffoli gates required (as these are the dominant factors in surface code-based resource estimates), and neglect the costs of Clifford gates. Our approach closely follows that of Ref. [14], with an improved use of Hamming weight phasing (HWP) [24].
We distribute the total budget for error in energy estimation (δ = ∆ P E + ∆ T S + ∆ syn ) roughly as follows: 33% to Trotter error ∆ T S ; 66% to phase estimation error ∆ P E ; and 1% to rotation synthesis error ∆ syn . In practice, we numerically optimise the error budget allocated to rotation synthesis error, but the optimal choice only differs slightly from 1%. With this split of the error budget, one finds [14,24]  and second (right) order commutator bounds for a 2D uniform electron gas system with rs = 5, and 49 electrons, as a function of the number of spin-orbitals used to resolve the system. Fermionic and Pauli commutator bounds could not be calculated for all datapoints, due to the large memory requirements of those approaches. The 'projected Pauli bounds' were obtained as described in Appendix F. For the Cholesky decomposition the chemical potential was shifted by the minimum value that ensured V was positive definite.
we outline in Appendix G, each Trotter step can be implemented withÕ(N 2 ) non-Clifford gates for an N spin-orbital problem. Therefore, the total algorithm complexity isÕ(N 2 W 1/2 2 /δ 3/2 ) where W 2 contains some dependence on N and η. The primary focus of our work has been to tighten the values of W 2 and we expect a factor C reduction in W 2 will lead to a corresponding factor C 1/2 runtime improvement.
We numerically count the non-Clifford resources for a range of different r s , η, and N values, as shown in Table II. We consider an architecture that distills T gates as its non-Clifford resource. We compare the gate counts obtained by our Trotter-based approach to those obtained using the Trotter-based approach of Ref. [14] (which used the fermionic commutator bound on the Trotter error) and those obtained using the qubitization-based method of Ref. [20]. We consider an extensive error bound δ = 1 mHartree per electron, consistent with leading classical approaches [35,49,50]. It is too memory intensive to calculate the fermionic commutator ('FC') bounds for 16 × 16 systems, showing the limitations of the prior art.
Comparing the gate counts obtained using our novel Trotter error bounds to those obtained using the existing fermionic commutator bound, we observe a reduction in T count by a factor of between 2.3−12.7×. This improvement is more pronounced at lower filling fractions, demonstrating the anticipated benefit of using the fermionic seminorm. The largest of these improvements stems from a reduction in Trotter error by a factor of 150 for η = 10, N = 288, r s = 10. In the high accuracy regime of 49 electrons in ∼ 1000 spin-orbitals, we would expect our bounds to provide an order-of-magnitude improvement over the prior art as this would be similar to the improvements showcased by our η = 10, N = 288 results.
Comparing our results to those of qubitization, we see that qubitization consistently (for r s = 5, 10) achieves a lower T count for the systems considered, by a factor of 2 − 11×. This comes at a cost of using 5 − 7× more ancilla qubits. We show in Appendix H that for 2D Jellium at large r s values, the cost of qubitization is roughly independent of r s (when N, η are fixed). In contrast, the cost of our Trotter-based approach scales as 1/r s . These scalings are evident in Table. II. As such, the Trotter-based approach will be the more suitable method for calculations probing the phase diagram of 2D Jellium, which target r s ≥ 20 [47]. In contrast, qubitization will likely perform better for the warm, dense phase (r s < 1 [54,55]). Our Trotter-based approach also scales less efficiently with target error than qubitization (δ −3/2 vs δ −1 ), and so the advantage of qubitization will also decrease if the target error in our calculations is loosened. A comparison of resource estimates for phase estimation of Jellium, using three different methods. We consider an energy error budget of δ = 1 mHa per electron. 'Our best' refers to the Trotter-based phase estimation considered in this work, using our best bound for the Trotter error. 'FC' refers to the Trotter-based phase estimation considered in Ref. [14], which uses the fermionic commutator bound for the Trotter error. The fermionic commutator (FC) bounds are too memory intensive to be calculated for the 16 × 16 systems. In both Trotter methods, 16 additional qubits are used (14 for Hamming weight phasing, one for phase estimation, and one for gate synthesis). 'Qubitization' refers to the post-Trotter approach considered in Ref. [20] (which we discuss in Appendix H). Four T gates can be used to implement a Toffoli gate, so the total aggregated T count for the algorithm is NT + 4N tof .
As a final caveat, this analysis assumes that we can prepare the main register in the desired energy eigenstate. If we are only able to prepare a state with overlap γ < 1, then the circuit depth required is increased by a factor of 1/γ. For the sake of comparison with prior art, we assume that γ = 1, but note that it is an open question whether an eigenstate with sufficient overlap can be prepared [11,58,59]. It will be necessary to repeat the phase estimation process a number of times, to ascertain that phase estimation has found the desired eigenstate. One can also consider other methods of phase estimation, such as that of Ref. [60], which requires an increased number of repetitions of the algorithm, but that has coherent circuit depth independent of γ.

VII. Discussion
We have demonstrated a substantial benefit of our approach to calculating Trotter errors, both in terms of tightness of the bound and the classical runtime and memory complexity. We have primarily focused on second-order Trotter in the plane wave dual basis, but our techniques naturally generalize. For more compact basis sets, fewer orbitals are required, but the Hamiltonian contains O(N 4 ) terms instead of O(N 2 ). In such a compact basis set, the spectral and Cholesky decompositions are still applicable [13], but it is unclear whether an analogue of the cosine decomposition could be used to obtain an even tighter bound in the low-filling fraction regime. Fourth-order Trotter may produce results competitive with those here [23,61], if 5 4 · W 4 < W 2 2 , and a similarly low-overhead compilation of the Trotter circuit can be found. While the methods introduced in this work apply straightforwardly to higher-order Trotter, calculating the fourth-order bounds would require time scaling as O(N 7 ), making it a potentially costly endeavour.
While this work has focused on the performance of Trotter methods, so-called post-Trotter methods [20,[62][63][64][65][66] are known to have superior asymptotic performance with respect to target error. These methods have also leveraged Hamiltonian factorizations to reduce costs [16][17][18]. Trotter methods often possess good constant prefactors in the runtime and require few additional ancilla qubits, compared to post-Trotter methods. As such, it has been proposed [14,24] that Trotter methods could perform better at some tasks in the pre-asymptotic regime. The gate counts presented in Sec. VI B show that our Trotter approach can be competitive with post-Trotter methods like qubitization, in some regimes of interest, and will even use fewer gates than qubitization for large enough Wigner-Seitz radius. It is currently unclear whether second quantized post-Trotter methods can similarly exploit low-filling fractions, which appears to strengthen the case for Trotter methods in this regime. Working in first quantization, one could certainly exploit low-filling fractions, but quantum algorithms would need to be substantially modified to work in this setting [67,68].

VIII. Acknowledgements
We thank Hsin-Yuan (Robert) Huang, Fernando Brandao, Mario Berta and Michael Kastoryano for discussions through this project. Yuan Su's contribution to this project was made while at Caltech. He was supported in part by the National Science Foundation RAISE-TAQS 1839204 and Amazon Web Services, AWS Quantum Program. The Institute for Quantum Information and Matter is an NSF Physics Frontiers Center PHY-1733907.
We consider simulating the following class of interacting electrons where a † j and a k are the fermionic creation and annihilation operators, n l are the occupation-number operators, T and V are coefficient matrices, and the summation is over N spin orbitals. We seek to bound the fermionic seminorm of the nested commutators We know from [Ref. [15], Eq. (60)] that (A2) Applying [Ref. [15], Eq. (77)], we get the following expansion which implies through [Ref. [15], Proposition 10] Similarly, we have from [Ref. [15], Eq. (78)] which implies

B. Spectral decompositions
Here we review how a general electronic structure Hamiltonian can be factorised using spectral decompositions, with slight modifications allowing Cholesky decompositions to be used. If the Hamiltonian is given in the form then we first use the fermionic anti-commutation rules to rewrite it in "chemist notation" as follows where in H we have changed variables so that s → q → r → s and introduced V pqrs := h prsq . Since H 0 is free fermionic, we aim to find a factorization of H . We define a matrix V (pq),(sr) = V pqrs with composite indices (pq) and (sr) so that Hermiticity of H entails we can always choose V (pq),(sr) to be Hermitian so that V (pq),(sr) = V * (sr),(pq) . Indeed, if V (pq),(sr) is not initially Hermitian, we can always map V (pq),(sr) → (1/2)(V (pq),(sr) + V * (sr),(pq) ) and confirm that this transformation results in the same Hermitian H . Therefore, we can diagonalize the matrix with elements V (pq),(sr) so that where λ are real eigenvalues and u i,j are matrix elements of a unitary U . Substituting this into the expressions for H we get We define L := pq u pq, a † p a q which is the first bracketed factor above. Notice that by changing dummy variables in the summation p → s and q → r we also have L = sr u sr, a † s a r . Taking the Hermitian conjugate, we have L † = sr u * sr, a † r a s , which corresponds to the second bracketed factor in Eq. (B8). Therefore, where L is free-fermionic with coefficient matrix X with matrix elements [X ] p,q = u pq, . Therefore, L = H(X ) and L † = H(X † ). Therefore, Computing the fermionic seminorm is significantly easier for H(A) with Hermitian A. In general, the individual factors H(X ) and H(X † ) might not be Hermitian, even though the full Eq. (B10) is Hermitian. There are two possible solutions to enforce Hermiticity of the factors. Following Sec IV. A of [10], we can always decompose X in terms of a Hermitian and skew-Hermitian part X = A + iB so that A and B are Hermitian. Then for each term we have The term −H(i[A , B ]) is free-fermionic and Hermitian and so can be added to the free-fermionic part H 0 . The full expression for H therefore will have 2L terms of the form H(A ) 2 or H(B ) 2 . The above approach is fully general, but results in a doubling of the number of terms in the summation. In some cases, we can directly ensure Hermiticity of H(X ) without any increase in the number of terms. Here we expand on the discussion given in [13,16] but warn the reader that [13] contains notational errors. When the basis set used for the fermionic orbitals is real-valued (such as for Gaussian basis sets or the plane wave dual basis), then the Hamiltonian constants V pqrs are real-valued and have an 8-fold symmetry [13] so that Recall that in the original decomposition for V pqrs we had Eq. (B7). Since the V (pq),(sr) matrix is real and Hermitian, it is therefore diagonalizable by an orthogonal transformation. Since matrix elements of orthogonal transforms are real, we have u * rs, = u rs, and so Next, we will show that we can always map X → (X + X † )/2 and verify that the new decomposition gives the same total Hamiltonian. The transformation X → (X + X † )/2 maps u pq, → (u pq, + u * qp, )/2 = (u pq, + u qp, )/2 and so Using the 8-fold symmetry of Eq. (B13), we have that the Hamiltonian is unchanged under this transform. Note that this approach is essentially a proof that for real-valued orbitals, the skew-Hermitian components can be made to vanish.
For a Hamiltonian spectrally decomposed as described above as we consider a Trotter decomposition where each term in the product formula implements evolution under one of the terms H(h) or λ H(X )H(X ). The first-order commutator bound on the Trotter error is given by (defining The first term can be expanded using Applying the triangle inequality and Hölder inequality, the bound on the fermionic seminorm of the first term is given by The second term can be expanded using Applying the triangle inequality and Hölder inequality, the bound on the fermionic seminorm of the second term is given by Thus This proves Eq. (21) in the main text. Lastly, we discuss how the above decompositions are related to the spectral decomposition used in section IV B. For the relevant special case of the plane wave dual basis (which is a basis of real orbitals) we have V (pq),(sr) = V ps δ p,q δ s,r . In other words, V ps is simply the nonzero sub-block of V (pq),(sr) . As such, it is equivalent to diagonalize the smaller matrix V ps .
And using the simplification of eq. (C15) we get D. Second-order commutator bounds for plane wave dual decompositions In this Appendix, we apply the formulae for the second-order commutator bound for a general decomposition to the decompositions introduced in the main text.
Comparing the spectral decomposition in the main text to Eq. (C3), Eq. (C16), we observe that The second-order commutator bounds are then given by Comparing the Cholesky decomposition in the main text to Eq. (C3), Eq. (C16), we observe that The second-order commutator bounds are then given by Comparing the cosine decomposition in the main text to Eq. (C3), Eq. (C16), we observe that H(X q )H(Y q ) := H(A ν )H(A ν ) such that A ν ∈ {C ν , S ν }. Because these are diagonal matrices, all of the commutators between different X matrices, different Y matrices, and between X and Y matrices vanish. The resulting commutator bounds are given by (D6)

E. Bound computation details
In this section we outline how the commutator bounds discussed in the main text were calculated. For the case of the first-order Fermionic commutator bounds, we must evaluate i,j p,q [T ij a † i a j , V pq n p n q ]. Commutators where p, q are both distinct from i, j trivially commute, leading to O(N 3 ) distinct commutators to store. As a result, the time cost for the algorithm is O(N 4 ) and the memory cost is O(N 3 ). For the second-order bounds, there are O(N 4 ) terms to store, and the algorithm has time cost O(N 6 ). We store the resulting commutators, collect like terms, and then apply the triangle inequality to the sum.
Pauli commutators can be evaluated in a similar manner. After evaluating the fermionic commutators to calculate the error operator, we apply the Jordan-Wigner transform to obtain the error operator written as a sum of tensor products of Pauli operators. Local fermionic operators are mapped to O(N )-local Jordan-Wigner operators, which increases the memory required by a factor of O(N ). The

F. Projected Pauli bounds
The high memory requirements of the Pauli commutator bound make it impractical to calculate W Pauli 2 for N ≥ 128. However, it is evident from Fig. 2 that the Pauli bounds appear to be only a constant factor better than the fermionic commutator bounds. As a result, we can estimate the Pauli commutator bounds for larger N values, using the available fermionic commutator datapoints. We wish to predict W Pauli 2 at N = 162, 200, 242, 288, which for 49 electrons, corresponds to a filling fraction of 0.302, 0.245, 0.202, 0.170, respectively. In Fig. 3, we plot the ratio between the second-order fermionic commutator bound W Ferm 2 , and W Pauli 2 , for a range of N values, varying the number of electrons such that the filling fraction is kept approximately constant. We observe that as the number of orbitals used increases, the ratio gradually increases. The ratio decreases as the filling fraction decreases. As a result, we assume that at N = 162, 200, 242, 288, the Pauli bound outperforms the fermionic bound by roughly a factor of 8.  3. The ratio between the second-order fermionic commutator bound W F 2 and the second-order Pauli commutator bound W P 2 , as a function of the number of spin-orbitals used, for a homogeneous electron gas system with rs = 5. The number of electrons in each calculation is varied, in order to match the specified filling fraction as closely as possible.

G. Phase estimation resource costs
We first discuss the resources used to implement a single Trotter step of time evolution. When the number of Trotter steps is large, the difference in gate count per Trotter step between implementing e i t 2 Hv e itHt e i t 2 Hv and e i t 2 Ht e itHv e i t 2 Ht is negligible. The final term of each Trotter step can be merged with the first term of the next, so that each Trotter step contains one implementation of e iHt and one of e iHv . The difference in total gate count between these two approaches is thus determined by the difference in Trotter error of the orderings.
For the uniform electron gas, H v = p =q V pq n p n q contains N (N − 1)/2 terms, and so can be implemented by an equivalent number of arbitrary angle Z rotations. However, the translational invariance of Jellium leads many of these rotations to be of the same angle. As discussed in Ref. [14], these rotations can be implemented via HWP in groups of size N/2. In practice, we counted the multiplicity of the terms in V pq , and used HWP to reduce the number of arbitrary rotations required. This contributes an O(N 2 log( −1 )) gate complexity to each Trotter step. Low and Wiebe [69] proposed an alternative approach that would need only O(N log(N ) log( −1 )) gates, but with a significant constant factor overhead that makes it more expensive in the regime considered here.
Changing from the plane wave dual to the plane wave basis can be accomplished using either the fermionic fast Fourier transform (FFFT, when the lattice sides are a power of two) [14,22,70,71], or using Givens rotation circuits [72,73]. These approaches have similar costs for Jellium [14]. The FFFT has a recursive structure, and requires L 2 log 2 (L) non-Clifford gates when applied to L qubits. The FFFT must be applied multiple times when changing the basis of a grid in multiple dimensions. For a L x × L y spinful lattice, we require 2L x applications of the FFFT on L y qubits, and 2L y applications of the FFFT on L x qubits (for a d-dimensional spinful lattice of side L, we require 2dL d−1 applications of the FFFT) [14]. Ref. [14] determined that implementing the FFFT requires 26 T gates for 8 qubits, and 81 T gates for 16 qubits. Givens rotations can be used to perform a single-particle orbital basis change, regardless of whether the number of orbitals considered is a power of 2. We follow the approach outlined in Ref. [14]. A single Givens rotation requires two non-Clifford gates, in the form of two arbitrary rotations (by the same angle). A basis change on M qubits requires M 2 Givens rotations. As with the FFFT, we perform the Givens rotations a number of times to change basis in multiple dimensions. For an L x × L y spinful lattice, we require 2L x implementations of the basis change on L y qubits, and 2L y implementations of the basis change on L x qubits. For the former case (with corresponding changes for the latter), we require 2L x × 2× Ly 2 arbitrary rotations. These can be parallelised into Ly 2 groups of size 4L x . The T/Toffoli cost of implementing these arbitrary rotations can be reduced using Hamming weight phasing.
Rotating into the plane wave basis diagonalises the kinetic operator, enabling us to implement it with N arbitrary rotations in the worst case. As it is efficient to classically diagonalise the kinetic coefficient matrix T pq , we can determine the multiplicity of each eigenvalue, and then use HWP to reduce the number of arbitrary rotations required. Overall, implementing this contributes a cost O(N log(N ) log( −1 )) per Trotter step.
To perform phase estimation we must implement not just a circuit approximating e iHt , but a circuit that approximates e iHt controlled on the state of an ancillary register. We can implement a controlled arbitrary rotation at double the cost of the un-controlled operation [74]. However, Ref. [72] introduced an approach known as directionally controlled phase estimation, that reduces the cost of controlled time evolution to be the same as the uncontrolled circuit, when implemented with symmetric product formulae (this approach was elaborated upon further in Refs. [12,14]).
The key insight is that one instance of U 2 (t) can be used to implement [|0 a 0| a ⊗ U 2 (−t) + |1 a 1| a ⊗ U 2 (t)], which for the purposes of phase estimation is equivalent to performing [|0 a 0| a ⊗ I + |1 a 1| a ⊗ U 2 (2t)]. In addition to halving the number of arbitrary rotations required, this optimization effectively doubles the time duration used for phase estimation. We use an adaptive variant of phase estimation that uses a single ancilla qubit [75]. As discussed in Ref. [14], this approach uses N P E applications of directionally controlled phase estimation to learn the energy eigenvalue to a root mean squared error of We note that this formula includes the reduction from t → 2t due to the use of directionally controlled phase estimation. The Trotter error contributes an error ∆ T S = W t 2 where W is the commutator bound constant. A third source of error of error are synthesis errors ∆ syn = O(log(N R / )) where N R is the number of arbitrary Z axis rotations in the algorithm. We distribute errors between these three sources using the approach outlined in Appendix F of Ref. [24].
for Jellium is given by Eq. (54) in Ref. [20] as where for a Hamiltonian written as H = a h a H a (with ||H a || = 1), λ = a |h a |, N is the number of spin-orbitals, and δ is the target energy error. The number of logical ancilla qubits required is given by Eq.(55) of Ref. [20] It is interesting to consider how the gate count scales as a function of r s . We have that for Jellium, λ = λ t + λ v . For a dim-d system, we can see directly from the Hamiltonian coefficients in Eq.(23) (using that Ω ∝ ηr d s ) that In these expressions we have implicitly assumed that N is held constant. As a result, in 2D λ scales as O(η −1 r −2 s ) + O(1). This can be contrasted with our second-order Trotter approach. We have that For d = 2, W 2 = O(η −2 r −4 s ) + O(η −1 r −2 s ). We note that this bound on the Trotter error may be very loose (in terms of the scaling with the number of electrons), as it does not use commutativity of terms in the Hamiltonian or the fermionic seminorm (c.f. Eq. 13).
If we fix η, N and vary r s , we see that for d = 2 and large r s , the cost of qubitization is independent of r s , while our Trotter-based approach scales as O(r −1 s ). Thus, the cost of Trotter-based approaches in 2D reduce as the value of r s is increased, while the cost of qubitization is roughly independent of r s . This is evident in the results presented in Table II.