Quantum Simulation of Quantum Field Theory in the Light-Front Formulation

Quantum chromodynamics (QCD) describes the structure of hadrons such as the proton at a fundamental level. The precision of calculations in QCD limits the precision of the values of many physical parameters extracted from collider data. For example, uncertainty in the parton distribution function (PDF) is the dominant source of error in the $W$ mass measurement at the LHC. Improving the precision of such measurements is essential in the search for new physics. Quantum simulation offers an efficient way of studying quantum field theories (QFTs) such as QCD non-perturbatively. Previous quantum algorithms for simulating QFTs have qubit requirements that are well beyond the most ambitious experimental proposals for large-scale quantum computers. Can the qubit requirements for such algorithms be brought into range of quantum computation with several thousand logical qubits? We show how this can be achieved by using the light-front formulation of quantum field theory. This work was inspired by the similarity of the light-front formulation to quantum chemistry, first noted by Kenneth Wilson.

Quantum simulation of relativistic quantum field theory poses new challenges. Among these challenges are the absence of any fixed particle-number formulation of relativistic quantum theory, multiple particle types with varying statistics, complicated interactions and observables, nontrivial vacuum structure and gauge invariance. Nevertheless, quantum simulation is the only efficient approach to the study of general QFTs in the non-perturbative regime.
However, ab initio digital quantum simulation of general QFTs using existing techniques requires hundreds of thousands of logical qubits [62,63]. This is far beyond what is required for Shor's algorithm, which has been the subject of serious architectural studies estimating requirements of several thousand logical qubits [64,65].
Can the ideas explored for quantum simulation of quantum chemistry be used to enable simulation of QFT on quantum computers with several thousand logical qubits? Fortunately we can be guided by Wilson, who suggested [66] that the light-front formulation of QFT [67,68,69,70] is amenable to the orbital representations used in chemistry. The light-front formulation is now well-developed [71]. Among its notable advantages are a trivial vacuum and the absence of ghost fields. The linearity of equations of motion further reduces the number of independent variables. While the discretized light-front quantization (DLCQ) [72] provides a natural framework for simulating fundamental interactions ab initio [73], the basis light-front quantization (BLFQ) [73] approach is well suited for constructing effective theories; in the present work, we focus on the former.
The main goal of the current paper is to demonstrate that the light-front formulation is advantageous for digital quantum simulation. First, the second-quantized form of the lightfront Hamiltonian permits a highly efficient encoding scheme, with qubit requirements scaling logarithmically with the spacetime dimension. This reduces qubit requirements by several orders of magnitude: for example, the qubit numbers for the calculation in [62] are reduced from ∼ 400000 to ∼ 1360 qubits. Second, the Hamiltonian in this encoding is sparse 1 [74,72], so one can employ sparsity-based simulation algorithms that are almost optimal in all parameters [13,16]. Third, the light-front approach is well-adapted to calculation of static quantities such as parton distribution functions, hadronic tensors, form factors, and decay constants [75]. All of these can be calculated directly from the light-front wave function, within the Fock space sector with some fixed light-front momentum 2 . This leads to a simple form for the corresponding qubit measurement operators.
These advantages apply to any light-cone formulation of a relativistic field theory. In the present work we focus on the DLCQ approach, which amounts to solving the fundamental theory in the plane-wave basis [72]. Within a complementary approach, BLFQ, one instead chooses the set of basis functions and effective interactions that are most suitable for a particular problem of interest. Quantum simulation algorithms based on this latter technique will be investigated in subsequent work.
We focus on the quantum computation of static quantities, which in our context refers to single-particle properties such as parton distribution functions (PDFs) [77,72,79], form factors [76], and decay constants [78]. PDFs constitute the dominant source of uncertainty on multiple cross section predictions at the LHC [80,81], substantially affecting the reach of searches for new physics at high final state masses. They affect experimental results as they limit the accuracy to which precision electroweak observables can be extracted from LHC data. To give an example, the difference in the W + and W − masses as measured at the LHC is 29.2 ± 28 MeV, with the PDF uncertainty accounting for 23.9 MeV, or 85%, of the uncertainty [82].
Lattice calculation techniques on a classical computer have been very successful at calculating static properties of quantum chromodynamics (QCD) in non-perturbative regimes. For example, they provide the most precise way of determining heavy quark masses, improving the uncertainty on c-quark mass over non-lattice predictions by more than a factor of two [83,84]. Similar improvement has been obtained on the estimate of the strong coupling constant α S [83]. This has a direct impact on the high-energy collider physics programs, as the parametric uncertainty on heavy-quark masses is the dominant uncertainty on the determination of the branching ratio of most of the Higgs boson decay channels. However, the static property that has the largest impact on physics at the current energy frontier is PDF.
Exploiting the factorization of short distance physics from universal large distance phenomena, the PDFs used at the LHC are obtained from a parameterized asymptotic form at low resolution (Q 2 ), perturbatively evolved to higher resolution at which cross sections in the partonic center-of-mass system are calculated. This low-Q 2 parametric form is largely responsible for the uncertainty in the knowledge of the PDF. A more precise prediction of the PDF in such a regime, and eventually at higher Q 2 , would therefore significantly improve the precision of theoretical predictions and of many experimental measurement results at LHC energies.
Currently, the dominant approach for performing ab initio QCD calculations in the strong coupling regime is lattice QCD (LQCD) [80,81]. Within the traditional approach to LQCD, one evaluates the PDFs indirectly, by calculating the matrix elements of local twist-two operators [80,81]. From a sufficient number of these operators the Mellin moments of PDFs can be reconstructed. In practice, one is limited to the first three moments because powerdivergent mixing between the operators occurs due to the reduced symmetries of the spacetime lattice. Considerable progress has been made recently by applying large momentum effective theory techniques [85]. Quasi-distributions [86,87,88,89,90,91,92,93,94,95] allow for the equivalent of higher moment calculations, by matching higher moment calculations to the effective field theory. More recent approaches include finding PDFs from the hadronic tensor [96,97,98,99,100] and Compton amplitude [101,102,103,104,105], using pseudo-PDFs [106,107,108,109,110,111,112], and calculating good lattice crosssections [113,114,115]. As noted above, these calculations are at present not sufficient to reduce the theoretical uncertainty due to the PDF in many high energy physics measurements such as [82].
Lattice QCD is based on path integral quantization, and thus requires Wilson gluon lines and loops to maintain the color gauge invariance. Finite size lattices with periodic boundary conditions constrain the simulations (as do the prescriptions for fermion sources), causing the fermion doubling problem. The numerical sign problem severely complicates Monte-Carlo sampling in strongly interacting fermionic systems. On the other hand, the second-quantized approach in light-front formulation avoids Wilson loops and gauge group discretization, but eventually must treat the gluon fields on an equal footing with the quark fields and their interactions.
Previous work on digital simulation of QFT has mainly focused on dynamic quantities like scattering cross sections [9,56,58,59]. The possibility of studying parton physics on quantum computers was first explored in [62]. These authors proposed an algorithm for calculating PDFs and the hadronic tensor of the massive Thirring model, based on equaltime quantization of the lattice Hamiltonian. However, because these approaches are based on an equal time lattice formulation of QFT they lead to daunting qubit requirements.
We study the computation of static quantities by digital quantum simulation in the light front formulation. In Sec. 1, we review the light-front treatment of a simple 1 + 1D model containing coupled fermion and scalar fields [116,117,118]. In the front form, the Hamiltonian matrix of the field theory quantized in a box is block-diagonal. Each finite size block approximates the Hilbert space of the theory with a certain precision, the so-called harmonic resolution. The harmonic resolution therefore plays the same role as the number operator in for example, simulations of quantum chemistry. For this model, we introduce an analogue of the QCD parton distribution function, and propose an algorithm for its calculation. In Sec. 2, we present the algorithm for calculating this quantity on a quantum computer, and provide detailed resource estimates for qubit and gate count. In Sec. 3, we discuss the generalization to 3 + 1D QCD, estimate the required qubit resources, and describe how to calculate PDFs, decay constants, and form factors using our algorithm.

Light-front Quantization of the + 1D Yukawa Model
We now consider a simple 1 + 1D model with Lagrangian [117,118] Here φ and ψ are mutually interacting real boson and fermion fields. As in QCD, due to confinement emerging at low energies, the eigenstates of the interacting theory can be thought of as composite particles -bound states which are made of quanta of fields φ and ψ, the partons. We will introduce analogues of the QCD parton distribution functions in this model. In [9,56,62], authors studied algorithms based on equal-time quantization and spatial discretization of the wave function. We instead use light-cone coordinates x ± and work with the second-quantized formulation of the theory known as the discretized light-cone quantization (DLCQ) [72].
For a system quantized in a box x − ∈ (−L, L), the plane wave expansions of the free fields are where Λ is the momentum cutoff and u is a momentum-independent spinor normalized to unity (unlike in equal-time quantization, where u n depends on the momentum quantum number n). Following [117,118], in (3) we impose periodic boundary conditions. The discretized momenta and energies of the free particles are where m is either the boson or fermion bare mass. The creation and annihilation operators obey canonical commutation relations: [a m , a † n ] = δ mn , {b m , b † n } = δ mn , {d m , d † n } = δ mn . When quantizing in equal-time (x 0 , x 1 ) coordinates, a complete set of commuting observables (CSCO) for the theory is given by the charge Q, momentum P , and energy E. Under the transformation (2) to LF coordinates, these become P ± = E ± P . The charge Q, P + , and P − form a CSCO in the light-cone coordinates [119].
The dimensionless operators K (the so-called harmonic resolution) and H (which we shall call the Hamiltonian) are defined by P + = 2π L K, P − = L 2π H. In terms of these light-front operators, the invariant mass operator of the theory can be expressed as A study of bound state masses and their renormalization was performed in [118]. We defer this discussion until Sec. 2.6. Note that as one switches from P + and P − to the dimensionless operators K and H, the particular value of L may only become important at the stage of converting from light-cone coordinates to equal-time quantities. As it follows from eq. (5), the value of L is irrelevant for calculating the mass spectrum. As we shall see later, neither will it enter the expression for parton distribution functions (which is to be expected, since the latter describe the relative parton momentum distributions within the bound state).
The second-quantized expressions for H, K, and Q in terms of the creation and annihilation operators are obtained from Lagrangian (1) by means of the Noether procedure [117]. The charge and harmonic resolution are: The Hamiltonian H is a sum of four types of terms: H M is a (diagonal) mass term, while H V , H S and H F contain a number of interaction terms qubic and quartic in creation and annihilation operators (see App. A). The elements of the Fock space are labeled by orbital occupancies for the fermionic, antifermionic, and bosonic degrees of freedom: In eq. (8) we only list modes with non-zero occupancies, the hat is used to collectively denote all the particle species.
The crucial fact is that the spectrum of the operator P + is bounded from below, unlike that of the equal-time momentum P . In the equal time formulation, in an inertial reference frame, the Fock space sector of any fixed total momentum contains an infinite number of multi-particle states with that momentum, since those states can contain arbitrary numbers of left-and right-moving particles whose momenta cancel each other. Therefore, in order to obtain a Hilbert space of a finite dimension, one has not only to introduce a momentum cutoff, but also to limit the number of bosonic quanta in a Fock state.
To see how this changes in the light-front, consider an observer moving at the speed of light to the left in equal-time coordinates. In the light-front formulation, this observer has constant light-front coordinate (i.e., is stationary), so to the observer all massive particles appear to be moving to the right, and have positive light-front momentum. Therefore, there can be no cancellation of momenta due to left-and right-moving particles. This implies that in a theory quantized in a box there exists a finite number of states with a given value of K. Thus, by restricting to a particular value of K one naturally obtains a finite-dimensional Hilbert space without the need to cut off the dimension of the Hilbert space by hand. For a fixed eigenvalue of Q, it turns out that the blocks of the Hamiltonian H corresponding to larger eigenvalues of K represent the Hilbert space of the system with a higher resolution.
Within a block of fixed harmonic resolution, the Hamiltonian is proportional to the mass matrix, eq. (5). Diagonalization of a fixed-K block of M 2 gives a set of bound states |Ψ K, s with masses M K, s -these are the physical states of the interacting theory: Increasing K results in considering more bound states, with higher resolution. Each state s = s * first appears at some K s * , and is also contained in all the blocks with K > K s * . By diagonalizing Hamiltonian blocks of relatively small K one can get a good idea of the general form of the spectrum (see Fig. 1 in [118] and the accompanying discussion). For a fixed K, the lowest eigenvalues of the mass matrix in the Q = 0 and Q = 1 sectors correspond to the physical (renormalized) boson and fermion masses. This gives a constraint (the so-called renormalization condition), obtained by insisting that these physical masses match their known empirical values. From this constraint we determine the bare masses appearing in the Hamiltonian matrix, which produces the rest of the physical spectrum upon diagonalization; see the discussion in Sec. 2.6.
Although the momentum cutoff Λ in (8) is not used to truncate the Hamiltonian matrix it corresponds to, it explicitly appears in the Hamiltonian due to the presence of the so-called self-induced inertias (see App. A). These play an important role in the mass renormalization [118], and are related to vacuum polarization and self-energy terms in the equal-time quantization [68]. In the next section we will show how the wave functions of mass eigenstates can be used to calculate the analogues of QCD PDFs in this model.

Parton distribution functions
The light-front approach to QCD is appealing because numerous quantities of practical interest, such as PDFs, elastic form factors, and decay constants can be calculated directly from the light-front wave function [72,73]. The PDF, f (x), represents the probability of finding a parton of type carrying a certain momentum fraction x = p + n /P + = n/K, where 0 < x ≤ 1, inside a bound state (hadron) with total light-front momentum P + = 2πK/L. Given a bound state of the interacting theory, the PDF can be calculated as an expectation value of the single mode number operator summed over all the quantum numbers other than the longitudinal momentum [75,120,121] (see also Sec. 3.2). However, since in our model the longitudinal momentum is the only quantum number, the PDFs of a particular bound state |Ψ K can be calculated simply as with the number operators of different parton species given by These define the number of partons carrying momentum fraction x = n/K inside a hadron studied at harmonic resolution K. Measuring the expectation values as in (10) results in evaluating PDFs at K points: x = 1/K, 2/K, . . . , 1. For a properly normalized state (with Ψ K |Ψ K = 1) of total charge Q, the normalization of PDFs is given by which reflects that fact that momenta and charges of partons should add up to those of the hadron. The PDF is also a function of the probing scale Q 2 , which is the magnitude of momentum exchanged in a scattering process. The probing scale Q 2 can be introduced by imposing a cut-off on bound states [72]. A particular way of doing this is achieved by only considering Fock states |{ p j , w j } of invariant momentum squared below Q 2 in the expansion of the bound state |Ψ (Q) K . In the absence of spin and transverse directions this constraint is: where the sums go over all the excited parton modes. 3 As follows from eq. (2), in the LF formalism the states of large equal-time momentum are those with either p + → ∞ or p + → 0. While the former option is automatically excluded in a block of fixed K, condition (13) ensures that the light-front momenta are not too small. In terms of truncated bound states, one calculates the PDFs at probing scale Q 2 as: This quantity is simply an expectation value of an unintegrated number operator, which may be calculated using a quantum computation that we discuss in the Sec. 2. For a fixed K, there exists an upper bound on the free invariant mass squared (the lefthand side of eq. (13)). This sets the maximum energy scale Q 2 max (K) up to which one can calculate PDFs at the given harmonic resolution. PDFs calculfvated at a particular value of Q 2 can be evolved according to the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations [122,123,124,79,75,125], including beyond Q 2 max (K) if needed: this is illustrated in

Quantum Simulation in the Light Front
In this Section we present an algorithm for simulation of QFT in the front form on a digital quantum computer. In Sec. 2.1 we describe the scaling with harmonic resolution of the dimension of the Hilbert space. In Sec. 2.2 we present three encodings, and explain the trade-offs between efficiency of encoding and simplicity of encoded operations. In Sec. 2.3 we discuss the cost of time evolution. In Sec. 2.4 we discuss the preparation of bound statesthe eigenstates of the interacting Hamiltonian. In Sec. 2.5, we discuss the measurement procedure used to obtain PDFS which, as was explained in Sec. 1.2, reduces in the DLCQ formalism to evaluating the expectation value of the number operator (14). The methods we describe in this Section are optimal in both qubits and gates up to logarithmic factors.
In what follows all the resource requirements will be given in terms of the harmonic resolution K (the dimensionless light-cone momentum), for two reasons. First, numerous quantities of physical interest can be calculated within a single Fock space sector with fixed K. In 1 + 1D, one example is given by the PDF, our focus in this paper; another example is the hadronic tensor [62]. A straightforward generalization to 3 + 1D would allow one to calculate electromagnetic form factors and decay constants: this is discussed in Sec. 3. In such calculations, K will define the number of points at which these quantities are evaluated. The calculation of dynamical quantities like cross-sections requires wave packets expanded in a basis of states having different total light-front momenta. In this case one would consider all the blocks of size up to some maximum K. In both cases K controls the resolving power of the theory in the light-front, and so is a natural quantity with which to express computational cost.

Hilbert space dimension
In 1 + 1D, for each harmonic resolution K we have a finite-dimensional Hilbert space D K , which can be further split into blocks D K,Q of fixed charge Q. A lower bound on the dimension of D K is given by considering bosonic configurations only, which belong to the D K,0 subblock. These are labeled by integer partitions of K, where the momenta n j are the parts of the partition, and the occupancies w j are the multiplicities: The number of partitions of K is denoted p(K): its asymptotic behavior is log 2 p(K) = Θ( √ K) (see for example Ch. 5 of [126]). Therefore dim D K ≥ dim D K,0 ≥ p(K).
Adding fermions and antifermions gives a subleading correction to the dimension of D K , since their occupancies are w j , w j ∈ {0, 1} due to the Pauli exclusion principle. Restricting to a particular value of Q does not change the aymptotic behavior either: since all purely bosonic configurations of the D K−Q(Q+1)/2,0 sector can be turned into those of D K,Q by adding fermions with momenta ranging from 0 to Q. Therefore, independent of whether or not we restrict to a fixed value of Q, the asymptotic behavior of the Hilbert space dimension is Θ(exp( √ K)).

State encoding
The light-front representation of the theory has given us a formulation in terms of orbitals representing fermions, antifermions or bosons with given momentum. We can use an analogue of the direct mapping as in quantum chemistry [7,127], which amounts to assigning a particular qubit register to each momentum mode. For fermions, this means having a single qubit for each fermionic degree of freedom. The anticommuting creation and annihilation operators can be defined with the aid of the Jordan-Wigner, Bravyi-Kitaev, or other related mappings [128,129,25,23,130,131]. The mapping of bosonic degrees of freedom has been previously studied in [132,133,134,135]. We consider two variations of the direct scheme which differ in how the bosons are encoded. The resulting encoding schemes will be referred to as the direct-direct or direct-compact mappings.
Within a block of harmonic resolution K, the occupancy w n of the bosonic mode of momentum n is bounded by the requirement that the maximum momentum carried by that mode is at most the total light-front momentum: w j ≤ r n , where r n = K/n . The direct-direct mapping, first introduced in [132], uses a unary encoding requiring r n qubits for storing r n + 1 levels of each bosonic mode. This results in a total of O(K log K) qubits. The bosonic creation and annihilation operators acting on the n-th mode are represented by a sum of r n 2-local terms. However, since the locality of the fermionic operators is at least logarithmic in K, one may naturally want to trade locality of bosonic operators for a reduced number of qubits.
We therefore describe the direct-compact mapping, which uses a binary encoding of the occupation number of the bosonic modes and requires log 2 r n qubits for encoding r n levels, giving a total of O(K) qubits. In this case, the creation and annihilation operators contain a sum of r n terms, each of which is log 2 r n -local. This encoding was recently used to describe molecular vibrations [134,133,135] and is described in detail in [134]. A related mapping is described in [133].
The optimal encoding in terms of qubit resources is the compact encoding scheme. This was first described for chemistry in [7] and efficient algorithms were given in [21,24]. For our model the compact mapping stores only the momentum modes with nonzero occupancies: For such an encoding, the number of qubits scales as O( √ K log K); by comparing this to the Hilbert space dimension (see Sec. 2.1) we can see that it is indeed optimal up to logarithmic factors. The compact encoding is discussed in detail in App. B. Use of this fully compact scheme requires simulation algorithms which depend on the sparsity of the Hamiltonian in the chosen basis. Fortunately, methods based on sparsity scale optimally with almost all simulation parameters [17,13]. The sparsity of the Hamiltonian of our model is shown in Fig. 3 and discussed in detail in Sec. 2.3. We focus our efforts on quantifying the simulation complexity of the compact mapping because of its optimality. The properties of the three mappings are summarized in Table 1.

Mapping
Qubit number, Q Hamiltonian locality Hamiltonian sparsity Table 1: Dependence on K of properties of the three encodings of the Fock space in 1 + 1D. The direct mappings (which store the occupancies of all momentum modes) require O(K) qubits to encode K modes. In these schemes the Hamiltonian is a sum of Pauli terms of locality O(log K). The compact mapping stores the occupancies of only nonempty momentum modes, giving an asymptotic scaling of O( √ K) qubits (which is optimal up to logarithmic factors -see Sec. 2.1). However, the Hamiltonian is no longer local in this encoding, so we must instead use sparsity-based techniques for simulation (discussed in Sec. 2.4 and App. C).

Time evolution at constant harmonic resolution
The goal of our simulation algorithm is first to prepare the eigenstates of the interacting quantum field theory described by Lagrangian given in eq. (1). In each sector of fixed harmonic resolution K and charge Q, the lowest mass-energy particle is a physical particle of the theory. We then aim to perform measurements on the state to determine properties of these composite particles such as PDFs and form factors.
State preparation is a basic element of any quantum simulation algorithm. In this section we give bounds on the cost in terms of quantum gates required to evolve a state in a subspace of fixed harmonic resolution K for time t, to precision . We use the methods of [17,13,136], which are optimal in all relevant parameters.
Sparse Hamiltonians may be specified efficiently by two oracles: functions that can be called to give the defining information for the Hamiltonian. In App. C we give details of implementing two oracles needed by the methods of [17,13]. The first is O F -an oracle that enumerates the positions of non-zero entries of the Hamiltonian in a given row.
and an additional gates are required.
To simulate time evolution in a subspace of constant harmonic resolution K for time t in the compact mapping we have n = O(

State preparation
State preparation by any of the standard schemes requires simulation of time evolution. These schemes include phase estimation, adiabatic state preparation [137,138,9] as well as variational approaches [40] and quantum imaginary time evolution [139]. Adiabatic state preparation performs time dependent evolution under a Hamiltonian varying along a path connecting the target Hamiltonian to some initial, simple Hamiltonian. The minimum time this evolution can take while preserving the system in the ground state is determined by the minimum gap along the path. To determine the cost of adiabatic state preparation we must bound the spectral gap along a chosen adiabatic path, either rigorously or by invoking physical arguments. Assuming a gapped adiabatic path can be found, one must quantify the cost of simulation of evolution under a time-dependent Hamiltonian to perform adiabatic state preparation.
In our system, K controls the precision with which the theory describes the field theory in the front-form. We consider adiabatic paths such that the max norm of the Hamiltonian is everywhere upper bounded by the max norm of the Hamiltonian with the final K. We conjecture that amongst such paths a gapped adiabatic path exists that connects the theory at low K to the theory at high K. The property of the space of paths that we shall use in the analysis below is that the max norm of the Hamiltonian varies as O(tK/T ) for t ∈ [0, T ], where T is the length of the adiabatic evolution.
We will use the results of [136] to bound the cost of adiabatic state preparation. Specifically, we use Theorem 10 of [136] which, given a Hamiltonian on n qubits, with sparsity d and a bound on the integral of the maximum matrix element of ||H|| max along the path, ||H|| max,1 , gives the number of queries to H F and H O required as: and a number of additional gates scaling as O(d||H|| max,1 n). This method requires two additional oracles: one to compute a scaling factor, and one to compute the time-dependent max norm. Many paths obeying our max norm condition can by realized with O(log K) gates and so only change the scaling by logarithmic factors. The cost of adiabatic evolution for time T by this method is given by setting d = O(K 2 ), ||H|| max,1 = KT , and using the costs of the oracles O H and O F above to obtain a total cost of O(K 4 T ), which is the same as that for time-independent simulation. It remains a matter of future work to provide rigorous results on the efficiency of this, and other, adiabatic state preparation procedures.
Our Hamiltonian commutes with both K and Q, and we are interested in preparing states of specific charge in a sector of fixed K. Our compact encodings do not restrict to states of fixed charge. Exact evolution under the Hamiltonian will preserve the expectation values of these quantities given by the initial state. However, time dependent evolution or approximate evolution under the Hamiltonian may cause leakage to states of different K and Q for the direct mappings, and states of different Q in the case of compact mappings. We can therefore improve our state preparation by using phase estimation of K and Q after adiabatic evolution to project back to the desired sector. If the leakage to sectors of incorrect K or Q is small, then with high probability phase estimation of those operators will project us to the desired value. In the low-probability case that phase estimation results in the incorrect value of K or Q, we simply discard the result and start the state preparation procedure again.
The existence of small example problems for small K makes the implementation of such calculations on Noisy Intermediate Scale Quantum (NISQ) computers a possibility. Such calculations would attempt to variationally minimize the expectation of the invariant mass in a given sector of K and Q, and then perform measurements to estimate the PDFs in this variational ansatz [40]. An alternative to variational optimization of an ansatz is to use another heuristic such as quantum imaginary time evolution (QITE) [139,140]. The BLFQ formulation [73] is particularly useful when one considers NISQ implementations [141].

Measurement
One of the benefits of the LF formulation of QFT for quantum simulation is the simple form of measurement operators. This is due to the fact that one can calculate various observables directly from the LF wave function [72]. The determination of PDFs values at a fixed total light-front momentum K can be accomplished by estimation of single-mode operators,as can be seen from eq. (10). 4 This task can be performed efficiently on a quantum computer.
For the direct mappings, eq. (11) could be written simply as a sum over projectors on the qubit registers corresponding to the occupancies of modes of momentum n. For the compact mappings, we also wish to construct an operator whose eigenstates are the compact Fock states, and whose eigenvalue for a particular Fock state is the occupation of a particular mode n. The task is more complicated for the compact encoding, because a particular momentum mode is not always encoded in the same register of qubits, but it is still efficiently tractable. In brief, if we wish to extract the occupancy of a momentum mode n, then on each register that encodes a mode n j and its occupation w j we must perform a locally-controlled operation that adds w j to a fixed ancillary register if and only if n j = n. After performing this task on the fermions, antifermions, and bosons, the ancillary register will encode the total occupation of n, which we can then simply read out. This method is efficient. By employing a slightly more complicated scheme, we can further improve the efficiency of this operation; details are given in App. C.3, and the resulting number of CNOTs and single-qubit operations required is O( √ K + p) for p bits of precision. The dependence of the PDF on the probing scale Q 2 is introduced by imposing a momentum cut-off as in equation (13). This is also easily accomplished within the second-quantized formalism. Indeed, since the Fock states are the eigenstates of the free Hamiltonian, calculating the quantity on the LHS of (13) can be achieved by running the phase estimation algorithm for the free Hamiltonian. It only remains to introduce an ancillary register storing the information on whether the particular Fock state has to be kept in the expansion.

Mass renormalization
Up to this point we have not discussed the issue of mass renormalization in our model, which was studied in detail in [118]. Thus, we have implicitly assumed that the parameters in the Lagrangian are also those which would be measured in a thought experiment. This is only approximately correct for weak couplings and fails in the strong coupling regime. In order to correctly determine the eigenvalues and eigenstates of the M 2 operator for arbitrary values of coupling we must proceed as follows. Given as input the finite renormalized (i.e., physically observed) masses m B and m F and the coupling λ, we first determine the values of the bare constants m B and m F appearing in the Lagrangian. At a given harmonic resolution K, this is achieved by varying the bare masses m B , m F for fixed λ to satisfy the condition: which implies that the lowest eigenvalues in the Q = 0 and Q = 1 sectors of the mass matrix are associated with the physical boson and fermion masses. Having thus determined the bare couplings m B and m F , the M 2 operator now reproduces the spectrum of the theory at harmonic resolution K.
To solve (21) one performs a gradient search in the two-dimensional space (m B , m F ) [118]. The starting value m (0) B can be analytically calculated from the K = 2 sector of the model: where α 2 is a function of Λ defined as in (56)  B into the second line of (21) and performing gradient search in m F : Each iteration of the gradient search corresponds to a run of an algorithm described in the previous section. When conditions (18) are satisfied with the desired precision, the wave functions can be used for calculating the PDFs. The renormalizability of the theory guarantees the convergence of the method in the Λ → ∞ limit: the physical masses do not depend on the cutoff Λ [118]. Moreover, the fact that the Hamiltonian norm only depends on Λ logarithmically (see App. A.3) means that choosing sufficiently large Λ to obtain convergence does not require exorbitant resources.
If one moves further to calculating dynamic quantities, such as cross-sections, one would similarly have to perform the renormalization of the coupling constant λ, which would amount to performing a gradient search in the three-dimensional space of bare couplings with an additional condition on a certain amplitude. 5

Ab Initio Simulation of QCD in the Light Front
We now discuss how the machinery developed in Secs. 1 and 2 can be generalized to QCD in 3 + 1D. We briefly review the notations of QCD, discuss the qubit requirements, and present the expressions for observables in a form suitable for quantum simulation. We will see that our method gives asymptotic improvement in the scaling of qubit resources with cutoffs over previous simulation methods based on equal time formulation [63]. This results in several orders of magnitude fewer qubits for the smallest physically-meaningful cutoffs [62].

Light-front QCD in 3 + 1D
QCD is a field theory of Dirac fermions (quarks) interacting via an SU (3) 'color' gauge field. Due to the non-abelian nature of the gauge group, the mediators of the color interaction (gluons) carry the color charges themselves and can directly interact with each other. The fermionic field Ψ c,α transforms under the fundamental representation of the gauge group, c = 1, 2, 3 being the index in the color space, and α = 1, 2, 3, 4 being the Dirac spinor index. The spinor Ψ c,α of color c and its adjoint Ψ c,α = Ψ † c,β (γ 0 ) βα each have four complex components. The gauge field vector potentials A µ cc transform under the adjoint representation of the gauge group, and can be expanded as A µ cc = A µ a T a cc , where A µ a are the eight real vector fields, while T a cc are the generators of the gauge group obeying 6 f rsa being the structure constants of the su (3) algebra.
For N f flavors of quarks, the QCD Lagrangian acquires the following form: where is the color-electromagnetic field tensor, m f are the masses of different quark flavors, and is the covariant derivative.
In DLCQ, we shall use the collective label ξ containing the following degrees of freedom for gluons and quarks: where a is the color index in the adjoint representation, c is the color index in the fundamental representation, f is the flavor index, and λ is polarization or helicity. The discretized light-front momentum n = 2πk z /L is analogous to that in 1 + 1D, while n ⊥ = (n x , n y ) is the dimensionless internal momentum, defined by k ⊥ = (k x , k y ) = 2π n ⊥ /L ⊥ , which is introduced in order to separate the center-of-mass motion of the composite state. For a Fock state |{ξ j , w j } , where the sum goes over all the partons. In 3 + 1D, one immediately benefits from using the light-front formulation of non-Abelian gauge theories because of the vacuum triviality and the absence of ghost fields [72]. However, the presence of the transverse directions necessitates an additional momentum cut-off Λ ⊥ . The Hamiltonian matrix remains sparse [74], allowing one to use the algorithms discussed above.
Furthermore, in the DLCQ all the momentum modes, including those of massless bosons, necessarily carry a nonzero light-front momentum [72], i.e., n > 0 in eq. (28a). 7 Hence, although the qubit requirements in arbitrary dimension increase relative to the 1 + 1D case, their scaling with harmonic resolution K only increases to O(K), since in the worst case the state may be composed of K modes with light-front momentum one, all having distinct transverse momenta. Note that using the compact encoding (in the sense of only storing the occupied modes) is crucial: the number of unoccupied modes scales as the product of the momentum cutoffs over all dimensions. 8 In order to simulate the full QCD Lagrangian with harmonic resolution K and transverse momentum cut-off Λ ⊥ , an upper bound on the number of required qubits to store the lightfront wavefunction for QCD in 3 + 1 dimensions is: (see App. C for a more detailed analysis in the 1 + 1D case). The helicity is encoded by a single qubit because the LF Dirac spinor has two 'good' (independent) components [149]. The number of flavors n f taken into consideration depends on the probing scale Q 2 . Evaluation of eqn. 30 for the computation of [62], which requires 400000 qubits for an equal time calculation on a 20 3 lattice, yields Q = 1360 qubits (after additionally including ancillas required for the computations -see App. C). This reduction in qubit numbers will become even more dramatic with increasing lattice size and cutoffs.

Parton distribution functions
In QCD, all the information about the hadronic part of a scattering process is encoded within the so-called hadronic tensor W µν , also known as the forward Compton amplitude [125], where |P is a hadronic state of four-momentum P (averaging over spins is implied unless |P is spinless), and J µ ( is the quark current operator (q f being the quark charges; the sum is taken over all the quark flavors).
In the case of deep inelastic scattering l(k) + p(P ) → l(k − q) + X(P + q), where l is the lepton of momentum k, p is the proton of momentum P , and X is the final hadronic state of momentum P + q, one can write W µν in terms of two scalar structure functions W 1,2 as where q 2 = −Q 2 , and x = Q 2 /(2P · q) [125]. According to the optical theorem, the crosssection of such an inclusive process is given by the imaginary part of W µν , and hence of W 1 and W 2 . Within the parton model (i.e., to the zeroth order in the strong coupling constant α S ), Im W 1 and Im W 2 can be expressed as where f f (x, Q 2 ) are the parton distribution functions (PDFs), and the sum is taken over all the quark flavors contributing at the energy scale Q 2 . 9 In the light-front formalism, PDFs represent the probability of finding a parton of a given type carrying the longitudinal momentum fraction inside a hadron of a total longitudinal momentum P + . Formally, quark PDFs are defined as matrix elements of the quark field operators 10 [75]: where we suppress the Q 2 -dependence and assume that P |P = 1. The gluon PDF emerges as one further evaluates eq. (33) to the first order in α S . It is defined as The normalization of PDFs is given by which reflects the fact that that the individual momenta and charges of partons sum up to those of the hadron. Upon substituting the LF free field expansion, (35) and (36) acquire a simple form in terms of the number operators [120,121,75]: where In (39), the operator b † ξ creates a quark with quantum numbers ξ = {n, n ⊥ , λ, c, f}, while a † ξ creates a gluon with quantum numbers ξ = {n, n ⊥ , λ, a}.
As in 1 + 1D, within the discretized light-front formalism, introducing the momentum transfer Q 2 amounts to cutting off the total four-momentum of the Fock states in the expansion of the hadronic state [72]: For a Fock state |{ξ j , w j } whose total four-momentum is given by a particular Lorentz-invariant cut-off is provided by only considering Fock states of total invariant mass below Q 2 [72]: where m j , x j = n j /K and w j are the parton masses, light-front momentum fractions and occupancies, respectively; the intrinsic momenta k ⊥ are defined as in (29), and the sum goes over all the occupied parton modes. As we discussed in Sec. 1.2 and illustrated in Fig. 1, for fixed harmonic resolution (and transverse cutoff in 3 + 1D) the left-hand side of eq. (42) is bounded above by some energy scale Q 2 max (K, Λ ⊥ ). Once calculated at some scale Q 2 ≤ Q 2 max (K, Λ ⊥ ), the PDFs can be evolved according to the DGLAP equations [122,123,124,79,75,125].
Expression (39) is appealing from the quantum computational perspective because the number operator can be measured efficiently -see App. C.3. This remains true if one wants to exclude certain Fock states from consideration according to (42). In App. C we illustrate this by providing an explicit realization of these measurements for the model (1).
As we mentioned above, the RHS of eq. (33) is the zero-order term in the perturbative expansion of the hadronic tensor in the powers of the strong coupling constant. It is obtained from eq. (31) by replacing the full Heisenberg currents J µ (x) with the currents j µ (x) of the non-interacting theory. Within the traditional approach, one calculates the hadronic tensor by obtaining higher-order perturbative corrections to eq. (33). The paradigm of quantum simulation naturally suggests a different way to proceed: by switching to the interaction picture, we can move all of the complexity of the interacting theory into the state preparation stage. Doing so will have a minor effect on computational resources while allowing us to keep the measurement operators unchanged. Most importantly, such a calculation would be nonperturbative.

Form factors and decay constants
In this subsection we derive expressions for the electromagnetic form factor of a hadronic state (similar to (39) for the PDF) and for the decay constant. For a spinless state, such as a meson, the electromagnetic form factor F (Q 2 ) is defined as [150,72]: Switching to the Drell-Yan frame implies directing the incident hadron along the z-axis, and setting photon's momentum q µ transverse to this direction: where M is the hadron's mass.
In the LF the full Heisenberg current J µ (0) in eq. (43) can be set equal to the free quark current j µ (0) [72]. Similarly to (43), one can define the electroweak form factor by replacing J µ (x) with the chiral current J 5 In the LF, the expression for the form factor then takes the following form [76]: where the first sum goes over all the Fock states, while ξ st indicates the choice of the struck quark in |{ξ j , w j } and varies over all the quark modes (with charges q ξ st ). The state |{ξ j , w j } ξ st differs from |{ξ j , w j } only in its transverse momenta: Note that the final expression for the form factor in (45) involves state |P , but not |P [72,76]. To describe the corresponding measurement, we can formally rewrite (45) as Note that F(Q 2 ) is not a Hermitian operator, which should not surprise us since the form factor is generally allowed to take complex values. Nonetheless, the real and imaginary parts of F (Q 2 ) can be obtained by measuring the following Hermitian operators: Such measurements can be performed efficiently (see the discussion in App. C for the analogous case of PDFs). This would amount to constructing a circuit identifying the Fock state and transforming it according to (46). Importantly, the final state of each mode in (46) depends only on its initial state, and is not conditioned on the rest of the modes. Moreover, since each initial Fock state is mapped onto a linear combination of a small number (at most √ 2K ) of final states, the matrix F(Q 2 ) is sparse. The LF wave functions can also be used to calculate meson decay constants. For scalar (s) and pseudoscalar (p) mesons, those can be written in terms of the vector and axial quark current operators as [151,78,72]: Since decay constants are linear in the wave function, their measurement has to be designed somewhat differently from that of PDFs and form factors, which are bilinear in the wave function (see eqs. (38) and (43)). Up to a constant coefficient dependent on a particular state, the decay constant is the integral of the wave function over all the two-particle Fock states. For example, for π ± one has [72]: Since only the magnitude of the decay constant is physically significant, its calculation reduces to evaluating v|P for some fixed vector |v . This measurement can be performed efficiently as long as one can efficiently prepare the state |v from a computational basis state. Since, as we noted above, the decay constant only requires integrating over two-particle Fock states, the vector |v turns out to be itself a linear combination of two-particle Fock states. Thus, it is indeed efficiently preparable.

Discussion and Perspectives
We have demonstrated several advantages of the light-front formulation for quantum simulation of quantum field theory. The qubit requirements in the light-front approach are greatly reduced as compared with those in equal-time quantization. This is due to the smaller number of physical degrees of freedom, and the fact that the sum of occupancies in a Fock state is upper bounded by K for fixed harmonic resolution K.
The Hamiltonian matrix at fixed harmonic resolution in the LF formalism is sparse [74,72], enabling us to make use of optimal simulation algorithms [13,136]. These algorithms require O(tK 4 ) gates to simulate time evolution for time t, with logarithmic dependence on error. For state preparation by simulation of adiabatic evolution, in the case of adiabatic paths whose max norm is bounded by K, we require O(T K 4 ) gates. Proving that such paths in fact obey the adiabatic theorem is a topic for future work.
The LF formalism allows one to calculate various measurable quantities directly from the bound state wave functions. We demonstrated how such observables can be efficiently calculated on a quantum computer, using as our main example the analogue of the parton distribution function. Quantum computation of these observables has been considered by other authors prior to this work in [152,62]. Some of the advantages of the light-front discussed in detail here were presented in [153]. We hope future work will further develop all approaches to these problems.
In (1 + 1)D for harmonic resolution K the qubit requirements scale as O( √ K log K) in the compact encoding, which is optimal up to logarithmic factors. The compact encoding of light-front Fock states was shown to be extendable to higher-dimensional field theories. The qubit scaling increases to O(K), which is a significant improvement compared to equal-time quantization. For a 20 3 grid in momentum space with n f = 5 and n c = 3, equation (30) gives 1360 qubits -much less than 4 × 10 5 qubits on the grid of the same size in equal time quantization estimated in [62]. This is comparable to the number of logical qubits required to factor a 1024 bit RSA key using Shor's algorithm [64].
In higher dimensions, more observables can be calculated within a Hamiltonian block of a fixed longitudinal momentum. Those include decay constants, form-factors, generalized parton distributions functions, transverse-momentum-dependent distributions. As a possible direction of future work one could consider direct evaluation of the hadronic tensor. Instead of calculating it perturbatively (with the zeroth order being the parton-model approximation considered in the present paper), one could switch to the interaction picture, thus keeping the measurement operators unchanged while slightly complicating the state preparation.
Further development and optimization of the simulation techniques is warranted. Encoding schemes that restrict to a particular block of both K and Q would not change the asymptotic scaling of the qubit requirements but might be practically useful. Similarly, improvements to the implementation schemes given herein could yield significant reduction in gate numbers even if scaling improvements cannot be achieved. Such improvements likely require the development of software allowing the simulation of this algorithm, as has been developed for quantum algorithms for quantum chemistry [154].
An important issue arising within the DLCQ approach to QCD, which we have not addressed in the present work, is the effect of gluon zero-modes, whose absence in the freefield expansion is critical to our ability to use the compact encoding scheme. Zero-modes play an important role in the light-front formulation, incorporating all the complexity of the theory related to the non-triviality of the vacuum in the equal-time formulation [142,143,144,145,146,147]. In the context of DLCQ, taking zero-modes into account may result in the appearance of new, non-canonical interactions [72]. As noted in [72], while the longitudinal confinement is immanent in the light-cone quanitization of QCD due to the linear growth of effective potential in the x − direction, the interactions arising from zero-modes may be responsible for the transverse confinement.
As an alternative to the fundamental QCD Lagrangian, one can use effective low-energy theories. At the level of simulation, this amounts to changing the set of basis states to the one better resembling the bound state wave function. The latter approach seems to be particularly appealing due to the recent success of the so-called basis light-front quantization (BLFQ) technique [73]. Within this method, the effective Lagrangian, respecting all the symmetries of the full QCD, is solved in the basis provided by an exactly solvable model emerging from the AdS/QCD [155,156,157,72,158,159,160,161,162,163,164,165,166]. We leave these, and other further details of the application to 3 + 1D, to future work.
Following [117,118], we write the Hamiltonian as where where c n = a n / √ n. The expressions in (52)- (55) are called the mass, vertex, seagull and fork parts of the Hamiltonian, respectively.

A.1 Self-induced inertias
The mass term contains the so-called self-induced inertias α n , β n , γ n , the cutoff-dependent quantities whose appearance is a general phenomenon in the LF framework not specific to the particular theory under consideration. Those are defined as where We must upper bound these quantities as they contribute to the norm of the Hamiltonian, which in turn determines the simulation complexity. We can first evaluate the sums as follows: where H n is the nth harmonic number. Using the well-known bounds on the harmonic numbers we can upper and lower bound the self induced inertias as a function of the cutoff [118]: Because n ≤ K thes bounds show that all self-induced inertias scale with K and Λ as Θ(log(K/Λ)).

A.2 Hamiltonian sparsity
The sparsity of the Hamiltonian in the Fock basis is the maximum number of final states onto which a single initial Fock state is mapped under the action of the Hamiltonian. The mass terms in the Hamiltonian, H M (52), are proportional to number operators, hence are diagonal in the Fock basis, and so do not contribute to the sparsity. We analyze the sparsity of the terms in the Hamiltonian, H V (53), H S (54), and H F (55), by separately finding the numbers of nonzero images of each set of terms with a given form, for a generic initial state, then summing these and maximizing the resulting expression. We call the state whose number of nonzero images is maximal the sparsity-determining state (SDS).
Each Hamiltonian term contains annihilation operators that will map a state to zero unless they act on an occupied mode. Thus the number of nonzero images is largest for a state in which all bosons have distinct momenta, because this maximizes the number of Hamiltonian terms in which the annihilation operators do not map the state to zero. Therefore, all occupation numbers in the SDS are zero or one, so the SDS is determined by the sets F, F, and B of occupied fermionic, antifermionic, and bosonic momenta (respectively).
To obtain an upper bound on the sparsity, we assume that every term whose annihilation operators act on occupied modes maps the initial state to a distinct nonzero image. This is a relaxation of the actual condition in two respects: first, some of the nonzero images thus obtained may not be distinct, and second, fermionic and antifermionic creation operators acting on occupied modes will map the state to 0 rather than to a nonzero image. However, we will see that the upper bound we obtain by ignoring these reductions to the sparsity will nonetheless turn out to be asymptotically tight.
We consider sets of terms of a fixed form, e.g., {b † k b m c † l | k, l, m ∈ {1, 2, ..., K}, k +l = m}. These sets are represented in eqs. (60a)-(60p) by a characteristic element, e.g., b † k b m c † l . In each set, the modes for each ladder operator vary over {1, 2, ..., K}, under the constraint that total momentum is conserved, i.e., the sum of the momenta of the annihilation operators is equal to the sum of the momenta of the creation operators, which is equal to the transferred momentum. Each term is thus associated to a particular value of transferred momentum.
For each set of terms, the transferred momentum will be a sum over the possible sets of occupied modes corresponding to the annihilation operators in the set of terms. The summand will be the sum over the possible assignments of the transferred momentum to the creation operators in the set of terms. We can thus tabulate the numbers of nonzero images for each set of terms. We will use the following facts repeatedly: if some transferred momentum m is to be divided between two outgoing modes, this may be accomplished in m − 1 ways (since each mode must have nonzero momentum), obtaining m − 1 nonzero images. If there is only one outgoing mode, then clearly it must possess all of the transferred momentum, giving only one nonzero image. The numbers of nonzero images for each set of terms are given by: Here the term on the left is the representative element of an entire set of terms of that type. Our upper bound on the total number of nonzero images of the full Hamiltonian is the sum of these: To simplify the above expression, let denote the total momenta possessed by fermions, antifermions, and bosons (respectively) in the initial state. The sum of these must be the total momentum, i.e., K F + K A + K B = K. Furthermore, by (74), the constraints on the sizes of the sets of momenta are 1 ≤ |F| = I F ≤ √ 2K F , 1 ≤ |F| ≤ √ 2K A , and 1 ≤ | B| ≤ √ 2K B . Thus our upper bound on the number of nonzero images of the Hamiltonian becomes: Since |F|, |F|, | B| scale at most as the square root of K, only the first three terms in this expression grow as K 2 : all others grow at most as K 3/2 . Thus for large K, the sparsity is maximized by maximizing the first three terms in (63): But this expression is clearly maximized when all of the initial momentum is carried by a single particle, i.e., one of F, F, or B is {K}, and the other two are empty. Which we should choose is determined by the remaining terms in (63): the maximizing choice is B = {K}, |F| = |F| = 0. Substituting these assignments into (63) gives which is thus our upper bound for the sparsity. Direct evaluation of the sparsity of the Hamiltonian for small K shows that this bound holds for all K. The results are plotted in Fig. 3.
To obtain a lower bound, note that out of all contributions to the sparsity in (60a)-(60p), the largest is in (60m); the maximizing term type is b † k d † n c † l c m . We choose this set of terms rather than that in (60k) or (60l), because for the terms b † k d † n c † l c m all creation operators act on different kinds of particles, which removes the possibility of double-counting nonzero images. The SDS that maximizes (60m) is the same as the SDS for the full Hamiltonian: a single initial boson with momentum K. Therefore, the fermion and antifermion creation operators b † k , d † n can create a nonzero image for each allowed value of (k, n), since all initial fermion and antifermion occupation numbers are zero. Thus, the sparsity of the set of terms b † k d † n c † l c m is exactly the value of (60m) when B = {K}, and forms our lower bound on the sparsity of the full Hamiltonian: The sparsity of the Hamiltonian is therefore tightly bounded by eq. (65) and eq. (67), is 1 2 K 2 at leading order in K, and is Θ(K 2 ).

A.3 Hamiltonian norm
We define H max as the largest matrix element of H in absolute value: As follows from eq. (52), the Fock states having largest eigenvalues (considered as the eigenstates of the free Hamiltonian) are the ones with high bosonic occupancy, such as |; ; 1 K−2 , 2 or |1; 1; 1 K−2 11 . When acting on those, the number operator a † n a n produces a factor of order O(K).
The same logic is applicable to the interaction terms (53)- (55), none of which can result in more than a linear dependence on K. Lastly, the self-induced inertias scale asymptotically with K and Λ as O(log K/Λ) (59). Altogether, this results in the following bound for the Hamiltonian norm: We may use the relations amongst the various Hamiltonian norms in [167] and the sparsity of the Hamiltonian from A.2 to bound the spectral norm of H as:

Appendix B Compact Mapping
In this appendix, we describe how to efficiently encode an occupation number state (i.e., a partition of momentum) in 1 + 1D at fixed Harmonic resolution as a qubit state. This construction may be applied to bosons (with occupancies in [0, K]), or fermions, or antifermions (with occupancies in {0, 1}). We encode an occupation number state as a list of pairs: {(n 1 , w 1 ), (n 2 , w 2 ), ..., (n I , w I )}, where the n i are the distinct momenta (part sizes) that appear, and the w i = [1, K] are the corresponding occupancies numbers (we do not store modes with w i = 0). Each n i , w i is an integer in [1, K], so we encode each pair (n i , m i ) in a register of size 2 log 2 K qubits. The total number of qubits required is then since we use I registers.
Label the ith register X i : then the complete qubit states used to encode the momentum partitions may be written |X 1 , X 2 , ..., X I . I should be equal to the maximum number of distinct part lengths in any partition of K, but most partitions will not have I distinct part lengths, so in general the states will not use all of the registers X I . In order to uniquely associate an encoding state to each partition, we choose the following conventions: 1. Occupied momenta will be arranged in decreasing order of size.
2. If I ≤ I momenta are occupied, they are encoded in the first I of the X i .
How do we determine I, the maximum number of distinct part lengths in any partition of K? Simply adding up the momenta gives so we may obtain a tight bound on I by noting that the least value of K for a given I is obtained by setting n i = i and w i = 1 for i ∈ {1, 2, ..., I}. In this case, K is a triangular number, i.e., The bound (74) is satisfactory for analyzing the asymptotic qubit requirements, but to set the number of qubits we want to choose the minimum possible integer I such that a partition of K contains at most I distinct parts: this turns out to be To see that (75) gives the minimal I, let K be the largest triangular number less than or equal to K; then the minimal I exactly satisfies since for K , (73) applies exactly. From this we obtain i.e., 8K + 1 is an odd square. This implication reverses: if 8K + 1 = J 2 for odd J, then we can choose I = 1 2 (J − 1) and we obtain K = I(I+1)

2
, so K is triangular. The largest odd integer less than or equal to some arbitrary x is 2 x−1 2 + 1, so for arbitrary K the largest odd J whose square is less than or equal to 8K + 1 is: Thus 8K + 1 = J 2 for J determined in this way, so which is (75).
For example, suppose the momentum is K = 6 (chosen to be a triangular number, for convenience). Then I = 3, and the possible partitions are encoded as where each X i = (n i , w i ) is encoded in 2 log 2 (6) = 6 qubits, 3 to encode each of n i and w i . In the case of our algorithm the momentum is partitioned among fermions, antifermions, and bosons. Let K continue to denote the total momentum, summed over the fermions, antifermions, and bosons. For bosons, we use exactly the mapping of Fock states to qubit states described above; we still require I bosons = I as given in (75), since in some states all of the momentum K is possessed by bosons. For fermions and antifermions, we use the mapping described above, but with the occupation numbers restricted to be 0 or 1. Since only momenta that are present are represented in the state, this means that for all momenta that are present, w i restricted to be 1. Thus we may drop the occupation numbers w i entirely, and simply keep a list of the fermion and antifermion momenta that are present. We still require I fermions = I antifermions = I as given in (75), since in some states all of the momentum K is possessed by fermions or antifermions. Thus our complete Fock states are stored as {n 1 , n 2 , ..., n I ; n 1 , n 2 , ..., n I ; ( n 1 , w 1 ), ( n 2 , w 2 ), ..., ( n I , w I )} , where n i , n i , n i ∈ {1, ..., K} denote the fermion, antifermion, and boson momenta that are present in the state, and w i ∈ {1, ..., K} denote the occupation numbers of the occupied boson momenta. Thus the total number of qubits that this encoding requires is

Appendix C Implementation
In this section we describe the details of the implementation of two oracles necessary for the sparse simulation algorithm in 1 + 1D. These oracles are ubiquitous in methods for simulation of sparse Hamiltonians and were defined for the elctronic structure problem in [21,24]. There are two differences in the definition of these oracles for the model defined in Sec. 1. Firstly, we do not rely on any analogue of the Slater rules that define the nonzero matrix elements of the configuration-interaction (CI) matrix. Instead we use the second quantized representation of the Hamiltonian to enumerate nonzero elements of a row or column. Secondly, for the electronic structure Hamiltonian the matrix elements are defined in terms of integrals over basis functions whereas ours are simple functions of the momenta.
One could adapt the methods of [21,24] to the enumeration of nonzero matrix elements. This would required the computation of the analog of the Slater rules and their implementation in the Hamiltonian oracle. This analysis would more complex than that for electronic structure due to the presence of bosons and fermions and the more complex form of the second quantized Hamiltonian. The analysis would also need to be repeated for each model considered, whereas the results we give here can be generalized more directly to any Hamiltonian in second quantized form.
We describe three quantum subroutines: 1. A subroutine that enumerates the positions of the nonzero matrix elements of given row of the Hamiltonian in the Fock basis. This subroutine is described in App. C.1, and requires O( √ K) local gates.
2. A subroutine that, given a pair of Fock states each with total momentum K, computes the first p bits of the matrix element connecting the states in an ancilla register. This subroutine is described in App. C.2, and requires O(p + K) local gates.
3. A subroutine that permits evaluation of the number operator for a given momentum mode (see also Sec. 2.5). This subroutine is described in App. C.3, and requires O( √ K + p) local gates.

C.1 Matrix element enumeration
The Hamiltonian connects a pair of Fock states only if they are the same or the difference between them is exactly two fermions or antifermions (i.e., either two fermions, two antifermions, or a fermion and an antifermion) and either one or two bosons. In other words, given a Fock state |ψ = |n 1 , n 2 , ..., n I ⊗ |n 1 , n 2 , ..., n I ⊗ |( n 1 , w 1 ), ( n 2 , w 2 ), ..., ( n I , w I ) , we may generate the states |ψ is connected to by listing the possible changes the various terms in the Hamiltonian may make to |ψ . We represent these possible changes as lists, which we will denote ∆. The nonzero elements of a row or column in the Hamiltonian will be indexed by ordering the set of changes giving rise to the nonzero elements starting from the Fock state labeling the row. The ith nonzero element of the row will be labeled by ∆ i . Each ∆ has the form: where k + 1 is the momentum of the first momentum state whose occupancy will be increased by one (by a creation operator), and t 1 is indicates what type of particle it is (fermion, antifermion, or boson). Similarly, k ± 2 and k ± 3 are the momenta of the second and third momentum states whose occupancy is changed, which may be added or removed (since a term in the Hamiltonian possesses between one and three creation operators); t 2 and t 3 indicate their types. Finally, k − 4 is the momentum of the fourth momentum state whose occupancy changes, which if present, must be lowered, since no term in the Hamiltonian contains four creation operators; t 4 indicates its type. Note that the ordering of increases and decreases in occupancy here is not that given by order of the creation and annihilation operators in the second quantized representation of Hamiltonian terms.
If one or more of these is not needed (because the change being described involves fewer than four particles), then the corresponding k ± j is set to zero. Thus either all four k ± j are nonzero (describing changes due to terms containing four ladder operators), only first three are nonzero (describing changes due to terms containing three ladder operators), or all k ± j are zero (describing the connection of |ψ to itself via the number operators in (52)). Finally, in order to ensure that the ∆ associated to any particular change is unique, we require that particles added appear first, followed by particles removed, and subject to that rule, types are ordered as fermions in increasing order of momentum, then antifermions in increasing order of momentum, then bosons in increasing order of momentum. This induces an ordering on the ∆, which lets us enumerate the ∆ i , and hence the nonzero matrix elements in a row.
In addition to these rules, the change must conserve momentum, i.e., Since k ± i ∈ {1, 2, ..., K} (together with a bit encoding the sign) and t i ∈ {0, 1, 2} for each i, the number of distinct ∆ appears to scale like K 4 . However, the momentum conservation constraint (85) means that one of the k ± i is determined by the other three, so in fact there are fewer than 3 4 K 3 2 2 = 324K 3 distinct ∆s (3 4 possible valuations for the t i , K 3 possible valuations for the k i , and 2 2 possible valuations for the two ±s). This still overcounts the nonzero elements because the sparsity of the Hamiltonian in the Fock basis is Θ(K 2 ). This is because the full set of ∆s considered here can act on certain states to produce unphysical occupancies -either fermion occupancies greater than one, boson occupancies above cutoff, or negative occupancies. Hence not every ∆ returns a matrix element. We denote the number of ∆ by L = O(K 3 ), and index them as ∆ i for 0 ≤ i ≤ L − 1.

Operations:
We enumerate the ∆ i as follows. First we note that the Hamiltonian is a sum of O(1) types of term labelled by tuples of momentum orbitals subject to momentum conservation constraints. For example, consider vertex terms of type b † k b m c † l , where k = m + l and 3 ≤ k ≤ K.
Given the tuple k, m l we can construct the ∆ corresponding to this term with O(log K) operations. It only remains to show how to enumerate all tuples (k, l, m). In fact, we only need to enumerate (k, l) with k > l and compute m = k − l. The first few tuples (k, l) are (2, 1), (3,1), (3,2), (4,1). Let n(k, l) be the number of the tuple (k, l), with n(2, 1) = 1. Then n(k, l) = (l − 1) + n(k, 1) and: from which we obtain: Hence k, l, m can be computed from n by O(1) elementary arithmetic operations, the most costly of which is the square root. Therefore, the enumeration of the ∆ i corresponding to terms of type b † k b m c † l requires O((log K) 2 log log K) elementary gates. Similar arguments apply to seagull and fork terms.
We can now implement the oracle O F with action: Where i enumerates the fock states |φ i such that ψ| H |φ i = 0. The index i runs over all types of terms in the Hamiltonian and all labellings of each type of term by momentum tuples as discussed above. The mapping (88) may be implemented in three steps, using an additional ancillary register capable of encoding a ∆, also initialized to 0: The first step computes ∆ i from i and copies the Fock state |ψ to an ancilla register. Note that this Fock state ψ is a computational basis state of the qubits and so the no-cloning theorem does not forbid this operation. The Fock state |φ i is obtained by changing |ψ according to ∆ i if the resulting state |φ i is physical, or |φ i = |0 if changing |ψ according to ∆ i results in an unphysical state (for example, if ∆ i would remove a particle from a mode that is unoccupied in |ψ ). Finally, we invert the first mapping.
We now describe the second step in this mapping in detail. There are two substeps. First, we must check whether |ψ can be changed according to ∆ i . We check that for each fermion and antifermion added by ∆ i , the corresponding mode in |ψ is empty, and for each particle removed by ∆ i , the corresponding mode in |ψ is nonempty. To determine this by a reversible computation that can be made coherent, we append to ∆ i = (k + 1 , t 1 ; k ± 2 , t 2 ; k ± 3 , t 3 ; k − 4 , t 4 ) four ancillary qubits |c 1 , c 2 , c 3 , c 4 , initially all 0, and for each particle change (k ± i , t i ), flip the corresponding bit c i if the particle change cannot be performed on |ψ . Thus if c 1 = c 2 = c 3 = c 4 = 0, then |ψ F can be changed according to ∆ j . If any one of the c i is nonzero, then |ψ F cannot be changed according to ∆ i , so we set φ i = 0. We first consider adding particles to |ψ , then consider removing them.
For each (k + j , t j ) ∈ ∆ i , the mode is either present or absent in |ψ . This can be determined by O( √ K) gates. If the mode is present and t j is 2 (indicating that the added particle is a boson) then we simply increase its occupation by one. If the mode is not present, we must add it to |ψ F . Because the modes must appear in order of increasing momentum, adding a new mode requires shifting all modes with particle type t j and momentum greater than k j over by O(log K) qubits. We append the new mode to the register containing particles of type t j above the highest momentum mode present, k . We check if k j > k : if so, we have updated |ψ . If not, we exchange the new mode and mode k , requiring O( √ K) gates, and compare k j with the next smallest momentum mode. In the worst case the new mode has the smallest momentum and so this operation requires O( √ K log K) gates. For each (k − j , t j ) ∈ ∆ i , we must remove a particle of type t j and momentum k j . To do this, we reverse the method we used to add a mode, thus requiring O( √ K log K) gates. Beginning from the mode k of type t j with lowest momentum in |ψ F , check whether k j = k : if it is, then if its initial occupation is greater than one, decrease its occupation by one. If its initial occupation is one, remove this mode by setting the state of the corresponding mode register to 0, and then swap it to the end of the register. This completes the implementation of the second step in (89), which in turn completes the full mapping (88).

Operations:
Given the input states |ψ ⊗ |ψ , we wish to evaluate the corresponding matrix element. To do this, we attach ancillary registers whose state has the form where the f i , f j , f i , f i , f j are each qubits initially set to |1 . The s, s , s, s, s are each registers capable of encoding I as a binary number. The k, l, m, n, w, w are each registers capable of encoding K as a binary number, initially set to 0. The g i , g j are each registers capable of encoding K as a binary number, initially set to 1. We then perform the following operations: 1. For each pair (i, j) ∈ {1, 2, ..., I} 2 , perform the following mapping on the registers encoding |n i , n j , f i , f j (initially in the state |n i , n j , 1, 1 ): There are O(I 2 ) pairs and so the cost of this step is O(K log K) gates.
2. For each i = 1, 2, ..., I, perform the following mappings on the registers encoding |f i , s and |f i , s : The cost of this step is O( The cost of this step is O( √ K log K) gates.
4. For each pair (i, j) ∈ {1, 2, ..., I} 2 , perform the following mapping on the registers encoding |n i , n j , f i (initially in the state |n i , n j , 1 ): The cost of this step is O(K log K) gates.

5.
For each i = 1, 2, ..., I, perform the following mapping on the registers encoding |f i , s : The cost of this step is O( √ K log K) gates.
The cost of this step is O(K log K) gates.
The cost of this step is O(K log K) gates.
When step 9 has been completed, the states can be connected only when s = s = 1, i.e., if exactly one of the f i and one of the f j is 1. Thus when step 10 has been completed, l and n store the bosonic momenta whose occupations are changed in the case when the states may be connected; and when step 11 has been completed, w, w store the (larger) occupation numbers of the two bosonic momenta that change between the two states.
Having implemented all of the preceding operations, the matrix element ψ |H S,1 |ψ may be computed as follows: 1. ψ |H S,1 |ψ = 0 if and only if s = s = s = s = 1 and s = 0.
2. If the above condition holds, then ψ |H S,1 |ψ = w w ln 3. The matrix element (106) is a function of the numbers k, l, m, n, w, w , each of which has already been stored in its own register of log 2 K qubits. Thus we may compute the matrix element to any desired number of bits p and store it in a register of the same length, by an operation on the 6 log 2 K + p qubits involved. Computation of the square root requires two multiplications of two O(log K)-bit numbers, costing O(log K 2 ) gates. These two numbers are then divided, yielding a result with O(p) bits of precision, requiring O(p 2 ) operations. We then take the square root of this p-bit number, requiring O(p 2 log p) gates. To compute the second term in the case k = n we can either perform one addition and one subtraction of two O(log K)-bit numbers, followed by two divisions, or compute the common denominator and numerator, and perform one division. In either case the cost is O(log K + p 2 ). Thus calculating the matrix element requires O (log K) 2 + p 2 log p (107) gates.
The matrix elements due to other terms in the Hamiltonian may be evaluated using similar methods. Similar analyses will apply for each term, so the overall cost to evaluate a matrix element of the full Hamiltonian will be O(K log K + p 2 log p) = O(p 2 + K) (108) gates, where the dependence on p comes from the final calculation of the matrix element. We can also calculate the total number of qubits required for these operations. The input states in (91) each require I log 2 K qubits for fermions, I log 2 K qubits for antifermions, and 2I log 2 K qubits for bosons, for a total of 8I log 2 K ≤ 8 √ 2K log 2 K ≤ 12 qubits to encode the input states; as we would hope (since we have two input states) this is twice the number required to encode a single state in the compact mapping, as described in App. B, eq. (82). The ancillary registers in (94) require 5I qubits for the {|f i , |f j , |f i , | f i , | f j }, 5 log 2 I qubits for {|s , |s , |s , | s , | s }, 6 log 2 K qubits for {|k , |l , |m , |n , | w , | w }, and 2I log 2 K qubits for the { g i , g j }, for a total of 5I + 5 log 2 I + 6 log 2 K + 2I log 2 K ≤ 2 √ 2K log 2 K + 5 √ 2K + 11 log 2 K + p (110) qubits. Thus the total number of qubits required is upper bounded by 10 √ 2K log 2 K + 5 √ 2K + 11 log 2 K ≤ 15 √ K log 2 K + 8 √ K + 11 log 2 K + p (111) (where, again, p is the number of bits desired in the output matrix element).

C.3 Measurement of an occupation number
In order to evaluate PDFs, we need to be able to measure the number operator N , to estimate the expectation value (10). This means that for an encoded Fock state |n 1 , n 2 , ..., n I ⊗ |n 1 , n 2 , ..., n I ⊗ |( n 1 , w 1 ), ( n 2 , w 2 ), ..., ( n I , w I ) , we want to perform a measurement of the occupation number of some particular momentum mode n, summed over fermions, antifermions, and bosons.
To do this, we employ an ancillary register whose states have the form: |f 1 , f 2 , ..., f I ⊗ |f 1 , f 2 , ..., f I ⊗ | f 1 , f 2 , ..., f I ⊗ |s , where the f i , f i , f i are each qubits initially set to |0 , and s is a register of O(log K) qubits, initially set to 0. For the desired n, we iterate over the n i , n i , and n i , checking whether each is equal to n: if it is, then we set the corresponding f i , f i , or f i to 1. This requires 3I operations on O(log K) qubits, requiring O(I log K) gates. Now, we iterate over the f i and f i , adding their values to s. Each such operation is a binary addition on O(log K) qubits, and we implement 2I of them, so the total number of gates required is again O(I log K). After this step is complete, s will encode the total number of fermions and antifermions with momentum n (between 0 and 2).
Finally, we iterate over the pairs ( w i , f i ), adding the products of their values to s, i.e., Each such operation is a binary addition on O(log K) qubits, and we implement I of them, so the total number of gates required is again O(I log K). After this step is complete, s will encode the total number of fermions, antifermions, and bosons with momentum n. Once this routine is complete, we can sample the occupation number of mode n by measuring the qubits encoding s. The total cost is O(I log K) = O( √ K). The situation becomes only slightly more complicated when we impose the probing scale Q 2 as in (13). Now we wish to estimate the expectation value of the number operator N for the cutoff state |Ψ (Q) , as in (14). We use the following version of the cutoff condition (13): To calculate the left-hand of this expression, we employ an additional ancillary register whose states have the form: where |s is a register of p qubits (which we will use to store a floating point number, initially 0), and |t is a single qubit.
To evaluate the cutoff condition: 1. For each j = 1, 2, ..., I, perform the following mappings on the registers encoding |n j , s and |n j , s : |n j , s → |n j , s + m 2 j n j , |n j , s → |n j , s + m 2 j n j .
3. Perform the following mapping on the registers encoding |s , t (which will be in the state |s , 0 for some s set by the previous steps): When steps 1 and 2 have been completed, s will be an encoding (to the desired precision, set by its number of qubits p) of Thus the third step merely checks whether the value of s is bounded or not by the quantity Q 2 K , which is classically precomputed, and updates the qubit t accordingly. Then in order to get the expectation value of the number operator for the cutoff state |Ψ (Q) , as in (14), we compute the number operator as above for the non-cutoff state, but only for pairs (s, t) where t = 1; when t = 0 we throw away the sample.
This avoids sampling values corresponding to disallowed states. Note that if in some superposition of Fock states, those above the cutoff Q 2 K possess too much of the total probability, it may become inefficient to sample only from the allowed states. This situation can be avoided by keeping the imposed cutoff Q 2 not too far below the maximum energy scale Q 2 max (K); see Fig. 1 and the associated discussion. Each of the operations in step 1 above involves a division and an addition on p + log 2 K qubits: log 2 K for the n j or n j , and p for s . Each of the operations in step 2 above involves a division, a multiplication, and an addition on p + 2 log 2 K qubits: log 2 K for each of n j and w j , and p for s . We perform 2I of the first type, and I of the second; the final step is just a multiply-controlled NOT on p + 1 qubits, so the total number of CNOTs and single-qubit gates required is O(I + p) = O( √ K + p).