Parton Physics on a Quantum Computer

Parton distribution functions and hadronic tensors may be computed on a universal quantum computer without many of the complexities that apply to Euclidean lattice calculations. We detail algorithms for computing parton distribution functions and the hadronic tensor in the Thirring model. Their generalization to QCD is discussed, with the conclusion that the parton distribution function is best obtained by fitting the hadronic tensor, rather than direct calculation. As a side effect of this method, we find that lepton-hadron cross sections may be computed relatively cheaply. Finally, we estimate the computational cost of performing such a calculation on a digital quantum computer, including the cost of state preparation, for physically relevant parameters.


I. INTRODUCTION
Parton distribution functions (PDFs) and hadronic tensors have been widely studied as they compactly parameterize the hadronic structure as seen by high-energy probes. They provide the nonperturbative input to deep inelastic scattering cross sections, as well as the initial conditions for heavy ion experiments. As a result of this fundamental importance, significant effort has been expended on both the experimental and theoretical side [1].
The near-advent of small-scale universal quantum computers has stimulated research into their possible applications. Of particular interest here are the algorithms for achieving real-time evolution of a field theory, particularly gauge theories . Classically hard, real-time evolution is among the most natural operations available on a quantum computer. However, it is unclear exactly what quantities of physical interest may be obtained more easily by this method. Inclusive cross sections, for instance, are a natural target [45][46][47], but are already in practice obtainable classically [48]. In this paper we propose that PDFs and the hadronic tensor may be a physically relevant target on which "quantum supremacy" -that is, the ability of quantum computers to calculate quantities unobtainable classically -may be demonstrated. Additionally, this provides a cheaper way to compute lepton-hadron cross sections than that discussed in [45][46][47]. Finally, * hlamm@umd.edu † srl@umd.edu ‡ yyukari@umd.edu these methods could complement Euclidean lattice methods to obtain theoretical predictions in partonic physics with different systematics. The largest quantum computers currently available are limited to tens of qubits, and in the near future, computers with more than a few hundred, moderately noisy qubits cannot be reasonably expected. This situation has been termed the NISQ (Noisy, Intermediate-Scale Quantum) era. With an eye toward calculations achievable in the NISQ era, [49] proposed a method for computing the hadronic tensor, in the Regge limit, in the framework of color glass condensate EFT. In this work we discuss a more expensive prospect: first-principles calculations in the framework of Hamiltonian lattice field theory.
The method discussed here, as applied to 3 + 1 dimensional QCD, requires computers several orders of magnitude larger than anything expected in the NISQ era. In contrast, the 1 + 1 dimensional Thirring model we discuss requires substantially fewer resources and could plausibly be studied in the NISQ era as a toy model. Our purpose is not to present new calculations, but to describe how these observables can be computed and what sort of resources are required to obtain, on a quantum computer, these particular physical observables from first principles. To that end, we analyze the cost of adiabatic state preparation of the proton, as well as the practical cost of time-evolution on a sufficiently sized lattice.
For illustration purposes we work with the staggered Thirring model in one spatial dimension. The Hamiltonian of this model is where χ denotes a single-component fermion field, which becomes a two-component field in the continuum limit. The Hilbert space of this model is small enough to allow us to simulate the evolution of the system with classical resources on small lattices. The Thirring models differs importantly from QCD in not being a gauge theory; we will see that for the evaluation of the hadronic tensor this makes no significant difference to the method, but direct computations of PDFs are dramatically affected. The rest of this paper is structured as follows. First we describe the direct computation of the PDF in Sec. II. As we shall see, this is entirely impractical in a gauge theory, and so in Sec. III we show how the hadronic tensor is obtained instead. A theoretical prediction for the PDF may be extracted from the hadronic tensor via the same procedure used by experiments. In Sec. IV we discuss the preparation of the ground state of a proton -a prerequisite for any of these procedures. We conclude in Sec. V with a discussion of the future prospects for this and similar methods.

II. QUARK DISTRIBUTIONS
Central to hadronic physics are the quark and gluon distribution functions [50]. The distribution functions f (x) may be interpreted as giving the probability for a high-energy probe to see a parton with a given momentum xP within a hadron with momentum P . It is most natural to consider f (x) on the light cone; however, in this paper we will view it in equal-time quantization, as that is the framework in which we would simulated a field theory on a quantum computer.
We use the Thirring model for illustration, and by analogy with QCD we refer to the fundamental fermion as a 'quark' and the bosonic bound state as a 'meson'. Note that, because the Thirring model is not confining, the dressed fermion is itself an asymptotic state, of which a PDF may be computed. The quark distribution function in the Thirring model is simpler than that of QCD because no Wilson line is needed. It is given by where γ + ≡ 1 √ 2 (γ 0 + γ 1 ), and P + ≡ 1 √ 2 (P 0 + P 1 ) is the lightcone momentum of the incoming hadron. The expectation value is taken in the ground state of a hadron with momentum P . Here we see that the PDF is the Fourier transform of a time-separated correlator.
Our strategy will be to calculate the integrand on a quantum computer for many values of y, and then approximate the Fourier transform classically. In a form more readily obtained on a quantum computer, the quark correlator is The preparation of the hadron state |P is involved, and we postpone its discussion to Sec. IV. Assuming that |P is readily prepared, we must now translate this expression to the lattice, where ψ is staggered and y is discrete. It may be seen by taking the Fourier transforms of the fields ψ and χ that an appropriate "staggered PDF" is For brevity, we have used the notation |y| = 0 when y is even and 1 when y is odd. We must translate the operators χ † , χ, and e −iHt into bosonic qubit operators. The operators χ(y) and χ † (y) are anticommuting, and are constructed from bosonic qubit operators σ via the Jordan-Wigner transformation [51]. The construction of e −iHt follows the standard techniques of Trotterization. The Hamiltonian is split up into several mutually non-commuting terms H x , H y , and H z , each of which is easily diagonalized in isolation.
a yy (n)σ y (n)σ y (n + 1), and The operator that forms the matrix element in Eq. (4) is not Hermitian. In order to evaluate this non-Hermitian operator on a quantum computer, we decompose into a sum of unitary operators, each of which may be evaluatedà la [52] with the help of an ancillary qubit. The decomposition of the operator in Eq. (4) is: where χ x = χ+χ † and χ y = i(χ−χ † ). The coefficients C ij are determined from the Jordan-Wigner transformation: [52], each term in Eq. (6) is measured on a quantum computer by first preparing the state where a denotes the ancillary qubit, and then applying U ij controlled on the ancillary to |P . Measurements of σ x and σ y of the ancillary qubit on the resulting state give us the real and imaginary parts, respectively, of the terms in Eq. (6), P |U ij |P .
Two technical complications remain. First, because the correlators are evaluated along the light-cone the speed of light must be known precisely. Without the hypercubic symmetry of the Euclidean lattice preventing renormalization, the speed of light must be computed nonperturbatively. In principle, we could measure the speed of light on the quantum computer. This is a formidable task, entailing careful measurements of the dispersion relation near the continuum limit. But this is in fact unnecessary. The dispersion relation is reflected in the low energy portion of the spectrum of the Hamiltonian, which can be readily determined on an anisotropic Euclidean lattice [53,54]. Thus the speed of light on a quantum computer may be determined without any calculations being performed on a quantum computer, as long as the Hamiltonian limit is taken on both the classical and quantum machines. In the specific case at hand of the 1 + 1 dimensional Thirring model, the situation is even simpler: numerical experiments reveal that the speed of light in the continuum limit is 1 in lattice units.
The speed of light is only defined in the continuum limit. On the lattice, no exact light-cone exists, and 'space-like' separated fermionic operators need not exactly anticommute. As a result, the lattice PDF will not have the desired symmetry properties until the continuum limit is taken. Additionally, if periodic boundary conditions are used, care must be taken not to evaluate the quark correlator at separations larger than the spatial size of the lattice.
Finally we must take the Fourier transform. Only a finite number of values of the quark correlator may be computed, and naively taking the Fourier transform will show highly oscillatory artifacts (as in Euclidean lattice calculations [11]). In order to take the Fourier transform in a stable way, without these artifacts, we impose a Gaussian window, defining and first taking the limit L → ∞, and only then allowing → 0. These limits may be taken numerically. This completes the description of how to obtain the PDF of the Thirring model on a quantum computer, given an already-prepared hadronic state. In Fig. 1 is shown a calculation, by exact diagonalization of the Hamiltonian, of the fermion distribution function of the lowest-lying fermion state at vanishing and weak couplings.
The PDFs of asymptotic states in the Thirring model may be interesting in the context of quantum supremacy: the qubit cost of this calculation makes it potentially accessibly in the NISQ era, while classical algorithms struggle to obtain this observable. Ultimately, however, we would like to calculate the PDFs of mesons and baryons in QCD. Quark distribution functions in QCD are given by f i (x) = dx e ixP + y P |ψ i (y)γ + W (y; 0)ψ i (0) |P (10) where W (y; 0) is a light-like Wilson line connect y to the origin required to ensure gauge-invariance [50], and i enumerates the quark flavors. Several proposals have been advanced for how gauge field theories may be simulated on a quantum computer. Let us consider the scheme laid out in [41], which provides a procedure for computing Wilson loops. In this scheme, time-evolution is implicitly performed in A 0 = 0 gauge, so that a W (y; 0) is approximated by a sequence of spatial link operators applied at different points in time. (The time-like links are fixed to be the identity in this gauge.) W (y; 0) ≈ e iHy W (y; y − a)e −iHa · · · e −iHa W (a; 0) (11) In [41] it was shown that obtaining a time-separated correlator of two gauge links (i.e. a temporal Wilson loop) requires a second-order derivative to be taken numerically.
Here, perturbations to the Hamiltonian occur at every time slice between the two operators (necessarily, so that the Wilson line is approximately light-like). The order of the finite differencing needed is equal to the number of time slices affected. This high-order finite differencing is not practical even in the absence of quantum noise. Nevertheless, this is the only candidate we are aware of for directly computing correlators of the form Eq. (10) on a quantum computer.
Fortunately an alternative procedure can be constructed: one may compute an easier observable -the hadronic tensor -and extract the PDFs after the fact.

III. HADRONIC TENSOR
Unlike the PDF, the hadronic tensor is constructed of currents, J µ , which are each gauge-invariant, unlike the fermionic operators. Thus, the hadronic tensor does not require a Wilson line and the issue of high-order finite differencing is avoided. The hadronic tensor of a d-dimensional theory is given explicitly by for a given current J µ , where |P denotes a proton in the zero momentum frame. Here we will assume J µ =ψγ µ ψ, corresponding to the current coupling to the photon. In combination with the leptonic tensor L µν , W µν (q) may be directly related at leading order (in the QED coupling α) to the cross section of lepton-proton scattering via the exchange of a photon with momentum q by [1] In these equations, Q 2 = −q 2 , x = Q 2 /2P · q, y = P · q/P · k, and k = k − q.
With an eye toward implementation on a quantum computer, Eq. (12) has two critical features. First, it is gauge-invariant without the need for a Wilson line. In fact, because it involves only gauge-invariant operators, it can be defined and measured without reference to unphysical gauge-variant (Gauss-law-violating) states. Second, each operator in the correlator is individually Hermitian, so that the decomposition procedure of the previous section is not required.
Several options exist for measuring the integrand of the hadronic tensor. One may follow the procedure of the previous section closely, decomposing J µ (x) as a sum of N unitary operators. After this, the operator in the integrand becomes a sum of N 2 unitary operators, each of whose expectation values may be directly evaluatedà la [52]. However, as mentioned above, this is needlessly complicated for the hadronic tensor.
In this context, the procedure of [55] for studying linear response is more straightforward. Consider a unitary evolution operator The first derivative of the expectation value of this operator with respect to either vanishes. The second derivative gives the desired correlator.
Finally, a more sophisticated procedure for the calculation of linear response is given in [56].
After measuring Eq. (16) at many values of x, the Fourier transform may be taken classically, via Eq. (9). Alternatively, the (regulated) Fourier transform may be subsumed into the expectation value, yielding This expression requires only one quantum circuit for a desired value of q; however, when many values of q are to be obtained, it is no more efficient.
There is an important way in which these procedures are not analogous to those performed in the laboratory: on a quantum computer, one may introduce a current coupling only to particular flavors of fermions. This allows one to isolate a single-flavor distribution function or hadronic tensor in a straightforward way, without any fitting.
To apply this method to the Thirring model, one needs the staggered form of J µ : As mentioned, the leading order cross section for leptonhadron scattering may be computed once the hadronic tensor is in hand. This is not the first proposal for computing a scattering cross-section on a quantum computer; in [45][46][47] is detailed a procedure in which two asymptotic states are prepared adiabatically on a large lattice, and then allowed to propagate towards each other. When obtaining a cross section via the hadronic tensor, the need to prepare two asymptotic states is removed-reducing substantially the cost of state preparation. Instead, we prepare only a single asymptotic (zero-momentum, in fact) state, and probe it with arbitrary momentum. Additionally, this avoids the long-range nature of the QED interaction that complicates lattice calculations. This procedure is substantially simpler, but the tradeoff comes in that while [45][46][47] computes the full cross section, our procedure is perturbative: to obtain higher-order contributions in α, one must calculate multiple matrix elements defined by additional current insertions.
If one ultimately wants the PDF, one can extract it from W µν (q) via a procedure analogous to how the experimental determinations are done. To first review, there are a number of processes where collinear factorization can be proven (e.g. deep inelastic scattering, Drell-Yan, weak boson production, and inclusive jet production). Here we consider deep inelastic scattering, but similar expressions are derivable for the other processes. The cross section σ eP →eX can be schematically decomposed into where i and j run over all species of parton, f i are parton distributions, P i→j is the splitting function required to match all experimental data at a single scale and can be computed perturbatively, and σ ej→ej is the hard partonic cross section. Theoretical expressions like Eq. (20) are used to numerically fit parameterized PDFs to the experimental data over large ranges of kinematics [1]. The complicated nature of the perturbative splitting function and hard cross section in addition to the need to perform two convolutions prove to make this process highly nontrivial.
In the same spirit, the hadronic tensor that would be obtained by a quantum computer can be defined in terms of the PDFs as (21) whereŴ µν are partonic tensors that couple to external currents. Thus, if one desires the PDFs, they can be extracted by numerical fits to parameterized PDFs when the hadronic tensor is computed in the kinematic regime of collinear factorization's validity. It is important to note that our procedure can, by allowing different four-momentum for the hadronic states, be trivially generalized to computing Generalized Parton Distributions [57] and similar small changes to obtain other distributionssomething that is not trivially possible for Euclidean field theory.

IV. STATE PREPARATION
Thus far we have neglected to discuss the preparation of the state |P . This is not a trivial matter, and in this section we rectify the situation. For concreteness, let us assume that we are preparing a zero-momentum proton on a three-dimensional QCD lattice.
A great deal of literature discusses the problem of preparing a ground state of a given Hamiltonian [45-47, 58, 59]. Typically, formal analysis of the efficiency of ground state preparation methods is not available, and the in-practice performance cannot yet be measured. Here we consider adiabatic state preparation [45][46][47], where we will be able to make some crude estimates of the cost of preparing a proton. Note that this does not imply adiabatic state preparation is the most efficient method in practice, simply that it is the easiest to analyze without the need to test on a full-scale quantum computer. Other possibilities for ground state preparation, not to be discussed here, are the spectral comb [58] and hybrid methods [49,59].
When a method is cast as preparation of a ground state, it is typically still applicable to the preparation of a state like |P : the preparation of the proton may be translated to ground state preparation. To do this, consider only the sector of Hilbert space which is translation-invariant (thus zero-momentum) and has quantum numbers of the proton. As long as time evolution e −iHt maps this sector to itself (as standard Trotterized time evolution does), we may safely use the algorithm to prepare the "ground state" of this sector specifically, which is of course the zero-momentum proton. Note that this "trick" is exactly the same as what is needed to prepare a gauge-invariant ground state, where we must restrict ourselves to the physical subspace of the Hamiltonian.
We now estimate the costs of adiabatic state preparation in the context of QCD. We make several assumptions about the spectrum of the lattice model. In the restricted Hilbert space with baryon number 1, the lowest-lying state should be the zero-momentum proton. If simulating pure QCD (that is, in the absence of weak interactions), we may further restrict to the states with isospin of the proton, thus avoiding a small gap between the proton and the neutron. Finally, on a finite volume the gap between a particle and a slowly moving particle is O( 1 L ). If we restrict to the zero-momentum subspace, then the gap between the proton and the nearest-energy state is the pion mass, m π ≈ 135 MeV.
Adiabatic state preparation begins by preparing the ground state for some modified Hamiltonian, for which the ground state is known with great precision. In this case we chose the free theory. The ground state of the baryon-number-1 sector is three zero-momentum fermions in a box. As the gauge coupling is 0, the configuration of the gauge fields in the ground state is also Gaussian, and may be prepared efficiently.
Adiabatic state preparation proceeds by slowly deforming the Hamiltonian, over a time T , from the initial (in this case, free) Hamiltonian H 0 to the desired Hamiltonian H T . The simplest trajectory we can pick increases the coupling from 0 to its desired physical value at a constant rate. The deformation must be done slowly, and more slowly when the gap in H t is small. The adiabatic theorem guarantees we will remain in the ground state (with high probability) as long asḢ/∆ 2 1, where ∆ is the gap andḢ is the rate of change of the Hamiltonian [60]. Thus to estimate the performance of the adiabatic procedure we must estimate the gap along the trajectory. At the end of the trajectory, as mentioned, the gap is large: about 1/7 the mass of the proton. This part of the evolution can be done quickly. At vanishing coupling, the outlook is less rosy. The 'proton' fills the lattice, and excited states are simply back-to-back low-momentum states of two fermions. (The massless glue excitations can be removed with appropriate boundary conditions, or by using the 1080-element approximation to SU (3) [41,42].) The gap, therefore, is O(1/L), and we see that the adiabatic algorithm will require O(L 2 ) time-evolution steps to keep the ground state.
Not all ground-state preparation methods may have this scaling. As an example, spectral combing [58] does not require a trajectory starting from weak coupling (thus avoiding the 1/L gap), and appears in numerical studies to scale like ∆ −1 . However, in the absence of large-scale tests, it is unclear how faithfully the method actually prepares a ground state.

V. DISCUSSION
With the procedure completely described, we may estimate what manner of quantum resources are required for a practical calculation of the hadronic tensor of the proton. For the sake of specificity we will assume that the proposal of [41,42] is used; namely, that SU (3) gauge theory is approximated by the 1080-element discrete subgroup, and the calculation is done in A 0 = 0 gauge with no further gauge fixing. The qubit cost of other methods is expected to be similar. The applicability of the S(1080) approximation is limited by the lattice spacing at which the lattice theory undergoes a phase transition. It was shown in [42] that it remains a good approximation down to a lattice spacing of a = 0.08 fm, suggesting that it is sufficient for low-and intermediate-momentum probes of hadronic physics [13].
The qubit costs are easiest to count: representing a single element of the group S(1080) on a link must cost at least 11 qubits, and therefore an L 3 lattice requires, taking into account a 12-component (4 spinor by 3 color) Wilson fermion at each site, ∼ 50L 3 qubits to store. Timeevolution brings in the need for some number of ancillary qubits, but this need not scale with the volume of the system, and must therefore be negligible. Assuming a lattice spacing a = 0.1 fm and a 20 3 lattice -chosen to be well within the range of applicability of S(1080) and large enough to fit a proton with moderate finitevolume effects, ∼ 4 × 10 5 qubits are required to perform the calculation.
Estimating the gate cost associated with the calculation is difficult, not least because circuits have not yet been produced for any concrete proposal for simulating SU (3) gauge theory. Nevertheless, we can give the scaling of the algorithm with volume. The procedure is dominated by state preparation, which requires time evolution for O(L 2 ) steps. In computing the hadronic tensor we must evaluate L 3 × T matrix elements, where T gives the length of time evolution needed in approximating the Fourier transform. As the J µ (x) at a single time are mutually commuting, L 3 measurements can be performed simultaneously, without needing to re-evolve the system. Assuming the evolution time T to be proportional to L, we find that O(L 3 ) time evolution steps will be required. Each time evolution step, of course, scales with the volume of the lattice, so the total scaling of the procedure should be O(V 2 ). This is comparable to the scaling of calculations performed on the Euclidean lattice.
In this paper we have detailed a possible application for large-scale quantum computers beyond the NISQ era: first-principles calculations of parton physics in field theories, particularly QCD. As a side effect, this would allow the calculation of hadron-lepton cross sections in the standard model, more cheaply than existing proposals. Although not explored in this paper, a nearer-term target of this method is the physics of bound states in fewer dimensions.