Light mesons with one dynamical gluon on the light front

We obtain the light meson mass spectroscopy from the light-front quantum chromodynamics (QCD) Hamiltonian, determined for their constituent quark-antiquark and quark-antiquark-gluon Fock components, together with a three-dimensional confinement. The eigenvectors of the light-front effective Hamiltonian provide a good quality description of the pion electromagnetic form factor, decay constant, and the valence quark distribution functions following QCD scale evolution. We also show that the pion's gluon densities can be probed through the pion-nucleus induced $J/\psi$ production data. Our pion parton distribution functions provide excellent agreement with $J/\psi$ production data from widely different experimental conditions.

With the framework of BLFQ [26], we consider an effective LF Hamiltonian and solve for its mass eigenvalues and eigenstates at the scales suitable for low-resolution probes. With quarks and gluons being the explicit degrees of freedom for the strong interaction, our Hamiltonian incorporates LF QCD interactions [13] relevant to constituent quark-antiquark and quark-antiquark-gluon Fock components of the mesons with a complementary three-dimensional (3D) confinement [30]. By solving this Hamiltonian in the leading two Fock sectors and fitting the constituent parton masses and coupling constants as the model parameters, we obtain a good quality description of light meson mass spectroscopy. We evaluate the pion electromagnetic form factor (EMFF) and the parton distribution functions (PDFs) from our resulting light-front wave functions (LFWFs) obtained as eigenvectors of this Hamiltonian. The former characterizes the spatial extent of a hadron, while the latter describe the longitudinal momentum distribution of partons within the hadron. Since we interpret our model as appropriate to a low energy scale, we apply QCD evolution of the PDFs to higher momentum scales, and compare results with experiment.
One salient issue can be addressed with our approach-the important role of an initial gluon distribution at a low energy scale that contributes to all (sea and gluon) the parton distributions under QCD scale evolution. Experimentally, the pion PDFs are measured with the pion-nucleusinduced Drell-Yan process [40], where the differential cross section is sensitive to the valence quark densities at the parton's light-front momentum fraction x > 0.2. Sea and gluon distributions contribute over a wider range of x. Meanwhile, pion-induced charmonium production is dominated by the quark-antiquark and gluongluon fusion partonic processes so that the available J/ψ production data [41][42][43][44] are sensitive to both the quark and gluon distributions of the incident pion [45]. Within the framework of the color evaporation model (CEM) [46][47][48], we study the sensitivity of J/ψ production to our pion PDFs and confirm that consistency with the data depends sensitively on the shape and magnitude of the pion's gluon density. Thus, the pion's initial gluon distribution at the valence quark regime can contribute significantly to interpreting the J/ψ production data.

Mass spectra and LFWFs
The bound state problem on the LF is cast into an eigenvalue problem of the Hamiltonian: P − P + |Ψ = M 2 |Ψ , where P ± = P 0 ± P 3 employs the LF Hamitonian (P − ) and the longitudinal momentum (P + ) of the system, respectively, with the mass squared eigenvalue M 2 . At fixed LF time (x + = t + z), the meson state can be expressed in terms of various quark (q), antiquark (q) and gluon (g) Fock components, where the LFWFs ψ (... ) correspond to the probability amplitudes to find different parton configurations in the meson. At the initial scale where the mesons are described by |qq and |qqg , we adopt the LF Hamiltonian P − = P − QCD + P − C , where P − QCD and P − C are the LF QCD Hamiltonian and a model for the confining interaction. With one dynamical gluon in LF gauge [13] where ψ and A µ are the quark and gluon fields, respectively. T represents eight adjoint matrices of the SU (3) gauge group, and γ + = γ 0 + γ 3 , where γ µ are the Dirac matrices. The first two terms in Eq. (2) are the kinetic energies of the quark and gluon, where m 0 and m g are the masses of the bare quark and gluon. While the gluon mass is zero in QCD, we allow a phenomenological gluon mass to be fit to the spectra in our low-energy model. The last two terms describe their interactions with coupling g s . Following a renormalization procedure developed for positronium in a basis consisting of the |eē and |eēγ [49,50], we introduce a mass counter term, δm q = m 0 − m q that represents the quark mass correction due to the quantum fluctuations to the higher Fock sector, where m q is the renormalized quark mass.
Referring to Ref. [51], we allow an independent quark mass m f in the vertex interaction.
We adopt confinement in the leading Fock sector following [30] with κ being the strength of the confinement. The transverse confinement corresponds to the LF holographic potential, where the holographic variable is x(1 − x) r ⊥ [15]. The confining potential, Eq. (3), reproduces a symmetric 3D confinement in the nonrelativistic limit and it has been successfully applied to mesons [21,[29][30][31][32][33][34][35][36][37] and to the nucleon [38]. Confinement in the qqg sector relies solely on the cutoffs of the BLFQ basis functions introduced below.
Following BLFQ, for each Fock-particle we employ a plane-wave, e −ip + x − /2 , to describe its longitudinal motion and two dimensional harmonic oscillator ("2D-HO") wave function, Φ nm ( p ⊥ ; b) with scale parameter b, to describe its transverse degrees of freedom [27]. The longitudinal motion is confined to a box of length 2L with antiperiodic (periodic) boundary conditions for fermions (bosons). Therefore, the longitudinal momentum p + = 2πk/L, where k = 1 2 , 3 2 , 5 2 , ... for fermions and k = 1, 2, 3, ... for bosons. We neglect the zero mode for bosons. For all manybody basis states, we rescale the total longitudinal momentum P + = i p + i using K = i k i such that P + = 2π L K. For given parton i, the longitudinal momentum fraction x is then defined as x i = p + i /P + = k i /K. The 2D-HO wave function carries the radial quantum number n and angular quantum number m. With the additional quantum number λ for the helic-  [20], and the BSE (magenta-square) [53].
ity, each single-parton basis state is identified usingᾱ = {x, n, m, λ}. In the case of Fock sectors allowing for multiple color-singlet state (i.e. beyond Fock sectors retained here) we need an additional label to distinguish each color singlet state. Additionally, our many-body basis states have well defined total angular momentum projection We truncate the infinite basis space by introducing limits K and N max , such that i (2n i + |m i | + 1) ≤ N max in longitudinal and transverse directions, respectively. The N max truncation enables factorization of transverse center of mass motion [28] and manifests a natural ultra-violet (UV) regulator The LFWFs of the mesons in momentum space are then expressed as where ψ N =2 ({α i }) and ψ N =3 ({α i }) are the components of the eigenvectors relevant to the Fock sectors |qq and |qqg , respectively, in the BLFQ basis obtained from diagonalizing the full Hamiltonian matrix. Meanwhile, φ nm ( p ⊥ , b) represents the 2D-HO basis functions.
We select {N max , K} = {14, 15}, the HO scale parameter b = 0.29 GeV and set our model parameters {m q , m g , κ, m f , g s } = {0.39 GeV, 0.60 GeV, 0.65 GeV, 5.69 GeV, 1.92} to fit the known masses of π, ρ, a 0 , a 1 , π , and π 1 . Figure 1 shows the mass spectra for unflavored light mesons and compares with experimental data compiled by the particle data group (PDG) [52] and with the predictions of LF holographic QCD (LFHQCD) [15]. We also include results based on light-front holography [20]  and the Bethe Salpeter equation (BSE) [53] for comparison. For states not involved in the fit, b 1 , a 2 , and ρ , our results deviate more significantly from the experimental data. However, we especially note that we are able to fit the hybrid state π 1 , which may be viewed as a qq meson with a vibrating gluon flux tube [52].

Pion EMFF and decay constant
We first employ our resulting LFWFs for the pion EMFF F (Q 2 ) via where p = p + q, Q 2 = −q 2 and the electromagnetic current J µ EM (z) = f e fq (z)γ µ q(z) with f =d, u and the quark electric charge e f . Taking µ = +, the EMFF can be expressed in terms of the meson LFWFs using the Drell-Yan-West formula [62], δ 2 ( p ⊥j ) and for a struck parton, Our prediction for the EMFF of the charged pion is compared to the experimental data in Fig. 2. Note that our choice of N max = 14, implies the UV regulator Λ UV ∼ b √ N max ≈ 1 GeV. Our agreement with the precise low Q 2 EMFF data is consistent with our expectation that our predictions are most reliable in the low Q 2 regime.
The decay constant f P of a pseudoscalar meson is defined as the local vacuum-to-hadron matrix element: We obtain f Th. π = 138 MeV for the pion from the |qq Fock sector alone (the |qqg sector does not contribute) close to the experimental f Exp. π = 130.0 ± 0.2 MeV [52].
Pion PDFs.-The PDF is the probability of finding a collinear parton carrying momentum fraction x. With our LFWFs, the valence quark (antiquark) and the gluon PDFs in the pion are given by where i = q,q, g labels the valence quark, valence antiquark, and gluon, respectively. At our model scale the PDFs for the valence quark (antiquark) are normalized as 1 0 f q/q (x)dx = 1, and those PDFs together with the gluon PDF satisfy the momentum sum rule: We solve the next-to-next-to-leading order (NNLO) Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations [63][64][65] of QCD numerically using the higher order perturbative parton evolution toolkit [66] to evolve our PDFs from our model scale (µ 2 0 ) to a higher scale (µ 2 ). We determine µ 2 0 by requiring the result after evolution to produce the total first moments of the valence quark and the valence antiquark distributions from the global QCD analysis, x valence = 0.48 ± 0.01 at µ 2 = 5 GeV 2 [67]. This results in µ 2 0 = 0.34 ± 0.03 GeV 2 and we then evolve our initial PDFs to the relevant experimental scale µ 2 = 16 GeV 2 . While employing the DGLAP equations, we impose the condition that the running coupling α s (µ 2 ) saturates in the infrared at a cutoff value of max{α s } ∼ 1 [24,32,33,38,39]. Note that the sea quark distributions are absent in the initial scales of our model. Figure 3 shows our results for the pion PDFs, where we compare the valence quark distribution after QCD evolution with the data from the E615 experiment [40] as well as the reanalysis of the E615 experiment [69,70]. We also include the pion PDFs previously obtained in BLFQ-NJL model [31][32][33] based on a valence Fock representation for comparison. The error bands in our evolved PDFs are manifested from an adopted 10% uncertainty in our initial scale. We find a good agreement between our prediction for the pion valence quark distribution and The red lines correspond to BLFQ-NJL predictions [32]. Results are compared with the original analysis of the FNAL-E615 experiment [40] data and with its reanalysis (E615 Mod-data) [69]. Lower panel: Our result for the pion gluon PDF at µ 2 = 4 GeV 2 is compared with the global fits, JAM [67] and xFitter [68].
The gluon distribution significantly increases in our approach compared to that in the BLFQ-NJL model as well as to the global fits [67,68] as can be seen from Fig. 3. The BLFQ-NJL model is based on the pion valence Fock sector and gluons are generated solely from the scale evolution. However, our model includes a dynamical gluon at the initial scale and results in a larger gluon PDF at large-x (> 0.2) after scale evolution. Our result is also supported by a recent study from DSEs [12]. We find that gluons carry {39.5, 42.1, 43.9, 44.6, 45.1}% of pion momentum at {1.69, 4, 10, 16, 27} GeV 2 , respectively.

Pion-nucleus induced J/ψ production
Finally, we perform the next-to-leading order (NLO) calculations of the differential cross sec-tions for charmonium production by a pion beam to compare with available J/ψ production data with a special focus on the role of our model's gluon distribution. Among several theoretical approaches [46][47][48][75][76][77][78], we adopt the CEM [46][47][48], assuming a constant probability for cc pairs to hadronize into a given charmonium. The CEM has only one effective parameter (F ) that accounts for the probability and provides a good description of many features of fixed-target J/ψ cross section data with proton beams [79,80] and the collider data at the Relativistic Heavy Ion Collider and Large Hadron Collider [81,82].
The differential cross section for J/ψ production from the pion-nucleus collision in the CEM is given by [45,[83][84][85] The variables S and s denote the square of center-of-mass energy of the colliding πN and the interacting partons, respectively, and m c , m D , and M cc are the masses of the charm quark, D meson, and cc pair, respectively. Theσ ij is the short-distance differential cross section of heavy quark pair production at NLO, calculable as a perturbation series in the strong coupling [83]. The f π ± and f N are the PDFs of the pion and the target nuclei, respectively, evaluated at the factorization scale µ F . To complement our model for the pion PDFs, we adopt the nuclear PDFs from nCTEQ 2015 [86] and consider the factorization scale µ F = 2m c and the renormalization scale µ R = m c with m c = 1.5 GeV to perform the calculation for the cross section [85]. Figure 4 illustrates the cross section dσ/dx F for the J/ψ production as a function of x F and compares with various experimental data, including the E672 and E706 experiments with 515 GeV pions and beryllium target [41], the E705 experiment with 300 GeV pions and lithium target [42], the NA3 experiment with 200 GeV pions and hydrogen target [43], and the WA11 experiment with 190 GeV pions and beryllium target [44]. We determine the hadronization factor F , as an overall normalization parameter, by the best χ 2 fit to those data, shown as the black lines in Fig. 4. The values of F and the corresponding χ 2 /ndf values of all the best fits are also displayed in the plot ("ndf" represents "number of degrees-of-  Figure 4: dσ/dx F for the π − nucleus → J/ψX process as a function of x F . The data are taken from FNAL E672, E706, E705, CERN NA003 and CERN WA11 experiments [41][42][43][44]. The negative contribution from qg is presented here by multiplying with "(−1)".
freedom"). We obtain good χ 2 /ndf values for all four data sets. We also present the individual qq, gg and qg contributions to the total cross sections in the NLO calculation. Notably, the gg contribution dominates the cross section up to x F ∼ 0.5 and it decreases dramatically toward larger x F . The contribution from gg is increasing as the energy increases. Meanwhile, the qq contribution exhibits slower falloff and the relative importance of qq rises at the large-x F region. This is because the valence quark (antiquark) dominates the PDFs at large-x. The contribution from qg is negative (as signified by the "(−1)" in the labels) and relatively small. Overall, we observe that the J/ψ production data are sensitive to the pion PDFs, especially the large-x (x > 0.2) gluon distribution of pions. Our pion PDFs provide good agreement with data from widely different experimental conditions. Our finding also supports the study reported in Ref. [45].

Conclusion and outlook
We have solved the light-front QCD Hamiltonian for the light mesons by considering them within the constituent quark-antiquark and the quark-antiquark-gluon Fock spaces. Together with a three-dimensional confinement in the leading Fock sector, the eigenvalues of the Hamiltonian in Basis Light Front Quantization provide a good quality description of the light mesons' mass spectra. The LFWFs obtained as the eigenvectors of this Hamiltonian were then employed to produce the pion EMFF and the initial PDFs. We have obtained excellent agreement with the experimental data in the low-Q 2 regime for the pion EMFF. The PDFs at a higher experimental scale have been computed based on the NNLO DGLAP equations and we obtain reasonable agreement with the experimental data for the valence quark distribution.
We have also calculated the differential cross sections for pion-induced J/ψ production using the CEM framework at NLO and have obtained good agreement with available data [41][42][43][44]. Our study indicates that a proper description of J/ψ production data imposes a strong constraint on the pion's PDFs. We confirm the sensitivity of the J/ψ production data to the pion's gluon density in the valence-quark regime. The resulting LFWFs can be employed to study other quark and gluon distributions, such as the generalized parton distributions, the transverse momentum dependent parton distributions as well as the double parton correlations etc., in the light mesons. On the other hand, this work can be extended to higher Fock sectors to incorporate, for example, sea degrees of freedom as well.