Simulating Hadronic Physics on NISQ devices using Basis Light-Front Quantization

The analogy between quantum chemistry and light-front quantum field theory, first noted by Kenneth G. Wilson, serves as motivation to develop light-front quantum simulation of quantum field theory. We demonstrate how calculations of hadron structure can be performed on Noisy Intermediate-Scale Quantum devices within the Basis Light-Front Quantization framework. We calculate the light-front wave functions of pions using an effective light-front Hamiltonian in a basis representation on a current quantum processor. We use the Variational Quantum Eigensolver to find the ground state energy and wave function, which is subsequently used to calculate pion mass radius, decay constant, elastic form factor, and charge radius.

Quantum simulation of quantum field theory (QFT) is a promising application of quantum computing that has recently seen a surge in interest . In our previous work [33] we demonstrated that the light-front (LF) quantization of quantum field theory provides a natural framework for ab initio arXiv:2011.13443v1 [quant-ph] 26 Nov 2020 digital quantum simulation. The Discretized Light-Cone Quantization (DLCQ) technique allows one to significantly reduce the quantum-computational resources when compared to lattice approaches.
In the present paper, we continue our program of investigating quantum simulation in the light-front (LF) formulation, and develop an approach to simulating field theory based on the Basis Light-Front Quantization (BLFQ) [34,35] technique. The basisfunction expansion allows us to further reduce the need for computational resources, making calculations accessible for existing quantum devices.
In [33] we developed quantum algorithms based on simulating time evolution and adiabatic state preparation. In this work, we instead aim for near-term devices by adopting the Variational Quantum Eigensolver (VQE) paradigm: this allows us to implement a demonstration on the IBM Vigo quantum processor. VQE is a hybrid quantum-classical algorithm. The classical computer optimizes the expectation value of the Hamiltonian, which is repeatedly evaluated on a quantum computer. The classical computer/algorithm optimizes the parameters of this ansatz to minimize the ground state energy.
The resulting parametrized quantum circuit approximately prepares the ground state wave function of the Hamiltonian. Thus once the VQE procedure is complete, we can compute the expectation values of other observables in this approximate ground state.
The DLCQ and BLFQ paradigms provide alternative approaches to describing relativistic interactions. While both are, in principle, ab initio frameworks, DLCQ studies the system starting from the light-front Hamiltonian quantized in the traditional free-field basis placed on a discrete momentum grid. BLFQ starts from the same light-front Hamiltonian but quantizes it in terms of modes tailored to the symmetries of the system under consideration in order to construct an computationally efficient representation of the Hamiltonian. Since each represents a choice of basis spaces for the fields, they should yield the same results in fully converged calculations, i.e., in their respective continuum limits.
Having much in common with ab initio quantum chemistry and nuclear theory, the BLFQ formulation provides an ideal framework for benchmarking NISQ devices and testing existing algorithms on physically relevant problems such as the calculation of hadronic spectra [36][37][38][39] and parton distribution functions (PDFs) [40][41][42]. In essence, BLFQ amounts to (1) Choosing the effective field theory most efficiently describing the problem of interest, (2) Quantizing the system in the light-cone coordinates, (3) Nonperturbatively solving the theory in the most suit-able basis. This results in an efficient representation of the QFT problem under study. One typically starts with a fixed-particle-number formulation, effectively reducing the QFT setting to a relativistic quantum-mechanical many-body problem. In many cases, already at this level one can obtain results with suitable precision to make meaningful comparisons with experimental results [37][38][39][40][41][42][43].
In this article, we consider the dynamics of valence quarks for light mesons on the light front using the Hamiltonian from [43]. This Hamiltonian includes the kinetic energy, the confinement potential in both the longditudinal and the transverse directions [37], and the Nambu-Jona-Lasinio interaction [44] to account for the chiral interactions among quarks. We limit ourselves to the valence Fock sector of mesons while working with relative momentum variables. The dependence of the light-front wave functions for these valence quarks on the relative momentum is expanded in terms of orthonormal basis functions. After implementing finite cut-offs in this expansion, the light-front Hamiltonian becomes a hermitian matrix in the basis representation. We use the same scheme as in Ref. [43] to fix our model parameters at each choice of basis cut-offs. We illustrate our ideas by choosing an effective Hamiltonian for the light meson system, and running the VQE minimization on the IBM Vigo machine to calculate the squared pion mass. Using the resulting wave function, we calculate squared mass, decay constant, mass radius, electromagnetic form factor, and charge radius of the pion.
The two different ansätze we consider in this paper are based on different ways of encoding physical states on the quantum computer. Within the direct encoding one stores the occupancies of the secondquantized states in the unary form and uses the Unitary Coupled Cluster ansatz for state preparation. Within a more efficient compact encoding, one stores the occupancies in the binary form, which requires logarithmically fewer qubits and allows one to prepare an ansatz state using the arbitrary state preparation algorithms, given that the particle number is fixed and small.
In Sec. II, we provide a summary of the BLFQ formalism and a representation of basis functions, which we used throughout the paper. In Sec. III, we derive expressions for various observables in the chosen basis. In Sec. IV we describe two variations of the VQE algorithm, and show the results of running it on an existing quantum computer.

A. Overview
The light-front quantization approach specifies the commutation relation of fields at equal light-front time [45]. In contrast to the Lagrangian formulation of equal-time quantization, the field theory dynamics after light-front quantization is governed by a light-front Hamiltonian [45] responsible for the light-front time evolution of the system. Quantizing a QFT on the light front has a number of advantages: triviality of the vacuum, absence of ghost fields in the light-cone gauge, Hamiltonian sparsity, and the simple form of observables in terms of the wave functions. Within the light-front approach, the bound state masses and associated wave functions are solvable from the light-front-time-independent Schrödinger equation.
In [33] we developed a simulation algorithm based on the DLCQ, allowing one to reduce the resource requirements for the ab initio simulation of QFT by a few orders of magnitude. In this work, we further reduce the computational requirements into the range of the capabilities of existing quantum devices. We achieve this by employing the framework of BLFQ [34,35].
Within BLFQ, a field is expanded in terms of second-quantized Fock states representing occupancies of modes (first-quantized basis functions), and there is no a priori limit on the degrees of freedom [34,46]. Accordingly, our algorithms are designed to efficiently simulate QFT applications where particle number is not conserved. However, for QFTs at low resolution, BLFQ is often restricted to the valence degrees of freedom, allowing the adoption of this restriction in order to implement quantum simulations on an existing quantum chip. These experiments represent the first stage shown in Fig. I, which illustrates a progression of methods that scale towards fault-tolerant simulation of QFTs in the quantum supremacy regime. However, the methods we propose apply to the first three stages in Tab. I. The final stage was discussed in [33].
Previous development of BLFQ for the heavy mesons is partially based on the holographic confinement potential between the valence quark and antiquark in the holographic transverse directions [37][38][39]. This potential is supplemented by a longitudinal confinement potential to attain a 3-dimensional spherical confinement potential in the nonrelativistic limit. These potentials are constructed independent of the spins for the quark and the antiquark and they are governed by a single overall strength parameter.
In addition to the kinetic energy and the confinement potentials, they form the baseline Hamiltonian that is analytically solvable and defines our basis functions [37,43]. These basis functions possess desired spatial symmetries and boost invariances. The derived effective one-gluon exchange interaction based on the gauge dynamics serves as the spin-orbit interaction and incorporates a running coupling [37].
In this article we adopt the Hamiltonian in Ref. [43] for the light mesons. Specifically the same confinement potential forms as those in Ref. [37] are implemented. However, we do not include the one-gluon exchange because the interactions for light quarks manifest from the chiral symmetry, which is insufficiently accounted for by a perturbative expansion of the gauge interaction. Instead, we resort to the Nambu-Jona-Lasinio (NJL) model for the chiral interaction of these quarks [44,47,48]. Within our basis representation, the matrix elements of the NJL interaction can be calculated analytically [43]. We compute the lowest mass eigenvalue and its corresponding eigenvector, its light-front wave function, using the algorithm to be described in Section IV. We then calculate observables based on this eigenvector.

B. The effective Hamiltonian of the BLFQ-NJL model
The light-front wave functions (LFWFs) of the valence quarks for the π + meson and the K + meson have been solved from Ref. [43] in the basis light-front quantization (BLFQ) framework using Nambu-Jona-Lasinio interactions [44,[47][48][49] on a classical computer. Specifically, one first truncates the light-front wave-function for the mesons to the valence quark Fock sector such that the state vector is expressed as where P = k+p is the total momentum of the meson, x = k + /P + is the longitudinal momentum fraction carried by the valence quark, and κ ⊥ = k ⊥ − x P ⊥ is the relative transverse momentum. In order to solve for the LFWFs for the valence quarks inside light mesons, we adopt the effective Hamiltonian that can be represented as a basis- diagonal term and the NJL interaction: The basis-diagonal term H 0 contains the kinetic energy of the valence quarks, the transverse confinement potential, and the longitudinal confinement potential. In the valence Fock sector of mesons, this term takes the form of where x is the longitudinal momentum fraction carried by the valence quark and κ ⊥ is the relative transverse momentum of the valence quarks. The masses of the valence quark and the valence antiquark are given by m and m, respectively. In a addition, b specifies the strength of the confinment potentials. This part of the Hamiltonian has analytic solutions that constitute the basis states for the BLFQ approach as will be seen in detail in Subsection II C.
When quarks in the confinement region are the retained degrees of freedom, the strong interaction among them can be understood to arise from the global chiral symmetry, an approximate symmetry of quantum chromodynamics. To model this chiral interaction, we employ the interaction in the scalar-pseudoscalar channel of the color-singlet NJL model [44]. Specifically, we ignore both the instantaneous interaction and the self-energy correction from the NJL interaction to obtain the following term in the total Hamiltonian: Here ψ is the fermion field operator, G π is the NJL coupling constant, and P + is the total light-front longitudinal momentum of the system. We then expand eq. (4) into relevant combinations of ladder operators for the quark fields. In the basis representation, this term further takes the form of a hermitian matrix, the elements of which can be calculated analytically [43].
In this work, we solve the eigenvalue problem defined by eq. (2) in the total angular momentum J z = 0 block with the lowest eigenstates of H 0 forming the longitudinal and radial basis states for the interacting Hamiltonian. In this representation, the effective Hamiltonian takes the form of a 4-by-4 matrix indexed by the basis quantum number θ that specifies the angular and spin excitations. The explicit expressions for elements in this matrix are given in Appendix A 2.

C. The basis function representations of wave functions for valence quarks of mesons
We adopt the following expansion of the light-front wave function for the valence quarks given by eq. (1): where ψ nmlrs is the expansion coefficient, φ nm is a 2-dimensional (2D) harmonic oscillator (HO) eigenfunction, and χ l is the longitudinal basis function.
Here r and s are the spin indices of the quark and the anti-quark. Each term in eq. (5) is an eigenfunction of H 0 in eq. (3). Explicitly, φ nm is defined as with tan(ϕ) = q 2 /q 1 and L |m| n being the associated Laguerre function. The parameter b sets the scale of the harmonic oscillator eigenfunction, which we choose to be identical to the confining strength in the light-front Hamiltonian. Meanwhile, χ l (x) is given by with P (α,β) l (z) being the Jacobi polynomial and When we solve the eigenvalue problem defined by the BLFQ-NJL Hamiltonian, the following cutoffs on the basis quantum numbers following Ref. [43] are imposed: Because truncations on different basis quantum numbers are independent, we call this truncation scheme the orthogonal enumeration. Such a scheme allows us to solve simultaneously for eigenstates with different azimuthal angular momentum projection J z since it is a good quantum number in this basis. The size of the Hamiltonian in the basis representation with this orthogonal enumeration is n H -by-n H , with n H = 4(N max + 1)(2M max + 1)(L max + 1) .
However, the capacity of NISQ devices motivates further reduction in the dimension of the Hilbert space spanned by our basis representation. Specifically, because eigenfunctions of this Hamiltonian have fixed azimuthal angular momentum projection J z , the basis quantum number θ indexes specific combinations of the spin and orbital bases in the orthogonal enumeration as specified in Appendix A 1. Each basis state in the fixed J z block is then given by the basis quantum numbers n, l, and θ. In the limit of M max = 2, the unitary transformation that relates the bases in the fixed J z blocks to those in the orthogonal enumeration is given by Table V. The degeneracy in the basis quantum number θ in each J z block is apparent in Table V. For example when J z = 0, this degeneracy is d θ = 4. With a given set of (n, l, θ) in a given J z block, we take the convention such that the index of this basis is given by For a given index, the corresponding basis quantum numbers can be easily calculated. Consequently, the size of the Hamiltonian for a fixed J z in this new enumeration becomes which is smaller than n H given by eq. (10). This provides an example of how one may exploit the symmetries embedded in the chosen BLFQ to achieve gains in computational efficiency.

III. COMPUTING OBSERVABLES FROM THE VALENCE LFWF
One of the many advantages of the light-front approach to quantum field theories is that observables for bound states can be easily extracted from lightfront wave functions. Explicitly, measurement operators corresponding to physical observables usually take a simple form, resulting in efficient measurements on a quantum computer (see Sec. IV and App. D). In this section, we demonstrate how to calculate the decay constant, mass radius, valence parton distribution function, and elastic form factor.

A. The decay constant
The meson decay constants are defined as the matrix elements of current operators between the vacuum and the meson wavefunctions [37]. They correspond to amplitudes of the wavefunctions at the coordinate-space origion. Specifically, the decay constants for scalar mesons (f S ), pseudoscalar mesons (f P ), vector mesons (f V ), and axial vector mesons respectively. Here the polarization vector for the vector mesons is defined as with e ⊥ ± = (1, ±i)/ √ 2. In terms of the valence-sector light-front wavefunctions, expressions for these decay constants are reduced into [37] with the condition m J = m + s 1 + s 2 = 0 specifying that only the states with zero angular momentum projections are used in the calculation. Here N c = 3 is the number of colors.
In the basis representation, the integrals over the longitudinal momentum fraction and the relative transverse momenta in eq. (15) can be evaluated exactly. Details of this calculation can be found in Appendix B 1. Since the decay constant is linear in the wave function, we only need to calculate these integrals for each basis function. Subsequently, the decay constants in the basis representation are given by where the longitudinal integrals L l (a, b; α, β) are defined and given analytically in Appendix B 3. Because the overall phase of the LFWF remains undetermined by the Hamiltonian, only the absolute value of the decay constant carries physical significance. Once the LFWF |ψ in our basis representation is known on a quantum computer, the calculation of the corresponding decay constant can be thought of as computing | v|ψ | for some fixed |v .

B. The mass radius
The mass radius is the square root of the expectation value for the relative transverse separation of the valence quarks. It can be calculated from the valence two-body wave-function based on eq. (33) of Ref. [37]. Specifically for the pseudoscalar mesons, we have whereψ rs x, r ⊥ is the light-front wave-function depending on the longitudinal momentum fraction x and the relative transverse coordinate r ⊥ . It is related to the momentum-space wave-function by the Fourier transform in the transverse momenta κ ⊥ . Explicitly, we havẽ with and tan φ r = r 2 /r 1 .
To calculate the mass radius in terms of expansion coefficients ψ nmls1s2 , we first need to evaluate the fol-lowing dimensionless integrals of the basis functions: We then have the square of the radius given by The explicit expression for the matrix I m (n , m , l ; n, m, l) is available in Appendix B 2, which takes the form of a hermitian operator in our basis representation.

C. Parton distribution function of valence quarks
The probability of finding a quark inside a meson carrying momentum fraction x is given by which is interpreted as the PDF for the valence quark. The PDF for the valence antiquark is given by f (1 − x) [40][41][42]. We use the solutions of ψ nmlrs defined in eq. (5) to calculate the valence PDFs of mesons.
Notice that eq. (23) defines a bilinear of the lightfront wave-functions in the basis representation. To compute the PDF with the LFWFs obtained from a quantum computer, let us rewrite eq. (23) as Elements of the density matrix ρ l , l defined in eq. (24b) can be evaluated as the expectation value of the corresponding projection operators on a quantum computer, and subsequently used to calculate the PDF in eq. (24a).

D. The elastic form factor for pseudoscalar mesons
To calculate the elastic form factors from the lightfront wavefunctions within the impulse approximation where the photon interacts with the meson through the quark-photon vertex, we apply the following formula [50,51] within the Drell-Yan frame P + = P + : with q = P − P and Q 2 = −q 2 . The operator inside the Dirac bracket is the charge density operator on the light front, with e f being the charge carried by the quark of flavor f in units of the elementary charge and the summation running over all quark flavors. Additionally, e q is the charge of the quark (e u = +2/3 for an up quark). While e q is the charge of the antiquark (e d = −1/3 for an anti-down quark). Detailed derivation of eq. (25) is given in Appendix C 1.
In the basis representation, we apply the Talmi-Moshinsky transform to simplify the integrals in the transverse momentum, leaving the longitudinal integral to be evaluated numerically for each Q 2 . Following steps in Appendix C 2, we rewrite the electromagnetic form factors into a bilinear form of the valence wave-function: The operatorC is defined according to eq. (C19). At a given Q 2 , the form factor can be calculated using the LFWFs obtained from a quantum computer by taking the expectation value of the Hermitian oper-atorC(n , m , l ; n, m, l; Q 2 ).
Specifically, the elastic form factors of the pseudoscalar mesons are given by The charge radius is then specified by the first Taylor expansion coefficient of the elastic form factor at the origin:

IV. QUANTUM-COMPUTATIONAL METHODS
The standard approach to quantum simulation of Hamiltonian dynamics of QFTs is as follows [52]: a) initialize the system in a certain state of the free Hamiltonian, b) adiabatically turn on the interaction, c) if necessary, evolve the system with time, and d) measure the energy (or another observable) of the system using the phase estimation algorithm [2,7,52]. However, near-term devices cannot perform this procedure due to limits in qubit numbers and in gate fidelities. This motivates the variational quantum eigensolver (VQE), an approach to finding Hamiltonian eigenvalues in which a NISQ device is used as a part of a hybrid quantum-classical algorithm [53]. In VQE, a quantum comp uter prepares a given variational state and evaluates the Hamiltonian expectation value, which a classical computer performs a gradient search to minimize (see Fig. 1). To prepare the variational state, we adopt an ansatz specified by parameters θ, which are controlled by the classical minimization.
While in [33] we focused on ab initio simulations which are likely to become available in the faulttolerant regime, in this paper we investigate the use of NISQ devices for high-energy nuclear physics calculations on the light front. Therefore, unlike in [33], we formulate the problem as a VQE instance.

Classical Minimization Algorithm
Quantum Computer

Updated parameters θ
The value of the cost function E( θ) = ψ( θ)|H|ψ( θ) We begin by briefly reviewing the VQE method. For the VQE algorithm to be efficient and accurate, it is essential to come up with a parametrized ansatz state |ψ( θ) that is easy to prepare and is expected to have significant overlap with the true ground state. Below we shall consider different choices of ansatz state preparation prodecures, encodings of the physical states in a quantum computer, and classical optimization algorithms.
While using VQE for simulating a Hamiltonian problem, the major steps are: 1. Define the state/operator mapping, i.e., a correspondence between the physical states and the multi-qubit states of a quantum computer, as well as the mapping between the operators acting on these spaces.
2. Choose a parametrized ansatz state. One typically writes the ansatz state as where |ψ 0 is a fixed reference state, and U ( θ) is the VQE ansatz operator. One possibility is to choose the form of U ( θ) to resemble the form of the Hamiltonian evolution operator e iHt [54].
3. Once the state |ψ( θ) is prepared for a given set of parameters θ, one evaluates the cost function by measuring the expectation value of the multi-qubit Hamiltonian operator: The algorithm can only be considered efficient if the number of measurements grows polynomially with the problem complexity (discussed further below).
4. The value of the cost function is then sent to the classical optimizer, which either determines the set of parameters for the next iteration of the algorithm, or terminates the algorithm if the desired precision has been achieved.
We shall explore two approaches to simulating problems in the BLFQ formulation, based on two different encoding schemes. The first of these is the direct encoding, widely used in quantum chemistry [7,55]. In such an encoding scheme, one assigns a particular set of qubit registers to each physical (basis) degree of freedom. In application to purely fermionic systems, one may use one qubit to encode one fermionic second-quantized mode, which leads one to the Jordan-Wigner (JW) encoding [56]. Thus, one needs N qubits in order to encode N fermionic modes. The fermionic raising and lowering operators are represented by N -local multi-qubit operators, due to the need to enforce anticommutation relations. One can alternatively employ the Bravyi-Kitaev (BK) encoding [57][58][59][60] that uses N qubits to store N fermionic modes, with operators being only log N -local. Circuits implementing VQE ansatz operators are typically based on trotterization [54].
The second encoding we employ is compact encoding, as was explored in [33] for front-form physics, and in the context of quantum chemistry in [7]. The idea is to only store the occupied modes of multiparticle Fock states. With this encoding, one can simulate time evolution using sparsity-based techniques, which are optimal in all parameters [17,33,[61][62][63][64]. These methods, however, require too many gates to be used to prepare ansätze on NISQ devices. Instead, we can use arbitrary state preparation as long as we restrict to small fixed numbers of particles in the system, as discussed below.

A. Direct encoding
In order to run a simulation on a current quantum device, in this work we consider a scenario where the particle number is fixed. However, since quantum advantage is likely to be achieved only in the multiparticle regime, it is essential for our methods to be extendable to this more general scenario.
A natural way of formulating a multi-particle problem is by using the second-quantized formalism. Consider a Hamiltonian of the form where Here h ij represents the single-body interactions, while h ijkl and higher-order terms correspond to many-body interactions. For the first experimental implementation, we restrict ourselves to H = H 1 (with h ij being the meson valence sector BLFQ Hamiltonian matrix), for two reasons. First, owing to the efficiency of the BLFQ formulation, considering the single-body part of the Hamiltonian is oftentimes enough to give reasonably good results [36,40,42,43,65]. Second, this is suitable for benchmarking, paralleling state-of-the-art experimental results in quantum simulation of chemistry [66].
Within the JW encoding, the multi-qubit states | . . . f 2 f 1 f 0 mimic the second-quantized fermionic states: the qubit f i stores the occupancy of the (i + 1)-th orbital (see Table II). In order to enforce anticommutation relations, the fermionic creation and annihilation operators are represented by N -local multi-qubit operators [56]. We shall use this encoding for the rest of the section; simulation in the Bravyi-Kitaev encoding is discussed in App. E.

Basis state index Direct encoding Compact encoding
|1000 |11   Table II. The multi-qubit representation of the four physical one-particle states in the direct and compact encodings.
Since one typically solves the problem in a basis found by means of some classical approximation, the reference state |ψ 0 can be chosen to have a simple form in terms of basis vectors; in the simplest case, it may coincide with one of the basis vectors. Next, we would like to design an ansatz operator that acts on the reference state to prepare an ansatz state that ideally has large overlap with the exact ground state. An example of such an operator is the Unitary Coupled Cluster (UCC) [54]. Choosing the form of the ansatz operator to resemble the form of the Hamiltonian ensures that one can explore the regions of the Hilbert space that can be reached via the Hamiltonian time evolution, and also guarantees that the symmetries are preserved. For the Hamiltonian of the form (32) one writes the UCC as [54] U ( θ) = e T −T † , T = T 1 + T 2 + . . . , where occ and virt denote occupied and unoccupied orbitals in the reference state |ψ 0 . Physically, the action of the UCC operator allows one to transfer "some amplitude" from initially occupied orbitals to the unoccupied ones. For real Hermitian Hamiltonians, the coefficients in (33) are real.
We would now like to translate (33), which was written in terms of the fermionic operators, into its qubit representation. According to the JW transformation [56], the qubit operators are obtained as where X i , Y i , Z i are the Pauli matrices acting on qubit i. Substituting (34) into (33) generates a mapping: where P j are the Pauli operators, while α j are the corresponding real coefficients. Trotterization of the expression above leads to where ρ is the Trotter number, which can be typically chosen quite small in VQE [67], unlike the case of simulating time evolution.
The traditional approach to calculating the expectation value as in (30) amounts to expanding the Hamiltonian in the basis of Pauli operators using (34): The expectation values of individual Pauli terms on the RHS of (37) can be efficiently measured via sampling from the state |ψ( θ) [53]. The optimal parameters θ * , obtained upon successful termination of the VQE algorithm, allow one to prepare the VQE approximation to the ground state of the system. By analogy with (37), this can be used to calculate the expectation value of any observable bilinear in the wave function (i.e., of the form ψ( θ)| O|ψ( θ) ), such as mass radius, PDF, or elastic form factor. Observables linear in the wave function (i.e., of the form v|ψ( θ) , where |v is a constant vector), such as the decay constant, eq. (15b), can be calculated using the simple circuit shown in Fig. 2.

Efficiency analysis
When proposing new algorithms for quantum simulations on NISQ devices, it is essential to elicit their scaling properties in order to distinguish aspects of a particular simulation that may lead to quantum advantage from those that cannot. The question relevant to the present paper is which aspects of the few-qubit calculations we can scale up to several hundred qubits. Concomitantly, which aspects of the few-qubit calculations are amenable to efficient classical calculations and which are not. Figure 2. Estimating the magnitude of the inner product v|ψ( θ) for fixed |v . Up to the first dashed line, the circuit prepares the VQE ansatz state by applying the ansatz circuit U ( θ) to the first set of registers, resulting in the state |Φ1 = |ψ( θ) ⊗ |0 . The next rotation, R |v →|1 ⊗n , represents any unitary operator that maps the state |v to the state |1 ⊗n . Thus the state |Φ2 at the second dashed line is given by The quantity v|ψ( θ) = v|ψ( θ) 2 is found as the square root of the probability for the ancilla qubit to collapse into the state |1 .
The stages of one shot of a VQE calculation are ansatz preparation and measurement of all Hamiltonian terms. Many shots with fixed ansatz parameters are required to obtain one estimate of the expectation value of the Hamiltonian. Many estimations of the expectation value of the Hamiltonian are required to optimize the ansatz parameters. Typically the resources required to optimize a given VQE ansatz to a fixed precision cannot be bounded theoretically. This is what makes VQE a heuristic method. However we can determine the computational cost of each step in a single shot and ensure that the quantum gates and qubits required scale polynomially with the problem size. We can also be sure that no known efficient classical algorithm exists for large-scale versions of the problem.
As a prototypical example, consider the problem of finding a ground state in quantum chemistry using VQE. The parameters describing the complexity of the problem are the total number of orbitals, N , and the number of electrons in the system, M (i.e., the number of occupied orbitals). Using the direct mapping requires N qubits for encoding physical states. The second-quantized Hamiltonian operator can be written as a polynomial in ladder operators. Those, in turn, are each represented by a polynomial number of Pauli operators, each of which is at most log N -local. Therefore, the measurement of the Hamiltonian operator can be replaced with the measurement of a polynomial number of elementary operators. Each of those can be measured with precision using O( −2 ) samples [53,68].
All that remains is to quantify the operational resources for preparing the ansätze. In quantum chemistry in the direct mapping, these are typically prepared using a unitary coupled cluster (UCC) operator. The UCC ansatz operator including single and double excitations (UCCSD) contains O(N 2 M 2 ) free parameters. Application of the (trotterized) ansatz operator to the initial state is realized by a circuit containing a polynomial number of gates, since the action of fermionic ladder operators can be represented by a polynomial number of gates in the case of direct encoding.
Let us now see how these arguments can be naturally extended to the case of quantum field theory (QFT). First of all, we note that unlike in quantum chemistry, where the number of particles is conserved, QFT allows for processes of creation and annihilation of particles. Nevertheless, QFT does have an operator similar to the non-relativistic number operator -namely, the total momentum operator. Indeed, since this relativistic momentum operator commutes with the Hamiltonian, one can solve the problem within a Fock space sector of a fixed total momentum.
This analogy extends to the terms in the Hamiltonian.
In quantum chemistry, the Hamiltonian operator can be written as a polynomial of fermionic creation and annihilation operators, containing O(poly(N )) terms. In QFT, the secondquantized Hamiltonian operator can be written as a polynomial of ladder operators, containing O(poly(Λ)) terms, where Λ is the momentum cutoff. As in chemistry, in the direct encoding those can be represented by log Λ-local Pauli operators, whose total number consequently also scales as O(poly(Λ)).
To obtain a finite-dimensional Hilbert space in the equal time quantization, one would have to impose an additional cutoff on the number of excitations in each bosonic mode. However, in the LF formalism the maximum number of excitations is automatically limited by harmonic resolution K, the dimensionless light-cone momentum [33,69]. Within the BLFQ, the role of Λ and K is played by N max , M max , and L max cutoffs (introduced in Sec. II C). Therefore, all the resources for a single VQE estimation of the QFT Hamiltonian expectation value based on the direct encoding and UCC will grow polynomially in momentum cutoffs and precision.

B. Compact encoding
In our previous work [33] we explored the possibility of using the compact encoding for simulating physics on the light front. This amounts to only storing information about occupied modes in the Fock states. In the simplified setting considered in Sec. II B, due to the usage of relative coordinates, the only information we store is the index of the single occupied orbital. While in the direct mapping the index of the occupied orbital was stored in the unary form, requiring N qubits for N orbitals, in the compact mapping it is stored in the binary form, requiring log 2 N qubits for N orbitals. Therefore, in the case when the single-body Hamiltonian matrix h ij is of size N × N = 2 n × 2 n , one would use all the basis states of the 2 n -dimensional Hilbert space of n qubits (see Table II).
In the compact mapping, an equation analogous to (34) would contain an exponential number of terms on the RHS, thus making the usage of the UCC inefficient. Instead, one may employ any of existing arbitrary state preparation algorithms [70]. While their complexity is exponential in the number of qubits, in our case the number of qubits is itself logarithmic in the problem cutoffs.
For the direct encoding, in order to measure the expectation value of the Hamiltonian we can express the Hamiltonian in terms of Pauli operators. Any observable of size N × N = 2 n × 2 n can be expanded in the basis of 4 n = N 2 Pauli operators defined on n qubits: where c α are real. It should be emphasized that the logarithmic scaling of the number of qubits required as a function of the problem cutoffs implies that the Hilbert space dimension is polynomial in the cutoffs. This also implies that classical approaches to this problem are efficient. We are considering these specific initial problems as benchmarks where the results obtained can be compared to the known classical solution as an evaluation of the NISQ device itself.

Efficiency analysis
Within the VQE regime, the approach to quantum simulation based on the compact mapping is more efficient than the one based on the direct encoding when solving a two-body problem in the relative variable basis. As one starts to consider the problem in the multi-particle setting, the number of qubits required for storing physical states in the compact encoding is nearly optimal [33]. Despite that, one faces serious problems at the stages of state preparation and measurement. Since the complexity of arbitrary state preparation algorithms scales exponentially with the number of qubits, and the number of qubits grows linearly with the number of occupied modes, those algorithms can only be used if the number of particles is fixed and small. Of course, in principle, one could use sparsity-based techniques for state preparation, but this produces gate counts that are not feasible in the NISQ era. Therefore, coming up with a good ansatz for a multi-particle state in the compact encoding is an important task, which we leave for future work.
Another problem arises at the measurement stage: the number of Pauli terms in the expansion of the Hamiltonian grows exponentially with the number of qubits. This motivates the development of VQE techniques for sparse Hamiltonians.

V. RESULTS
In this section we describe numerical and experimental results of implementing VQE for a sample QFT problem, namely simulation of a pion in the minimal BLFQ representation. In order to run our simulation on an existing device, we shall use the 4 × 4 light meson BLFQ Hamiltonian from Sec. II B corresponding to J z = 0 sector in Table V ( We can analyze the VQE calculation in a few steps: a) Check that the classical optimizer is working correctly. To eliminate any errors arising due to sampling, we begin with evaluating the Hamiltonian expectation values exactly, using the statevector representation. b) Determine the number of steps required to reach the desired precision when evaluating the expectation value via sampling from the exact distribution. This gives the lower bound on the number of samples, and models the situation of using a "perfect quantum computer." c) Evaluate expectation values on the IBM Vigo quantum processor. d) Use error mitigation techniques to postpro-cess the results obtained on the quantum computer. We shall perform these steps using both the direct and compact encodings, evaluating the Hamiltonian eigenstate as well as other observables discussed in Sec. III.
The multi-qubit states representing the four physical basis states are shown in Table II. The states in the direct encoding can be thought of as JW-encoded states. Therefore, we use the JW transformation for calculating the corresponding multi-qubit Hamiltonian: where each term is a tensor product of single-qubit Pauli matrices I, X, Y, Z. (In what follows, we use this convention when expanding Hermitian matrices in Pauli terms acting on qubits.) As an ansatz operator, one could use the UCCS (no Doubles) operator. According to (33), T 1 will contain N terms, each of which is N -local in the JW encoding and log N -local in the BK encoding. In the former case, the circuit will contain O(N 2 ) gates, while in the latter only O (N log N ) gates. However, in order to further improve the gate count in the direct encoding-based algorithm, instead of the UCCS ansatz, we design a simple parametrized circuit of depth O(log N ) using O(N ) gates shown in Fig. 3a, which is capable of preparing an arbitrary superposition (with real amplitudes) of single-occupied states in the JW encoding. (This circuit is a generalization of the circuit proposed in [71] for preparing W N states.) The multi-qubit representation of the Hamiltonian in the compact encoding is obtained from (39), calculating the coefficients by applying (38): The ansatz state is prepared using the circuit shown in Fig. 3b, which prepares an arbitrary two-qubit state with real amplitudes. For both encodings, we need to estimate expectation values of Pauli operators. To do this, we first rotate to a basis in which the desired Pauli operator is diagonal, then measure single-qubit Z operators. The desired operator is a product of some set of single-qubit Z operators in this basis.   Figure 3. Ansatz circuits for preparing an arbitrary superposition of single-particle Fock states with real coefficients. For the direct encoding (a), we use a generalization of a circuit from [71] for preparation of WN states. For the binary encoding (b), we use arbitrary state preparation, with all single qubit rotations replaced by Ry(θ) gates, where Ry(θ) denotes a single-qubit rotation through an angle θ about the y-axis.
The classical optimization was performed using various algorithms from the python.scipy.optimize library.
In agreement with [54], the best convergence to the true ground state was achieved with L-BFGS-B [72] and COBYLA [73] methods. The latter showed better convergence, and it is used in all the following calculations. Depending on the choice of the initial guess state, the optimizer was typically reaching 4-digit precision after ∼ 10 1 -10 2 steps (for a good initial guess, such as (0, −1/ √ 2, 1/ √ 2, 0) T ) and up to a few hundred steps for randomly chosen initial state. In rare cases, the minimization was not converging.
Next, we determined the number of samples from the exact distribution required to reach the desired precision, which is expected to scale as O(1/ 2 ) [68]. To do so, we calculated the relative error for determining the Hamiltonian's expectation value in the true ground state using the classical simulation (the corresponding parameters of the circuits were obtained via the optimization at the previous stage). We performed 1000 experiments with a fixed number of samples, and calculated the RMS relative errors in determining the ground state expectation value over each set of experiments. The results on Fig. 4 indicate that on an ideal quantum computer we would need to generate ∼ 10 6 samples per Pauli term in order to reach 2% precision, and ∼ 4 · 10 6 samples to reach 1% precision. Fig. 7 shows the relative errors for the energy, decay constant, and mass radius, evaluated in the approximate ground state obtained via the VQE minimization procedure. The expressions for all observ-ORJ1XPEHURIVDPSOHV ORJ5HODWLYHHUURU 'LUHFWHQF &RPSDFWHQF ables are obtained from the corresponding BLFQ matrices in analogy with eqs. (40) and (41); the explicit expressions can be found in App. D. Note that all the observables have a dominant contribution from the unity term (IIII in the direct encoding and II in the compact encoding), whose expectation value is exactly 1. Therefore, in Fig. 7 we also show    Pion elastic form factor is used to calculate the charge radius, obtaining the values given in Tab. IV (charge radius is defined in eq. 28). Datapoints for the quantum simulation on the IBM Vigo processor used 8192 samples per term, with and without measurement error mitigation. The results measured on the quantum computer are in good agreement with the exact ones due to the strong contribution to the measurement operators from the identity term.
the expectation values for observables from which this term has been subtracted, which in certain cases improves the relative precision of results. The expectation values without the unit terms are the quantities actually measured on the quantum computer, while those including the unit terms are the physically relevant numbers, so the relative errors in both are of interest. In order to calculate the decay constant, one can use the circuit shown in Fig. 2 or Pauli measurements; we use the latter option to minimize the number of gates. The elastic form factors, eq. (27), are shown in Fig. 6, and the corresponding charge radii, eq. (28), are presented in the Table IV. In both cases, the results obtained on the quantum computer are in good agreement with the exact ones. This is to be expected, because the corresponding measurement operators have a large contribution from the identity operator.
With our choice of cutoffs, the calculation of PDFs in the compact encoding reduces to measurement of II, while in the direct encoding, it reduces to the projector onto the computational subspace (spanned by the single-occupancy Fock states).

VI. DISCUSSION
In our paper, for the first time we simulated high energy nuclear physics in the light front formulation on existing devices. We considered a detailed example in which we studied a relativistic analog of hydrogen, the pion. We studied this problem in a fixed particle number formulation as a benchmarking test for existing devices, and as a preparation for moving on to the mixed particle number formulation. Using the basis light front quantization (BLFQ) formalism, we demonstrated how small quantum computers can  be used for calculating baryonic spectra and various observables. Adopting an effective interaction (suggested by AdS/QCD correspondence [74]) and a set of basis functions (motivated by the particular problem of interest [34]) allows one to significantly reduce the computational resources, and to obtain reasonable results for realistic theories on devices having just a few qubits.
Within the VQE approach to quantum simulation, studied in the present work, we considered various encodings and state preparation procedures, some of which were naturally suggested by our experience in quantum chemistry. Together with our previous paper [33], this work defines a spectrum of methods for quantum simulation of quantum field theories. On this spectrum, one can move from restricted models, to be simulated on existing quantum devices, all the way to full ab initio simulation of QCD in 3+1D, to be simulated on future fault-tolerant quantum computers. Future work will expand and improve our methods on both ends of the spectrum. The next step in near-term simulation will be to switch to single-particle coordinates, providing the framework for mixed particle number simulations where quan-tum advantage is possible. Because the light-front Hamiltonian conserves the angular momentum in the z direction, the Hamiltonian in our basis representation can be diagonalized into blocks of fixed J z . This is equivalent to combining the spin quantum numbers s 1 and s 2 with the magnetic quantum number m to form a new quantum number θ. Specifically when M max = 2, the unitary transformation from the original BLFQ basis to this block-diagonal form is given by Table V. When N max = L max = 0, the light-front effective Hamiltonian in the J z = 0 block takes the form of a 4-by-4 matrix. The subscripts of the matrix then index the basis quantum number θ. The explicit expressions for these matrix elements are given in the following equations. Here κ is strength of the confining potential which we set identical to the basis scale b. The parameter G π is the coupling constant of the NJL interaction. Functions L (a, b) and L(a, b) both stand for L 0 (a, b; α, β) given in Appendix B 3.
Explicitly, the matrix elements of this Hamitlonian in our basis representation is given by Appendix B: Analytical expressions for integrals of basis functions

Integrals for the calculation of the decay constant
When calculating the decay constants using the valence LFWFs of mesons, we encounter the follow-ing integral: where L l (1/2, 1/2; α, β) is given by Eq. (B7). We have also used

Integrals for the mass radius
To evaluate the transverse integrals in Eq. (21), we first define ρ = b x(1 − x) r ⊥ . After this substitution of variables we obtain The integrals over the product of generalized Laguerre polynomials can be obtained by the orthonormality relations. In order to apply such relations, we convert L Here when n = 0, the second term drops out. We then have n !n! (n + |m|)!(n + |m|)!
(B5) with n ∈ N by default.
To evaluate L l (a, b; α, β) numerically, we first rewrite Eq. (B7) as L l (a, b; α, β) = 2l + α + β + 1 4π with C l,m ≡ (−1) l−m Γ(l + 1)Γ(l + α + β + 1) Γ(m + 1)Γ(l + α − m + 1) × Γ(l + α + 1)Γ(l + β + 1) Γ(l − m + 1)Γ(β + m + 1) We then obtain the following recurrence relations for C l,m : C l,0 C l−1,0 = − (l + β)(l + α + β) l(l + α) × α/2 + a + l β/2 + b + α/2 + a + l + 1 (for l ≥ 1) , The longitudinal integral L l (a, b; α, β) can then be calculated by first generating and then summing the following sequences: Let us first expand the quark current operator in terms of creation and annihilation operators. In agreement with the light-front quantization condition, the Dirac field operator at a given light-front time x + = 0 is expanded according to where the flavor indices are implicit. Here u s (p) and v s (p) are solutions of the Dirac equation for free fermions. Meanwhile, the creation and annihilation operators satisfy these anti-commutation relations: while other anti-commutation relations all vanish. We have defined the integral measure in the momentum space as dp = The reduced delta-function is defined according to These conventions ensure that one reduced deltafunction can be utilized to eliminate one momentumspace integration. With these definitions, the charge density operator becomes lim x→0 e f ψ(x)γ + ψ(x) = lim x→0 s s dp dp Here we have made use of u s (p )γ + u s (p) = 2 p + p + δ s s and v s (p )γ + v s (p) = 2 p + p + δ s s . We have only kept terms of relevance to the valence Fock sector of mesons. The form factors then becomes σ dp dp The anti-commutation relation for the creation and annihilation operator can be used to deduce σ dp dp d s ( σ dp dp d s ( The expression for the form factors is then reduced to where we have defined and Meanwhile, the following reductions of deltafunctions hold in the Drell-Yan frame: Therefore the expression for the form factors is reduced to that given on the right-hand side of Eq. (25) with q ⊥ = P ⊥ − P ⊥ .

The electromagnetic form factors in the basis representation
In the basis representation, Eq. (25) becomes = n m l nml rs ψ * n m l rs ψ nmlrs where we have applied shifts in the transverse momentum and φ * nm = φ n,−m . We then apply the Talmi-Moshinsky (TM) transform to simplify the integrals in the transverse momentum [65,76]. Specifically, the we have which corresponds to for the quark contribution and for the anti-quark contribution. The coefficient C(n , −m , n, m; N, M, n, m) can be computed with established procedures [65,76]. The following observations made specifically for the valence Fock sector of mesons will be helpful in enumerating terms after in the TM transform.
• Because the TM transform cannot change the total magnetic projection of the orbital angular momentum, we must have −m + m = M + m.
• The integral in k ⊥ will select the terms with m = 0, leaving other values of m not contributing to the integral.
• Because the mesons obtained from the lightfront Hamiltonian have fixed magnetic projection for the sum of the spin and orbital angular momenta, when the spins of the two wavefunction in a bilinear are identical, so are their magnetic quantum numbers m and m.
These observations further confine us to m = −m + m = 0, which is expected since the electromagnetic form factors have no angular dependence. The integrals over the momenta of the light-front wave functions then becomẽ C(n , m , l ; n, m, l; Q 2 ) With the aid of the TM transform, the electromagnetic form factors in the basis representation be- Below we provide the multi-qubit expressions for observables discussed in Sec. III for the J z = 0 sector of the BLFQ pion Hamiltonian.
In order to reduce the gate count, rather than using the circuit from and calculate the decay constant using v|ψ( θ) = ψ( θ) |v v| ψ( θ) . The mass radius matrix is given by eq. (22), which when expressed in terms of qubit operators is: The ρ l=0,l =0 density matrix of the parton distribution function is given by eqs. We calculate the elastic form factor matrix F P (Q 2 ), eq. (27), by discretizing Q 2 on the interval 0 ≤ Q 2 ≤ 5152900, evaluating the matrix C(n , m , l ; n, m, l; Q 2 ) for each value of Q 2 , and expanding it in terms of Pauli operators. For the sake of brevity, we do not include these explicit expressions for each point.
Appendix E: Bravyi-Kitaev encoding Both Jordan-Wigner and Bravyi-Kitaev encoding allow one to store the second-quantized fermionic states in a quantum computer. Within the Jordan-Wigner encoding, each qubit stores the occupancy of a particular orbital [56]. Within the Bravyi-Kitaev encoding, the information about parity is distributed equally between the operators and states [57]. In practice, one typically uses the more efficient BK encoding. While a single fermionic orbital in BK encoding is represented by up to O(log N ) qubits (instead of O(N ) in JW), the creation and annihilation operators are represented now by log N -local multi-qubit operators [57].
The BK-encoded basis states | . . . b 2 b 1 b 0 can be obtained from the JW-encoded states by means of the linear transformation b i = ij P ij f j [59], where the entries of P ij are {0, 1}, and multiplicaation modulo 2 is implied. Such a transformation of states can be implemented efficiently on a quantum computer. Since the matrix P ij is lower-triangular [59], multiplication modulo 2 can be performed on qubits using the CNOT gates, starting from the bottom row. For example, in the case of four qubits, the encoded matrix has the form of

(E1)
The corresponding circuit is shown on Fig. 8.
In order to perform the simulation in the BK encoding, one adjusts the procedure outlined in Sec. IV A as follows: a) After preparing the JWencoded initial state |ψ 0 , one appends to the circuit a block converting JW-encoded states to BKencoded ones as on Fig. 8; b) Eq. (34) is replaced with its BK version which changes the coefficients α j in (35) and h i in eq. (37).