# Light mesons with one dynamical gluon on the light front

Jiangshan Lan<EMAIL_ADDRESS>Kaiyu Fu<EMAIL_ADDRESS>Chandan
Mondal<EMAIL_ADDRESS>Xingbo Zhao<EMAIL_ADDRESS>James P. Vary
<EMAIL_ADDRESS>
(BLFQ Collaboration)  Lanzhou University, Lanzhou 730000, China  Institute for
Modern Physics, Chinese Academy of Sciences, Lanzhou-730000, China  School of
Nuclear Science and Technology, University of Chinese Academy of Sciences,
Beijing 100049, China  CAS Key Laboratory of High Precision Nuclear
Spectroscopy, Institute of Modern Physics, Chinese Academy of Sciences,
Lanzhou 730000, China  Department of Physics and Astronomy, Iowa State
University, Ames, Iowa 50011, USA

###### Abstract

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.

###### keywords:

Light mesons , Dynamical gluons , Form factors , Parton distribution functions
, Light-front quantization

††journal: Physics Letters B

##  1 Introduction

Quantum chromodynamics (QCD) is the well established theory for the strong
interactions [1], where hadrons are viewed as confined systems of partons
(quarks and gluons). However, understanding color confinement and chiral
symmetry breaking remains incomplete and it is not yet possible to forecast,
from QCD first principles, the extensive experimentally measured hadron
spectroscopy. Successful theoretical frameworks for predicting some aspects of
hadron spectra and illuminating partonic structures are discretized space-time
euclidean lattice [2, 3, 4, 5, 6, 7] and the Dyson–Schwinger equations (DSEs)
of QCD [8, 9, 11, 10, 12]. Much progress is being made within the Hamiltonian
formulation of QCD quantized on the light front (LF) [13, 14]. Complementary
insights into nonperturbative QCD can be achieved from LF holography [15, 16,
17, 18, 19, 20, 21, 22, 23, 25, 24]. Meanwhile, basis light-front quantization
(BLFQ) provides a computational framework for solving relativistic many-body
bound state problems in quantum field theories [26, 29, 27, 28, 30, 31, 32,
33, 34, 35, 38, 39, 36, 37].

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-nucleus-induced 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 gluon-gluon
fusion partonic processes so that the available $J/\psi$ 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/\psi$ 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/\psi$ production data.

##  2 Mass spectra and LFWFs

The bound state problem on the LF is cast into an eigenvalue problem of the
Hamiltonian: $P^{-}P^{+}|{\Psi}\rangle=M^{2}|{\Psi}\rangle$, where
$P^{\pm}=P^{0}\pm 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 ($\bar{q}$) and
gluon $(g)$ Fock components,

$\displaystyle|\Psi\rangle=\psi_{(q\bar{q})}|q\bar{q}\rangle+\psi_{(q\bar{q}g)}|q\bar{q}g\rangle+\dots\,,$
(1)

where the LFWFs $\psi_{(\dots)}$ correspond to the probability amplitudes to
find different parton configurations in the meson. At the initial scale where
the mesons are described by $|q\bar{q}\rangle$ and $|q\bar{q}g\rangle$, we
adopt the LF Hamiltonian $P^{-}=P^{-}_{\rm QCD}+P^{-}_{C}$, where $P^{-}_{\rm
QCD}$ and $P^{-}_{C}$ are the LF QCD Hamiltonian and a model for the confining
interaction. With one dynamical gluon in LF gauge [13]

$\displaystyle P_{\rm QCD}^{-}=\int
d^{2}x^{\perp}dx^{-}\Big{\\{}\frac{1}{2}\bar{\psi}\gamma^{+}\frac{m_{0}^{2}+(i\partial^{\perp})^{2}}{i\partial^{+}}\psi$
$\displaystyle-\frac{1}{2}A_{a}^{i}\left[m_{g}^{2}+(i\partial^{\perp})^{2}\right]A^{i}_{a}+g_{s}\bar{\psi}\gamma_{\mu}T^{a}A_{a}^{\mu}\psi$
$\displaystyle+\frac{1}{2}g_{s}^{2}\bar{\psi}\gamma^{+}T^{a}\psi\frac{1}{(i\partial^{+})^{2}}\bar{\psi}\gamma^{+}T^{a}\psi\Big{\\}},$
(2)

where $\psi$ and $A^{\mu}$ are the quark and gluon fields, respectively. $T$
represents eight adjoint matrices of the $SU(3)$ gauge group, and
$\gamma^{+}=\gamma^{0}+\gamma^{3}$, where $\gamma^{\mu}$ 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\bar{e}\rangle$ and $|e\bar{e}\gamma\rangle$ [49, 50], we introduce a mass
counter term, $\delta 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]

$\displaystyle P_{\rm C}^{-}P^{+}$
$\displaystyle=\kappa^{4}\left\\{x(1-x)\vec{r}_{\perp}^{2}-\frac{\partial_{x}[x(1-x)\partial_{x}]}{(m_{q}+m_{\bar{q}})^{2}}\right\\}$
(3)

with $\kappa$ being the strength of the confinement. The transverse
confinement corresponds to the LF holographic potential, where the holographic
variable is $\sqrt{x(1-x)}\vec{r}_{\perp}$ [15]. The confining potential, Eq.
(3), reproduces a symmetric 3D confinement in the nonrelativistic limit and it
has been successfully applied to mesons [29, 21, 30, 31, 32, 33, 34, 35, 36,
37] and to the nucleon [38]. Confinement in the $q\bar{q}g$ 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^{-i{p^{+}x^{-}/2}}$, to describe its longitudinal motion and two
dimensional harmonic oscillator (“2D-HO”) wave function,
$\Phi_{nm}(\vec{p}_{\perp};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\pi k/L$,
where $k=\frac{1}{2},\frac{3}{2},\frac{5}{2},...$ for fermions and
$k=1,2,3,...$ for bosons. We neglect the zero mode for bosons. For all many-
body basis states, we rescale the total longitudinal momentum
$P^{+}=\sum_{i}p_{i}^{+}$ using $K=\sum_{i}k_{i}$ such that
$P^{+}=\frac{2\pi}{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 $\lambda$ for the helicity, each
single-parton basis state is identified using
$\bar{\alpha}=\\{x,n,m,\lambda\\}$. 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
$M_{J}=\sum_{i}\left(m_{i}+\lambda_{i}\right).$ We truncate the infinite basis
space by introducing limits $K$ and $N_{\rm max}$, such that
$\sum_{i}(2n_{i}+|m_{i}|+1)\leq N_{\rm max}$ in longitudinal and transverse
directions, respectively. The $N_{\rm max}$ truncation enables factorization
of transverse center of mass motion [28] and manifests a natural ultra-violet
(UV) regulator $\Lambda_{\rm UV}\sim b\sqrt{N_{\rm max}}$.

The LFWFs of the mesons in momentum space are then expressed as

$\displaystyle\Psi^{\mathcal{N},\,M_{J}}_{\\{x_{i},\vec{p}_{\perp
i},\lambda_{i}\\}}$
$\displaystyle=\sum_{\\{n_{i}m_{i}\\}}\psi^{\mathcal{N}}({\\{\overline{\alpha}_{i}}\\})\prod_{i=1}^{\mathcal{N}}\phi_{n_{i}m_{i}}(\vec{p}_{\perp
i},b)\,,$ (4)

where $\psi^{\mathcal{N}=2}(\\{\overline{\alpha}_{i}\\})$ and
$\psi^{\mathcal{N}=3}(\\{\overline{\alpha}_{i}\\})$ are the components of the
eigenvectors relevant to the Fock sectors $|q\bar{q}\rangle$ and
$|q\bar{q}g\rangle$, respectively, in the BLFQ basis obtained from
diagonalizing the full Hamiltonian matrix. Meanwhile,
$\phi_{nm}(\vec{p}_{\perp},b)$ represents the 2D-HO basis functions.

Figure 1: The mass spectra of unflavored light mesons. Our results (Black-
dashed bars) are compared with the experimental data taken from PDG (red-solid
bars) [52], the LFHQCD (blue-circle) [15], Gutsche et. al. (cyan-triangle)
[20], and the BSE (magenta-square) [53].

We select $\\{N_{\rm max},K\\}\,$=$\,\\{14,15\\}$, the HO scale parameter
$b\,$=$\,0.29$ GeV and set our model parameters
$\\{m_{q},\,m_{g},\,\kappa,\,m_{f},\,g_{s}\\}\,=\,\\{{0.39\,\rm
GeV},\,{0.60\,\rm GeV},\,{0.65\,\rm GeV}$, $\,{5.69\,\rm GeV},\,1.92\\}$ to
fit the known masses of $\pi$, $\rho$, $a_{0}$, $a_{1}$, $\pi^{\prime}$, and
$\pi_{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 $\rho^{\prime}$, our results deviate more significantly from the
experimental data. However, we especially note that we are able to fit the
hybrid state $\pi_{1}$, which may be viewed as a $q\bar{q}$ meson with a
vibrating gluon flux tube [52].

##  3 Pion EMFF and decay constant

We first employ our resulting LFWFs for the pion EMFF $F(Q^{2})$ via

$\displaystyle\langle\Psi(p^{\prime})|J_{\text{EM}}^{\mu}(0)|\Psi(p)\rangle=(p+p^{\prime})^{\mu}F(Q^{2})\,,$
(5)

where $p^{\prime}=p+q$, $Q^{2}=-q^{2}$ and the electromagnetic current
$J_{\mathrm{EM}}^{\mu}(z)=\sum_{f}e_{f}\,\bar{q}(z)\gamma^{\mu}q(z)$ with
$f=\bar{d},u$ and the quark electric charge $e_{f}$. Taking $\mu=+$, the EMFF
can be expressed in terms of the meson LFWFs using the Drell-Yan-West formula
[62],

$\displaystyle F(Q^{2})=$
$\displaystyle\sum_{f}e_{f}\sum_{\mathcal{N},\,\lambda_{i}}$
$\displaystyle\int_{\mathcal{N}}\Psi^{\mathcal{N},\,M_{J}=0\,*}_{\\{x_{i},\vec{p}^{\,\prime}_{\perp
i},\lambda_{i}\\}}\,\Psi^{\mathcal{N},\,M_{J}=0}_{\\{x_{i},\vec{p}_{\perp
i},\lambda_{i}\\}}\,,$ (6)

where
$\int_{\mathcal{N}}\equiv\prod_{i=1}^{\mathcal{N}}\int\left[\frac{dx\,d^{2}\vec{p}_{\perp}}{16\pi^{3}}\right]_{i}16\pi^{3}\delta(1-\sum
x_{j})$ $\delta^{2}(\sum\vec{p}_{\perp j})$ and for a struck parton,
$\vec{p}^{\,\prime}_{\perp i}=\vec{p}_{\perp i}+(1-x_{i})\vec{q}_{\perp}$
while $\vec{p}^{\,\prime}_{\perp i}=\vec{p}_{\perp i}-x_{i}\vec{q}_{\perp}$
for the spectators. Considering the frame, where $q^{+}=0$,
$Q^{2}=-q^{2}=\vec{q}^{\,2}_{\perp}$. The electric charge
$e_{u}(e_{\bar{d}})=\frac{2}{3}(\frac{1}{3})$, while $e_{g}=0$.

Figure 2: The EMFF of the pion. The data are taken from Refs. [54, 55, 56, 57,
58, 59, 60, 61].

Our prediction for the EMFF of the charged pion is compared to the
experimental data in Fig. 2. Note that our choice of $N_{\rm max}=14$, implies
the UV regulator $\Lambda_{\rm UV}\sim b\sqrt{N_{\rm max}}\approx 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_{\rm P}$ of a pseudoscalar meson is defined as the local
vacuum-to-hadron matrix element:

$\displaystyle\langle
0|\,\bar{q}(0)\gamma^{\mu}\gamma_{5}q(0)\,|\Psi(p)\rangle=ip^{\mu}f_{\rm
P}\,.$ (7)

We obtain $f_{\rm\pi}^{\rm Th.}=138$ MeV for the pion from the
$|q\bar{q}\rangle$ Fock sector alone (the $|q\bar{q}g\rangle$ sector does not
contribute) close to the experimental $f_{\rm\pi}^{\rm Exp.}=130.0\pm 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

$\displaystyle f_{i}(x)=\sum_{\mathcal{N},\,\lambda_{i}}$
$\displaystyle\int_{\mathcal{N}}\Psi^{\mathcal{N},\,M_{J}=0\,*}_{\\{x_{i},\vec{p}_{\perp
i},\lambda_{i}\\}}\,\Psi^{\mathcal{N},\,M_{J}=0}_{\\{x_{i},\vec{p}_{\perp
i},\lambda_{i}\\}}\,\delta(x-x_{i})\,,$ (8)

where $i=q,\bar{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 $\int_{0}^{1}f_{q/\bar{q}}(x)dx=1$, and those PDFs together
with the gluon PDF satisfy the momentum sum rule:
$\int_{0}^{1}\sum_{i}xf_{i}(x)dx=1$.

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 ($\mu_{0}^{2}$) to a higher scale ($\mu^{2}$). We
determine $\mu_{0}^{2}$ 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, $\langle x\rangle_{\rm
valence}=0.48\pm 0.01$ at $\mu^{2}=5$ GeV2 [67]. This results in
$\mu_{0}^{2}=0.34\pm 0.03$ GeV2 and we then evolve our initial PDFs to the
relevant experimental scale $\mu^{2}=16~{}{\rm GeV}^{2}$. While employing the
DGLAP equations, we impose the condition that the running coupling
$\alpha_{s}({\mu^{2}})$ saturates in the infrared at a cutoff value of
max$\\{\alpha_{s}\\}\sim 1$ [24, 32, 33, 38, 39]. Note that the sea quark
distributions are absent in the initial scales of our model.

Figure 3: The PDFs of the pion. Upper panel: the black lines are our results
evolved from the initial scale $(0.34\pm 0.03~{}\mathrm{GeV}^{2})$ using the
NNLO DGLAP equations to the experimental scale of $16~{}\mathrm{GeV}^{2}$. 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 $\mu^{2}=4$ GeV2 is compared with the global fits, JAM [67] and xFitter
[68].

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 reanalyzed E615 data, while the BLFQ-NJL model
prefers the original E615 data. In our current treatment, the pion valence PDF
falls off as $(1-x)^{1.77}$, favoring a slightly slower falloff compared to
the $(1-x)^{2}$ predicted by perturbative QCD [71] and calculations using DSEs
[72]. Note that Refs. [69, 70, 73] support our finding. On the other hand, the
BLFQ-NJL at same high-$x$ exhibits $(1-x)^{1.44}$ [32], supported by Ref.
[74].

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\\}$ GeV2, respectively.

##  4 Pion-nucleus induced $J/\psi$ production

Finally, we perform the next-to-leading order (NLO) calculations of the
differential cross sections for charmonium production by a pion beam to
compare with available $J/\psi$ 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 $c\bar{c}$ 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/\psi$ 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/\psi$ production from the pion-nucleus
collision in the CEM is given by [83, 84, 85, 45]

$\displaystyle\frac{d{\bf\sigma}}{dx_{\mathrm{F}}}\Big{|}_{J/\psi}$ (9)
$\displaystyle=F\sum_{i,j=q,\bar{q},g}\int^{2m_{D}}_{2m_{c}}dM_{c\bar{c}}\frac{2M_{c\bar{c}}}{S\sqrt{x_{F}^{2}+4M^{2}_{c\bar{c}}/S}}$
$\displaystyle\times\hat{\bf\sigma}_{ij}(s,m^{2}_{c},\mu_{F}^{2},\mu_{R}^{2})\,f^{\pi^{\pm}}_{i}(x_{1},\mu_{F}^{2})\,{f}^{N}_{j}(x_{2},\mu^{2}_{F})\,$

with $x_{1,2}=(\sqrt{x_{F}^{2}+4M^{2}_{c\bar{c}}/S}\pm x_{F})/2$ and
$x_{F}=x_{1}-x_{2}$. The variables $S$ and $s$ denote the square of center-of-
mass energy of the colliding $\pi N$ and the interacting partons,
respectively, and $m_{c}$, $m_{D}$, and $M_{c\bar{c}}$ are the masses of the
charm quark, $D$ meson, and $c\bar{c}$ pair, respectively. The
$\hat{\sigma}_{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^{\pi^{\pm}}$ and $f^{N}$ are the PDFs of the pion
and the target nuclei, respectively, evaluated at the factorization scale
$\mu_{F}$. To complement our model for the pion PDFs, we adopt the nuclear
PDFs from nCTEQ 2015 [86] and consider the factorization scale
$\mu_{F}=2m_{c}$ and the renormalization scale $\mu_{R}=m_{c}$ with
$m_{c}=1.5$ GeV to perform the calculation for the cross section [85].

Figure 4: $d{\bf\sigma}/dx_{F}$ for the
${\pi^{-}\rm{nucleus}}\rightarrow{J/\psi 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)$”.

Figure 4 illustrates the cross section $d\sigma/dx_{F}$ for the $J/\psi$
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
$\chi^{2}$ fit to those data, shown as the black lines in Fig. 4. The values
of $F$ and the corresponding $\chi^{2}/$ndf values of all the best fits are
also displayed in the plot (“ndf” represents “number of degrees-of-freedom”).
We obtain good $\chi^{2}/$ndf values for all four data sets. We also present
the individual $q\bar{q}$, $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}\sim 0.5$ and it decreases dramatically toward
larger $x_{F}$. The contribution from $gg$ is increasing as the energy
increases. Meanwhile, the $q\bar{q}$ contribution exhibits slower falloff and
the relative importance of $q\bar{q}$ 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/\psi$ 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].

##  5 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/\psi$ 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/\psi$ production data imposes a strong constraint on
the pion’s PDFs. We confirm the sensitivity of the $J/\psi$ 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.

## Acknowledgements

We thank J. Wu, H. Zhao, S. Xu, S. Jia and S. Platchkov for useful
discussions. C. M. thanks the Chinese Academy of Sciences President’s
International Fellowship Initiative for the support via Grants No. 2021PM0023.
C. M. is supported by new faculty start up funding by the Institute of Modern
Physics, Chinese Academy of Sciences, Grant No. E129952YR0. X. Z. is supported
by new faculty startup funding by the Institute of Modern Physics, Chinese
Academy of Sciences, by Key Research Program of Frontier Sciences, Chinese
Academy of Sciences, Grant No. ZDB-SLY-7020, by the Natural Science Foundation
of Gansu Province, China, Grant No. 20JR10RA067 and by the Strategic Priority
Research Program of the Chinese Academy of Sciences, Grant No. XDB34000000. J.
P. V. is supported by the Department of Energy under Grants No. DE-
FG02-87ER40371, and No. DE-SC0018223 (SciDAC4/NUCLEI). A portion of the
computational resources were also provided by Gansu Computing Center.

## References

  * [1] C. G. Callan, Jr., R. F. Dashen and D. J. Gross, “Toward a Theory of the Strong Interactions,” Phys. Rev. D 17, 2717 (1978).
  * [2] B. Joó et al. [USQCD], “Status and Future Perspectives for Lattice Gauge Theory Calculations to the Exascale and Beyond,” Eur. Phys. J. A 55, no.11, 199 (2019) [arXiv:1904.09725 [hep-lat]].
  * [3] A. Bazavov et al. [MILC], “Nonperturbative QCD Simulations with 2+1 Flavors of Improved Staggered Quarks,” Rev. Mod. Phys. 82, 1349-1417 (2010) [arXiv:0903.3598 [hep-lat]].
  * [4] Z. Fodor and C. Hoelbling, “Light Hadron Masses from Lattice QCD,” Rev. Mod. Phys. 84, 449 (2012) [arXiv:1203.4789 [hep-lat]].
  * [5] R. A. Briceno, J. J. Dudek and R. D. Young, “Scattering processes and resonances from lattice QCD,” Rev. Mod. Phys. 90, no.2, 025001 (2018) [arXiv:1706.06223 [hep-lat]].
  * [6] S. Durr, Z. Fodor, J. Frison, C. Hoelbling, R. Hoffmann, S. D. Katz, S. Krieg, T. Kurth, L. Lellouch and T. Lippert, et al. “Ab-Initio Determination of Light Hadron Masses,” Science 322, 1224-1227 (2008) [arXiv:0906.3599 [hep-lat]].
  * [7] P. Hagler, “Hadron structure from lattice quantum chromodynamics,” Phys. Rept. 490, 49-175 (2010) [arXiv:0912.5483 [hep-lat]].
  * [8] P. Maris and C. D. Roberts, “Dyson-Schwinger equations: A Tool for hadron physics,” Int. J. Mod. Phys. E 12, 297-365 (2003) [arXiv:nucl-th/0301049 [nucl-th]].
  * [9] C. D. Roberts and A. G. Williams, “Dyson-Schwinger equations and their application to hadronic physics,” Prog. Part. Nucl. Phys. 33, 477-575 (1994) [arXiv:hep-ph/9403224 [hep-ph]].
  * [10] A. Bashir, L. Chang, I. C. Cloet, B. El-Bennich, Y. X. Liu, C. D. Roberts and P. C. Tandy, “Collective perspective on advances in Dyson-Schwinger Equation QCD,” Commun. Theor. Phys. 58, 79-134 (2012) [arXiv:1201.3366 [nucl-th]].
  * [11] R. Alkofer and L. von Smekal, “The Infrared behavior of QCD Green’s functions: Confinement dynamical symmetry breaking, and hadrons as relativistic bound states,” Phys. Rept. 353, 281 (2001) [arXiv:hep-ph/0007355 [hep-ph]].
  * [12] A. Freese, I. C. Cloët and P. C. Tandy, “Gluon PDF from Quark dressing in the Nucleon and Pion,” [arXiv:2103.05839 [hep-ph]].
  * [13] S. J. Brodsky, H. C. Pauli and S. S. Pinsky, “Quantum chromodynamics and other field theories on the light cone,” Phys. Rept. 301, 299-486 (1998) [arXiv:hep-ph/9705477 [hep-ph]].
  * [14] B. L. G. Bakker, A. Bassetto, S. J. Brodsky, W. Broniowski, S. Dalley, et al. “Light-Front Quantum Chromodynamics: A framework for the analysis of hadron physics,” Nucl. Phys. B Proc. Suppl. 251-252, 165-174 (2014) [arXiv:1309.6333 [hep-ph]].
  * [15] S. J. Brodsky, G. F. de Teramond, H. G. Dosch and J. Erlich, “Light-Front Holographic QCD and Emerging Confinement,” Phys. Rept. 584, 1-105 (2015) [arXiv:1407.8131 [hep-ph]].
  * [16] S. J. Brodsky and G. F. de Teramond, “Hadronic spectra and light-front wavefunctions in holographic QCD,” Phys. Rev. Lett. 96, 201601 (2006) [arXiv:hep-ph/0602252 [hep-ph]].
  * [17] G. F. de Teramond and S. J. Brodsky, “Hadronic spectrum of a holographic dual of QCD,” Phys. Rev. Lett. 94, 201601 (2005) [arXiv:hep-th/0501022 [hep-th]].
  * [18] G. F. de Teramond and S. J. Brodsky, “Light-Front Holography: A First Approximation to QCD,” Phys. Rev. Lett. 102, 081601 (2009) [arXiv:0809.4899 [hep-ph]].
  * [19] T. Branz, T. Gutsche, V. E. Lyubovitskij, I. Schmidt and A. Vega, “Light and heavy mesons in a soft-wall holographic approach,” Phys. Rev. D 82, 074022 (2010) [arXiv:1008.0268 [hep-ph]].
  * [20] T. Gutsche, V. E. Lyubovitskij, I. Schmidt and A. Vega, “Chiral Symmetry Breaking and Meson Wave Functions in Soft-Wall AdS/QCD,” Phys. Rev. D 87, no.5, 056001 (2013) [arXiv:1212.5196 [hep-ph]].
  * [21] G. F. De Téramond and S. J. Brodsky, “Longitudinal dynamics and chiral symmetry breaking in holographic light-front QCD,” [arXiv:2103.10950 [hep-ph]].
  * [22] M. Ahmady, H. Dahiya, S. Kaur, C. Mondal, R. Sandapen and N. Sharma, “Extending light-front holographic QCD using the ’t Hooft Equation,” Phys. Lett. B 823, 136754 (2021) [arXiv:2105.01018 [hep-ph]].
  * [23] M. Ahmady, S. Kaur, S. L. MacKay, C. Mondal and R. Sandapen, “Hadron spectroscopy using the light-front holographic Schrödinger equation and the ’t Hooft equation,” Phys. Rev. D 104, no.7, 074013 (2021) [arXiv:2108.03482 [hep-ph]].
  * [24] J. Lan and C. Mondal, “Pion-nucleus induced Drell-Yan cross section in models inspired by light-front holography,” Phys. Lett. B 807, 135613 (2020) [arXiv:2007.05858 [hep-ph]].
  * [25] M. Ahmady, C. Mondal and R. Sandapen, “Dynamical spin effects in the holographic light-front wavefunctions of light pseudoscalar mesons,” Phys. Rev. D 98, no.3, 034010 (2018) [arXiv:1805.08911 [hep-ph]].
  * [26] J. P. Vary, H. Honkanen, J. Li, P. Maris, S. J. Brodsky, A. Harindranath, G. F. de Teramond, P. Sternberg, E. G. Ng and C. Yang, “Hamiltonian light-front field theory in a basis function approach,” Phys. Rev. C 81, 035205 (2010) [arXiv:0905.1411 [nucl-th]].
  * [27] X. Zhao, H. Honkanen, P. Maris, J. P. Vary and S. J. Brodsky, “Electron g-2 in Light-Front Quantization,” Phys. Lett. B 737, 65-69 (2014) [arXiv:1402.4195 [nucl-th]].
  * [28] P. Wiecki, Y. Li, X. Zhao, P. Maris and J. P. Vary, “Basis Light-Front Quantization Approach to Positronium,” Phys. Rev. D 91, no.10, 105009 (2015) [arXiv:1404.6234 [nucl-th]].
  * [29] Y. Li and J. P. Vary, “Light-front holography with chiral symmetry breaking,” [arXiv:2103.09993 [hep-ph]].
  * [30] Y. Li, P. Maris, X. Zhao and J. P. Vary, “Heavy Quarkonium in a Holographic Basis,” Phys. Lett. B 758, 118-124 (2016) [arXiv:1509.07212 [hep-ph]].
  * [31] S. Jia and J. P. Vary, “Basis light front quantization for the charged light mesons with color singlet Nambu–Jona-Lasinio interactions,” Phys. Rev. C 99, no.3, 035206 (2019) [arXiv:1811.08512 [nucl-th]].
  * [32] J. Lan, C. Mondal, S. Jia, X. Zhao and J. P. Vary, “Parton Distribution Functions from a Light Front Hamiltonian and QCD Evolution for Light Mesons,” Phys. Rev. Lett. 122, no.17, 172001 (2019) [arXiv:1901.11430 [nucl-th]].
  * [33] J. Lan, C. Mondal, S. Jia, X. Zhao and J. P. Vary, “Pion and kaon parton distribution functions from basis light front quantization and QCD evolution,” Phys. Rev. D 101, no.3, 034024 (2020) [arXiv:1907.01509 [nucl-th]].
  * [34] S. Tang, Y. Li, P. Maris and J. P. Vary, “$B_{c}$ mesons and their properties on the light front,” Phys. Rev. D 98, no.11, 114038 (2018) [arXiv:1810.05971 [nucl-th]].
  * [35] S. Tang, Y. Li, P. Maris and J. P. Vary, “Heavy-light mesons on the light front,” Eur. Phys. J. C 80, no.6, 522 (2020) [arXiv:1912.02088 [nucl-th]].
  * [36] J. Lan, C. Mondal, M. Li, Y. Li, S. Tang, X. Zhao and J. P. Vary, “Parton Distribution Functions of Heavy Mesons on the Light Front,” Phys. Rev. D 102, no.1, 014020 (2020) [arXiv:1911.11676 [nucl-th]].
  * [37] W. Qian, S. Jia, Y. Li and J. P. Vary, “Light mesons within the basis light-front quantization framework,” Phys. Rev. C 102, no.5, 055207 (2020) [arXiv:2005.13806 [nucl-th]].
  * [38] C. Mondal, S. Xu, J. Lan, X. Zhao, Y. Li, D. Chakrabarti and J. P. Vary, “Proton structure from a light-front Hamiltonian,” Phys. Rev. D 102, no.1, 016008 (2020) [arXiv:1911.10913 [hep-ph]].
  * [39] S. Xu, C. Mondal, J. Lan, X. Zhao, Y. Li, J. P. Vary [BLFQ], “Nucleon structure from basis light-front quantization,” Phys. Rev. D 104, no.9, 094036 (2021) [arXiv:2108.03909 [hep-ph]].
  * [40] J. S. Conway, C. E. Adolphsen, J. P. Alexander, K. J. Anderson, J. G. Heinrich, J. E. Pilcher, A. Possoz, E. I. Rosenberg, C. Biino and J. F. Greenhalgh, et al. “Experimental Study of Muon Pairs Produced by 252-GeV Pions on Tungsten,” Phys. Rev. D 39, 92-122 (1989).
  * [41] A. Gribushin et al. [E672 and E706], “Production of $J/\psi$ and psi (2S) mesons in $\pi$\- Be collisions at 515-GeV/c,” Phys. Rev. D 53, 4723-4733 (1996).
  * [42] L. Antoniazzi et al. [E705], “A Measurement of $J/\psi$ and psi-prime production in 300-GeV/c proton, anti-proton and pi+- nucleon interactions,” Phys. Rev. D 46, 4828-4835 (1992).
  * [43] J. Badier et al. [NA3], “Experimental J/psi Hadronic Production from 150-GeV/c to 280-GeV/c,” Z. Phys. C 20, 101 (1983).
  * [44] J. G. McEwen, B. Pietrzyk, R. Barate, P. Bareyre, P. Bonamy, P. Borgeaud, M. David, F. X. Gentit, G. Laurens and Y. Lemoigne, et al. “Measurement of the Gluon Structure Function of the Pion,” Phys. Lett. B 121, 198-202 (1983).
  * [45] W. C. Chang, J. C. Peng, S. Platchkov and T. Sawada, “Constraining gluon density of pions at large $x$ by pion-induced $J/\psi$ production,” Phys. Rev. D 102, no.5, 054024 (2020) [arXiv:2006.06947 [hep-ph]].
  * [46] M. B. Einhorn and S. D. Ellis, “Hadronic Production of the New Resonances: Probing Gluon Distributions,” Phys. Rev. D 12, 2007 (1975).
  * [47] H. Fritzsch, “Producing Heavy Quark Flavors in Hadronic Collisions: A Test of Quantum Chromodynamics,” Phys. Lett. B 67, 217-221 (1977).
  * [48] F. Halzen, “Cvc for Gluons and Hadroproduction of Quark Flavors,” Phys. Lett. B 69, 105-108 (1977).
  * [49] X. Zhao, “Advances in Basis Light-front Quantization,” Few Body Syst. 56, no.6-9, 257-265 (2015) [arXiv:1411.7748 [nucl-th]].
  * [50] X. Zhao, K. Fu, H. Zhao and J. P. Vary, “Positronium: an illustration of nonperturbative renormalization in a basis light-front approach,” PoS LC2019, 090 (2020) [arXiv:2103.06719 [hep-ph]].
  * [51] S. D. Glazek and R. J. Perry, “Special example of relativistic Hamiltonian field theory,” Phys. Rev. D 45, 3740-3754 (1992).
  * [52] M. Tanabashi et al. [Particle Data Group], “Review of Particle Physics,” Phys. Rev. D 98, no.3, 030001 (2018).
  * [53] C. S. Fischer, S. Kubrak and R. Williams, “Mass spectra and Regge trajectories of light mesons in the Bethe-Salpeter approach,” Eur. Phys. J. A 50, 126 (2014) [arXiv:1406.4370 [hep-ph]].
  * [54] S. R. Amendolia et al. [NA7 Collaboration], “A Measurement of the Space - Like Pion Electromagnetic Form-Factor,” Nucl. Phys. B 277, 168 (1986).
  * [55] C. J. Bebek et al., “Further measurements of forward-charged-pion electroproduction at large $\kappa^{2}$,” Phys. Rev. D 9, 1229 (1974).
  * [56] C. J. Bebek, C. N. Brown, M. Herzlinger, S. D. Holmes, C. A. Lichtenstein, F. M. Pipkin, S. Raither and L. K. Sisterson, “Measurement of the pion form-factor up to $q^{2}$ = 4-GeV2,” Phys. Rev. D 13, 25 (1976).
  * [57] C. J. Bebek et al., “Electroproduction of single pions at low epsilon and a measurement of the pion form-factor up to $q^{2}$ = 10-GeV2,” Phys. Rev. D 17, 1693 (1978).
  * [58] J. Volmer et al. [Jefferson Lab F(pi) Collaboration], “Measurement of the Charged Pion Electromagnetic Form-Factor,” Phys. Rev. Lett. 86, 1713 (2001).
  * [59] T. Horn et al. [Jefferson Lab F(pi)-2 Collaboration], “Determination of the Charged Pion Form Factor at Q**2 = 1.60 and 2.45-(GeV/c)**2,” Phys. Rev. Lett. 97, 192001 (2006).
  * [60] V. Tadevosyan et al. [Jefferson Lab F(pi)], “Determination of the pion charge form-factor for $Q^{2}=0.60$ GeV2 – $1.60$ GeV2,” Phys. Rev. C 75, 055205 (2007) [arXiv:nucl-ex/0607007 [nucl-ex]].
  * [61] G. M. Huber et al. [Jefferson Lab], “Charged pion form-factor between $Q^{2}=0.60$ GeV2 and $2.45$ GeV2. II. Determination of, and results for, the pion form-factor,” Phys. Rev. C 78, 045203 (2008) [arXiv:0809.3052 [nucl-ex]].
  * [62] S. J. Brodsky and G. F. de Teramond, “Light-Front Dynamics and AdS/QCD Correspondence: The Pion Form Factor in the Space- and Time-Like Regions,” Phys. Rev. D 77, 056007 (2008) [arXiv:0707.3859 [hep-ph]].
  * [63] Y. L. Dokshitzer, “Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics.,” Sov. Phys. JETP 46, 641 (1977) [Zh. Eksp. Teor. Fiz. 73, 1216 (1977)].
  * [64] V. N. Gribov and L. N. Lipatov, “Deep inelastic e p scattering in perturbation theory,” Sov. J. Nucl. Phys. 15, 438 (1972).
  * [65] G. Altarelli and G. Parisi, “Asymptotic Freedom in Parton Language,” Nucl. Phys. B 126, 298 (1977).
  * [66] G. P. Salam and J. Rojo, “A Higher Order Perturbative Parton Evolution Toolkit (HOPPET),” Comput. Phys. Commun. 180, 120 (2009).
  * [67] P. C. Barry, N. Sato, W. Melnitchouk and C. R. Ji, “First Monte Carlo Global QCD Analysis of Pion Parton Distributions,” Phys. Rev. Lett. 121, no.15, 152001 (2018) [arXiv:1804.01965 [hep-ph]].
  * [68] I. Novikov, H. Abdolmaleki, D. Britzger, A. Cooper-Sarkar, F. Giuli, A. Glazov, A. Kusina, A. Luszczak, F. Olness and P. Starovoitov, et al. “Parton Distribution Functions of the Charged Pion Within The xFitter Framework,” Phys. Rev. D 102, no.1, 014040 (2020) [arXiv:2002.02902 [hep-ph]].
  * [69] C. Chen, L. Chang, C. D. Roberts, S. Wan and H. S. Zong, “Valence-quark distribution functions in the kaon and pion,” Phys. Rev. D 93, no.7, 074021 (2016) [arXiv:1602.01502 [nucl-th]].
  * [70] M. Aicher, A. Schafer and W. Vogelsang, “Soft-gluon resummation and the valence parton distribution function of the pion,” Phys. Rev. Lett. 105, 252003 (2010) [arXiv:1009.2481 [hep-ph]].
  * [71] E. L. Berger and S. J. Brodsky, “Quark Structure Functions of Mesons and the Drell-Yan Process,” Phys. Rev. Lett. 42, 940-944 (1979).
  * [72] M. B. Hecht, C. D. Roberts and S. M. Schmidt, “Valence quark distributions in the pion,” Phys. Rev. C 63, 025213 (2001) [arXiv:nucl-th/0008049 [nucl-th]].
  * [73] K. D. Bednar, I. C. Cloët and P. C. Tandy, “Distinguishing Quarks and Gluons in Pion and Kaon Parton Distribution Functions,” Phys. Rev. Lett. 124, no.4, 042002 (2020) [arXiv:1811.12310 [nucl-th]].
  * [74] G. F. de Teramond et al. [HLFHS], “Universality of Generalized Parton Distributions in Light-Front Holographic QCD,” Phys. Rev. Lett. 120, no.18, 182001 (2018) [arXiv:1801.09154 [hep-ph]].
  * [75] C. H. Chang, “Hadronic Production of $J/\psi$ Associated With a Gluon,” Nucl. Phys. B 172, 425-434 (1980).
  * [76] E. L. Berger and D. L. Jones, “Inelastic Photoproduction of J/psi and Upsilon by Gluons,” Phys. Rev. D 23, 1521-1530 (1981).
  * [77] R. Baier and R. Ruckl, “Hadronic Collisions: A Quarkonium Factory,” Z. Phys. C 19, 251 (1983).
  * [78] G. T. Bodwin, E. Braaten and G. P. Lepage, “Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium,” Phys. Rev. D 51, 1125-1171 (1995) [erratum: Phys. Rev. D 55, 5853 (1997)]. [arXiv:hep-ph/9407339 [hep-ph]].
  * [79] R. Gavai, D. Kharzeev, H. Satz, G. A. Schuler, K. Sridhar and R. Vogt, “Quarkonium production in hadronic collisions,” Int. J. Mod. Phys. A 10, 3043-3070 (1995) [arXiv:hep-ph/9502270 [hep-ph]].
  * [80] G. A. Schuler and R. Vogt, “Systematics of quarkonium production,” Phys. Lett. B 387, 181-186 (1996) [arXiv:hep-ph/9606410 [hep-ph]].
  * [81] R. E. Nelson, R. Vogt and A. D. Frawley, “Narrowing the uncertainty on the total charm cross section and its effect on the J/$\psi$ cross section,” Phys. Rev. C 87, no.1, 014908 (2013) [arXiv:1210.4610 [hep-ph]].
  * [82] J. P. Lansberg, H. S. Shao, N. Yamanaka, Y. J. Zhang and C. Noûs, “Complete NLO QCD study of single- and double-quarkonium hadroproduction in the colour-evaporation model at the Tevatron and the LHC,” Phys. Lett. B 807, 135559 (2020) [arXiv:2004.14345 [hep-ph]].
  * [83] P. Nason, S. Dawson and R. K. Ellis, “The Total Cross-Section for the Production of Heavy Quarks in Hadronic Collisions,” Nucl. Phys. B 303, 607-633 (1988).
  * [84] P. Nason, S. Dawson and R. K. Ellis, “The One Particle Inclusive Differential Cross-Section for Heavy Quark Production in Hadronic Collisions,” Nucl. Phys. B 327, 49-92 (1989) [erratum: Nucl. Phys. B 335, 260-260 (1990)].
  * [85] M. L. Mangano, P. Nason and G. Ridolfi, “Fixed target hadroproduction of heavy quarks,” Nucl. Phys. B 405, 507-535 (1993).
  * [86] K. Kovarik, A. Kusina, T. Jezo, D. B. Clark, C. Keppel, F. Lyonnet, J. G. Morfin, F. I. Olness, J. F. Owens and I. Schienbein, et al. “nCTEQ15 - Global analysis of nuclear parton distributions with uncertainties in the CTEQ framework,” Phys. Rev. D 93, no.8, 085037 (2016) [arXiv:1509.00792 [hep-ph]].


# Number Entity Recognition

Dhanasekar Sundararaman1, Vivek Subramanian2,
Guoyin Wang2, Liyan Xu3, Lawrence Carin1
1 Duke University
2 Amazon Alexa AI
3 Emory University
<EMAIL_ADDRESS>

###### Abstract

Numbers are essential components of text, like any other word tokens, from
which natural language processing (NLP) models are built and deployed. Though
numbers are typically not accounted for distinctly in most NLP tasks, there is
still an underlying amount of numeracy already exhibited by NLP models. In
this work, we attempt to tap this potential of state-of-the-art NLP models and
transfer their ability to boost performance in related tasks. Our proposed
classification of numbers into entities helps NLP models perform well on
several tasks, including a handcrafted Fill-In-The-Blank (FITB) task and on
question answering using joint embeddings, outperforming the BERT and RoBERTa
baseline classification.

##  1 Introduction

Named entity recognition (NER) is the task of classifying nouns within a
document into categories to which they belong Sang and De Meulder (2003). Deep
learning approaches for NER have been successful recently because of their
ability to understand the underlying semantics of sentences and incorporate
this knowledge into tagging. Deep architectures for NER Lample et al. (2016)
include recurrent neural networks and attention based models including the
Transformer Vaswani et al. (2017).

While NER has been performed on a number of datasets, tagging types of
numbers, or number entities, has not been a principal focus of NER. This is
despite the fact that numbers are prevalent in most texts, comprising roughly
5% of tokens in a sentence Naik et al. (2019). Numbers can take on several
categories, including years, ages, phone numbers, dates, etc. Figure 1 shows
sample sentences in which numbers are classified into entity types.
Distinguishing these different types of numbers can greatly enhance
performance on downstream tasks.

Figure 1: Sequence classification of numbers, or number entity recognition.

In this paper, we present several methods for addressing number entity
recognition (NuER). We begin by demonstrating that state-of-the-art NLP models
are able to classify types of numbers that appear in text Wallace et al.
(2019); Zhang et al. (2020). We study this by utilizing BERT Devlin et al.
(2018) to perform NuER on a custom version of the Stanford Question Answering
Dataset (SQuAD) Rajpurkar et al. (2016), pruned to contain roughly 10K
sentences for which the answers consist of numbers, and annotated in-house for
numbers of six different entity types (year, percentage, date, count, age, and
size) Sundararaman et al. (2020). We hereby refer to this as the SQuAD-Num
dataset and the BERT classifier used to predict entities as the NuER model. We
discover that BERT model performs well at this task on most of the categories
listed above. The success of this model is made possible by the highly
contextualized embeddings (e.g., the appearance of the name of a month, say
March, next to a year, say 2003) inferred by BERT.

We then verify the effectiveness of BERT fine-tuned on SQuAD-Num to tag
numbers that appear in a out-of-domain, human-annotated, 2K-sentence subset of
the Numeracy 600K dataset Chen et al. (2019), which we call HA-Numeracy 2K. We
chose this dataset because of the abundance and diversity of numbers present.
The subset was sampled in such a way as to maintain the distribution of number
magnitudes that are originally present. We find that contextual embeddings are
again beneficial in this setting, as surrounding tokens in both the source and
target domains tend to be related. After verifying the quality of the
annotations on the small sample, we annotated the full Numeracy 600K dataset
using this model, which we hereafter call ENT-Numeracy 600K. Table 1 shows the
number of samples in each entity type in SQuAD-Num. We share these annotations
with the NLP community to further enhance research in this direction.

C | S | Y | P | D | A
---|---|---|---|---|---
1800 | 447 | 4355 | 291 | 418 | 82
Table 1: Number of entities in our annotation for each type in SQuAD-Num. The
table headers represent the following entity types respectively: Count, Size,
Year, Percentage, Date, Age.

Next, we jointly train embeddings for both tokens and number entities on the
SQuAD-Num dataset. Answers to questions involving numbers are predicted by
combining together token embeddings with trainable number entity embeddings.
We find that this addition enables BERT and RoBERTa models to better answer
questions in SQuAD involving numbers. We later analyze the effectiveness of
NuER through a handcrafted Fill-In-The-Blank (FITB) classification task. Our
proposed NuER model helps the models classify numbers better when presented
with the entity type information.

In summary, our contributions are as follows:

  1. 1.

We develop a BERT-based model for NuER which is used to tag numbers that
appear in the SQuAD-Num dataset and in the out-of-domain HA-Numeracy 2K.

  2. 2.

We release datasets containing number entity annotations for the NLP
community.

  3. 3.

We observe the effectiveness of NuER through a handcrafted Fill-In-The-Blank
(FITB) task to predict a masked number in a sentence.

  4. 4.

We then jointly train embeddings for tokens and number entities, and show that
numerical question answering performance improves when both are fed into the
model.

##  2 Related Work

Classification of numbers in plain text to entities has not been addressed
methodologically and deserves considerable attention. Naik et al. (2019)
highlighted the importance of developing specialized embeddings for numbers
such that basic properties including comparison and summation would be
preserved. Since then In addition, Sundararaman et al. (2020) designed a novel
loss function that captured these properties and demonstrated improved
performance on a variety of NLP tasks. Unlike their method which uses static
embeddings, we propose to supplement the entity information of numbers in a
contextual fashion to enable models to perform well on a number of tasks Geva
et al. (2020). There are existing methods on numerical word embeddings Berg-
Kirkpatrick and Spokoyny (2020); Spithourakis and Riedel (2018), but they
didn’t consider numbers as entities for the task of NER. NLP models have shown
to use supplemental knowledge to its advantage Tanwar et al. (2022);
Sundararaman et al. (2019); Subramanian and Sundararaman (2021). Existing
works have also shown that NLP models are intrinsically capable of numerical
reasoning. This potential can be tapped using simple modifications like those
proposed by Andor et al. (2019), in which the BERT model was used to perform
arithmetic operations. Like the above-mentioned works, we seek to leverage the
potential of existing NLP models to classify numbers into different entity
types and transfer such ability across related tasks to boost performance.

##  3 Preliminaries

For each corpus, we index sentences by $n=1,2,\ldots,N$. Each token in a given
vocabulary $\mathcal{V}$ is assigned a one-hot vector embedding
$\mathbf{v}\in\left\\{0,1\right\\}^{|\mathcal{V}|}$. Tokens in a sequence are
indexed by $i=1,2,\ldots,T_{n}$, where $T_{n}$ is the length of sequence $n$
(after padding). During embedding lookup, tokens are mapped to embeddings
$\mathbf{e}_{ni}\in\mathbb{R}^{d}$ ($d$ refers to the dimension of embedding).
For all tasks, we use pre-trained models including BERT and RoBERTa to form
latent representations of inputs; token embeddings are summed with position
embeddings $\mathbf{p}_{ni}\in\mathbb{R}^{d}$ and segment embeddings
$\mathbf{s}_{ni}\in\mathbb{R}^{d}$ prior to being fed into the Transformer
encoder. The output consists of a classification token
$\mathbf{c}_{n}\in\mathbb{R}^{d}$ and a sequence of contextualized embeddings
$\mathbf{o}_{ni}\in\mathbb{R}^{d}$, one for each token. For each model, we use
an Adam optimizer with $\eta=0.001$, $\beta_{1}=0.9$, $\beta_{2}=0.999$, and
$\epsilon=10^{-7}$ to minimize the empirical loss across samples in each
batch. More experimental settings are detailed in their respective sections.
We train models on Tesla K80 GPUs and TPUs.

##  4 Methods

###  4.1 Number entity recognition

BERT, introduced by Devlin et al. (2018), has been shown to exhibit state-of-
the-art performance for NER. Token embeddings are passed through a Transformer
encoder, generating contextualized representations for each subword token,
which are classified by the attention-based BERT classifier. We fine-tune the
BERT model on the SQuAD-Num dataset. Both questions and context for each
sequence are input to the NuER model, and the model outputs a sequence of tags
denoting whether each token is one of the six types of numbers or “other”
(seven total classes). The outputs of the NuER model for each element in the
sequence, given by the $\mathbf{o}_{i}$, are mapped through a linear layer (we
drop the $n$ subscript for clarity), and the resultant scalars are fed into a
softmax function:

$\displaystyle\mathbf{\hat{y}}_{i}=\frac{\exp(\mathbf{W}\mathbf{o}_{i})}{\sum_{k=1}^{K}\exp(\mathbf{w}_{k}^{\top}\mathbf{o}_{i})}$
(1)

where $\mathbf{w}\in\mathbb{R}^{d}$,
$\mathbf{W}=\left[\mathbf{w}_{1},\ldots,\mathbf{w}_{k},\ldots,\mathbf{w}_{K}\right]^{\top}\in\mathbb{R}^{K\times
d}$ are the parameters of the linear layer, $k=1,\ldots,K$ indexes the entity
labels, and the $\exp$ operator is applied elementwise. The model is trained
with categorical cross-entropy loss:

$\displaystyle\mathcal{L}=\sum_{i=1}^{T}\sum_{k=1}^{K}-y_{ik}\log(\hat{y}_{ik})$
(2)

where $y_{ik}\in\left\\{0,1\right\\}$ denotes whether the $i$th element
belongs to class $k$. We divide SQuAD-Num into 75% training, 10% validation,
and 15% testing. We fine-tune the NuER model for 20 epochs on SQuAD-Num and
evaluate our model on the out-of-domain HA-Numeracy 2K.

###  4.2 Jointly trained embeddings

We experiment with training embeddings for tokens jointly with number
entities, referred as Jointly trained Embeddings (JEM) Wang et al. (2018).
Specifically, we compare the performance of numerical question answering on
SQuAD-Num, with and without entity embeddings. As input to the pretrained
model, we feed in pairs $(\mathbf{v}_{i},\mathbf{z}_{i})$ of tokens
$\mathbf{v}_{i}$ and number entity labels $\mathbf{z}_{i}$.
$\mathbf{z}_{i}\in\left\\{0,1\right\\}^{K+1}$ is a one-hot vector, where $K=6$
indicates the six entities of numbers in the corpus. During embedding lookup,
the pairs of embeddings for each token are summed and fed as inputs:

$\displaystyle\mathbf{\tilde{e}}_{i}=\mathbf{e}_{i}+\mathbf{h}_{i}$ (3)

where $\mathbf{h}_{i}$ denotes the trainable number entity embedding for token
$i$. We compare the performance of this model after fine-tuning for two epochs
against the baseline, which is only fed the token identities as input. As is
standard for SQuAD, the model is trained with cross-entropy loss in order to
predict the probability of any token in the context being the start/end of a
span that answers the numerical question:

$\displaystyle\mathcal{L}=\sum_{i=1}^{T}-y_{i}\log(\hat{y}_{i})$ (4)

where here, $y_{i}\in\\{0,1\\}$ identifies the start/end of a span and
$\hat{y}_{i}$ is the corresponding probability predicted by our model.

|  | Year | Count | Percentage | Age | Size | Date | Total
---|---|---|---|---|---|---|---|---
SQuAD-Num | Precision | 98.61 | 88.66 | 95.08 | 94.12 | 87.84 | 100.00 | 95.38
Recall | 98.99 | 95.91 | 96.67 | 88.89 | 84.42 | 100.00 | 97.24
F1 | 98.8 | 92.15 | 95.87 | 91.43 | 86.09 | 100.00 | 96.3
ENT-Numeracy 600K | Precision | 97.96 | 92.00 | 91.67 | 75.00 | 54.55 | 96.88 | 94.26
Recall | 96.48 | 81.01 | 95.65 | 7.14 | 30.00 | 85.94 | 84.14
F1 | 97.21 | 86.16 | 93.62 | 13.04 | 38.71 | 91.09 | 88.91
Table 2: Performance metrics for NuER on SQuAD-Num and ENT-Numeracy 600K for
each of the six entity types.

###  4.3 FITB

In this subsection, we measure the effectiveness of NuER on classification
tasks. We handcraft a Fill-In-The-Blank (FITB) task of predicting a masked
number in a given input sentence. This is challenging as there are several
thousand unique numbers in the corpus to which a masked number could be
classified. It measures the true capability of BERT to fill in a number given
the context. Our model which we refer to as $BERT_{NuER}$ consists of BERT
tokenizer with six special tokens namely <size>, <age>, <count>, <date>,
<year>, and <percentage>. These special tokens are added next to the CLS as an
indication to the model about their respective number entity types. On the
classification side, a linear projection layer of size $M\times D$ is created
to learn about the six different entity types ($M=6$, $D=768$). These
projections are then concatenated with the BERT hidden states of 768
dimensions to obtain final set of dimensions $(768*2)\times M$ which is
projected again to a classification layer of $V$ output vocabulary items.
Here, the $V$ output vocabulary items represent the unique numbers that are
found in the corpus.

|  | Year | Count
---|---|---|---
Few-shot | Precision | 87.09 | 76.27
Recall | 97.32 | 89.59
F1 | 91.92 | 82.39
Full finetuning | Precision | 91.15 | 80.53
Recall | 98.49 | 95.33
F1 | 94.68 | 87.31
Table 3: Few-shot transfer vs fine-tuning performance on NuER entity types.

##  5 Results and Discussion

###  5.1 Number entity recognition

Table 2 shows the performance of our NuER model on each dataset. We find that
the model performs well across all categories on SQuAD-Num and in most
categories of ENT-Numeracy 600K. Performance drops are observed for the “Age”
and “Size” categories, partially because there are not enough training
instances of these two categories in SQuAD-Num. For “Age,” precision is high
while recall is low. This is a limitation of our work and could be because the
numbers that are classified as “Age” can typically be confused with other
entities including “Size,” suggesting that the model is only able to capture
age accurately when it is quite obvious in the dataset. For “Size,” both
precision and recall are low, perhaps because, unlike the other types of
numbers, size may have many different units (e.g., feet, meters, acres, etc.),
which make it more difficult to associate with the numerical value. Table 3
shows that the classification performance when only few of the samples are
used in fine-tuning is still comparable to the full fine-tuning performance on
two classes. Though it works well for few classes, the F1 scores drop
significantly for classes like Age and Size.

####  5.1.1 Training settings

For Table 2 experiments, the learning rate was set to $2\times 10^{-5}$. For
few-shot used in Table 3, only first 20 samples (using grid search, fewer
samples to retain maximum accuracy compared to full finetuning) from Year and
Count, while full fine-tuning uses all the samples.

Model | Version | Exact match | F1
---|---|---|---
$BERT_{B}$ | Baseline | 90.66 | 90.71
JEM | 91.17 | 91.33
$RoBERTa_{B}$ | Baseline | 95.83 | 95.88
JEM | 96.84 | 96.95
$BERT_{L}$ | Baseline | 93.19 | 93.28
JEM | 94.83 | 95.01
Table 4: Improvements on SQuAD-Num dataset using Jointly trained Embedding
(JEM) technique with $BERT_{B}$ (base), $BERT_{L}$ (large), and $RoBERTa_{B}$
(base). All numbers are tested for significance and have a standard deviation
of 0.05.

###  5.2 Jointly trained embeddings

Table 4 compares the performance of our model trained on numerical question
answering against the BERT baseline. Our model demonstrates gains in both the
exact match and F1 scores, suggesting that the additional information imparted
by knowledge of the type of numbers that are present in the question and
context aids the model in localizing the number within the context. A variety
of models including $BERT_{B}$ (base), $BERT_{L}$ (large) and $RoBERTa_{B}$
(base) Devlin et al. (2018); Liu et al. (2019) are used to understand if these
improvements are consistent. We found that in terms of overall performance,
$RoBERTa_{B}>BERT_{L}>BERT_{B}$, with $RoBERTa_{B}$ being the highest, while
in terms of performance improvements, $BERT_{L}>RoBERTa_{B}>BERT_{B}$ with
$BERT_{L}$ JEM obtaining an improvement of 1.64 exact match and 1.7 F1 points
respectively.

####  5.2.1 Training settings

BERT for QA consists of a pretrained model along with classification heads for
start and ends of spans. The model was trained for 2 epochs with a learning
rate of $2\times 10^{-5}$. During training, embeddings for numbers and their
entities are learned in parallel.

Figure 2: Qualitative comparison of predictions of numbers on the FITB task
between the BERT base and $BERT_{NuER}$. The left column shows BERT base
predictions, while the right columns shows $BERT_{NuER}$ predictions.
$BERT_{NuER}$ predictions align with a number entity more than BERT base.

###  5.3 FITB

Table 5 shows the performance improvements of $BERT_{NuER}$ over the BERT base
for FITB. The top-$k$ performance tells us that $BERT_{NuER}$ is able to
predict the correct answer within the top-$k$ entries. We see that the gap in
performance improvement increases as $k$ increases. This tells us that the
enriched model is more accurate in its predictions. The dist metric quantifies
the average absolute distance between pairs of actual value and the predicted
values. This ensures the quality of numeral predictions. For example, as shown
in the first column of Figure 2, when the ground truth is a Year, we expect
all the other predictions to also be years (as opposed to Count or other
classes). The dist metric calculates the spread of numbers while top-$k$
calculates the accuracy itself. We see that the dist is significantly lower
than the baseline model. Figure 2 shows that $BERT_{NuER}$ generates top-$k$
predictions whose elements are more consistent with the desired class label.

####  5.3.1 Training settings

For the FITB experiment, DistilBERT Sanh et al. (2019) with an Adam optimizer
and a learning rate of $2\times 10^{-4}$ is used. The model is trained for 3
epochs with a batch size of 64. On the architectural side, adding a linear
layer to the BERT hidden states yielded marginal improvements. But,
concatenating the linear layer responsible for learning number entity onto the
BERT hidden states followed by classification lead to significant
improvements.

Top-1 Top-2 Top-5 Top-10 Baseline 36.88 52.51 71.57 79.84 Dist 624 1279.06
3275.22 8240.56 Ours 37.69 53.9 73.1 81.47 Dist 554 1110 2829.73 6435.14

Table 5: Predictive performance on the FITB task. All numbers are tested for
significance and have a standard deviation of 0.01

##  6 Conclusions

We have defined number entities and find that this aids multiple NLP tasks,
via joint training, and classification. We also introduce a dataset
representing number entities and a methodology to annotate unseen numbers. In
the classification objective, we find that the addition of entities provides a
considerable boost in FITB evident through the quality of predictions. Through
our datasets and methods, we hope to incite further research at the
intersection of numeracy and NER.

## Acknowledgments

The authors would like to thank Vishwanath Seshagiri for helping with the
annotations on the SQuAD-Num dataset.

## References

  * Andor et al. (2019) Daniel Andor, Luheng He, Kenton Lee, and Emily Pitler. 2019.  Giving bert a calculator: Finding operations and arguments with reading comprehension.  _arXiv preprint arXiv:1909.00109_.
  * Berg-Kirkpatrick and Spokoyny (2020) Taylor Berg-Kirkpatrick and Daniel Spokoyny. 2020.  An empirical investigation of contextualized number prediction.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 4754–4764.
  * Chen et al. (2019) Chung-Chi Chen, Hen-Hsen Huang, Hiroya Takamura, and Hsin-Hsi Chen. 2019.  Numeracy-600k: learning numeracy for detecting exaggerated information in market comments.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 6307–6313.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_.
  * Geva et al. (2020) Mor Geva, Ankit Gupta, and Jonathan Berant. 2020.  Injecting numerical reasoning skills into language models.  _arXiv preprint arXiv:2004.04487_.
  * Lample et al. (2016) Guillaume Lample, Miguel Ballesteros, Sandeep Subramanian, Kazuya Kawakami, and Chris Dyer. 2016.  Neural architectures for named entity recognition.  _arXiv preprint arXiv:1603.01360_.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_.
  * Naik et al. (2019) Aakanksha Naik, Abhilasha Ravichander, Carolyn Rose, and Eduard Hovy. 2019.  Exploring numeracy in word embeddings.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 3374–3380.
  * Rajpurkar et al. (2016) Pranav Rajpurkar, Jian Zhang, Konstantin Lopyrev, and Percy Liang. 2016.  Squad: 100,000+ questions for machine comprehension of text.  _arXiv preprint arXiv:1606.05250_.
  * Sang and De Meulder (2003) Erik F Sang and Fien De Meulder. 2003.  Introduction to the conll-2003 shared task: Language-independent named entity recognition.  _arXiv preprint cs/0306050_.
  * Sanh et al. (2019) Victor Sanh, Lysandre Debut, Julien Chaumond, and Thomas Wolf. 2019.  Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter.  _arXiv preprint arXiv:1910.01108_.
  * Spithourakis and Riedel (2018) Georgios P Spithourakis and Sebastian Riedel. 2018.  Numeracy for language models: Evaluating and improving their ability to predict numbers.  _arXiv preprint arXiv:1805.08154_.
  * Subramanian and Sundararaman (2021) Vivek Subramanian and Dhanasekar Sundararaman. 2021.  How do lexical semantics affect translation? an empirical study.  _arXiv preprint arXiv:2201.00075_.
  * Sundararaman et al. (2020) Dhanasekar Sundararaman, Shijing Si, Vivek Subramanian, Guoyin Wang, Devamanyu Hazarika, and Lawrence Carin. 2020.  Methods for numeracy-preserving word embeddings.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 4742–4753.
  * Sundararaman et al. (2019) Dhanasekar Sundararaman, Vivek Subramanian, Guoyin Wang, Shijing Si, Dinghan Shen, Dong Wang, and Lawrence Carin. 2019.  Syntax-infused transformer and bert models for machine translation and natural language understanding.  _arXiv preprint arXiv:1911.06156_.
  * Tanwar et al. (2022) Ashwani Tanwar, Jingqing Zhang, Julia Ive, Vibhor Gupta, and Yike Guo. 2022.  Unsupervised numerical reasoning to extract phenotypes from clinical text by leveraging external knowledge.  _arXiv preprint arXiv:2204.10202_.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Advances in neural information processing systems_ , pages 5998–6008.
  * Wallace et al. (2019) Eric Wallace, Yizhong Wang, Sujian Li, Sameer Singh, and Matt Gardner. 2019.  Do nlp models know numbers? probing numeracy in embeddings.  _arXiv preprint arXiv:1909.07940_.
  * Wang et al. (2018) Guoyin Wang, Chunyuan Li, Wenlin Wang, Yizhe Zhang, Dinghan Shen, Xinyuan Zhang, Ricardo Henao, and Lawrence Carin. 2018.  Joint embedding of words and labels for text classification.  _arXiv preprint arXiv:1805.04174_.
  * Zhang et al. (2020) Xikun Zhang, Deepak Ramachandran, Ian Tenney, Yanai Elazar, and Dan Roth. 2020.  Do language embeddings capture scales?  _arXiv preprint arXiv:2010.05345_.