Feynman first proposed using one quantum system to simulate another [1]. A decade later, the first general quantum algorithms appeared [2][3][4][5][6], with applications to quantum chemistry [7], condensed matter [8], and high-energy physics [9]. Quantum simulation is now recognized as a significant future application of quantum computation [10][11][12], especially in the context of near-term devices. Quantum algorithms for quantum simulation with almost optimal scaling are now known [13][14][15][16][17]. Applications of these methods to condensed matter and quantum chemistry are well developed theoretically [18][19][20][21][22][23][24][25][26][27], and experiments have been performed on many different quantum architectures [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45].
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 quantum field theories (QFTs) in the nonperturbative regime.
Quantum simulation of high-energy physics can be approached via analog simulation of lattice gauge theories in cold atoms or ions [46][47][48][49][50][51][52], analog simulation using continuous variable quantum systems [53,54], or digital simulation of quantum field theories using conventional qubits and gates [9,[55][56][57]. Theoretical proposals for digital quantum simulation of quantum field theory [9,56,[58][59][60] were followed by experimental implementations of simple models such as the Schwinger model in 1 + 1 dimensions [(1 + 1)D] [37,61].
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 lightfront 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 this work, we focus on the former.
The main goal of this 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 space-time dimension. This reduces qubit requirements by several orders of magnitude: for example, the qubit numbers for the calculation in [62] are reduced from ∼400 000 to ∼1360 qubits. Second, the Hamiltonian in this encoding is sparse 1 [72,74], 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 this 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) [72,77,79], form factors [76], and decay constants [78]. PDFs constitute the dominant source of uncertainty on multiple cross-section predictions at the Large Hadron Collider (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 nonperturbative regimes. For example, they provide the most precise way of determining heavy-quark masses, improving the uncertainty on c-quark mass over nonlattice predictions by more than a factor of 2 [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 parametrized 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 power-divergent mixing between the operators occurs due to the reduced symmetries of the space-time lattice. Considerable progress has been made recently by applying large-momentum effective theory techniques [85]. Quasidistributions [85][86][87][88][89][90][91][92][93][94] 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 [95][96][97][98][99] and Compton amplitude [100][101][102][103][104], using pseudo-PDFs [105][106][107][108][109][110][111], and calculating good lattice cross sections [112][113][114]. 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 equal-time 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. II, we review the light-front treatment of a simple (1 + 1)D model containing coupled fermion and scalar fields [115][116][117]. 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 analog of the QCD parton distribution function, and propose an algorithm for its calculation. In Sec. III, we present the algorithm for calculating this quantity on a quantum computer, and provide detailed resource estimates for qubit and gate count. In Sec. IV, we discuss the generalization to (3 + 1)-dimensional [(3 + 1)D] QCD, estimate the required qubit resources, and describe how to calculate PDFs, decay constants, and form factors using our algorithm.
We now consider a simple (1 + 1)D model with Lagrangian [116,117]
(1) 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 analogs of the QCD parton distribution functions in this model.
In [9,56,62], authors studied algorithms based on equaltime 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].
Equal-time coordinates describe the Minkowskian spacetime as seen by a massive observer. The (1 + 1)D metric and gamma matrices can be chosen as g 00 = -g 11 = 1, g 01 = g 10 = 0, γ 0 = σ 3 , γ 1 = iσ 2 where σ 3 and σ 2 are the Pauli matrices in the standard basis. Light-front (LF) coordinates are obtained by performing the following coordinate transformation:
thus switching to the so-called light-front time (x + ) and distance (x -). Physically, we may think of this as describing the experience of a massless observer. The metric is g ++ = g --= 0, g +-= g -+ = 2, and the gamma matrices are defined as γ ± = γ 0 ± γ 1 . The only independent variables in the LF formulation of the theory (1) are the fields φ and ψ (+) , where
Here (±) act on the spinor field as projectors since ( (±) ) 2 = (±) and (+) + (-) = 1 [116].
For a system quantized in a box x -∈ (-L, L), the planewave expansions of the free fields are
where is the momentum cutoff and u is a momentumindependent spinor normalized to unity (unlike in equal-time quantization, where u n depends on the momentum quantum number n). Following [116,117], 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:
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 [115].
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 lightfront operators, the invariant mass operator of the theory can be expressed as
(
A study of bound state masses and their renormalization was performed in [117]. We defer this discussion until Sec. III F. 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 [116]. 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 Appendix 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 nonzero occupancies, and 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 multiparticle 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 [117] and the accompanying discussion). For a fixed K, the lowest eigenvalues of the mass matrix in the Q = 0 and 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. III F.
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 Appendix A). These play an important role in the mass renormalization [117], 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 analogs of QCD PDFs in this model.
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,118,119] (see also Sec. IV B). 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 cutoff 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 FIG. 1. At fixed harmonic resolution K, one can calculate PDFs up to the energy scale Q 2 max (K ). Once calculated at some energy scale, the PDFs can be evolved according to the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations. 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 Sec. III. For a fixed K, there exists an upper bound on the free invariant mass squared [the left-hand 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 calculated at a particular value of Q 2 can be evolved according to the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations [75,79,[120][121][122][123], including beyond Q 2 max (K ) if needed: this is illustrated in Fig. 1. We show bosonic and fermionic parton distribution functions for the massive Yukawa model defined in Sec. II evaluated at harmonic resolution K = 14 in Fig. 2.
In this section we present an algorithm for simulation of QFT in the front form on a digital quantum computer. In Sec. III A we describe the scaling with harmonic resolution FIG. 2. Bosonic and fermionic parton distribution functions, as defined in Eq. (10), for the massive Yukawa model defined by (1) evaluated for harmonic resolution K = 14. The values of parameters are chosen as in [117]: m B = 6.7, m F = 1, λ = 1, = 2048. Shown for the M = 18.96 eigenstate with different values of momentum cutoff: 2 , where Q 2 max = 40.2 2 . The choice
max corresponds to taking all the Fock states from the K = 14 sector into account.
of the dimension of the Hilbert space. In Sec. III B we present three encodings, and explain the tradeoffs between efficiency of encoding and simplicity of encoded operations. In Sec. III C we discuss the cost of time evolution. In Sec. III D we discuss the preparation of bound states: the eigenstates of the interacting Hamiltonian. In Sec. III E, we discuss the measurement procedure used to obtain PDFS which, as was explained in Sec. II B, reduce 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 + 1)D, one example is given by the PDF, our focus in this paper; another example is the hadronic tensor [62]. A straightforward generalization to (3 + 1)D would allow one to calculate electromagnetic form factors and decay constants: this is discussed in Sec. IV. 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.
In (1 + 1)D, 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, Chap. 5 of [124]). 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 asymptotic behavior either: for nonzero Q we have
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 )].
TABLE I. Dependence on K of properties of the three encodings of the Fock space in (1 + 1)D. 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. III A)]. However, the Hamiltonian is no longer local in this encoding, so we must instead use sparsity-based techniques for simulation (discussed in Sec. III D and Appendix C).
Qubit number Q Hamiltonian locality Hamiltonian sparsity
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 analog of the direct mapping as in quantum chemistry [7,125], 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 [23,25,[126][127][128][129]. The mapping of bosonic degrees of freedom has been previously studied in [130][131][132][133]. 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 [130], 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 nth 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 [131][132][133] and is described in detail in [132]. A related mapping is described in [131].
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. III A) we can see that it is indeed optimal up to logarithmic factors. The compact encoding is discussed in detail in Appendix 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 [13,17]. The sparsity of the Hamiltonian of our model is shown in Fig. 3 and discussed in detail in Sec. III C. 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 I.
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 [13,17,134], 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 infor-FIG. 3. Hamiltonian sparsity vs K. The curves label the upper and lower bounds on the sparsity, while the data points mark the exact sparsities for K = 3, 4, . . . , 19. The upper and lower bounds are given by
mation for the Hamiltonian. In Appendix C we give details of implementing two oracles needed by the methods of [13,17]. The first is O F , an oracle that enumerates the positions of nonzero entries of the Hamiltonian in a given row. O F is defined in Appendix C 1 where we show that the cost of O F for the compact mapping is O( √ K log K ). The second is O H , an oracle that computes the value of a nonzero entry to p bits of precision given its indices. O H is defined in Appendix C 2 where we show that the cost of O H for the compact mapping is O(K log K + p 2 log p).
Using Theorem 1 from [13], simulation of time evolution for time t under a Hamiltonian on n qubits of sparsity d and maximum matrix element ||H|| max to precision is given in terms of the parameter τ = d||H|| max t. The number of calls to O H and O F is
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
The number of oracle calls required is then O(tK 3 ), and the number of gates required for this number of calls is O(tK 4 ) if p is polylogarithmic in K. The number of additional gates required is O(tK 7/2 ) and so the overall simulation cost up to logarithmic factors is O(tK 4 ).
State preparation by any of the standard schemes requires simulation of time evolution. These schemes include phase estimation, adiabatic state preparation [9,135,136], as well as variational approaches [40] and quantum imaginary-time evolution [137]. Adiabatic state preparation performs timedependent 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 [134] to bound the cost of adiabatic state preparation. Specifically, we use Theorem 10 of [134] 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) [137,138]. The BLFQ formulation [73] is particularly useful when one considers NISQ implementations [139].
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 singlemode 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 Appendix C 3, and the resulting number of controlled-NOT (CNOT) 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 cutoff as in Eq. ( 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 left-hand side 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.
Up to this point we have not discussed the issue of mass renormalization in our model, which was studied in detail in [117]. 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 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 twodimensional space (m B , m F ) [117]. 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 (A6) below. The starting value m (0) F is then found by substituting m (0) 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 [117]. Moreover, the fact that the Hamiltonian norm only depends on logarithmically (see Appendix 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 threedimensional space of bare couplings with an additional condition on a certain amplitude.foot_4
We now discuss how the machinery developed in Secs. II and III can be generalized to QCD in (3 + 1)D. 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]. 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[T r ,
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 + 1)D, 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 + 1)D, 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 cutoff ⊥ . 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 lightfront momentum [72], i.e., n > 0 in Eq. (28a). 7 Hence, although the qubit requirements in arbitrary dimension increase relative to the (1 + 1)D 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. 8In order to simulate the full QCD Lagrangian with harmonic resolution K and transverse momentum cutoff ⊥ , an upper bound on the number of required qubits to store the light-front wave function for QCD in 3
[see Appendix C for a more detailed analysis in the (1
The helicity is encoded by a single qubit because the LF Dirac spinor has two "good" (independent) components [147].
The number of flavors n f taken into consideration depends on the probing scale Q 2 . Evaluation of Eq. ( 30) for the computation of [62], which requires 400 000 qubits for an equal-time calculation on a 20 3 lattice, yields Q = 1360 qubits [after additionally including ancillas required for the computations (see Appendix C)]. This reduction in qubit numbers will become even more dramatic with increasing lattice size and cutoffs.
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 [123]
where |P is a hadronic state of four-momentum P (averaging over spins is implied unless |P is spinless), and
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 (kq) + 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 /(2Pq) [123]. According to the optical theorem, the cross section 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 ex-pressed 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. 9 For example, n f = 3 at Q 2 = (1 GeV) 2 and n f = 5 at Q 2 = (90 GeV) 2 .
10 Equations ( 35) and (36) are written in the light-cone gauge A + a = 0, and, therefore, do not contain the Wilson line operator [75].
Upon substituting the LF free-field expansion, (35) and ( 36) acquire a simple form in terms of the number operators [75,118,119]
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 + 1)D, 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]:
(40) For a Fock state |{ξ j , w j } whose total four-momentum is given by
a particular Lorentz-invariant cutoff 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, lightfront 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. II B and illustrated in Fig. 1, for fixed harmonic resolution [and transverse cutoff in (3 + 1)D] 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 [75,79,[120][121][122][123].
Expression (39) is appealing from the quantum computational perspective because the number operator can be measured efficiently (see Appendix C 3). This remains true if one wants to exclude certain Fock states from consideration according to (42). In Appendix C we illustrate this by providing an explicit realization of these measurements for the model (1).
As we mentioned above, the right-hand side 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 noninteracting 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.
In this section 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 [72,148]
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 μ (x) = f q f f (x)γ μ γ 5 f (x). 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:
l ⊥ j = k ⊥ jx j q ⊥ + q ⊥ for the struck quark (ξ j = ξ st ), k ⊥ jx j q ⊥ for all other partons (ξ j = ξ st ). ( 46) 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 Appendix 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 [72,78,149] 0|J μ (0
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.
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 [72,74], enabling us to make use of optimal simulation algorithms [13,134]. 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 analog of the parton distribution function. Quantum computation of these observables has been considered by other authors prior to this work in [62,150]. Some of the advantages of the light-front discussed in detail here were presented in [151]. 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, Eq. ( 30) gives 1360 qubits, much less than 4 × 10 5 qubits on the grid of the same size in equaltime 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 distribution functions, transverse-momentumdependent 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 this 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 [152].
An important issue arising within the DLCQ approach to QCD, which we have not addressed in this work, is the effect of gluon zero modes, whose absence in the free-field 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 nontriviality of the vacuum in the equal-time formulation [140][141][142][143][144][145]. In the context of DLCQ, taking zero modes into account may result in the appearance of new, noncanonical 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 [72,[153][154][155][156][157][158][159][160][161][162][163][164]. We leave these, and other further details of the application to (3 + 1)D, to future work.
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 [117]:
Because n K these bounds show that all self-induced inertias scale with K and as [ln(K/ )].
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 [Eq. (A2)] 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 [Eq. (A3)], H S [Eq. (A4)], and H F [Eq. (A5)], 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., 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 (B4), the constraints on the sizes of the sets of momenta are 1
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 (A13):
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 (A13): the maximizing choice is B = {K}, |F | = |F| = 0. Substituting these assignments into (A13) 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 (A10a)-(A10p), the largest is
We choose this set of terms rather than that in (A10k) or (A10l) 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 (A10m) 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 (A10m) when B = {K}, and forms our lower bound on the sparsity of the full Hamiltonian:
Figure 3 (in Sec. III B) plots the upper and lower bounds together with exact values for the sparsity, which we calculated directly from the Hamiltonians with Q = 0 and K = 3, 4, . . . , 19. The sparsity of the Hamiltonian is therefore tightly bounded by Eqs. (A15) and (A17), is 1 2 K 2 at leading order in K, and is (K 2 ).
We define H max as the largest matrix element of H in absolute value:
As follows from Eq. (A2), the Fock states having largest eigenvalues (considered as the eigenstates of the free Hamiltonian) are the ones with high bosonic occupancy, such as 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 (A3)-(A5), 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/ ) [Eq. (A9)]. Altogether, this results in the following bound for the Hamiltonian norm:
We may use the relations amongst the various Hamiltonian norms in [165] and the sparsity of the Hamiltonian from Appendix A 2 to bound the spectral norm of H as
The bound (B4) 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 (B5) 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 , (B3) 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 (B5). 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 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 (B5), 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 (B5) since in some states all of the momentum K is possessed by fermions or antifermions. Thus, our complete Fock states are stored as
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
In this Appendix we describe the details of the implementation of two oracles necessary for the sparse simulation algorithm in (1 + 1)D. These oracles are ubiquitous in methods for simulation of sparse Hamiltonians and were defined for the electronic structure problem in [21,24]. There are two differences in the definition of these oracles for the model defined in Sec. II. First, we do not rely on any analog 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. Second, 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 require the computation of the analog of the Slater rules and their implementation in the Hamiltonian oracle. This analysis would be 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 Appendix 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 Appendix 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. III E). This subroutine is described in Appendix C 3, and requires O( √ K + p) local gates.
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
we may generate the state |ψ 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 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 occupancies are 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 (A2)]. 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 increas-ing 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.
Let the t i encode fermions, antifermions, and bosons as 0, 1, 2, respectively. Allowed changes in occupancy are then (1) 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 (C3) 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 labeled 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 = kl. 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 labelings of each type of term by momentum tuples as discussed above. The mapping (C6) 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. 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
) 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 (C7), which in turn completes the full mapping (C6).
We take the first set of terms in H S [Eq. (A4)] as our example: we call this set of terms H S,1 . The matrix elements due to the remaining sets of terms in the Hamiltonian may be computed using similar methods. Substituting the explicit expressions for c n and {•, •} into the first line of (A4) gives
and( n j , w j ) = ( n i , w i -1), (C10) for some set of indices i, j, i , j ( denotes symmetric difference).
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( √ K log K ) gates. (3) For each i = 1, 2, . . . , I, perform the following mappings on the registers encoding
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. (6) For each pair (i, j) ∈ {1, 2, . . . , I} 2 , perform the following mapping on the registers encoding |( n i , w i ), ( n j , w j ), f i , g i [initially in the state |( n i , w i ), ( n j , w j ), 1, 1 ]:
The cost of this step is O(K log K ) gates. t (7) For each pair (i, j) ∈ {1, 2, . . . , I} 2 , perform the following mapping on the registers encoding |( n i , w i ), ( n j , w j ), f j , g j [initially in the state |( n i , w i ), ( n j , w j ), 1, 1 ]:
The cost of this step is O(K log K ) gates. (8) For each i = 1, 2, . . . , I, perform the following mappings on the registers encoding | f i , g i , s and | f i , g i , s :
The cost of this step is O( √ K log K ) gates. (9) 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( √ K log K ) gates. (10) For each i = 1, 2, . . . , I, perform the following mappings on the registers encoding | f i , n i , l and | f i , n i , n :
The cost of this step is O( √ K log K ) gates. (11) For each i = 1, 2, . . . , I, perform the following mappings on the registers encoding | f i , w i , w and | f i , w i , w :
The cost of this step is O( √ K log K ) gates.
When steps 2 and 3 have been completed, the fermionic part of the condition for the states to be connected is satisfied if and only if s = s = 1, i.e., if exactly one of the f i and one of the f j is 1. k and m store the fermionic momenta whose occupations are changed in the case when the states may be connected.
When steps 4 and 5 have been completed, the antifermionic part of the condition for the states to be connected is satisfied if and only if s = 0. When steps 6 and 7 have been completed, if the final case in (C18) or (C19) holds for any i (or j), then the two states are not connected. Step 8 therefore adds at least 1 to s if and only if the final case in (C18) or (C19) holds for some i j. Thus, after these operations are complete, the states can only be connected if s = 0.
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
(C24)
(3) The matrix element (C24) 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] (C25)
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
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 (C9) each require I log 2 K qubits for fermions, (where, again, p is the number of bits desired in the output matrix element).
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 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
) 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, 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 side 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, note the following:
(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 .
(C35) (2) For each j = 1, 2, . . . , I, perform the following mappings on the registers encoding |( n j , w j ), s : |( n j , w j ), s → ( n j , w j ), s + w j m 2 j n j . (C36)
(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): 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 noncutoff 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).
Here, sparsity refers to the maximum number of off-diagonal elements in any column (or row, since the Hamiltonian is Hermitian).
This corresponds to switching to the Drell-Yan-West frame[72,73,[76][77][78][79].
Note that in Eq.(13) for the (1 + 1)D theory, the only dimensionful quantities on the left-hand side are the masses (but not the box size L). In higher dimensions, the left-hand side of Eq. (42) also depends on ⊥ and L ⊥ .
As we mentioned above, in a more general case one needs to measure the sum of single-mode operators over the transverse directions and additional quantum numbers [Eqs.(38) and(39) below]. This only changes the complexity polynomially.
Strictly speaking, since the coupling constant also has to be determined from an experiment (similar to the masses), one needs to implement this procedure even to calculate static quantities. Calculating cross sections amounts to expanding the wave packets within Hamiltonian blocks of different K. However, due to the exponential growth of the Hilbert space size with √ K, the qubit asymptotics will remain unchanged.
Note that one only has to be careful when raising and lowering space-time indices because their metric is nontrivial.
The light-front zero mode a 0 requires special treatment; in particular, it carries the information about the equal-time vacuum of the theory[72,[140][141][142][143][144][145]. By imposing antiperiodic boundary conditions on the LF fields, one may by able to completely eliminate the effect of zero modes[146].
Technically, the scaling including dimension in the light-front formulation is O(dK ), whereas the scaling including dimension in equal time is O(K d-1 ⊥ ). The factor of d in the light-front scaling is due to the necessity of encoding the value of each component of momentum for each occupied mode. However, for a fixed theory, d is a constant, so we may ignore it in the light-front scaling.
We do not use |; ; 1 K for it is the so-called angel state which decouples from the rest of the spectrum[117].