Nearly tight Trotterization of interacting electrons

We consider simulating quantum systems on digital quantum computers. We show that the performance of quantum simulation can be improved by simultaneously exploiting commutativity of the target Hamiltonian, sparsity of interactions, and prior knowledge of the initial state. We achieve this using Trotterization for a class of interacting electrons that encompasses various physical systems, including the plane-wave-basis electronic structure and the Fermi-Hubbard model. We estimate the simulation error by taking the transition amplitude of nested commutators of the Hamiltonian terms within the $\eta$-electron manifold. We develop multiple techniques for bounding the transition amplitude and expectation of general fermionic operators, which may be of independent interest. We show that it suffices to use $\left(\frac{n^{5/3}}{\eta^{2/3}}+n^{4/3}\eta^{2/3}\right)n^{o(1)}$ gates to simulate electronic structure in the plane-wave basis with $n$ spin orbitals and $\eta$ electrons, improving the best previous result in second quantization up to a negligible factor while outperforming the first-quantized simulation when $n=\eta^{2-o(1)}$. We also obtain an improvement for simulating the Fermi-Hubbard model. We construct concrete examples for which our bounds are almost saturated, giving a nearly tight Trotterization of interacting electrons.

We consider simulating quantum systems on digital quantum computers. We show that the performance of quantum simulation can be improved by simultaneously exploiting commutativity of the target Hamiltonian, sparsity of interactions, and prior knowledge of the initial state. We achieve this using Trotterization for a class of interacting electrons that encompasses various physical systems, including the plane-wave-basis electronic structure and the Fermi-Hubbard model. We estimate the simulation error by taking the transition amplitude of nested commutators of the Hamiltonian terms within the η-electron manifold. We develop multiple techniques for bounding the transition amplitude and expectation of general fermionic operators, which may be of independent interest. We show that it suffices to use n 5/3 η 2/3 + n 4/3 η 2/3 n o(1) gates to simulate electronic structure in the plane-wave basis with n spin orbitals and η electrons, improving the best previous result in second quantization up to a negligible factor while outperforming the first-quantized simulation when n = η 2−o (1) . We also obtain an improvement for simulating the Fermi-Hubbard model. We construct concrete examples for which our bounds are almost saturated, giving a nearly tight Trotterization of interacting electrons.
Lloyd's original work considered the simulation of k-local Hamiltonians. This was subsequently extended to the study of d-sparse Hamiltonians [1,10], which provides a framework that abstracts the design of quantum algorithms from the underlying physical settings. However, despite their theoretical value, algorithms for sparse Hamiltonian simulation do not always provide the fastest approach for simulating concrete physical systems. Hamiltonians arising in practice often have additional features beyond sparseness, such as locality [33,79], commutativity [22,27,74], and symmetry [35,80], that can be used to improve the performance of simulation. Besides, prior knowledge of the initial state [6,28,69,73] and the norm distribution of Hamiltonian terms [18,21,34,49,58] have also been proven useful for digital quantum simulation.
We show that a number of these features, in particular the sparsity, commutativity, and initial-state information, can be combined to give an even faster simulation. We achieve this improvement for a class of interacting-electronic Hamiltonians, which includes many physically relevant systems such as the plane-wave-basis electronic-structure Hamiltonian and the Fermi-Hubbard model. Our approach uses Trotterization-a method widely applied in digital quantum simulation.
Our analysis proceeds by computing the transition amplitude of simulation error within the η-electron manifold. To this end, we develop multiple techniques for bounding the transition amplitude/expectation of a general fermionic operator, which may be of independent interest. For an n-spin-orbital electronic-structure problem in the plane-wave basis, our result improves the best previous result in second quantization [6,27,54] up to a negligible factor while outperforming the first-quantized result [7] when n = η 2−o (1) . We also obtain an improvement for simulating the Fermi-Hubbard model. We construct concrete examples for which our bounds are almost saturated, giving a nearly tight Trotterization of interacting electrons.

Combining interaction sparsity, commutativity, and initial-state knowledge
Sparsity can be used to improve digital quantum simulation in multiple ways. A common notion of d-sparsity concerns the target Hamiltonian itself, where each row and column of the Hamiltonian contains d nonzero elements accessed by querying quantum oracles. As aforementioned, this provides an abstract framework for designing efficient simulation algorithms and is versatile for establishing lower bounds [10], although it sometimes ignores other important properties of the target system, such as locality, commutativity, and symmetry. Another notion of sparsity, closely related to our paper, considers the interactions between the underlying qubits or modes [14,57,63,87]. The sparsity of interactions does not in general imply the underlying Hamiltonian is sparse, but it provides a tighter bound on the number of terms in the Hamiltonian and may thus be favorable to digital quantum simulation.
Trotterization (and its alternative variants [25,27,34,55,62]) provides a simple approach to digital quantum simulation and is so far the only known approach that can exploit the commutativity of the Hamiltonian. Indeed, in the extreme case where all the terms in the Hamiltonian commute, we can simultaneously diagonalize them and apply the first-order Lie-Trotter formula S 1 (t) without error. Previous studies have also established commutator error bounds for certain low-order formulas [76] and specific systems [22,74]. An analysis of a general formula S p (t) is, however, considerably more difficult and has remained elusive until the recent proof of the commutator scaling of Trotter error [27].
A different direction to speeding up digital quantum simulation is to exploit information about the initial state. The error of digital quantum simulation is commonly quantified in previous work by the spectral-norm distance, which considers all possible states in the underlying Hilbert space. But if the state is known to be within some subspace throughout the simulation, then in principle this knowledge could be used to improve the algorithm. For instance, digital quantum simulation in practice often starts with an initial state in the low-energy subspace of the Hamiltonian, so a worst-case spectral-norm analysis will inevitably overestimate the error. To address this, recent studies have considered a lowenergy projection on the simulation error and provided improved approaches, using either Trotterization [6,28,69,73] or more advanced quantum algorithms [51], that can be advantageous when the energy of the initial state is sufficiently small.
Ideally, the sparsity of interactions, commutativity of the Hamiltonian, and prior knowledge about the initial state can all be combined to yield an even faster digital quantum simulation. This combination, however, appears to be technically challenging to achieve. Indeed, the state-of-the-art analysis of Trotterization represents the simulation error in terms of nested commutators of Hamiltonian terms with exponential conjugations [27,Theorem 5]. This representation is versatile for computing the commutator scaling of Trotter error, but it yields little information about the energy of the initial state. To the best of our knowledge, the only previous attempt to address this problem was made by Somma for simulating bosonic Hamiltonians [73], whose solution seems to have a divergence issue. 1 Instead, we combine the sparsity, the commutativity and the initial-state information to give an improved simulation of a class of interacting electrons.
Here, we consider simulating the following class of interacting electrons by Trotterization: where A † j , A k are the fermionic creation and annihilation operators, N l are the occupationnumber operators, τ , ν are coefficient matrices, and the summation is over n spin orbitals. The specific definitions of these fermionic operators are given in Section 2.2. We say the interactions are d-sparse if there are at most d nonzero elements within each row/column of τ and ν. This model represents various systems arising in physics and chemistry, including the electronic-structure Hamiltonians in the plane-wave basis [6] and the Fermi-Hubbard model [28,42].
To apply Trotterization, we need to express the Hamiltonian as a sum of elementary terms, each of which can be directly exponentiated on a quantum computer; see Section 2.1 for a review of this algorithm. For the electronic Hamiltonian (1), it suffices to consider the two-term decomposition H = T + V , as the exponentials of T and V can be directly implemented using various quantum circuits. For instance, all the terms in V commute with each other so e −itV = l,m e −itν l,m N l Nm , where each e −itν l,m N l Nm corresponds to a two-qubit operation under the Jordan-Wigner transformation. On the other hand, the exponential e −itT can be implemented by diagonalization, i.e., e −itT = U e −i λ N U † , where U can be efficiently implemented using Givens rotations [41,65]. In cases where τ and ν are translationally invariant τ j,k = τ j+q,k+q , ν l,m = ν l+q,m+q , we can implement e −itT using the fast fermionic Fourier transform [6] and a related circuit implementation exists for e −itV [54].
We now apply a pth-order Trotterization S p (t) to approximate the evolution of the electronic Hamiltonian (1) for time t. We prove the following bound on the error of this approximation. (1) with n spin orbitals, which we simulate using a pth-order formula S p (t). Then,

Theorem 1 (Fermionic seminorm of Trotter error). Let
Furthermore, if the interactions are d-sparse, Here, · is the spectral norm, · max is the max-norm denoting the largest matrix element in absolute value, and is the fermionic η-seminorm for number-preserving operator X, where |ψ η , |φ η are quantum states with η electrons.
This theorem follows from an inductive estimate of the fermionic seminorm of nested commutators of Hamiltonian terms, and will be formally proved in Section 3 and Section 4. Note that in order to use the prior knowledge of the initial state, we have considered the fermionic seminorm · η of Trotter error with respect to the η-electron manifold. This seminorm is closely related to other metrics used to quantify the impact of initial-state information to digital quantum simulation [6,28,69,73]; see Section 2.3 for a detailed discussion. The resulting bound depends on the number of electrons η, as well as the spectral norm τ , the max-norm τ max , ν max , and the sparsity of interactions d, but there is no dependence on the total number of spin orbitals n. This improves over previous work [27,Theorem F.5] where an explicit n-scaling seems unavoidable. Meanwhile, other prior estimates of the fermionic seminorm [6, Appendix G] [28, Theorem 13] did not exploit commutativity of the Hamiltonian and would introduce an additional factor of η p in the Trotter error bound. Our result thus improves the performance of digital quantum simulation by combining the initial-state information, the interaction sparsity, and commutativity of the Hamiltonian.
A common issue with the Trotterization algorithm is that existing analyses can be very loose for simulating certain physical systems. However, we address this with the following theorem, which shows that the asymptotic scaling of our bound is nearly tight. Theorem 2 (Tightness). For s, w > 0 and positive integer η ≤ n 2 , there exists an interacting-electronic Hamiltonian In addition, for u, w > 0 and positive integer d ≤ η ≤ n 2 , there exists a d-sparse interactingelectronic Hamiltonian (1) with n spin orbitals such that τ max = u, ν max = w, We prove the above theorem by choosing T = n−1 j,k=0 A † j A k and V = contribute to the Trotter error, as well as other types of nested commutators which do not dominate the error scaling (Proposition 1). Modulo an application of the triangle inequality, Theorem 2 then shows that our bound (2) overestimates the Trotter error by a factor of nη in the worst case, whereas (3) overestimates a factor of at most η. For p sufficiently large, this only contributes n o (1) and η o (1) to the gate complexity, respectively. In this sense, we have given a nearly tight Trotterization of interacting-electronic Hamiltonians (1). 2

Main techniques
The proof of Theorem 1 relies on multiple approaches we develop for bounding the fermionic seminorm, which may be of independent interest. Recall from (4) that the fermionic seminorm X η of a number-preserving operator X is the maximum transition amplitude of X within the η-electron manifold.
Our first approach is based on the observation that the fermionic seminorm of X can be alternatively represented using the expectation of X † X, i.e., We then upper bound X † X in terms of the particle-number operator N = j N j , so that the expectation scales with the number of electrons η = ψ η |N |ψ η instead of the total number of spin orbitals. Assuming X is a sum of product of fermionic operators, we contract the summation indices in X † X by using either diagonalization (Lemma 2) or an operator Cauchy-Schwarz inequality (Lemma 1) [61]. To extend this argument to general fermionic operators, we prove a Hölder-type inequality (Lemma 3) and apply it recursively to bound X † X. We detail this recursive approach in Section 3 and use it to prove (2). Our second approach starts by bounding the fermionic seminorm in terms of the maximum expectation value max |ψη | ψ η |X|ψ η |. We then expand X and |ψ η and give a combinatorial argument to count the number of "paths" which have nonzero contribution to the expectation (Proposition 11). We discuss this path-counting approach in more detail in Section 4 and use it to prove (3). It is worth mentioning that the path-counting technique can also be used to prove the following alternative bound for the pth-order Trotterization This is slightly weaker than (2) since τ ≤ n τ max always holds but not necessarily saturates, but in our application it yields the same gate complexity for the electronicstructure simulation in the second-quantized plane-wave basis. We discuss this further in Appendix B.
Note that the expectation of fermionic operators, when taken with respect to the computational-basis states, can be exactly computed using the so-called Wick's theorem [64,86]. However, this approach would introduce unnecessary term reordering which actually complicates our proof. In contrast, the underlying idea of path counting is conceptually simpler and may have potential applications in other contexts beyond the analysis of Trotter error.

Applications
The nearly tight Trotterization of electronic Hamiltonian (1) gives improved simulations of many systems arising in condensed matter physics and quantum chemistry, including the plane-wave-basis electronic-structure Hamiltonian and the Fermi-Hubbard model.
The electronic-structure problem considers electrons interacting with each other and some fixed nuclei. An efficient simulation of such systems could help understand chemical reactions, and provide insight into material properties. Here, we consider representing the electronic-structure Hamiltonian in the plane-wave basis [6]: where ω is the volume of the computational cell, κ µ = 2πµ/ω 1/3 are n vectors of planewave frequencies, µ are three-dimensional vectors of integers with elements in the interval [−n 1/3 , n 1/3 ], r j are the positions of electrons; ζ ι are nuclear charges; and r ι are the nuclear coordinates. We can represent this Hamiltonian in the form (1) with Assuming a constant system density η = O (ω), Theorem 1 then implies that This approximation is accurate for a short-time evolution. To simulate for a longer time, we divide the evolution into r steps and apply S p (t/r) within each step, obtaining Simulation Algorithm n, η η = Θ(n) Interaction-picture (Ref. [7], first quantization) Qubitization (Ref. [7], first quantization) Interaction-picture (Ref. [54], second quantization) O n 8/3 Trotterization (Ref. [6], second quantization) n 5/3 η 1/3 + n 4/3 η 5/3 n o(1) n 3+o (1) Trotterization (Ref. [27], second quantization) (1) Trotterization (Theorem 1, second quantization) Therefore, steps suffices to simulate for a constant time and accuracy with a pth-order Trotterization. Implementing each step using the approach of [54, Sect. 5] and choosing p sufficiently large, we obtain the gate complexity Up to the negligible factor n o (1) , this improves the best previous result in second quantization while outperforming the first-quantized simulation when n = η 2−o(1) . See Table 1 for a gate-count comparison. We discuss this in detail in Section 6.1. We also consider applications to the Fermi-Hubbard model, which is believed to capture the physics of some high temperature superconductors. This model is classically challenging to simulate [44,88], but is a potential candidate for near-term quantum simulation [16,17,48,66]. We have where j, k denotes a summation over nearest-neighbor lattice sites and σ ∈ {0, 1} labels the spin degree of freedom. The Fermi-Hubbard model represents a lattice system with nearest-neighbor interactions and, according to [22], can be simulated with O n 1+1/p gates using a pth-order Trotterization for a constant time and accuracy. On the other hand, recent work [28] shows that the Trotterization algorithm has gate complexity O nη 1+1/p when restricted to the η-electron manifold. By simultaneously using the sparsity of interactions, commutativity of the Hamiltonian and information about the initial state, we show in Section 6.2 that O nη 1/p gates suffices, improving both results for the Fermi-Hubbard model. 3 We conclude the paper in Section 7 with a discussion of the results and some open questions.

Preliminaries
In this section, we summarize preliminaries of this paper, including a discussion of the Trotterization algorithm and its error analysis in Section 2.1, a brief summary of the second-quantization representation in Section 2.2, and an introduction to the fermionic seminorm and its properties in Section 2.3.

Trotterization and Trotter error
The Trotterization algorithm approximates the evolution of a sum of Hamiltonian terms using exponentials of the individual terms. For the interacting-electronic Hamiltonian (1), it suffices to consider a two-term decomposition H = T + V , as the exponentials of T and V can be directly implemented on a quantum computer. Then, the ideal evolution under H for time t is given by e −itH = e −it(T +V ) , which can be approximated by a pth-order product formula S p (t), such as the first-order Lie-Trotter formula (17) and (2k)th-order Suzuki formulas [77] where u k := 1/(4 − 4 1/(2k−1) ). This approximation is accurate when t is small. To simulate for a longer time, we divide the evolution into r Trotter steps and apply S p (t/r) with Trotter error at most /r. We choose r sufficiently large so that the simulation error, as quantified by the spectral norm S r p (t/r) − e −itH , is at most . Trotterization (and its alternative variants) provides a simple approach to digital quantum simulation and is so far the only known approach that can exploit commutativity of the Hamiltonian. Indeed, in the extreme case where all the Hamiltonian terms commute, Trotterization can implement the exact evolution without error. Previous studies have also established commutator analysis of Trotter error for systems with geometrical locality [22] and Lie-algebraic structures [74], as well as certain low-order formulas [76], including the first-order Lie-Trotter formula and the second-order Suzuki formula An analysis of the general case is, however, considerably more difficult and has remained elusive until the recent proof of commutator scaling of Trotter error [27]. Here, we introduce a stronger version of that result which can be proved as in [26, Appendix C] by combining Theorem 3, 4, and 5 of [27] without invoking the triangle inequality. 4 Proposition 1 (Commutator representation of Trotter error). Let H = T + V be a twoterm Hamiltonian and S p (t) be a pth-order formula. Define H 0 = V and H 1 = T . Then, where γ γ γ ∈ {0, 1} p+1 are binary vectors 5 and j goes through a constant range of numbers (depending on the order p). Here, U γ γ γ,j (τ 1 , τ 2 ) and W γ γ γ,j (τ 1 , τ 2 ) are products of evolutions of T and V with time variables τ 1 and τ 2 and a γ γ γ (τ 1 , τ 2 ) are coefficients such that As an immediate application, we find that the spectral norm of the Trotter error scales with nested commutators of the Hamiltonian terms, i.e., Note that the use of max γ γ γ in place of γ γ γ does not change the scaling as γ γ γ only ranges over a constant number of binary vectors. We then divide the evolution into r steps and apply the triangle inequality to obtain It thus suffices to choose to ensure that the error of simulation is no more than . The above analysis is versatile for computing the commutator dependence of Trotter error. Unfortunately, the resulting bound does not use prior knowledge of the initial state and will in particular be loose if the initial state lies within a low-energy subspace. On the other hand, recent work of Şahinoğlu and Somma proposed a Trotterization approach for simulating low-energy initial states but the commutativity of the Hamiltonian was ignored in their analysis [69]. Here, we address this by simultaneously using commutativity of the Hamiltonian and prior knowledge of the initial state to improve the simulation of a class of interacting electrons. We obtain further improvement when the electronic Hamiltonian has sparse interactions. In the following, we introduce preliminaries about the second-quantization representation (Section 2.2) and the notion of fermionic seminorm (Section 2.3), on which our analysis will be based.

Second-quantization representation
In this section, we review several facts about the second-quantization representation that are relevant to our analysis. We refer the reader to the book of Helgaker, Jørgensen, and Olsen [37] for a detailed discussion of this topic.
We use the abstract Fock space to represent electronic Hamiltonians. Specifically, for a system of n spin orbitals, we construct a 2 n -dimensional space span{|c 0 , c 1 , . . . , c n−1 } spanned by the basis vectors |c 0 , c 1 , . . . , c n−1 , where c j = 1 represents that mode j is occupied and c j = 0 otherwise. General vectors in the Fock space, denoted by |ψ or |φ , are then given by linear combinations of these orthonormal basis vectors. We define the η-electron subspace span{|c 0 , c 1 , . . . , c n−1 , j c j = η}. By considering all 0 ≤ η ≤ n, we obtain the decomposition where k denotes the orthogonal direct sum. Using bold symbol c c c to represent an arbitrary fermionic configuration and |c c c| = j c j to denote the Hamming weight, we have We say that normalized vectors in the η-electron subspace form the η-electron manifold and denote an arbitrary such vector by |ψ η or |φ η . The n elementary fermionic creation operators are defined through the relations whereas the fermionic annihilation operators are defined by The use of † is justified by the fact that A † j is indeed the Hermitian adjoint of A j with respect to the inner product in the Fock space. We also introduce the occupation-number operators N j = A † j A j and add them together to get the particle-number operator N = n−1 j=0 N j . Fermionic creation and annihilation operators satisfy the canonical anticommutation relations where the Kronecker-delta function δ j,k is one if j = k and zero otherwise. Applying these, we obtain the following commutation relations of second-quantized fermionic operators.

Proposition 2 (Commutation relations of fermionic operators).
The following commutation relations hold for second-quantized fermionic operators: We say a fermionic operator is number-preserving if every η-electron subspace is invariant under the action of this operator. Equivalently, operator X is number-preserving if and only if it commutes with the particle-number operator, i.e., [N, X] = 0. Yet another equivalent definition is based on the notion of η-electron projections: letting Π η be orthogonal projections onto the η-electron subspaces, then X is number-preserving if and only if it commutes with every Π η , namely, [Π η , X] = 0. In the matrix representation, X is block-diagonalized by the set of η-electron projections {Π η }.
A special example of number-preserving operator is the particle-number operator N , which acts as a scalar multiplication by η within the η-electron subspace. Other examples include excitation operators A † j A k , occupation-number operators N l , and elementary exponentials in the Trotterization algorithm e The fermionic Fourier transform [29] as given below is also number-preserving: 6 In fact, the set of number-preserving operators contains identity and is closed under linear combination, multiplication, Hermitian conjugation, and taking limit.

Fermionic seminorm
We now introduce the notion of fermionic seminorm, which we use to quantify the error of the Trotterization algorithm that takes the prior knowledge of the initial state into consideration. For any number-preserving operator X and 0 ≤ η ≤ n, we define the fermionic ηseminorm as the maximum transition amplitude within the η-electron manifold: where |ψ η , |φ η are quantum states containing η electrons. 7 When there is no ambiguity, we drop the dependence on η and call X η the fermionic seminorm of X. As the name suggests and the following proposition confirms, the fermionic seminorm is indeed a seminorm defined on the closed unital †-subalgebra of number-preserving operators.
Proposition 4 (Seminorm properties). The following properties hold for the fermionic seminorm: Proof. We will only prove the third statement, as the remaining follow directly from the definition of the fermionic seminorm. We consider where Π η is the orthogonal projection onto the η-electron subspace and the last step follows from the Cauchy-Schwarz inequality. To proceed, we optimize over an arbitrary state |ϕ to get assuming Π η |ϕ = 0, as the case Π η |ϕ = 0 never leads to maximality. But on the other hand, This completes the proof of the third statement.
The fermionic seminorm, as defined in (33) by the maximum transition amplitude within the η-electron manifold, provides a reasonable metric for quantifying the error of digital quantum simulation with initial-state constraints. Indeed, a seminorm similar to our definition was used by Somma [73] for analyzing quantum simulation of bosonic Hamiltonians. However, we point out that this is not the only error metric that takes the prior knowledge of the initial state into account. Recent work [69] analyzed the lowenergy simulation of k-local frustration-free Hamiltonians by computing the spectral norm of Trotter error projected on the low-energy subspace. However, the following proposition shows that these two error metrics are the same for fermionic systems.
Proposition 5 (Fermionic seminorm as a projected spectral norm). For any numberpreserving operator X, it holds that Proof. The underlying idea behind this proposition is already hinted in the proof of Proposition 4. We have But on the other hand, assuming Π η |φ = 0 and Π η |ψ = 0, as the zero vector will not lead to maximality. The proposition then follows since number-preserving operator X commutes with the η-electron projection Π η .
Another common approach to quantify the error of digital quantum simulation is to take the maximum expectation max |ψη | ψ η | · |ψ η | within the η-electron manifold. This approach is used by previous work [6,28,65] and appears to give a natural metric when digital quantum simulation is used as a subroutine in phase estimation. We show that this only differs from our definition (33) by at most a constant factor, reaffirming the fermionic seminorm as a proper error metric for simulating fermionic systems. Proposition 6 (Transition amplitude and expectation). For any number-preserving operator X, the following statements hold: Note that the second inequality is tight by considering X = A † 0 A1 on a fermionic system with two spin orbitals and one electron.
Proof. The first statement follows from the fact that Π η XΠ η is Hermitian and that the spectral norm of a Hermitian operator is its largest eigenvalue in absolute value. For the second statement, The first inequality of Statement 3 is trivial. For the second inequality, we apply the polarization identity from which the claimed inequality follows by maximizing over states |ψ η and |φ η .
We now apply Proposition 1 to compute the fermionic seminorm of Trotter error, obtaining We find that the resulting error bound depends on the fermionic seminorm of nested commutators, and the performance of digital quantum simulation can thus be potentially improved by simultaneously exploiting commutativity of the Hamiltonian and prior knowledge of the initial state. However, the main difficulty here is to give a tight estimate of H γ p+1 , · · · [H γ 2 , H γ 1 ] η , which seems technically challenging to address. To this end, we develop two approaches for bounding the expectation/transition amplitude of general fermionic operators in Section 3 and Section 4 and prove our main result Theorem 1, establish the tightness of our bound in Section 5, and discuss applications and further implications of our result in Section 6 and Section 7.

Recursive bound on the expectation of fermionic operators
In this section, we present our first approach for bounding the expectation of fermionic operators, and thereby bounding the fermionic seminorm of Trotter error. We introduce in Section 3.1 the main techniques used in our approach, including an operator Cauchy-Schwarz inequality, a diagonalization procedure, and a Hölder-type inequality for the expectation value. We then describe our approach in detail and apply it to prove Eq. (2) of our main result Theorem 1. The proof is based on induction: we analyze the base case in Section 3.2 and the inductive step in Section 3.3, respectively.

Main techniques
Recall that the main technical challenge to estimate the simulation error of the electronic Hamiltonian (1) is to bound the fermionic seminorm H γ p+1 , · · · [H γ 2 , H γ 1 ] η , where γ j = 0, 1, H 0 = V and H 1 = T . Applying the commutation relations in Proposition 2, we see that we need to analyze a general fermionic operator of the form X = j j j,k k k,l l l w j j j,k k k,l l l · · · A † jx · · · N lz · · · A ky · · · (44) Our first approach starts by reexpressing the fermionic seminorm of X using the expectation of X † X: We note that X † X is a positive semidefinite operator, and an upper bound of it with respect to the partial ordering of positive semidefiniteness will therefore give a bound on the expectation value. We achieve this by contracting the corresponding indices in X and X † , using either an operator Cauchy-Schwarz inequality (Lemma 1) or diagonalization (Lemma 2).
where Hermitian operators are partially ordered according to the positive semidefiniteness.
This implies ± j,k from which the claimed inequality follows.

Lemma 2 (Diagonalization). For any finite list of operators {B j } and Hermitian coeffi-
where Hermitian operators are partially ordered according to the positive semidefiniteness.
Proof. Since µ is Hermitian, we may diagonalize it to µ by unitary transformation w as where µ is a diagonal matrix with all eigenvalues of µ as the diagonal elements. We then define B l := k w l,k B k so that But l B † l B l has identity as the coefficient matrix which is invariant under a change of basis: This completes the proof.
By applying Lemma 1 or Lemma 2, we can get a bound of X † X with respect to the partial ordering of positive semidefiniteness, with one pair of the corresponding indices in X and X † contracted. Indeed, these techniques were used by Otte to establish the boundedness of quadratic fermionic operators in infinite-dimensional Hilbert spaces [61]. However, the difficulty here is that we need to handle more complex products of fermionic operators in the Trotter error estimate. To this end, we prove a Hölder-type inequality for the expectation value, which allows us to apply Lemma 1 and Lemma 2 recursively to get a desired bound.

Lemma 3 (Hölder-type inequality for expectation). For any finite lists of fermionic operators {B
where we assume B j map the η-electron subspace to the ξ-electron subspace and C j are number-preserving. In terms of the fermionic seminorm, we have Proof. The claimed inequality follows from Using the above lemmas, we can now prove Eq. (2) of our main result Theorem 1 by induction. We analyze the base case in Section 3.2 and the inductive step in Section 3.3.

Single-layer commutator
We now prove Eq. (2) of our main result Theorem 1 by induction. In the base case, we consider simulating the interacting-electronic Hamiltonian (1) using the first-order formula S 1 (t). We know from (19) that where To this end, we apply Proposition 2 to expand the single-layer commutator [T, V ] into linear combinations of fermionic creation, annihilation, and occupation-number operators. We have At this stage, it is possible to directly bound the terms in the last equality using Lemma 1, Lemma 2, and Lemma 3 from the previous subsection. However, we will further commute the occupation-number operator in between the creation and annihilation operators, obtaining This additional commutation leads to an error bound with the same asymptotic scaling but a slightly larger prefactor. The benefit is that the analysis can be directly extended to handle the inductive step in the next subsection.
It is worth noting that the above six terms from the expansion of [T, V ] share a similar structure. Specifically, they all consist of creation operator A † j , annihilation operator A k , and (possibly) occupation-number operator N l , with one coefficient matrix τ and one Here, vertices in the graph denote the indices in the summation and edges represent the coefficients. Note that the graphs can be made directional so that they are one-to-one corresponding to fermionic operators, although this is not needed in our analysis and will not be further pursued here.
matrix ν. The main difference between these terms is that the coefficient matrix ν is acting on different indices. See Figure 1 for a graph illustration of this structure. We now bound the asymptotic scaling of the fermionic seminorm for each of the six terms in the commutator expansion.

Proposition 8 (Fermionic seminorm of single-layer commutator). Let
Proof. We present the proof of the first two statements here. The remaining justifications proceed in a similar way and are left to Appendix A.
whereτ j 1 ,k 1 is the complex conjugate of τ j 1 ,k 1 and we have used the anti-commutation rela- and apply the operator Cauchy-Schwarz inequality (Lemma 1): This implies Note that j 1τ j 1 ,k 1 τ j 1 ,k 2 gives the (k 1 , k 2 ) matrix element of τ † τ . Then, we define Now that the indices k 1 and k 2 are contracted, we can apply the Hölder-type inequality for the expectation value (Lemma 3). To this end, we let D † k 1 = m 1ν k 1 ,m 1 N m 1 and compute The first factor can be directly bounded as For the second factor, we have which implies Combining (64), (65), (66), (67), and (69) establishes the first statement.
For the second statement, we let X = j,k τ j,k ν k,k A † j A k and compute (70) Applying Lemma 1, (71) Performing diagonalization using Lemma 2, we have Note that we could directly bound the above operators as τ 2 ν 2 max N 2 and thereby complete the proof. But we choose to instead apply Lemma 3 so that the analysis can then be directly extended to analyze multilayer nested commutators. We have The proof of the second statement is now completed. See Appendix A for the proof of the remaining statements.

Multilayer nested commutators
We now analyze the error of simulating the interacting-electronic Hamiltonian (1) using a general pth-order formula S p (t). We know from (43) that where for each multilayer nested commutator To this end, we assume that H γ p+1 , · · · [H γ 2 , H γ 1 ] is expressed as a fermionic operator of the form j j j,k k k,l l l w j j j,k k k,l l l · · · A † jx · · · N lz · · · A ky · · · and analyze its commutator with either T or V . For the commutator with T , we have from Proposition 2 To develop some intuitions about these commutations, we introduce the notion of fermionic chain, which refers to a product of fermionic operators that has a creation operator on the left and an annihilation operator on the right. Then, the above commutations either extend an existing fermionic chain (in the case where commutator is taken with A † jx or A ky ), or create a new chain (in the case where commutator is taken with N lz ). On the other hand, we also apply Proposition 2 to compute the commutator with V : (78) Unlike the commutator with T , these commutations do not extend an existing chain or create a new chain. Rather, their effect is to append occupation-number operators to an existing chain.
We now apply (77) and (78) iteratively to expand a general multilayer nested commutator H γ p+1 , · · · [H γ 2 , H γ 1 ] . We summarize the structure of the resulting operator in the following proposition.
Proposition 9 (Structure of multilayer nested commutators). Let H = T +V = j,k τ j,k A † j A k + l,m ν l,m N l N m be an interacting-electronic Hamiltonian as in (1). Then, each nested commutator H γ p+1 , · · · [H γ 2 , H γ 1 ] where H 0 = V and H 1 = T is a linear combination of fermionic chains: for some integer q ≤ p. Here, we have B x,y = l ν l,jx N l , m ν jx,m N m , ν jx,jx I and Vertices in the graph denote the indices in the summation, whereas edges represent the coefficients. We color a vertex if there is no fermionic operator corresponding to this index (due to taking nested commutators). C x,z = l ν l,kx N l , m ν kx,m N m , ν kx,kx I, or they define fermionic subchains: The definition of fermionic subchain is similar to that of the fermionic chain, except we for some x 0 ≤ q and χ x, for some x 0 ≤ q . 9 See Figure 2 for a graph illustration of this structure.
Proof. We will analyze the structure of multilayer nested commutators by induction. and In each case, we see from (77) and (78) that the result is again a fermionic chain. Specifically, commutators T, A † jq and [T, A k 1 ] increase the "length" of the current fermionic chain from q to q + 1; commutators [T, B x,y ] and [T, C x,z ] either create a fermionic subchain or give the zero operator, or they can be computed recursively when B x,y and C x,z are fermionic subchains; commutators V, A † jq and [V, A k 1 ] do not increase the length q of the current fermionic chain, but they increase the number of B x,y and C x,z by one; commutators [V, B x,y ] and [V, C x,z ] either give the zero operator, or they can be computed recursively if B x,y and C x,z are fermionic subchains.
Each application of commutation rules (77) and (78) increases the number of terms by a factor of at most 3. The nested commutator H γ p+1 , · · · [H γ 2 , H γ 1 ] contains products of at most 2(p + 1) elementary fermionic operators, giving at most 6 p+1 (p + 1)! terms in total. Meanwhile, the number of τ or ν increases by one depending on whether H γ p+2 = T or H γ p+2 = V . The claim about the number preservation can be verified directly. This completes the inductive step.
Proposition 10 (Fermionic seminorm of fermionic chain and subchain). Let H = T +V = j,k τ j,k A † j A k + l,m ν l,m N l N m be an interacting-electronic Hamiltonian (1). Then, we can bound the fermionic seminorm of the fermionic chain X in (79) as whereas fermionic subchains B x,y , C x,z in (80) can be similarly bounded as Proof. We will prove this bound using Lemma 1, Lemma 2, and Lemma 3 in a similar way as in Proposition 8. Specifically, we write X = jq A † jq D jq , where Then, Applying the operator Cauchy-Schwarz inequality (Lemma 1), we obtain Next, we write D jq = y B q,y E jq , where Invoking the Hölder-type inequality for the expectation value (Lemma 3), we get Then, We perform diagonalization using Lemma 2, obtaining Next, we write F kq = z C q,z G kq , where Invoking again the Hölder-type inequality for the expectation value (Lemma 3), we get Note that we can write We can now iterate this procedure q times to get This completes the proof of (83). Essentially the same argument can be applied to bound the fermionic seminorm of fermionic subchains. The only difference is that we have additional coefficients β x,x 0 in B x,y and respectively χ x,x 0 in C x,z . But their indices will be contracted in the x 0 th and x 0 th iteraction of the above analysis and the coefficients can then be bounded by ν max , which completes the proof of (84).
We now apply Proposition 9 to expand each nested commutator H γ p+1 , · · · [H γ 2 , H γ 1 ] into fermionic chains and use Proposition 10 to bound their fermionic seminorm. The τ factors are already bounded by their spectral-norm τ in Proposition 10. To proceed, we need to further bound each B x,y η−1 and C x,z η−1 separately. We have 11 for B x,y and similar estimates hold for C x,z . In the case where B x,y or C x,z creates a fermionic subchain, we can estimate recursively using Proposition 10. In particular, we will introduce a factor of ν max η each time a subchain is created. We know from Proposition 9 that the number of τ factors in each chain agrees with the number of H 1 = T in the nested commutator, whereas the number of ν factors coincides with the number of H 0 = V . Since the number of fermionic chains is at most 6 p p!, we obtain the bound Here

Path-counting bound on the expectation of fermionic operators
We now present the second strategy for bounding the expectation of fermionic operators, and apply it to estimate the fermionic seminorm of Trotter error. Recall from (43) that where Hence to analyze the Trotter error, it suffices to bound the fermionic seminorm H γ p+1 , · · · [H γ 2 , H γ 1 ] η . We develop a general bound on this quantity in Section 4.1 based on a path-counting technique. We then use it to analyze the simulation of d-sparse interacting electrons in Section 4.2, proving Eq. (3) of our main result Theorem 1.
It is worth noting that our approach can also be adapted to establish (9), a bound slightly weaker than our main result (2) but sufficient for our applications. See Appendix B for details.

Path-counting bound
We start by bounding the transition amplitude between any two states in terms of the expectation value. Since H γ p+1 , · · · [H γ 2 , H γ 1 ] is antihermitian, we have We now aim to bound the expectation for any |ψ η . To this end, we first expand |ψ η as: where c c c is a configuration with η electrons, and the number of ones in c c c is given by the Hamming weight |c c c| = n−1 j=0 c j . Using the notation we expand everything to get where c c c 1 , c c c 2 are configurations with η electrons, and j, k only sum over indices such that the corresponding µ γ j,k = 0 (either τ or ν depending on γ). Using the commutation relations in Equation (77) and (78), we know that the nested commutator H γ p+1 can be written as a sum of for some a ∈ {0, 1} and a sequence of elementary fermionic operators. We call each term P a fermionic path and write P H γ p+1 j p+1 k p+1 , . . . , H γ 1 j 1 k 1 to mean P is one of the terms in the expansion of the nested commutator. If the nested commutator evaluates to zero, then we consider the set to be an empty set. One possible expansion of the nested commutator is presented in Proposition 9. This allows us to make a further expansion to yield where c τ ν = τ |γ γ γ| max ν p+1−|γ γ γ| max . We use the following proposition to characterize | c c c 1 |P |c c c 2 |.

Lemma 4. For any computational basis state |c c c where c c c is a fermionic configuration, and fermionic path
we have either P |c c c is a computational basis state with some phase ±1 or P |c c c = 0.
Proof. The proof follows from a simple induction. For the base case, we have P = (−1) a without any fermionic operator, so P |c c c is a computational basis state with some phase ±1. Now we consider the three cases: P = N l P , P = A k P , or P = A † j P . By induction, we have P |c c c is a computational basis state |c c c with some phase ±1 or P |c c c = 0. The latter is trivial. For the former case, we go through the following three cases.
• If N l is applied on |c c c , we check if site-l has an electron in configuration c c c . If site-l has an electron, then N l |c c c = |c c c ; otherwise, N l |c c c = 0.
• If A k is applied on |c c c , we check if site-k has an electron in configuration c c c . If sitek has an electron, then A k |c c c will remove the site-k electron and add some phase according to the rule in Equation (29); otherwise, A k |c c c = 0.
• If A † j is applied on |c c c , we check if site-j has an electron in configuration c c c . If site-j does not have an electron, then A † j |c c c will create an electron at site-j and add some phase according to the rule in Equation (28); otherwise, A † j |c c c = 0. Therefore, P |c c c is either a computational basis state with some phase ±1 or P |c c c = 0.

Corollary 1.
We have that | c c c 1 |P |c c c 2 | is either 0 or 1. Furthermore, c c c 1 ∈S | c c c 1 |P |c c c 2 | ≤ P |c c c 2 for any set S of configurations.
Next, we define a graph G = (V, E) where the vertices V are the second-quantized configurations with η electrons, and the weighted adjacency matrix for the edges E is defined as The . . .
which is equivalent to counting the number of fermionic paths that evaluate nonzero on the initial state |c c c 2 . The last inequality follows from Corollary 1. We now introduce the following lemma which relates the maximum degree and the quadratic form (109) we wish to bound.

Lemma 5. For any real symmetric matrix w ∈ R k×k with nonnegative entries and normalized real vectors
Proof. Let u 1 be an eigenvector corresponding to the largest eigenvalue λ 1 of w with u 1 = 1. By the Rayleigh quotient theorem [38,Theorem 4.2.2], for any v ∈ R k with v = 1, where v T denotes the vector transpose of v. Consider i * = argmax j (u 1 ) j . We assume (u 1 ) i * > 0 without loss of generality, for otherwise we multiply u 1 by −1. Then, we have This concludes the proof.
Using Lemma 5, we obtain an upper bound of | X | . . .
in terms of the maximum degree of the graph G. Finally, we arrive at the following proposition by combining the above bound with Equation (115).

Counting fermionic paths for d-sparse interactions
As an illustrative example, let us consider an upper bound of max c c c deg(c c c) for electronic Hamiltonians with d-sparse interactions. We will use the commutation relations in Equation (77) and (78), restated below We start with an intuitive argument. For every q = 2, . . . , p + 1, we have where P only contains fermionic operator acting on sites j 1 , k 1 , . . . , j q−1 , k q−1 . From the commutation relations (124) and (125), we see that at least one of j q , k q must match one of the indices j 1 , k 1 , . . . , j q−1 , k q−1 . Furthermore, for every j q , there are at most d k q 's that have non-zero coefficient in τ jq,kq (for γ q = 1) or ν jq,kq (for γ q = 0). Hence, we have the following bound . . .
The n factor follows from the fact that only one index can be freely choosen between 0, . . . , n − 1. And for any pair of indices j q , k q , one of them has to match the previous indices, while the other one can only choose from the d indices under the sparsity constraint. Hence we have the d p factor in the asymptotic bound. However, this analysis can be further improved using certain properties of P . Specifically, we will show that the rightmost fermionic operator in P can be either an annihilation operator A or an occupation-number operator N . This means that, for P |c c c η to be nonzero, the rightmost fermionic operator of P must act on the η occupied sites in the configuration c c c η . Therefore, we have the bound Similarly, we have . . .
Combining with the Trotter error bound (43) and the path-counting bound (122), we obtain Finally, since 1 ≤ |γ γ γ| = p+1 q=1 γ q ≤ p, we have This sketches the proof of the scaling in Eq. (3) of Theorem 1. A rigorous proof using induction is given in Proposition 12.
• All fermionic paths P start with either N or A, but not A † (we refer to the rightmost operator as the starting point).
• All fermionic paths P have at most q + 1 elementary fermionic operators.
• The number of fermionic paths P that start with a fermionic operator acting on a specific site i is at most (2d) q q!/2. For every q > 2, we use the induction hypothesis for q − 1 to prove the desired result. If γ q = 1, then we will take another commutator with T = jq,kq τ jq,kq A † jq A kq . We can see that all fermionic paths P H γq jqkq , . . . , H γ 1 j 1 k 1 come from the expansion of Using the commutation rule [X, Furthermore, since P has at most (q − 1) + 1 = q elementary fermionic operators, the expansion of [A † jq A kq , P ] will have at most q + 1 elementary fermionic operators.
We now prove an upper bound for the number of fermionic paths starting from a specific site i. If we take the commutator of A † jq A kq with a fermionic operator that is not the starting operator in P , then the starting operator is not affected. Because of the sparsity constraint and the δ function created by the commutation relation in Eq. (124), we have created at most 2d(q − 1)× more fermionic paths starting with site i. Now if we take the commutator of A † jq A kq with the starting operator A ky (for some index k y ) in the fermionic path P , then the starting operator becomes A kq and we have an additional δ ky,jq . In this case, k q can start from any site, but there will be at most d choices of j q , hence d choices of k y . This means we have created at most 2d× more fermionic paths starting with each site. The case where N lz is the starting operator can be analyzed in a similar way. Together, we have created at most 2dq× more fermionic paths starting with each site. This leads to an upper bound of fermionic paths for each fixed starting site. The analysis for γ q = 0 proceeds in a similar way using Eq. (125). This completes the inductive step for q.
Performing the induction over q from 2 to p + 1 shows that the number of fermionic paths starting with site i is at most Because P starts with either A or N , P |c c c η would be nonzero only if the starting fermionic operator acts on one of the η occupied sites in the configuration c c c η . Hence there are at most ηO(d p+1 ) fermionic paths with non-zero P |c c c η . Finally, recall from Lemma 4 that P |c c c η is either 0 or 1. Therefore, . . .
A similar argument can be used to prove the other claimed bound . . .
The argument uses the property that all fermionic paths P end with either N or A † but not A, which again follows from the commutation relations (124) and (125). The proof is now completed.
It is worth mentioning that the path-counting approach can also be used to analyze the simulation of non-sparse electronic Hamiltonians. The resulting bound, as given by (9), is slightly weaker than Eq. (2) of Theorem 1, but suffices for our applications to be discussed in Section 6.1. See Appendix B for details and proofs.

Tightness
We have already established multiple bounds in Theorem 1 on the fermionic seminorm of the Trotter error. However, a common issue with the Trotterization algorithm is that its error estimate can be very loose for simulating specific systems. Here, we prove Theorem 2 that demonstrates the tightness of our analysis for the interacting-electronic Hamiltonian Note that we may without loss of generality assume that n is even, for otherwise we restrict to the first n − 1 spin orbitals. Comparing with the definition of the interacting-electronic model (1), we see that the coefficient matrix τ is an all-ones matrix with spectral norm τ = n, whereas ν contains an all-ones submatrix on the top left corner with max-norm ν max = 1. Our goal is to lower-bound the fermionic seminorm [T, . . .
Due to the complicated commutation relations between T and V , a direct computation of [T, . . . [T, V ]] seems technically challenging. Instead, we perform a change of basis by applying the fermionic Fourier transform This gives We also define the η-electron states for η ≤ n 2 : (142) The following proposition shows that the above choice of operators and states almost saturates the fermionic seminorm of nested commutators. (141) and | ψ η , | φ η as in (142). Then,

Proposition 13. Define T , V as in
A proof of this proposition is given in Appendix C. By rescaling the Hamiltonian constructed in (141), we can demonstrate the tightness of our bound as follows. For any s, w > 0, we define the rescaled Hamiltonian Comparing with the definition of the interacting-electronic model (1), we see that τ = s and ν max = w. The above proposition then shows that where we have used the unitary invariance of the fermionic seminorm in the first equality. This establishes the first claimed bound in (5) of Theorem 2. Note that a similar example can be constructed to demonstrate the tightness of our bound for simulating sparse electronic Hamiltonians. Specifically, suppose we have u, w > 0 and positive integer 2 ≤ d ≤ η ≤ n 2 . 13 We may assume without loss of generality that d is even, for otherwise we use d − 1. We then define 13 The special case d = 1 can be handled separately by choosing T = uA † 0 A1 + uA † 1 A0 and V = wN0.
Comparing with the definition of the interacting-electronic model (1), we see that τ max = u and ν max = w. We also perform a fermionic Fourier transform to define T and V , but only to the first d spin orbitals Then, a similar calculation shows that This establishes the first claimed bound in (6) of Theorem 2.

Lower-bounding
Recall from the previous section that we have constructed the electronic Hamiltonian (139) to prove the tightness of our bound. Comparing to the definition of the interactingelectronic model (1), we see that the coefficient matrix τ has spectral norm τ = n, whereas coefficient matrix ν has max-norm ν max = 1. Our goal in this subsection is to lower-bound the fermionic seminorm To this end, we define the η-electron states for η ≤ n 2 : Similar to the previous subsection, we may assume that n is even. We have the following proposition showing that the fermionic seminorm of nested commutators is nearly attained.

Proposition 14.
Define T , V as in (139) and |ψ η , |φ η as in (149). Then, A proof of this proposition is given in Appendix D. By rescaling the Hamiltonian constructed in (141), we can demonstrate the tightness of our bound as follows. For any s, w > 0, we define the rescaled Hamiltonian as in (144). Comparing with the definition of the interacting-electronic model (1), we see that τ = s and ν max = w. The above proposition then shows that This establishes the second claimed bound in (5) of Theorem 2. Note that a similar example can be constructed to demonstrate the tightness of our bound for simulating sparse electronic Hamiltonians. Specifically, for u, w > 0 and integer 2 ≤ d ≤ η ≤ n 2 , 14 we define the electronic Hamiltonian as in (146). Comparing with the definition of the interacting-electronic model (1), we see that τ max = u and ν max = w. A similar calculation then shows that This proves the second claimed bound in (6) of Theorem 2.

Applications
The class of interacting-electronic Hamiltonians (1) encompasses various quantum systems arising in physics and chemistry, for which the performance of digital quantum simulation can be improved using our result. As for illustration, we consider improving quantum simulation of the plane-wave-basis electronic structure in Section 6.1 and the Fermi-Hubbard model in Section 6.2.

Plane-wave-basis electronic structure
Simulating the electronic-structure Hamiltonians is one of the most promising applications of digital quantum computers. Recall that in the second-quantized plane-wave basis, such a Hamiltonian takes the form where ω is the volume of the computational cell, κ µ = 2πµ/ω 1/3 are n vectors of planewave frequencies, µ are three-dimensional vectors of integers with elements in the interval [−n 1/3 , n 1/3 ], r j are the positions of electrons, ζ ι are nuclear charges, and r ι are the nuclear coordinates. We further rewrite the second term as which is valid since we estimate the simulation error within the η-electron manifold. Comparing with the definition of interacting-electronic model (1), we see that (155) 14 Similar to above, the special case d = 1 can be handled using T = uA † 0 A1 + uA † 1 A0 and V = wN0.
To proceed, we need to bound the spectral norm τ and the max-norm ν max of the coefficient matrices. We have where the first equality follows from [6, Eq. (F10)], the second equality follows from [6, Eq. (F11)-(F13)] and the third equality follows from [6, Eq. (F7) and (F9)]. We also consider a constant system density η = O (ω) following the setting of [6]. Applying Theorem 1, we find that a pth-order formula S p (t) can approximate the evolution of electronic-structure Hamiltonian with Trotter error This approximation is accurate for sufficiently small t. To evolve for a longer time, we divide the evolution into r steps and use S p (t/r) within each step, which gives an approximation with error (158) To simulate with accuracy , it suffices to choose Note that this can also be achieved using the weaker bound (9) from path counting, since both τ and n τ max have the same asymptotic scaling.
To simplify our discussion, we consider digital quantum simulation with constant time and accuracy, obtaining We further implement each Trotter step using the approach of [54,Sect. 5], and obtain a quantum circuit with gate complexity which implies by choosing the order p sufficiently large. Up to a negligible factor n o (1) , this gate complexity improves the best previous result of the electronic-structure simulation in the second-quantized plane-wave basis. This is because our approach improves the performance of digital quantum simulation by simultaneously exploiting commutativity of the Hamiltonian and prior knowledge of the initial state, whereas previous results were only able to employ at most one of these information. Indeed, previous work [6, Appendix G] gave a Trotterization with error bound Their approach used the initial-state information by computing the Trotter error within the η-electron manifold, but the commutativity of the Hamiltonian was ignored, giving a simulation with gate count n 5/3 η 1/3 + n 4/3 η 5/3 n o(1) worse than our result. On the other hand, the work [27,Proposition F.4] used commutativity of the Hamiltonian to show and gave a simulation with complexity n 7/3 η 1/3 n o(1) , whereas Ref. [54] gave an interactionpicture approach with cost O n 8/3 η 2/3 polylog(n) . Our new result matches these when η and n are comparable to each other, but can be much more efficient in the regime where η is much smaller than n.
Interestingly, our asymptotic result remains conditionally advantageous even when compared with the first-quantized simulations. There, the best previous approach is the interaction-picture approach [7] with gate complexity O n 1/3 η 8/3 polylog(n) , larger than our new complexity when n = O η 2−2/(p+1) polylog(n) . A related approach was described in [7] based on qubitization, which has a similar performance comparison with our result. 15 See Table 1 for details.
We mention however that there is one caveat when ignoring the factor n o(1) in our above discussion. This is achieved by choosing the order p of Trotterization sufficiently large, which can result in a gate complexity with an unrealistically large prefactor depending on the definition of higher-order formulas. 16 Nevertheless, recent work suggests that Trotterization remains advantageous for simulating the plane-wave-basis electronic structure even with a low-order formula [42], to which our paper provides new theoretical insights.

Fermi-Hubbard model
We also consider applications of our result to simulations of the Fermi-Hubbard Hamiltonian, which models many important properties of interacting electrons. This Hamiltonian is defined as where j, k denotes a summation over nearest-neighbor lattice sites and σ ∈ {0, 1}. We note that this Hamiltonian can be represented in terms of a sparse interacted Hamiltonian. Indeed, in the one-dimensional case, we have where j = 0, 1, . . . , n − 1 and σ = 0, 1. Comparing with the definition of interactingelectronic model (1), we see that so the coefficient matrices τ and ν are 2-sparse. Similar analysis holds for the higherdimensional Fermi-Hubbard model, with the sparsity d = 2 m where m is the dimensionality of the lattice.
We can therefore apply Theorem 1 to conclude that a pth-order formula S p (t) approximates the evolution of Fermi-Hubbard Hamiltonian with Trotter error assuming s, v, and m are constant. For r steps of Trotterization, we apply the triangle inequality to get To simulate with constant time and accuracy, it thus suffices to choose giving gate complexity 17 The Fermi-Hubbard model only contains nearest-neighbor interactions and, according to [22], can be near optimally simulated with O n 1+1/p gates. On the other hand, recent work [28] shows that Trotterization algorithm has gate complexity O nη 1+1/p when restricted to the η-electron manifold. Our result improves over those previous work by combining the sparsity of interactions, commutativity of the Hamiltonian and information about the initial state.

Discussion
We have given improved quantum simulations of a class of interacting electrons using Trotterization, by simultaneously exploiting commutativity of the Hamiltonian, sparsity of interactions, and prior knowledge of the initial state. We identified applications to simulating the plane-wave-basis electronic structure, improving the best previous result in second quantization up to a negligible factor while conditionally outperforming the first-quantized simulation. We obtained further speedups when the electronic Hamiltonian has d-sparse interactions, which gave faster Trotterization of the Fermi-Hubbard model. We constructed concrete electronic systems for which our bounds are almost saturated, providing a provable guarantee on the tightness of our analysis. 17 By exploiting locality, we can perform e −itV with O (n) gates. An additional logarithmic factor will be introduced if we implement e −itT using the fermionic Fourier transform. However, this can be avoided by further decomposing T = T odd + Teven as in [22], so that H = T odd + Teven + V . The analysis of Trotter error proceeds in a similar way as in Section 4, which establishes the claimed gate complexity of O nη 1/p . Our focus has been on the asymptotic performance of digital quantum simulation throughout this paper. However, we believe that the techniques we have developed can also be used to give quantum simulations with low constant-prefactor overhead, for instance, through a more careful application of our Proposition 9, Proposition 10, Proposition 11 and Proposition 12. Such improvements would especially benefit the simulation of planewave-basis electronic structure, where many pairs of Hamiltonian terms commute and the number of electrons can be significantly smaller than the spin-orbital number. Existing numerical studies almost exclusively used the second-order Suzuki formula [5,19,42,65,84] and their results did not fully leverage the commutativity of the Hamiltonian and the initial-state knowledge, which may then be improved by the techniques presented here using general Trotterization schemes.
Our analysis is applicable to a class of electronic Hamiltonians of the form H = j,k τ j,k A † j A k + l,m ν l,m N l N m . By imposing further constraints on the coefficients, we may somewhat sacrifice this generality but instead get further improvement on the simulation performance. One possibility is to consider the subclass of systems that are translationinvariant, i.e., τ j,k = τ j+q,k+q and ν l,m = ν l+q,m+q . This translational invariance is used in the circuit implementations for both our applications (electronic-structure Hamiltonians and Fermi-Hubbard model), but is nevertheless ignored in the proof of our upper bounds (Theorem 1) and tightness result (Theorem 2). By incorporating additional features of the Hamiltonian such as translational invariance, it is plausible that our current complexity estimate can be further improved.
A natural problem that has yet to be addressed here is the simulation of electronicstructure Hamiltonians in a more compact molecular basis. Such Hamiltonians typically take the form more complex than the electronic model (1) considered here. In this case, the exponentials of the two-body terms j,k,l,m h j,k,l,m A † j A k A † l A m do not have a convenient circuit implementation and our current approach is not directly applicable. This motivates further developments of hybrid quantum simulation, in which Trotterization is combined with more advanced quantum algorithms to speed up digital quantum simulation. We leave a detailed study of such problems as a subject for future work.
More generally, we could consider digital quantum simulations of other types of physical systems, such as bosonic systems [70] or fermion-boson interacting systems [72]. We hope our techniques could offer insights to such problems and find further applications in digital quantum simulation beyond what have been discussed here.
For the third statement of Proposition 8, we let X = j,k,l τ j,k ν l,k A † j N l A k and compute (172) Applying the operator Cauchy-Schwarz inequality (Lemma 1) similarly as in (63), We now perform diagonalization using Lemma 2, obtaining Using the Hölder-type inequality for the expectation value (Lemma 3), we have This completes the proof of the third statement of Proposition 8.
For the fourth statement, we let X = j,k,m τ j,k ν j,m A † j N m A k and compute Applying the operator Cauchy-Schwarz inequality (Lemma 1), We now use the Hölder-type inequality for the expectation value (Lemma 3) to get The second fermionic seminorm can be directly bounded as whereas the first seminorm can be bounded using diagonalization (Lemma 2) This completes the proof of the fourth statement of Proposition 8. For the fifth statement, we let X = j,k τ j,k ν j,j A † j A k and compute (180) Applying the operator Cauchy-Schwarz inequality (Lemma 1), (181) We now use the Hölder-type inequality for the expectation value (Lemma 3) to get The second fermionic seminorm can be directly bounded by ν 2 max , whereas we perform diagonalization to the first seminorm (Lemma 2): This completes the proof of the fifth statement of Proposition 8. For the sixth statement, we let X = j,k,l τ j,k ν l,j A † j N l A k and compute Applying the operator Cauchy-Schwarz inequality (Lemma 1), We now use the Hölder-type inequality for the expectation value (Lemma 3) to get The second fermionic seminorm can be directly bounded by ν 2 max η 2 , whereas we perform diagonalization to the first seminorm (Lemma 2): This completes the proof of the sixth statement of Proposition 8.

B Counting fermionic paths for non-sparse interactions
In this appendix, we use the path-counting technique to prove (9) for non-sparse interacting electrons. We will make use of the following commutation relations which are slightly different from the ones used before. These relations can be derived in a similar way as in Equation (77) and (78) Proof. We will prove the following claims by induction on q = 2, . . . , p + 1.
• All fermionic paths P are products of A † i A j and N k .
• All fermionic paths P have at most q + 1 elementary fermionic operators.
• The number of fermionic paths P that start with a fermionic operator acting on a specific site i is at most The base case q = 2 can be easily verified by noting that we only need to consider [T, V ] or [V, T ]. For every site i, there are at most 6nη fermionic paths starting with this site, all of which are products of A † i A j and N k . This is because there are at most three summation indices. The rightmost index must be equal to i and the indices for N k , A j have at most η choices, while the remaining index has n possible choices, giving a total of nη choices. The additional factor of 6 comes from the number of different expansion terms in Equation (188), (189), (190). Furthermore, every fermionic path consists of at most 3 fermionic operators. These established the claims for the base case q = 2.
For every q > 2, we now use the induction hypothesis for q − 1 to prove the claims for q. If γ q = 1, then we take another commutator with T = jq,kq τ jq,kq A † jq A kq . We can see that all fermionic paths P H γq jqkq , . . . , H γ 1 j 1 k 1 come from the expansion of Using the commutation rule [X, . . Y κ , we show that the claims hold for q as follows. When we take the commutation of A † jq A kq with A † j or A k , we know from (189) that one free index will be introduced, resulting in an additional factor of n. When we take the commutation of A † jq A kq with N l , we will remove the fermionic operator N l and replace it with A † jq A kq , which removes a factor of η and adds an additional factor of nη. Additionally, there are at most (q − 1) + 1 fermionic operators in P . Hence the number of fermionic paths P that start with a fermionic operator acting on site i is at most Furthermore, in both cases, we add at most one additional fermionic operator. Therefore, all fermionic paths will have at most (q − 1) + 1 + 1 = q + 1 fermionic operators. And [A † jq A kq , P ] remains a product of A † i A j and N k . The inductive step for γ q = 0 follows from a similar argument. The first two claims can be directly verified. For the last claim, we proceed in a slightly different way as follows. When we take the commutator of N † jq N kq with A † j or A k , we will add N jq or N kq to the fermionic path, which results in an additional factor of η. When we take the commutator of N † jq N kq with N k , the commutator is equal to zero. Hence the number of fermionic paths that start with a fermionic operator acting on site i is at most We have thus shown that the claims hold for q. Performing the induction on q from 2 to p + 1 shows that the number of fermionic paths starting with site i is at most Because each fermionic path P is a product of A † i A j and N k , P |c c c η would be nonzero only if the rightmost fermionic operator acts on one of the η occupied sites in the configuration c c c η . Hence there are at most ηO(n p+1 q=1 γq η p+1 q=1 (1−γq) ) fermionic paths with non-zero P |c c c η . Finally, recall from Lemma 4 that P |c c c η is either 0 or 1. Therefore, we have The other bound . . .
can be similarly proved using the fact that the leftmost fermionic operator of P † must act on one of the η occupied sites in any fixed configuration.
(206) We will see that this is the dominant contribution to the effective commutator that is at least Ω (nη).

The second group contains terms
which does not dominate the result scaling.