Quantum simulation costs for Suzuki-Trotter decomposition of quantum many-body lattice models

Quantum computers offer the potential to efficiently simulate the dynamics of quantum systems, a task whose difficulty scales exponentially with system size on classical devices. To assess the potential for near-term quantum computers to simulate many-body systems we develop a formalism to straightforwardly compute bounds on the number of Trotter steps needed to accurately simulate the time evolution of fermionic lattice models based on the first-order commutator scaling. We apply this formalism to two closely related many-body models prominent in condensed matter physics, the Hubbard and t-J models. We find that, while a naive comparison of the Trotter depth first seems to favor the Hubbard model, careful consideration of the model parameters and the allowable error for accurate simulation leads to a substantial advantage in favor of the t-J model. These results and formalism set the stage for significant improvements in quantum simulation costs.


I. INTRODUCTION
Despite rapid growth in accessible computational resources, simulating quantum systems on classical computers has remained an elusive challenge. This difficulty arises from the exponential scaling of complexity with system size necessary to accurately model quantum systems. Quantum computers, at least in theory, present a solution to this problem. By relying on hardware that is itself quantum in nature, quantum computers can be used to efficiently simulate the dynamics of other quantum systems [1][2][3]. Quantum simulation is poised to be a potent tool, with potential applications in a diverse array of areas, including high-energy physics and biology [4]. In particular, the quantum simulation of electronic structure promises to have particular significance for quantum chemistry [5][6][7] and materials science [8].
A critical step in any quantum simulation problem is to translate the Hamiltonian that generates the dynamics for the simulated system into a qubit Hamiltonian that can be implemented on a quantum computer. For the case of electronic structure, this translation can be accomplished using the celebrated Jordan-Wigner transformation [9], which maps the fermionic creation and annihilation operators to the Pauli operators. The cost of this mapping in terms of circuit design can be quantified by the "Pauli depth," the maximum length of Pauli operators needed to implement any given term of the qubit-mapped Hamiltonian. A notable challenge in maintaining a low Pauli depth for fermionic Hamiltonains is the need to account for the underlying anticommutivity of the fermionic operators. To do so, the parity for each orbital must be stored using a string of Pauli Z operators [10]. A consequence of the parity information being stored nonlocally in the "Jordan-Wigner * myersn1@vt.edu string" is that each fermionic operator is, in the most general case, mapped to O(N ) Pauli operators, where N is the total qubit number [7] (unless non-local gates are used [11]). This non-locality is an inevitable consequence of the fermion sign problem, which is responsible for the classification of such fermionic lattice models as NP-hard [12]. With this in mind, in this work we focus on the Jordan-Wigner mapping as our primary tool for encoding fermionic Hamiltonians in qubit form.
Along with encoding the Hamiltonian in qubit form, quantum simulation also requires a method of implementing the time-evolution of the system. This can be done by decomposing the time evolution operator into a series of local gate operations by means of the Suzuki-Trotter expansion [13,14], a method known as "Trotterization" [15]. As the expansion is only exact in the limit of an infinite number of expansion steps, practical applications require truncating at a finite number of steps, and consequently introducing a truncation error, . This error, along with the desired evolution time, τ , are the model independent parameters that determine the required number of Trotter steps for the simulation, r. Thus, , τ , and r are the three key quantities for benchmarking Trotterized quantum simulation. In general, the choice of and τ will depend on the simulation observable of interest and the model energy scale. For example, an accurate quantum simulation to determine the energy gap [16] of a particular many-body model will require an significantly smaller than the typical energy spacing. We therefore use the dimensionless quantity r /τ 2 as a model independent measure of the Trotter cost.
Together, the Pauli depth and maximum number of Trotter steps, r, provide two significant benchmarks for the computational cost of quantum simulation. This, in turn, has led to significant study in how to optimize these costs, such as by circuit modifications that lead to the cancellation of Jordan-Wigner strings and reordering the terms in the Suzuki-Trotter expansion to reduce the error in each step [11].
In this work, rather than seeking to optimize these costs for a particular model, we instead compare the costs between two prominent models for quantum condensed matter, namely the Hubbard and t-J models. We avoid performing detailed gate counts or circuit depth estimations, as has been previously explored [11,[17][18][19]. Our goal is to first answer the question of which condensed matter models are most amenable to simulation independent of any specific implementation or choice of simulation parameters, before proceeding with any detailed minimization of computational costs. In this way we can first demonstrate that any advantage either model shows in computational cost is inherent to the structure of the model itself.
To date, the Hubbard model has served as one of the foremost models targeted for quantum simulation. Numerous approaches, have been pursued, including experimental analog simulations using using atoms confined in optical lattices [20][21][22][23][24][25] as well as quantum dots [26] and NMR systems [27]. The Hubbard model has also received focus as test-bed model for simulation on nearterm quantum computers [28][29][30][31]. With this in mind, it is important to identify whether or not the Hubbard model is truly the most efficient model, in terms of simulation resources, that captures the physics of interest.
The t-J model, which arises out of a perturbative treatment of the Hubbard model in the limit where the on-site interaction dominates over the hopping energy [32][33][34][35], is of particular interest due to its potential to model hightemperature superconductivity [35][36][37][38][39]. An important feature of the t-J model Hamiltonian is the presence of Gutzwiller projection operators that eliminate the possibility of any site being doubly-occupied [38,[40][41][42]. This significantly reduces the Hilbert space of the t-J model in comparison to the more general Hubbard model, which motivates the supposition that the t-J model should be more tractable to simulate on near-term quantum computers. To confirm whether this hypothesis bears out, we apply the Jordan-Wigner transformation along with optimized bounds on the Trotter error scaling [43] to benchmark the simulation costs for two-dimensional Hubbard and t-J models on a square lattice in terms of the Pauli and Trotter depths.
There remains debate about the presence of superconducting behavior in the ground state of the doped t-J model. Unbiased classical simulation using exact diagonlization techniques to probe this issue have been limited to 20 sites for the case of the 2D t-J model [44]. Quantum simulation could help settle the debate. With the growing, but still very limited, qubit resources available to existing quantum computers, determining the model that makes the most efficient use of those resources while maintaining the potential to reveal impactful new physics is of significant importance to both the quantum computing and condensed matter communities. To address this question, we develop a formalism for applying the optimized commutator bound [43] to fermionic lattice models that yields analytical expressions for the Trotter depth in terms of the model parameters. We then apply this formalism to demonstrate that the bound on the Trotter depth is significantly lower for the t-J model in comparison to the Hubbard model in the parameter regime of validity. These results suggest that the t-J model should be a prominent candidate for simulation on near-term quantum computers.
In Section II we begin by reformulating the upper bound on Trotter depth given in Ref. [43] for the case of square lattice models with both open and periodic boundary conditions. Then, in Section III we review the 2D Hubbard and t-J models including a mapping of the spinful models onto bipartite spinless lattices, before we apply the Jordan-Wigner transformation in Section IV. From the Jordan-Wigner transformed Hamiltonians we can immediately evaluate the minimum Pauli depth for each model. In Section V we apply our expanded form of the Trotter bound in order to compare the maximum Trotter depth for the t-J and Hubbard models in both one and two dimensions. Notably, we show that the Trotter depth for both models scales linearly with the system size in both 1D and 2D. Furthermore, we find that, within the parameter regimes where a valid comparison can be made, the t-J model is significantly less costly to simulate. Finally, in Section VI we summarize our results and offer some perspectives on their significance for near-term quantum simulation.

II. BOUNDING TROTTER DEPTH
For a Hamiltonian that can be decomposed as H = Γ γ H γ it is well established that the bound on the Trotter depth scales with Γ γ ||H γ || [45,46], where ||A|| denotes the spectral norm of operator A. However, this bound is typically a significant overestimate, as it neglects to account for any commutativity within the terms of the Hamiltonian. In Ref. [43] it was shown that this bound can be improved to, This bound is a powerful improvement and can be orders of magnitude smaller than bounds arising from just the norm of the Hamiltonian terms. Note that Eq. (1) is specific to a first order product formula. The bound derived in Ref. [43] is further generalized to any order, but for the purposes of this manuscript we consider only the first order result, and leave generalizations to higher orders to be explored in future work. Before applying this bound to any specific models, let us first consider it in more detail. First, we note that the bound is primarily dependent on the spectral norms of commutators of the individual Hamiltonian terms. For Jordan-Wigner transformed fermionic Hamiltonians, the Hamiltonian will always consist of a sum over products of Pauli matrices, along with some prefactors. The spectral norm of any product of Pauli matrices is always one, as demonstrated in Appendix A. Thus, if each H γ consists of only products of Pauli matrices, || [H γ1 , H γ2 ] || will simply be given by the absolute value of the product of the respective prefactors of H γ1 and H γ2 .
Fermionic lattice Hamiltonians are typically expressed as a sum over lattice sites, H = i,j H i,j . In the following we assume translationally invariant lattices, up to boundary terms. As stated above, after applying the Jordan-Wigner transformation each single site term, H i,j , can be decomposed into D subterms, Let us consider some general properties of each A i1,j1 i2,j2 . i2,j1 The first property follows from We can verify the second is true by expanding the summations as, Ny,j1 .
(6) We see that for each term A i1,j1 i2,j2 in Eq. (5) there will be a corresponding A i2,j2 i1,j1 term in Eq. (6). Finally, we also note that, since each single-site term H i,j in a lattice Hamiltonian is identical to each other single-site term, up to a change in indices, we have A i,j i+a,j+b = A k,l k+a,l+b . This property is in essence a statement that the rectangular lattice is invariant under horizontal and vertical translations, up to boundary terms. Using these properties we can expand Eq. (1) as, This expression constitutes the centerpiece of our formalism. We pause here for a clarification of notation.
Using the condition A i,j i+a,j+b = A k,l k+a,l+b we see that all A i,j i,j are equal, A 1,1 1,1 = A 1,2 1,2 = A 2,1 2,1 = .... Thus we have replaced the summation that would appear in Eq. (7) with a product over the total number of A i,j i,j terms, i.e.
For simplicity's sake, we have chosen i = j = 1, but we stress that Eq. (7) is true for any choice of i and j. However, to account for the lattice boundaries, for any choice other than i = j = 1, all j indices should be considered modulo N x and all i indices should be considered modulo N y . We emphasize further that Eq. (7) is completely general in regards to hopping, including all-to-all interactions. As we will demonstrate in subsequent sections, modeldependent rules and specification of boundary conditions can considerably simplify this expression.
To ensure that we have expanded Eq. (1) properly, we can check to make sure Eq. (7) has the same number of terms. We see that Eq. (1) has N 2 x N 2 y total terms. Using the summation identity, we find that Eq. (7) has, total terms. Therefore, as expected, Eqs. (1) and (7) have an equal number of terms.

III. THE HUBBARD AND t-J MODELS
The 2D Hubbard model is [41], where t is the hopping energy, U is the on-site interaction energy, and c † j,s and c j,s are the fermionic creation and annihilation operators obeying the anticommutation relations, Note that n ij,s = c † ij,s c ij,s and n ij = n ij,↑ +n ij,↓ . We consider a rectangular lattice with N = N x N y sites. Each pair of indices (i, k = 1, . . . , N y and j, l = 1, . . . , N x ) denote the x, y position of the site in the lattice, respectively, as illustrated in Fig. 2. We also assume a closed system with a fixed number of particles and thus do not include a chemical potential term in the Hamiltonian.
The t-J model can be derived from the Hubbard model by working in the regime where U/t 1 and considering only the lower energy subspace of unoccupied and singly-occupied sites [38,41]. The Hamiltonian for the t-J model is, where J ≡ 4t 2 /U . Here the P are the Gutzwiller projection operators that ensure the system does not admit any doubly-occupied sites. In Fig. 1 we provide a schematic representation demonstrating the conditions under which the t-J model emerges from a perturbative treatment of the Hubbard model.
For later convenience when applying the Jordan-Wigner transformation it will be useful to express the t-J Hamiltonian entirely in terms of the creation, annihilation, and number operators, The full derivation for this expression is provided in Appendix B.
Examining Eq. (13) we see that the local projection operators leave the Heisenberg terms unchanged. This is expected, as the spin-spin interaction only arises under the condition of single occupancy. A natural question is then whether the projection operators need to be applied to the J terms at all. Assuming an initial state with no double occupancy, the existence of the projection operators on the hopping terms is sufficient to keep the state restricted to the subspace consisting of only unoccupied and singly-occupied sites. However, simulating the dynamics of the t-J model on a real quantum device will inevitably introduce errors, some of which may lead to the system straying out of the restricted Hilbert space. With this in mind, we leave the projection operators in place on the J terms as a limited form of protection against this class of errors.
A deeper degree of protection that would also guard against initial states with doubly occupied sites can be implemented using global projection operators that also project the implicit identity operators that act on every other site in the lattice for each term in Eq. (13). However, these global projection operators are highly nonlocal and come at a prohibitive cost in terms of Pauli depth. Thus, for the purpose of this analysis, we assume the initial state used for any simulation of the t-J model does not contain any double occupancy.

IV. THE JORDAN-WIGNER TRANSFORMATION
It is well established that fermionic creation and annihilation operators can be mapped to Pauli (qubit) operators by means of the Jordan-Wigner transformation. In 1D this transformation takes the form of [9,10], Here σ + ≡ (σ x + iσ y )/2, σ − ≡ (σ x − iσ y )/2 with σ x , σ y , and σ z being the usual Pauli matrices and σ 0 being the identity matrix.
In order to capture the full Jordan-Wigner string in two dimensions, as illustrated in Fig. 2, each fermionic operator must transform as follows, Note that this implementation of the 2D Jordan-Wigner transformation matches with the result derived in Ref. [47]. Using Eq. (15) we have, Without loss of generality we assume i < k and j < l. Equation (16) then becomes, We note that Eq. (17) has three distinct "segments" of Jordan-Wigner strings. We will refer to these segments as string A, and string C, respectively. From a quick inspection of the bounds on the summations for each string we see that together the three segments account for all the sites between (i, j) and (k, l). This is illustrated in Fig. 3.
We are now ready to apply the Jordan-Wigner trans-formation to the Hubbard and t-J Hamiltonians. The first step in doing so is to remove the spin dependence from the Hamiltonians. We can do this by mapping from a spinful lattice to a bipartite spinless lattice. This technique is illustrated graphically in Fig. 4.
Making the substitutions, we arrive at for the Hubbard model and, for the t-J model. Applying the Jordan-Wigner transformation as given in Eq. (14) and simplifying using as well as, we arrive at, For the Hubbard model. Similarly, for the t-J model we find, The t-J model arises from the low energy Hubbard model in the limit U/t 1. As a check for the consistency of the Jordan-Wigner transformations, we numerically calculate the energy spectra of a four-site, fourparticle 1D Hubbard and t-J model using the Pauli representation we derived in section IV. We plot the resulting spectra for U = 10, t = 0.1 and J = 4t 2 /U = 0.004 in Fig. 5. We see that, as expected, the t-J energy spectrum overlaps with that of the low energy states of the Hubbard spectrum.
As an additional verification, we also calculate the energy eigenstates for the Hubbard and t-J models under the same parameters using the ALPS (Algorithms and Libraries for Physics Simulations) software package [48]. To within numerical precision, we obtain the same spectrum as from the exact diagonlization of our Jordan-Wigner transformed Hamiltonians.
Comparing Eqs. (26) and (27) we see that the maximum Pauli depth of the t-dependent terms is identical between the two models. Both contain at least one term that consists of a product of N x Pauli operators. In contrast, the maximum Pauli depth of the U -dependent term in the Hubbard model is only two, half that of the Jdependent term in the t-J model.

V. t-J AND HUBBARD TROTTER DEPTH
With the full Jordan-Wigner transformed Hamiltonians, we can now apply Eq. (7) to bound the number of Trotter steps necessary to simulate the time evolution of the Hubbard and t-J models. Note that, for completeness, we examine the 1-norm scaling in Appendix C. This comparison emphasizes the improvement offered by the commutator scaling, particularly for the Hubbard model, which can be over an order of magnitude lower, even for small system sizes. In the following analysis we also consider specifically the case of open boundary conditions. In Appendix D we extend these results to the case of periodic boundary conditions. We begin by considering some common features of both models. As both the Hubbard and t-J models are restricted to only nearest neighbor hopping, we have A i1,j1 i2,j2 = 0 for all (i 2 , j 2 ) = (i 1 + 1, 1 ≤ j 2 ≤ j 1 ) and for all (i 2 , j 2 ) = (i 1 , j 1 ≤ j 2 ≤ N x ). These conditions account for the segments of Jordan-Wigner string that are not shared between site (i, j) and its vertical neighbor (i + 1, j), as discussed in Section IV. The first set of nonzero A i1,j1 i2,j2 accounts for the segment of Jordan-Wigner string stretching from site (i, j) rightward to site (i, N x ), and the second accounts for the segment stretching leftward from site (i + 1, 1) to site (i + 1, j).
Furthermore, by inspecting Eqs. (26) and (27) we find that all A i,j i,j+p>1 and A i+1,j i,j+p>1 vanish. This is due to the fact that the only Pauli operators with the same indices H i,j will share with H i,j+p>1 arise from the Jordan-Wigner strings, which consist of only σ z operators. Thus, each commutator will take the form of either, all of which are zero by the properties of the Pauli operators.
Thus, for the case of the Hubbard and t-J models, Eq. (7) simplifies to, We reiterate here that, like Eq. (7), this equation is valid for any choice of i and j, and that for any choice other than i = j = 1 all j indices should be considered modulo N x and all i indices modulo N y .

A. Hubbard
In order to get an expression for Eq. (29) in terms of our model parameters (t and U ) we need to compute each A i2,j2 i1,j1 . In Table I we break up each single site term in the Hubbard Hamiltonian, Eq. (26), into terms that consist only of products of Pauli operators, along with some common prefactors. Then using, where each H δ i,j is listed in Table I, we calculate each A i2,j2 i1,j1 . This process is considerably simplified since each H δ i,j consists of a tensor product of Pauli operators. Since the norm of a tensor product of Pauli operators is always one, we immediately know each A i2,j2 i1,j1 will be of the form, where a H , b H and c H are positive integers determined by the number of nonzero commutators in A i2,j2 i1,j1 with the corresponding prefactor, and the factor of 2 arises from the Pauli commutation relations.
Inspecting the terms in Table I We can repeat the same process for A i,j i,j+1 and A i,j i+1,j where in both cases we find 16 non-zero commutators, 8 with the prefactor t 2 /4 and 8 with |tU |/8. Thus, Finally, we see that for A i+1,j i,j+1 there are only 4 nonzero commutators, all corresponding to t 2 /4. Thus, We can now combine all these results with Eq. (29) to yield, Here we note several important features of this result. First, we see that the bound scales quadratically with t and linearly with U . However, even more notably, it scales linearly with the total number of lattice sites, N ≡ N x N y , despite the presence of the lattice spanning Jordan-Wigner strings in Eq. (26). This linear scaling arises due to the fact that all commutators between nonnearest-neighbor sites will be of the form given in Eq. (28) and will thus vanish. In Fig. 6 we plot the functional form of the Trotter depth from Eq. (35) as a function of t, U , and N for the case of a 6 by 6 lattice. As the choices for values of and τ will generally depend on the observables of interest for the simulation, we quantify the overall Trotter cost in terms of the problem-independent quantity r /τ 2 . For comparison, we also include the Trotter depth for the case of periodic boundary conditions (Appendix D) and for a 1D Hubbard model with an equivalent number of lattice sites (Appendix F). Here we see the quadratic scaling with t and linear scaling with U and N illustrated for both 1D and 2D models. Significantly, the difference in cost between 1D and 2D simulations is less than a factor of ten for small system sizes, but still larger than what can be efficiently simulated classically.

B. t-J
We can follow an exactly analogous process for the t-J model. In Table II we break up each single site term in the t-J Hamiltonian, Eq. (27), into terms that consist only of products of Pauli operators, along with some common prefactors. Then using, where each H δ i,j is listed in Table II, we calculate each A i2,j2 i1,j1 in Eq. (29). For the t-J model, each A i2,j2 i1,j1 will be of the form, Unlike the Hubbard model, neither a tJ , b tJ , or c tJ will always be zero, since the J terms arise from a second-order perturbation that mixes both hopping and interaction. Let us again consider each A i1,j1 i2,j2 individually. For A i,j i,j , of the 4096 commutators in (37), 1600 are non-zero. Of these 1600 commutators, 384 of them correspond to t 2 /32, 1024 to |tJ|/128, and 192 to J 2 /256. Thus we have a tJ = 384, b tJ = 1024, and c tJ = 192 forA i,j i,j . Combining these coefficients with Eq. (37) we arrive at, For A i,j i,j+1 and A i,j i+1,j we find 960 nonzero commutators, with 256 corresponding to t 2 /32, 512 to |tJ|/128, and 192 to J 2 /256. Thus, For A i+1,j i,j+1 there are only 480 nonzero commutators, with 128 corresponding to t 2 /32, 256 to |tJ|/128, and 96 to J 2 /256. Thus, Combining these results with Eq. (29), we find that the functional form of the commutator bound on the Trotter depth for the 2D t-J model is, As in the case of the Hubbard model, we see that the bound on Trotter depth for the 2D t-J model scales linearly with N ≡ N x N y and quadratically with t. However, while the Hubbard model scaled linearly with U , the t-J model scales quadratically with J. This is expected, as the J term in the t-J model arises from a second-order perturbation that mixes both the hopping and interaction.
In Fig. 7 we plot the functional form of the Trotter depth from Eq. (41) as a function of t, J, and N , again for a 6 by 6 lattice. As before, we include the Trotter depth for the case of periodic boundary conditions (and for a 1D t-J model with an equivalent number of lattice sites. As we observed for the Hubbard model, we see that the Trotter cost for both the 1D and 2D models scales linearly with N , with less than a factor of ten difference between them at N = 36.

C. Comparison: Hubbard vs t-J
With the functional form of the Trotter bound for both Hubbard and t-J, we can now make a direct comparison between the models. As we do so, we must be careful in making sure our comparison is fair, as the relevant parameters in the Hubbard model are t and U while in the t-J model the relevant parameters are t and J. In Fig. 8a we compare the two bounds treating t as fixed and U and J as independent parameters that are varied over the same range of values. In this case, we see that the bound on the Hubbard model is significantly lower than the t-J model.
However, this comparison is somewhat disingenuous, as we know that U and J are not independent, but are related by J ≡ 4t 2 /U . Furthermore, the t-J model is a valid approximation to the Hubbard model only under the condition that U/t 1. To account for both of these factors, in Fig. 8b we replace J in the t-J model bound with 4t 2 /U and fix U/t = 100. In this case, we see that the bound for the t-J model is over an order of magnitude lower than for the Hubbard model. The reason for this behavior can be intuitively seen when considering Fig. 5. For the Hubbard model, U sets the energy scale, with the gaps between relevant energy levels being on this or-der. For the t-J model, J sets the energy scale. Thus, a Hubbard simulation that wants to capture the behavior of the model at the energy scale of t-J must be able to resolve energy gaps on the order of J, which for the t and U parameters used in Figs. 5 and 8, are two orders of magnitude smaller than U . For the same accuracy, the Hubbard simulation must have much greater fidelity, and therefore requires a larger Trotter depth 1 .
It is important to note that in Fig. 8, as in our other comparison plots, we have plotted the bound in terms of the dimensionless quantity r /τ 2 , which is the maximum number of Trotter steps multiplied by the allowed error divided by the square of the dimensionless time parameter. When comparing the Hubbard and t-J model, the implicit assumption of this approach is that the allowed error, , is the same for both models. If we wish to treat U and J as independent parameters, as in Fig. 8a, an alternative approach that maintains the fairness of the comparison is to realize that must be different for each model. For the Hubbard model it is logical to set as a fraction of U while for the t-J model it is logical to set it as a fraction of J. Since the Hubbard model must resolve energy spacing on the order of J to have accuracy comparable to the t-J model, and since U J, this means H tJ . As an illustration of this behavior, let us consider a specific example for parameters of t = 0.1 and U = 10. In order to resolve the energy spectra of both models with equal accuracy, we need to pick an significantly lower than the smallest characteristic energy scale, which in this case is given by J = 4t 2 /U . Let us choose = 0.1J. Thus we see that for the Hubbard model we have an that is five orders of magnitude smaller than the characteristic model energy scale, while for the t-J model we have an that is only one order of magnitude smaller than the characteristic energy scale. Assuming an evolution time of τ = 10t and plugging all our parameters into Eqs. (35) and (41) we arrive at r H com ≈ 1.34 × 10 6 and r tJ com ≈ 4.15 × 10 4 , indicating that the Hubbard model is over 30 times more costly to simulate than the t-J model at these parameters.
From these results, it is clear that the best choice of model to simulate depends on the desired parameter range, namely the ratio of U/t. With this in mind, in Fig. 8c we plot the bound on the Trotter depth as a function of this ratio. We see that as U/t → 0 the bound for the Hubbard model vanishes. This is expected, as in the limit of vanishing U with finite t all terms in the Hubbard Hamiltonian commute. Conversely, as U/t → 0 the bound for the t-J model blows up. This occurs since the t-J bound grows with J and J is inversely related to U/t. Physically, this can be seen as corollary to the behavior discussed in the previous two paragraphs. As we move beyond the parameter range satisfying U t, the energy scale grows much larger than the typical spacing of the t-J model, and eventually the assumptions going into the derivation of the model break down.
As the number of Trotter steps is directly proportional to the circuit depth, this advantage translates to a significant reduction in overall gate count for simulating the t-J model in comparison to Hubbard. These results indicate that, on NISQ devices where gate errors and decoherence are limiting factors on circuit design, the t-J model is a more amenable target for near-term simulation.

VI. CONCLUDING REMARKS
In this work, we constructed a lattice-fermion formalism for studing the Trotter depth of important condensed models. We applied our approach to the Hubbard and t-J models, motivated by the supposition that the reduced Hilbert space of the t-J model due to the elimination of any doubly-occupied states would naturally lead to lower computational costs. Our results show that the situation is significantly more nuanced than this straightforward hypothesis.
In terms of the Pauli depth, we find that we gain no advantage from the reduced Hilbert space. While it is true that the t-J model has only three basis states for each orbital, |0 , |↑ , and |↓ , in comparison to the four basis states of the Hubbard model, which also allows |↑↓ , both cases still require two qubits to represent their respective bases. In this case, the t-J model is actually disadvantaged, as extra qubit overhead is required to implement the projection operators that ensure the system never strays out of the reduced Hilbert space. We note, however, that this is naturally a consequence of using a fermion-qubit mapping. If we instead applied a qutrit based mapping for the t-J model, the additional overhead would be mitigated (an example of a fermion-qutrit mapping for the 2D t-J model can be found in Ref. [49]). While the difference in Pauli depth is not significant in the 2D case, as the Jordan-Wigner transformed Hamiltonians for both the Hubbard and t-J models contain terms that span the full width of the lattice, it plays a larger role in the 1D case where these non-local Jordan-Wigner strings cancel out (see Appendix F).
In terms of the bound on Trotter depth, we find the situation to be more in-line with our initial hypothesis. If we naively treat U and J on the same footing, the bound first appears to favor the Hubbard model by around two orders of magnitude. However, if we properly account for the fact that the t-J model is a valid approximation to the Hubbard model only under the condition U >> t, and that J is in reality a function of U and t, then we see that the bound significantly favors the t-J model. Thus, the most efficient choice of model depends on the energy scale that the simulation needs to resolve.
We provide the full form of the Jordan-Wigner transformed Hamiltonians for the Hubbard and t-J models. In particular, we include an explicit qubit representation of the 2D t-J model without relying on alternative approaches, such as the auxiliary particle scheme [50]. While this work provides a first benchmark of comparison between the Hubbard and t-J models, there remain numerous avenues for future exploration. Chief among these is the extension to higher order product formulas, using the generalized bounds provided in Ref. [43]. More compact fermion-qubit mappings that reduce the Pauli depth through the use of ancilla qubits [51][52][53][54][55][56][57][58] or via circuit design that leads to the cancellation of Jordan-Wigner strings [11] could also be considered. However, we note that, as the Pauli depth for both the Hubbard and t-J models is equivalent outside of 1D, any compact mapping based on reducing the Jordan-Wigner string will benefit both models equally and will not impact that advantage in Trotter depth that the t-J model possesses.
A more detailed accounting of the resource overhead could also be carried out in order to provide specific gate count estimates for each model, which will serve as a necessary prerequisite to implementing simulations on realworld devices. Furthermore, it would be of interest to determine if the t-J model shares a similar advantage under other simulation schemes, such as hybrid quantumclassical variational algorithms [59,60]. Finally, we note that the general nature of the formalism developed in Section II makes it straightforward to extend this approach to other lattice models with different geometries or to orbital models for application to quantum chemistry simulations. In this appendix we prove the spectral norm of a Kronecker product of any number of Pauli matrices is always unity. Let us consider the Kronecker product of N Pauli matrices, where α = 0, x, y, z with σ 0 being the identity matrix. As transposition is distributive over the Kronecker product, and as the Pauli Matrices are Hermitian, the matrix in Eq. (A1) must be Hermitian. Thus, the spectral norm will be given by the absolute value of its largest magnitude eigenvalue.
The eigenvalues of σ x , σ y , and σ z are each {−1, 1} while the eigenvalues of σ 0 are {1, 1}. As the eigenvalues of the Kronecker product of two square matrices are given by the products of the eigenvalues of the individual matrices, we can see that the eigenvalues of Eq. (A1) will be −1 and 1, each with multiplicity 2 N/2 . Thus, we note, and Thus we have, Now using [41], we can express Eq. (B4) as, Thus, The projection operators can be explicitly accounted for by replacing the fermionic operators in Eq. (B1) with the projected operators [38], Note that the projected operators in Eq. (B8) provide a local implementation of the the Gutzwiller projection through the elimination of any terms that would ever produce a doubly-occupied site. However, the projectors in Eq. (B8) will not protect against an initial state of the system that contains doubly-occupied sites. In other words, the implicit identity operators that act on every site pq = ij, kl contained in each term in Eq. (B9) are not projected.
Appendix C: 1-norm bound on Trotter depth In this appendix we determine a looser bound on the Trotter depth for the Hubbard and t-J models by applying the operator norm scaling derived in Ref. [43]. For the case of the 1-norm, the bound for a first order product formula is given by, As in the calculation of the commutator bound, we begin by decomposing the full Hamiltonian into terms that consist only of products of Pauli matrices, For the purposes of Eq. (C1), H γ = H δ ij . Thus we have, where N ≡ N x N y is the total number of sites in the lattice. Since each H δ ij consists of a product of Pauli matrices, along with some prefactor, its norm will simply be the absolute value of the prefactor.
For the 2D Hubbard model there are 12 distinct H δ ij terms, listed in Appendix E. Eight of these terms have a prefactor of t/2 and four have a prefactor of U/4. However, one of the U/4 terms is simply an identity term and can be disregarded. Thus for the Hubbard model we have, For the 2D t-J model there are 64 distinct H δ ij terms, also tabulated in Appendix E. Half of these terms have a prefactor of t/16 and the other half have a prefactor of J/16. However, two of the J/16 terms are identity terms and can be disregarded. Thus for the t-J model we have, In order to compare the bounds given by the 1-norm scaling to those given by the commutator scaling, we consider the ratio, In Fig. 9 and 10 we plot this ratio for the 2D Hubbard and t-J models, respectively. We see that, at fixed N , the commutator bound is significantly tighter for both models. Notably, the commutator bound grows linearly with N while the 1-norm bound grows quadratically with N . Comparing Eqs. (35) and (C4) we see that for any fixed U and t there is no lattice size, N , where the 1-norm bound will be tighter. However, this is not true for the case of the t-J model. Comparing Eqs. (41) and (C5) we see that at small lattice sizes the 1-norm bound will be tighter due to the large polynomial of t and J in the commutator bound. This is illustrated in Fig. 10c where at small N with fixed t and J the ratio of the bounds is less than one, indicating that the commutator bound is larger. Furthermore, we also see that the ratio of the bounds is monotonic as a function of t and J for the t-J model, but that the Hubbard model displays a clear minimum as a function of both t and U . While a minimum in the ratio as a function of N also occurs between N = 1 and N = 2 for both models, this does have physical significance, as in practice N is a discrete quantity.

Appendix D: Periodic boundary conditions
In this appendix we determine the Trotter depth for the 2D Hubbard and t-J models for the case of periodic boundary conditions. We begin by simplifying Eq. (7) to include only nearest neighbor hopping, as is the case in both the Hubbard and t-J models, where b.t. denotes the contribution from the boundary terms. In Fig. 11 we illustrate the boundary terms that arise for a rectangular lattice of dimension N x by N y . Let us first consider the terms shown in Fig. 11a. We can see that we have N x terms that take the form A Similarly, the Hamiltonian for the 1D t-J model is, − (1 − n 2j−1 )n 2j n 2j+1 (1 − n 2j+2 ) . (F2) Note that we have expressed both these Hamiltonians in spinless form, which can be derived in much the same manner as in the 2D case, by mapping a single spinful chain to two spinless chains.
Applying the standard Jordan-Wigner transformation to the above Hamiltonians we arrive at, For the Hubbard model. Similarly, for the tJ model we find, Note that, unlike the 2D case, for the 1D models we find that the Jordan-Wigner strings cancel out, leading to significantly reduced Pauli depth. In this case, the presence of the projection operators in the t-J model plays a more significant role, leading to a maximum Pauli depth of four for the t-J model in comparison to three for Hubbard. Taking the case N y = 1 in Eq. (29) we find the general form, Noting that for both the 1D Hubbard and t-J models A j,j+p>1 = 0, this simplifies to, In the case of the 1D Hubbard model we have A 1,1 = 4|tU | and A 1,2 = 2t 2 + 2|tU |. Thus, r H com,1D = τ 2 4N |tU | + 2(N − 1)(2t 2 + 2|tU |) .