# Stress out of charmonia

Siqi Xu  Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou
730000, China  School of Nuclear Science and Technology, University of Chinese
Academy of Sciences, Beijing 100049, China  Xianghui Cao  Department of Modern
Physics, University of Science and Technology of China, Hefei, Anhui 230026,
China  Tianyang Hu  Department of Modern Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China  Yang Li  Department of Modern
Physics, University of Science and Technology of China, Hefei, Anhui 230026,
China  Department of Physics and Astronomy, Iowa State University, Ames, Iowa
50010, U.S.  Xingbo Zhao  Institute of 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  James P. Vary
Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50010,
U.S.

###### Abstract

We investigate the gravitational form factors of charmonium. Our method is
based on a Hamiltonian formalism on the light front known as basis light-front
quantization. The charmonium mass spectrum and light-front wave functions were
obtained from diagonalizing an effective Hamiltonian that incorporates
confinement from holographic QCD and one-gluon exchange interaction from
light-front QCD. We proposed a quantum many-body approach to construct the
hadronic matrix elements of the energy momentum tensor $T^{++}$ and $T^{+-}$,
which are used to extract the gravitational form factors $A(Q^{2})$ and
$D(Q^{2})$. The obtained form factors satisfy the known constraints, e.g. von
Laue condition. From these quantities, we also extract the energy, pressure
and invariant mass squared distributions of the system. We find that hadrons
are multi-layer systems. We also find that there is a tachyonic core within
charmonium.

##  I Introduction

The year 2023 marked the 50 year anniversary of quantum chromodynamics (QCD).
However, our understanding of the strong interaction remains incomplete [1].
One of the challenges is the strong-coupling nature of QCD. As such, the
strong force within hadrons is dramatically different from what the QCD
Lagrangian naïvely tells us, in contrast to weakly coupled systems, e.g. the
hydrogen atom. The primary operators that describe the distribution of stress
within a system are the energy-momentum tensor (EMT) $T^{\mu\nu}$ [2]. These
operators also dictate how the system gravitates, and conversely, how it
responds to the gravitational field [3, 175]. The access to the hadronic
matrix elements (HME) of the EMT has been a hot topic in recent literature.
Kumano et al. extracted the EMT of the pion from the process
$\gamma^{*}\gamma\to\pi^{0}\pi^{0}$ measured on Belle [5]. Burkert et al.
obtained the stress inside the proton using deeply virtual Compton scattering
data and deeply vector meson production data collected from Jefferson Lab [6,
7, 8, 9, 10]. Duran et. al. accessed the gluon contributions for the proton
using near-threshold electroproduction of $J/\psi$ [11]. The precision of
these experiments is expected to dramatically improve in the forthcoming era
of electron-ion colliders [12, 13, 14]. Lattice simulation of the nucleon and
the pion EMT are recently available with a pion mass
$M_{\pi}=170\,\mathrm{MeV}$ [15, 16] (cf. [17, 18, 19, 20]). The hadronic EMT
were also investigated in various theories and phenomenological models, e.g.,
[21, 22, 23, 24, 25, 26, 27, 28, 29, 172, 31, 32, 33, 34, 35, 36, 37, 38, 176,
40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 45, 179, 54, 55]. See
Refs. [56, 57] for recent reviews.

One of the challenges to access the hadronic EMT is that $T^{\mu\nu}$ contains
the interaction, which is strong at low energy. In this work, we present the
first calculation of the charmonium EMT. Charmonium, first discovered 50 years
ago in the November revolution, is an ideal system to probe the properties of
the strong force. This system consists of a pair of charm and anti-charm
quarks. The relatively large charm quark mass $m_{c}\gg\Lambda_{\text{QCD}}$
makes the perturbative QCD description of the quark-antiquark interaction
possible. On the other hand, with the binding energy
$\alpha_{s}m_{c}\gtrsim\Lambda_{\text{QCD}}$, a non-perturbative treatment is
required to describe the properties of the QCD bound states including the
stress within. Hence, this bound system involves the interplay of the
perturbative and non-perturbative physics of QCD. It is sometimes dubbed as
the “hydrogen atom of QCD” [58].

Our method is based on basis light-front quantization (BLFQ) [59]. This method
stems from the light-front Hamiltonian formalism, which finds its roots from
two lineages [60, 61, 62, 63, 64, 65, 66, 67, 68, 69]. One is Dirac [70] and
Weinberg’s [71] efforts to simplify relativistic dynamics, which leads to the
discovery of the front form of quantum field theory [72, 73, 74, 75, 76, 77,
78, 79, 80, 81, 82, 83, 84, 85]. The other is the light-cone current algebra
[86, 87, 88, 89, 90, 91] and Feynman’s parton model [92, 93, 94, 95, 96]
emerging in deep inelastic scattering [97, 98]. These two complementary
perspectives are unified in the light-front Hamiltonian formalism, which
allows us to solve QCD as a quantum many-body problem, as described by
discretized light-cone quantization (DLCQ) [99] and BLFQ [59].

BLFQ was successfully applied to a number of systems including charmonium, and
to a number of hadronic observables including mass spectra, radiative
transitions, parton distributions and transverse momentum distributions [101,
102, 103, 104, 105, 106, 107, 108, 109, 52, 110, 111, 112, 113, 114, 115, 116,
117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131,
132, 133, 134, 129, 135, 136, 137, 138, 139, 140, 141, 142]. For charmonium,
the effective Hamiltonian combines confinement from holographic light-front
QCD [143] for the long-distance physics as well as one-gluon exchange
interaction from light-front QCD [103] for the short-distance physics [110].
The computed mass spectrum is shown to be within 40 MeV deviation with the PDG
compilation [111, 115, 117, 118]. The obtained light-front wave functions
(LFWFs) were used to investigate a variety of observables, including the decay
constants [110, 111], electromagnetic form factors [114], parton distribution
functions [116], and GPDs [114]. In particular, the parameter-free predictions
of the radiative transition widths and radiative transition form factors are
in remarkable agreement with the experimental measurements whenever available
[119, 120, 121].

The rest of the article is organized as follows. Sect. II briefly introduces
the light-front Hamiltonian formalism including the basis function
representation chosen for BLFQ. Sect. III describes the energy-momentum tensor
and the light-front densities associated with these operators. The light-front
wave function representation of the gravitational form factors $A(Q^{2})$ and
$D(Q^{2})$ are presented in Sect. IV. We present the numerical results in
Sect. V. Finally, we conclude in Sect. VI.

##  II Hamiltonian field theory on the light front

In light-front QCD, the quantum state of a hadron $|\psi_{h}\rangle$ can be
obtained by solving the light-front Schrödinger equation [63],

$i\frac{\partial}{\partial
x^{+}}|\psi_{h}(x^{+})\rangle=\frac{1}{2}P^{-}|\psi_{h}(x^{+})\rangle\,.$ (1)

Here, we adopt the light-front coordinates, $x^{\pm}=x^{0}\pm x^{3}$,
$\vec{x}_{\perp}=(x^{1},x^{2})$. The hadronic state vector
$|\psi_{h}(p,j,m_{j})\rangle$ can be further classified by its momentum
$p^{\mu}$, total spin $j$ and spin magnetic projection $m_{j}$, viz. [144,
63],

$\displaystyle
P^{\mu}|\psi_{h}(p,j,m_{j})\rangle=p^{\mu}|\psi_{h}(p,j,m_{j})\rangle,$ (2)
$\displaystyle\vec{S}^{2}|\psi_{h}(p,j,m_{j})\rangle=j(j+1)|\psi_{h}(p,j,m_{j})\rangle,$
(3) $\displaystyle
S^{+}|\psi_{h}(p,j,m_{j})\rangle=m_{j}|\psi_{h}(p,j,m_{j})\rangle,$ (4)

where, $P^{\mu}$ is the 4-momentum operator, $\vec{S}^{2}=W^{2}/P^{2}$ is the
total spin operator, constructed from the Pauli-Lubanski vector
$W^{\mu}=\frac{1}{2}\varepsilon^{\mu\nu\rho\sigma}P_{\nu}J_{\rho\sigma}$,
where $J^{\alpha\beta}$ is the generalized angular momentum operator, i.e.,
the generator of the Lorentz symmetry, and $S^{+}=W^{+}/P^{+}$ is the magnetic
projection of the total spin operator, also known as the light-front helicity
operator. The light-front energy for a free hadron is
$p^{-}=(p^{2}_{\perp}+M^{2})/p^{+}$, where $M$ is the invariant mass of the
particle.

These operators are not all kinematical, i.e. some of these operators contain
interactions. The dynamical nature of operators depends on the choice of the
initial surface, i.e., the quantization scheme. Light-front quantization with
initial surface $x^{+}=0$ turns out to be equipped with the maximal number (7)
of kinematical operators out of the 10 Poincaré generators [70]. The remaining
3 generators originate from the same local interaction in quantum field theory
[64],

$\displaystyle P^{\mu}_{\text{int}}=\,$
$\displaystyle\omega^{\mu}\int{\mathrm{d}}^{4}x\,\delta(\omega\cdot
x)\mathcal{H}_{\text{int}}(x),$ (5) $\displaystyle J^{\mu\nu}_{\text{int}}=\,$
$\displaystyle\int{\mathrm{d}}^{4}x\,\delta(\omega\cdot
x)(x^{\mu}\omega^{\nu}-\omega^{\mu}x^{\nu})\mathcal{H}_{\text{int}}(x)\,.$ (6)

Here,
$\omega^{\mu}=(\omega^{+},\omega^{-},\vec{\omega}_{\perp})=(0,2,\vec{0}_{\perp})$
is the null vector that indicates the orientation of the light-front
quantization surface. On the other hand, the Poincaré generators stem from the
Noether currents, e.g.,

$P^{\mu}=\int{\mathrm{d}}^{3}x\,T^{+\mu}(x)\,,$ (7)

where, $T^{\mu\nu}$ is the EMT. In particular, the density of the light-front
interaction energy $\mathcal{H}_{\text{int}}$ is related to the $+-$ component
of the EMT, $\mathcal{H}_{\text{int}}=\frac{1}{2}T^{+-}_{\text{int}}$.
Momentum conservation (2) imposes constraints on the forward HMEs of the EMT.
It is straightforward to show,

$\langle\psi_{h}(p,j,m_{j}^{\prime})|T^{+\mu}(0)|\psi_{h}(p,j,m_{j})\rangle=2p^{+}p^{\mu}\delta_{m_{j},m^{\prime}_{j}}\,.$
(8)

Similar constraints can be obtained from the conservation of the angular
momentum (3–4) as shown by Lorcé et. al. [145].

In QCD, the force is strong and we need a non-perturbative approach to solve
the eigenvalue problems (2–4). Inspired by the success of ab initio nuclear
structure theory, especially the no-core shell model (NCSM) [146], Vary et.
al. proposed BLFQ as a non-perturbative computational framework to solve
light-front QCD [59]. In BLFQ, the light-front Schrödinger equation is cast
into a matrix eigenvalue equation,

$\sum_{j}H_{ij}\psi_{j}=M^{2}\psi_{i}\,,$ (9)

where $H_{ij}=\langle i|H|j\rangle$ and $\psi_{i}=\langle i|\psi\rangle$. Here
$|i\rangle$ is the basis state and $|\psi\rangle$ is the state vector of a
hadron, and $H=P^{+}P^{-}-\vec{P}_{\perp}^{2}$ is equivalent to the light-
front Hamiltonian operator $P^{-}$ since $P^{+}$ and $\vec{P}^{2}_{\perp}$ are
kinematical. The basis states are chosen to retain the maximal kinematical
symmetries of the system. A convenient choice is the 2D harmonic oscillator
function in the transverse direction and Jacobi polynomials in the
longitudinal direction. The many-body basis is then constructed from the
tensor product of the single-particle basis. Other choices are also available,
such as the discretized momentum basis adopted in DLCQ [147], conformal basis
functions employed in Hamiltonian conformal truncation method [148], and
coherent basis in light-front coupled cluster method [149]. These bases have
been shown successful in applications to low-dimensional or simpler quantum
field theories [68].

In BLFQ, one adopts an effective Hamiltonian suitable for the finite basis
space with low resolution. Two approaches are adopted. In the top-down
approach, one starts from QCD and perform a similarity renormalization group
evolution [150, 100, 151, 152, 153]. In the bottom-up approach, one starts
with a phenomenological effective Hamiltonian that embodies the key physics of
the system [110, 111]. The latter is similar to the nuclear shell model in
nuclear structure calculations, where the role of the first approximation, the
nuclear shell model, is played by holographic light-front QCD, as shown by
Brodsky and de Téramond [143].

In this work, we adopt the phenomenological effective Hamiltonian with a
confining interaction from light-front holographic QCD. The form is remarkably
simple, $V_{q\bar{q}}\propto x(1-x)r^{2}_{\perp}$, where,
$x=p^{+}_{q}/P^{+}_{h}$ is the longitudinal momentum fraction of the quark,
and $\vec{r}_{\perp}$ is the relative transverse separation of the quark and
antiquark. This potential generates a 2D harmonic oscillator wave functions as
the eigenfunction, which coincides with the basis function in BLFQ in the
transverse direction which provides computational advantages such as enabling
factorization of the center-of-mass motion and facilitation of the
transformation between relative and single-particle coordinates. A
longitudinal interaction of the ’t Hooft type is supplemented to generate
confinement in the longitudinal direction [154, 107, 155, 108]. We further
supplement confinement with a one-gluon exchange interaction obtained from
light-front QCD. This piece of interaction is essential for describing the
short-distance physics as well as the spin structure [103].

The basis functions are the eigen-functions of the confining interactions.
Their analytic forms are,

$\displaystyle\phi_{nm}(\vec{v}_{\perp})=\,$ $\displaystyle
N_{nm}\kappa^{-1}\Big{(}\frac{v_{\perp}}{\kappa}\Big{)}^{|m|}e^{im\theta-\frac{v^{2}_{\perp}}{2\kappa^{2}}}L_{n}^{|m|}\Big{(}\frac{v^{2}_{\perp}}{\kappa^{2}}\Big{)},$
(10) $\displaystyle\chi_{l}(x)=\,$ $\displaystyle
N_{l}x^{\frac{\beta}{2}}(1-x)^{\frac{\alpha}{2}}P_{l}^{(\alpha,\beta)}(2x-1),$
(11)

where,
$N_{l}=\sqrt{4\pi(2l+\alpha+\beta+1)}\sqrt{\frac{\Gamma(l+1)\Gamma(l+\alpha+\beta+1)}{\Gamma(l+\alpha+1)\Gamma(l+\beta+1)}}$,
$N_{nm}=\sqrt{\frac{4\pi n!}{(n+|m|)!}}$,
$\vec{v}_{\perp}=\vec{p}_{\perp}/\sqrt{x(1-x)}$ is the holographic momentum,
$v_{\perp}=|\vec{v}_{\perp}|$, and $\theta=\arg\vec{v}_{\perp}$,
$x=p^{+}/P^{+}$ is the longitudinal momentum fraction. The basis parameters
$\alpha={2m_{\bar{q}}(m_{q}+m_{\bar{q}})}/{\kappa^{2}}$,
$\beta={2m_{q}(m_{q}+m_{\bar{q}})}/{\kappa^{2}}$, where $\kappa$ is the
strength of the confining interaction and $m_{q}=m_{\bar{q}}=m_{c}$ is the
effective quark mass. Then effective Hamiltonian for the $c\bar{c}$ system is,

$H_{\text{eff}}=\sum_{I}h_{I}a^{\dagger}_{I}a_{I}+\frac{1}{2}\sum_{I,J,I^{\prime},J^{\prime}}V_{IJI^{\prime}J^{\prime}}a^{\dagger}_{I}a^{\dagger}_{J}a_{J^{\prime}}a_{I^{\prime}}\,.$
(12)

Here, $I,J,I^{\prime},J^{\prime}$ each represent a collection of single-
particle quantum numbers: radial quantum number $n$, angular quantum number
$m$, longitudinal quantum number $l$, spin $s$, color $i$ and parton flavor
$f$, e.g., $I=\\{n,m,l,s,i,f\\}$. The one-body effective Hamiltonian
$h_{I}=4m_{c}^{2}+2\kappa^{2}(2n+|m|+l+\frac{3}{2})+({\kappa^{4}}/{4m_{c}^{2}})l(l+1)$
incorporates the kinetic energy as well as the long-distance interactions,
i.e. confinement. We adopted the one-gluon exchange interaction as the two-
body effective interaction, which is the dominant physics at short distance
for charmonium. For charmonium, we further impose a valence approximation,
i.e. keeping only the $c\bar{c}$ Fock sector.

BLFQ adopts the $N_{\text{max}}$ basis truncation scheme to render the basis
space finite. Specifically for charmonium, we impose the conditions,

$2n+|m|+1\leq N_{\text{max}},\quad 0\leq l\leq L_{\text{max}}.$ (13)

$N_{\text{max}}$ defines the UV and IR regulators of the basis space as
$\Lambda_{\textsc{uv}}=\sqrt{N_{\text{max}}}\kappa$, and
$\lambda_{\textsc{ir}}=\kappa/\sqrt{N_{\text{max}}}$. $L_{\text{max}}$ defines
the resolution of the longitudinal momentum fraction $\delta
x=L^{-1}_{\text{max}}$. For simplicity, we have tied the basis scale parameter
to the confining strength $\kappa$.

We diagonalized the Hamiltonian with different $N_{\text{max}},L_{\text{max}}$
to obtain the mass squared eigenvalues and the LFWFs. The model parameters
$\kappa$, and $m_{c}$ were fit to the charmonium mass spectrum below the
$D\overline{D}$ threshold, which is a good approximation for states with
narrow resonance. Since radiative corrections are neglected, we adopt
$N_{\text{max}}=L_{\text{max}}=8$ for observables that are sensitive to the
short-distance quantum fluctuations, which corresponds to a UV resolution
$\Lambda_{\textsc{uv}}\approx M_{c\bar{c}}$. To estimate the uncertainty
associated with the basis resolution, we quote the difference between the
$N_{\text{max}}=L_{\text{max}}=8$ results and the
$N_{\text{max}}=L_{\text{max}}=16$ results. The details of the model can be
found in Ref. [111].

From the basis functions, we can reconstruct the valence LFWFs,

$\psi^{(m_{j})}_{s\bar{s}/h}(x,\vec{k}_{\perp})=\sum_{n,m,l}\psi^{(m_{j})}_{s\bar{s}/h}(n,m,l)\\\
\times\phi_{nm}(\vec{k}_{\perp}/\sqrt{x(1-x)})\chi_{l}(x)\,,$ (14)

where, $\psi^{(m_{j})}_{s\bar{s}/h}(n,m,l)$ are the basis coefficients
obtained from the diagonalization. The state vector can be represented by the
LFWFs as,

$|\psi_{h}(P,j,m_{j})\rangle=\sum_{s,\bar{s}}\int_{0}^{1}\frac{{\mathrm{d}}x}{2x(1-x)}\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{3}}\\\
\times\psi_{s\bar{s}/h}^{(m_{j})}(x,\vec{k}_{\perp})\frac{1}{\sqrt{N_{c}}}\sum_{i=1}^{N_{c}}b^{\dagger}_{si}\big{(}xP^{+},\vec{k}_{\perp}+x\vec{P}_{\perp}\big{)}\\\
\times
d^{\dagger}_{\bar{s}i}\big{(}(1-x)P^{+},-\vec{k}_{\perp}+(1-x)\vec{P}_{\perp}\big{)}|0\rangle.$
(15)

Here, $P^{\mu}=(P^{-},P^{+},\vec{P}_{\perp})$ is the 4-momentum of the bound
state $h$. $\vec{k}_{\perp}=\vec{p}_{\perp}-x\vec{P}_{\perp}$ is the relative
transverse momentum, and $x=p^{+}/P^{+}$ is the longitudinal momentum
fraction. The charmonium LFWFs are available in Mendeley Data repository
[156].

##  III Energy-momentum tensor

By virtue of Poincaré symmetry, the HMEs of the EMT for a spin-0 particle,
such as $\eta_{c}$ and $\chi_{c0}$, can be parametrized in terms of covariant
tensors [56],

$\langle P+{\frac{1}{2}}q|T^{\mu\nu}(0)|P-{\frac{1}{2}}q\rangle=\\\
2P^{\mu}P^{\mu}A(-q^{2})+\frac{1}{2}(q^{\mu}q^{\nu}-q^{2}g^{\mu\nu})D(-q^{2})\,,$
(16)

where Lorentz scalar functions $A(-q^{2})$ and $D(-q^{2})$ are known as the
gravitational form factors (GFFs).

The 2D Fourier transform of the HMEs in the Drell-Yan frame $q^{+}=0$ can be
interpreted as the hadronic EMT on the transverse plane [157, 158, 159, 160,
161, 49, 162],

$\mathcal{T}^{\alpha\beta}(\vec{r}_{\perp};P)=\int\frac{{\mathrm{d}}^{2}q_{\perp}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{r}_{\perp}}t^{\alpha\beta}(\vec{q}_{\perp};P)\,,$
(17)

where, $t^{\alpha\beta}(\vec{q}_{\perp};P)$ is the normalized hadronic matrix
elements,

$t^{\alpha\beta}(\vec{q}_{\perp};P)=\frac{1}{2P^{+}}\langle
P+\frac{1}{2}q|T^{\alpha\beta}(0)|P-\frac{1}{2}q\rangle\,.$ (18)

The factor $1/(2P^{+})$ is due to the normalization of the state vector:

$\langle
p^{\prime}|p\rangle=2p^{+}(2\pi)^{3}\delta(p^{+}-p^{\prime+})\delta^{2}(p_{\perp}-p^{\prime}_{\perp})\,.$
(19)

Note that $t^{\alpha\beta}$ hence $\mathcal{T}^{\alpha\beta}$ depends on the
c.m. momentum $P^{\mu}$ of the hadron [159]. The density of the light-front
longitudinal momentum $P^{+}$ is,

$\mathcal{T}^{++}(\vec{r}_{\perp};P)=P^{+}\int\frac{{\mathrm{d}}^{2}q_{\perp}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{r}_{\perp}}A(q_{\perp}^{2})\,.$
(20)

Similarly, the density of the light-front transverse momentum
$\vec{P}_{\perp}$ is,

$\mathcal{T}^{+i}(\vec{r}_{\perp};P)=P^{i}_{\perp}\int\frac{{\mathrm{d}}^{2}q_{\perp}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{r}_{\perp}}A(q_{\perp}^{2})\,.$
(21)

From these expressions, the Fourier transform of GFF $A(q_{\perp}^{2})$,

$\mathcal{A}(r_{\perp})=\int\frac{{\mathrm{d}}^{2}q_{\perp}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{r}_{\perp}}A(q_{\perp}^{2}),$
(22)

can be interpreted as the matter density. In Ref. [161], this density is
called the internal light-front longitudinal momentum density. We emphasize
that $\mathcal{A}(r_{\perp})$ describes the convective distribution of all
three momenta. Note that for spin-$\frac{1}{2}$ particles, such as the
nucleon, there is an extra contribution from spin,

$\mathcal{T}^{+i}_{ss}(\vec{r}_{\perp};P)=P^{i}_{\perp}\mathcal{A}(r_{\perp})+\big{(}\nabla\times\vec{\mathcal{S}}(r_{\perp})\big{)}^{i}_{\perp}\,.$
(23)

Here, $\vec{\mathcal{S}}(r_{\perp})=2\mathcal{J}(r_{\perp})\vec{s}$ is the
spin density, $\vec{s}=s\hat{z}$, and $\mathcal{J}(r_{\perp})$ is the 2D
Fourier transform of the GFF $J(Q^{2})$, which exists for spin-$\frac{1}{2}$
particles. N.B. $\mathcal{J}(r_{\perp})$ is normalized to 1/2. The spin
current $\nabla\times\vec{\mathcal{S}}$ is non-convective.

Momentum conservation (8) implies the conservation of matter,

$A(0)=1\;\Leftrightarrow\;\int{\mathrm{d}}^{2}r_{\perp}\,\mathcal{A}(r_{\perp})=1.$
(24)

The r.m.s radius of the matter distribution is,

$r^{2}_{A}\equiv\frac{3}{2}\int{\mathrm{d}}^{2}r_{\perp}\,r^{2}_{\perp}\,\mathcal{A}(r_{\perp})=-6A^{\prime}(0).$
(25)

The distribution of the light-front energy $P^{-}=(P^{2}_{\perp}+M^{2})/P^{+}$
in a hadron state with mass $M$ reads [49],

$\mathcal{T}^{+-}(r_{\perp};P)=\frac{P^{2}_{\perp}\mathcal{A}(r_{\perp})+\mathcal{M}^{2}(r_{\perp})}{P^{+}}\,$
(26)

where,

$\displaystyle\mathcal{M}^{2}(r_{\perp})=\,$
$\displaystyle\int\frac{{\mathrm{d}}^{2}q_{\perp}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{r}_{\perp}}M^{2}(q^{2}_{\perp}),$
(27) $\displaystyle M^{2}(q^{2}_{\perp})=\,$
$\displaystyle\Big{(}M^{2}+\frac{1}{4}q^{2}_{\perp}\Big{)}A(q^{2}_{\perp})+\frac{1}{2}q_{\perp}^{2}D(q^{2}_{\perp}).$
(28)

The first term ${\mathcal{A}(r_{\perp})P^{2}_{\perp}}/{P^{+}}$ represents the
light-front energy due to the motion of the particle. The second term in (26)
can be interpreted as the distribution of the invariant mass squared $M^{2}$
[49]. Hence, we have abused the notation and introduced the mass squared form
factor $M^{2}(q^{2}_{\perp})$ and its Fourier transform, the mass squared
density $\mathcal{M}^{2}(r_{\perp})$. In Ref. [161], this density [Eq. (27)]
is termed the internal light-front energy density. For spin-$\frac{1}{2}$
particles, the spin term also contributes to the light-front energy density
$\mathcal{T}^{+-}(r_{\perp})$,

$\mathcal{T}^{+-}_{ss}(r_{\perp};P)=\\\
\frac{P^{2}_{\perp}\mathcal{A}(r_{\perp})+\mathcal{M}^{2}(r_{\perp})+\vec{P}_{\perp}\cdot(\nabla\times\vec{\mathcal{S}}(r_{\perp}))}{P^{+}}\,.$
(29)

The last term vanishes upon integration.

Light-front energy conservation (8) implies [49],

$M^{2}(0)=M^{2}\;\Leftrightarrow\;\int{\mathrm{d}}^{2}r_{\perp}\,\mathcal{M}^{2}(r_{\perp})=M^{2}.$
(30)

From (28), the above condition is equivalent to

$\lim_{q_{\perp}\to 0}q_{\perp}^{2}D(q_{\perp}^{2})=0.$ (31)

This condition is known as the von Laue condition [56]. It is useful to define
the r.m.s. radius of the invariant mass squared distribution as,

$\displaystyle r^{2}_{M^{2}}=\,$
$\displaystyle\frac{3}{2M^{2}}\int{\mathrm{d}}^{2}r_{\perp}\,r_{\perp}^{2}\,\mathcal{M}^{2}(r_{\perp})$
(32) $\displaystyle=\,$ $\displaystyle
r_{A}^{2}-\frac{3}{2}\lambda^{2}_{C}(1+2D)\,.$ (33)

Here, $\lambda_{C}=1/M$ is the Compton wavelength of the hadron. $D=D(0)$ is
called the $D$-term, also called the Druck term, and is dubbed as the last
global unknown of hadrons [56].

The above definitions show the advantage of light-front formalism in
extracting frame-independent hadronic densities. In general, the hadronic EMT
can be parametrized by (the proper) energy density $\mathcal{E}(x)$, pressure
$\mathcal{P}(x)$ and (traceless) shear $\Pi^{\alpha\beta}(x)$ as,

$\mathcal{T}^{\alpha\beta}=(\mathcal{E}+\mathcal{P})u^{\alpha}u^{\beta}-\mathcal{P}g^{\alpha\beta}+\Pi^{\alpha\beta}.$
(34)

Here, $u^{\alpha}\equiv P^{\alpha}/\sqrt{P^{2}}$ is the 4-velocity vector. The
energy density can be extracted as,

$\mathcal{E}(r_{\perp})=\int\frac{{\mathrm{d}}^{2}q_{\perp}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{r}_{\perp}}E(q^{2}_{\perp}),$
(35)

where, the energy form factor $E(q^{2}_{\perp})$ is defined as,

$E(q^{2}_{\perp})=\Big{(}M+\frac{q^{2}_{\perp}}{4M}\Big{)}A(q^{2}_{\perp})+\frac{q^{2}_{\perp}}{4M}D(q^{2}_{\perp})\,.$
(36)

For spin-$\frac{1}{2}$ particles, the energy form factor should include a spin
contribution,

$E(q^{2}_{\perp})=\Big{(}M+\frac{q^{2}_{\perp}}{4M}\Big{)}A(q^{2}_{\perp})\\\
+\frac{q^{2}_{\perp}}{4M}\Big{[}D(q^{2}_{\perp})-2J(q^{2}_{\perp})\Big{]}\,.$
(37)

This is in contrast to $M^{2}(q^{2}_{\perp})$, the invariant mass squared form
factor, in which, the spin contribution is absent. The energy-momentum
conservation (24) & (30) [or (31)] implies

$E(0)=M\;\Leftrightarrow\;\int{\mathrm{d}}^{2}r_{\perp}\,\mathcal{E}(r_{\perp})=M,$
(38)

i.e. energy density is normalized to the total energy in local rest frame,
i.e. the mass. The r.m.s. radius of the energy distribution is,

$r^{2}_{E}=\frac{3}{2M}\int{\mathrm{d}}^{2}r_{\perp}\,r^{2}_{\perp}\mathcal{E}(r_{\perp})=r^{2}_{A}-\frac{3}{2}\lambda^{2}_{C}\big{(}1+D\big{)}.$
(39)

The pressure can be extracted as,

$\mathcal{P}(r_{\perp})=-\frac{1}{6M}\int\frac{{\mathrm{d}}^{2}q_{\perp}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{r}_{\perp}}q^{2}_{\perp}D(q_{\perp}^{2}).$
(40)

The form factor $q^{2}_{\perp}D(q^{2}_{\perp})$ is also known as the hadronic
cosmological constant [163], in analogy to the quantum expectation value of
the EMT with respect to the vacuum, $\langle
0|T^{\mu\nu}|0\rangle=g^{\mu\nu}\Lambda$. One can readily identify
$\Lambda=\frac{1}{2}q^{2}D(q^{2})$. The von Laue condition (31) implies that
hadrons are in mechanical equilibrium,

$\int{\mathrm{d}}^{2}r_{\perp}\,\mathcal{P}(r_{\perp})=0.$ (41)

The shear is also related to the GFF $D(q^{2}_{\perp})$,

$\displaystyle\Pi^{ij}(\vec{r}_{\perp})=\,$
$\displaystyle\frac{1}{4M}\int\frac{{\mathrm{d}}^{2}q_{\perp}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{r}_{\perp}}$
(42)
$\displaystyle\times\big{(}q^{i}_{\perp}q^{j}_{\perp}-\frac{1}{3}\delta^{ij}q^{2}_{\perp}\big{)}D(q^{2}_{\perp}),$
(43) $\displaystyle=\,$
$\displaystyle\frac{1}{2}\big{(}\delta^{ij}-3\hat{r}^{i}_{\perp}\hat{r}^{j}_{\perp}\big{)}\mathcal{P}(r_{\perp}),$
(44)

where, $\hat{r}_{\perp}^{i}=r^{i}_{\perp}/r_{\perp}$. In fact, pressure and
shear are the trace and traceless parts of the Cauchy stress tensor
$\mathcal{C}^{\alpha\beta}=\Pi^{\alpha\beta}-\mathcal{P}\Delta^{\alpha\beta}$.

Note that the energy density $\mathcal{E}(r_{\perp})$ is related to from the
invariant mass squared density $\mathcal{M}^{2}(r_{\perp})$ through,

$\mathcal{M}^{2}(r_{\perp})=M\Big{[}\mathcal{E}(r_{\perp})-\frac{3}{2}\mathcal{P}(r_{\perp})\Big{]}.$
(45)

Physically, the proper energy density $\mathcal{E}$ is generally assumed
positive, a constraint known as the weak energy condition, which provides a
constraint on $D$:

$D\leq\frac{2}{3}M^{2}r^{2}_{A}.$ (46)

Similarly, a related condition, known as the strong energy condition
$\mathcal{E}\geq 3\mathcal{P}$ implies, $D\leq({2}/{5})M^{2}r^{2}_{A}$.

A related quantity of interest is the scalar density,
$\theta=\mathcal{T}^{\alpha}_{\;\;\alpha}$, which obeys the classical relation
[164],

$\theta(r_{\perp})=\mathcal{E}(r_{\perp})-3\mathcal{P}(r_{\perp}).$ (47)

The scalar radius is
$r^{2}_{\theta}=r^{2}_{A}-\frac{3}{2}\lambda_{C}^{2}(1+3D)$.

The $D$-term is related to the second moment of the pressure,

$D\equiv
D(0)=M\int{\mathrm{d}}^{2}r_{\perp}\,r^{2}_{\perp}\mathcal{P}(r_{\perp}).$
(48)

Polyakov et al. conjectured that $D<0$ for a mechanically stable system [56].
Intuitively, a negative $D$ represents a repulsive core along with an
attractive periphery. A negative $D$-term suggests a chain of inequalities,

$r_{\theta}>r_{M^{2}}>r_{E}.$ (49)

For charmonium, as we will show below, $D\sim-5$. Then, $r_{E}>r_{A}$ for this
system. Then, we obtain a layered picture of charmonium as shown in Fig. 1.

Figure 1: A 3D picture of the normalized matter density
$\mathcal{A}(r_{\perp})$, energy density $\mathcal{E}(r_{\perp})$, invariant
mass squared density $\mathcal{M}^{2}(r_{\perp})$, scalar density
$\theta(r_{\perp})$ and pressure $\mathcal{P}(r_{\perp})$ on the light front.
All of these densities, except the pressure, are normalized to unity upon
integration over the transverse plane. The pressure $\mathcal{P}$ is
normalized to the hadron mass $M$.

The layered picture is consistent with the physics of QCD. From the LFWF
representation in Sect. IV, partons with large-$x$ contribute most to the
matter density $\mathcal{A}(r_{\perp})$. Therefore, at the core of a
charmonium are the valence quarks. For the mass squared density
$\mathcal{M}^{2}(r_{\perp})$, small-$x$ partons, known as “wee partons”, also
have a significant contribution since the light-front energy
$p^{-}\propto(p^{2}_{\perp}+m^{2})/xP^{+}$. The fact that $r_{M^{2}}>r_{A}$
implies that the distribution of the wee partons is broader than the valence
quarks. The scalar trace $T^{\alpha}_{\;\;\alpha}$ contains the anomalous
gluon fields, whose radius is largest. Therefore, at the outermost of a
charmonium, there are glueballs and meson clouds. Note that the quark and the
antiquark in our model are effective degrees of freedom. Even though there is
no explicit gluons, some effects of them have been incorporated.

In the literature, it is popular to consider the decomposition of the hadron
mass into contributions from different species (quarks & gluon), with each
contribution defined from a gauge-invariant operator. Since mass in relativity
is not additive, what actually gets decomposed is the proper energy $E$. To
remove the kinetic energy contributions, each species is assumed to be
stationary. Physically, this decomposition is in analogy to the hydrostatics
of coupled multifluids [164, 165]. In principle, the invariant mass squared
can be decomposed in a similar fashion, except that there is no need to assume
that the hadron or its components are at rest since the invariant mass squared
$M^{2}$ is a Lorentz scalar. In the previous work [49], some of us proposed to
decompose the invariant mass squared $M^{2}$ into a free (or kinetic) part and
an interaction (or potential energy) part. This decomposition is more familiar
in classical and quantum mechanics although it is not gauge invariant as is
the case of the Coulomb potential. Nevertheless, it provides insights into the
strong force within charmonium, as we will show in Sect. V.

##  IV Light-front wave function representation

###  IV.1 Gravitational form factor $A$

In principle, the GFFs $A(-q^{2})$ and $D(-q^{2})$ can be extracted from any
two independent HMEs $t^{\alpha\beta}=\langle
P+\frac{1}{2}q|T^{\alpha\beta}(0)|P-\frac{1}{2}q\rangle/(2P^{+})$. We adopt
the energy-momentum density $T^{+\mu}$ to extract these form factors. These
components provide the correct value of the GFFs in the forward limit as long
as the 4-momentum is conserved [49]. Furthermore, covariant light-front
analysis shows that the hadronic matrix elements of these components are free
of spurious zero-mode contributions [49]. In the Drell-Yan frame $q^{+}=0$,

$\displaystyle t^{++}(\vec{q}_{\perp};P)=\,$ $\displaystyle
P^{+}A(q^{2}_{\perp}),$ (50) $\displaystyle t^{+i}(\vec{q}_{\perp};P)=\,$
$\displaystyle P^{i}_{\perp}A(q^{2}_{\perp}),$ (51) $\displaystyle
t^{+-}(\vec{q}_{\perp};P)=\,$
$\displaystyle\frac{P^{2}_{\perp}A(q^{2}_{\perp})+M^{2}(q^{2}_{\perp})}{P^{+}},$
(52)

where, recall,
$M^{2}(q^{2}_{\perp})=(M^{2}+\frac{1}{4}q^{2}_{\perp})A(q^{2}_{\perp})+\frac{1}{2}q_{\perp}^{2}D(q^{2}_{\perp})$.

The GFF $A$ can be extracted from either $T^{++}$ or $T^{+i}$. Indeed, the
former is the traditional “good current”. Since $T^{++}$ and $T^{+i}$ do not
contain the interaction, $A(Q^{2})$ can be represented in terms of the LFWFs
as [166],

$A(q^{2}_{\perp})=\sum_{n}\sum_{\\{s_{i}\\}}\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}k_{i\perp}\big{]}_{n}\\\
\times\sum_{j}x_{j}\psi_{n}(\\{x_{i},\vec{k}_{i\perp},s_{i}\\})\psi_{n}(\\{x_{i},\vec{k}_{i,j\perp},s_{i}\\}),$
(53)

where,

$\vec{k}_{i,j\perp}=\begin{cases}\vec{k}_{i\perp}-x_{i}\vec{q}_{\perp},&\text{spectator:
}i\neq j\\\ \vec{k}_{i\perp}+(1-x_{i})\vec{q}_{\perp},&\text{struck parton:
}i=j\\\ \end{cases}$ (54)

And the $n$-body integration measure is defined as,

$\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}k_{i\perp}\big{]}_{n}=\frac{1}{S_{n}}\prod_{i=1}^{n}\int\frac{{\mathrm{d}}x_{i}}{2x_{i}}2\delta(\sum_{i}x_{i}-1)\\\
\times\int\frac{{\mathrm{d}}^{2}k_{i\perp}}{(2\pi)^{3}}(2\pi)^{3}\delta^{2}(\sum_{i}k_{i\perp})$
(55)

where, $S_{n}$ is the symmetry factor. Using the coordinate space
representation developed in our previous work [49], $A(q^{2}_{\perp})$
becomes,

$A(q_{\perp}^{2})=\sum_{n}\sum_{\\{s_{i}\\}}\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}r_{i\perp}\big{]}_{n}\\\
\times\Big{|}\widetilde{\psi}_{n}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})\Big{|}^{2}\sum_{j}x_{j}e^{i\vec{r}_{j\perp}\cdot\vec{q}_{\perp}}\,,$
(56)

where,

$\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}r_{i\perp}\big{]}_{n}=\frac{1}{S_{n}}\prod_{i=1}^{n}\int\frac{{\mathrm{d}}x_{i}}{4\pi
x_{i}}4\pi\delta(\sum_{i}x_{i}-1)\\\
\times\int{\mathrm{d}}^{2}r_{i\perp}\delta^{2}(\sum_{i}x_{i}\vec{r}_{i\perp})\,.$
(57)

Then, the matter distribution becomes a one-body density (OBD),

$\mathcal{A}(r_{\perp})=\sum_{n}\sum_{\\{s_{i}\\}}\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}r_{i\perp}\big{]}_{n}\Big{|}\widetilde{\psi}_{n}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})\Big{|}^{2}\\\
\times\sum_{j}x_{j}\delta^{2}(\vec{r}_{\perp}-\vec{r}_{j\perp})\,.$ (58)

And the corresponding r.m.s. matter radius is,

$r^{2}_{A}\equiv-6\frac{1}{A(0)}\frac{{\mathrm{d}}A(Q^{2})}{{\mathrm{d}}Q^{2}}\Big{|}_{Q^{2}\to
0},\\\
=\frac{3}{2}\sum_{n}\sum_{\\{s_{i}\\}}\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}r_{i\perp}\big{]}_{n}\Big{|}\widetilde{\psi}_{n}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})\Big{|}^{2}\\\
\times\sum_{j}x_{j}r_{j\perp}^{2}.$ (59)

As we mentioned, the matter density mainly samples contributions from the
large-$x$ partons, i.e. valence quarks.

In this work, we only consider the valence Fock sector contribution. The LFWF
representation in this ansatz can be written as,

$A(q_{\perp}^{2})=\sum_{s,\bar{s}}\int_{0}^{1}\frac{{\mathrm{d}}x}{2x(1-x)}\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{3}}\psi_{s\bar{s}}(x,\vec{k}_{\perp})\Big{\\{}(1-x)\\\
\times\psi^{*}_{s\bar{s}}(x,\vec{k}_{\perp}-x\vec{q}_{\perp})+x\psi^{*}_{s\bar{s}}(x,\vec{k}_{\perp}+(1-x)\vec{q}_{\perp})\Big{\\}},$
(60)

Similarly, the matter density reads,

$\mathcal{A}(\vec{r}_{\perp})=\sum_{s,\bar{s}}\int_{0}^{1}\frac{{\mathrm{d}}x}{4\pi
x(1-x)}\Big{\\{}(1-x)\Big{|}\widetilde{\psi}_{s\bar{s}}\Big{(}x,\frac{\vec{r}_{\perp}}{x}\Big{)}\Big{|}^{2}\\\
+x\Big{|}\widetilde{\psi}_{s\bar{s}}\Big{(}x,-\frac{\vec{r}_{\perp}}{1-x}\Big{)}\Big{|}^{2}\Big{\\}}.$
(61)

In basis representation,

$A(q^{2}_{\perp})=\sum_{s,\bar{s},n,l,n^{\prime},l^{\prime}}\psi^{*}_{s\bar{s}}(n^{\prime},m^{\prime},l^{\prime})\psi_{s\bar{s}}(n,m,l)\int_{0}^{1}\frac{{\mathrm{d}}x}{4\pi}\\\
\times\chi_{l^{\prime}}(x)\chi_{l}(x)\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{2}}\Big{[}(1-x)\phi^{*}_{n^{\prime}m^{\prime}}\big{(}\vec{k}_{\perp}-\sqrt{\frac{x}{1-x}}\vec{q}_{\perp}\big{)}\\\
+x\phi^{*}_{n^{\prime}m^{\prime}}\big{(}\vec{k}_{\perp}+\sqrt{\frac{1-x}{x}}\vec{q}_{\perp}\big{)}\Big{]}\phi_{nm}(\vec{k}_{\perp})\,.$
(62)

Note that $m=m^{\prime}=m_{j}-s-\bar{s}$ according to angular momentum
conservation, but is otherwise arbitrary though limited by the total angular
momentum of the state. Similarly, the basis representation of the matter
density $\mathcal{A}(\vec{r}_{\perp})$ reads,

$\mathcal{A}(\vec{r}_{\perp})=\sum_{s,\bar{s},n,l,n^{\prime},l^{\prime}}\psi^{*}_{s\bar{s}}(n^{\prime},m^{\prime},l^{\prime})\psi_{s\bar{s}}(n,m,l)\int_{0}^{1}\frac{{\mathrm{d}}x}{4\pi}\\\
\times\chi_{l^{\prime}}(x)\chi_{l}(x)\Big{[}\frac{x^{2}}{1-x}\widetilde{\phi}^{*}_{n^{\prime}m^{\prime}}\Big{(}\sqrt{\frac{x}{1-x}}\vec{r}_{\perp}\Big{)}\widetilde{\phi}_{nm}\Big{(}\sqrt{\frac{x}{1-x}}\vec{r}_{\perp}\Big{)}\\\
+\frac{(1-x)^{2}}{x}\widetilde{\phi}^{*}_{n^{\prime}m^{\prime}}\Big{(}\sqrt{\frac{1-x}{x}}\vec{r}_{\perp}\Big{)}\widetilde{\phi}_{nm}\Big{(}\sqrt{\frac{1-x}{x}}\vec{r}_{\perp}\Big{)}\Big{]}.$
(63)

###  IV.2 Gravitational form factor $D$: kinetic part

The GFF $D(q^{2}_{\perp})$ and the mass squared form factor
$M^{2}(q^{2}_{\perp})$ can be extracted from either the light-front energy
density $T^{+-}$, or the transverse stress $T^{11}+T^{22}$ and $T^{12}$. From
the covariant light-front dynamics analysis [49], $T^{11}+T^{22}$ and $T^{12}$
are contaminated by spurious form factors and may violate the von Laue
condition111We recently discovered a new symmetry that may protect $T^{12}$
from the contamination of the spurious form factor [167]. It is not clear
whether this symmetry is applicable to the present system.. $T^{+-}$ is the
light-front energy density and it contains the interaction. While this
operator is more involved, it is constrained by energy conservation which
ensures the von Laue condition even after model truncation. We therefore use
$T^{+-}$ to extract GFF $D$. Our method is general enough and can be applied
to other bound-state systems. The HME of $T^{+-}$ consists of two parts:

$t^{+-}=\frac{P^{2}_{\perp}A(q^{2}_{\perp})+M^{2}(q^{2}_{\perp})}{P^{+}}\,.$
(64)

In the Breit frame $P_{\perp}=0$, only the invariant mass squared form factor
$M^{2}(q^{2}_{\perp})$ survives. GFF $D$ is related to $M^{2}(q^{2}_{\perp})$
as,

$q_{\perp}^{2}D(q_{\perp}^{2})=2M^{2}(q^{2}_{\perp})-2(M^{2}+\frac{1}{4}q^{2}_{\perp})A(q^{2}_{\perp})\,.$
(65)

In Ref. [49], we propose to split $T^{+-}$ into two parts
$T^{+-}=T^{+-}_{0}+T^{+-}_{\text{int}}$: a free part $T^{+-}_{0}$ and an
interaction part $T^{+-}_{\text{int}}$. The corresponding mass squared form
factor is also split into two pieces,
$M^{2}(q^{2}_{\perp})=M^{2}_{0}(q^{2}_{\perp})+M^{2}_{\text{int}}(q^{2}_{\perp})$.
The free energy density $T^{+-}_{0}$ is diagonal in Fock space, and the
corresponding mass squared form factor admits an exact LFWF representation,

$M^{2}_{0}(q^{2}_{\perp})=\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}k_{i\perp}\big{]}_{n}\sum_{j}\psi_{n}^{*}(\\{x_{i},\vec{k}_{i,j\perp}^{+},s_{i}\\})\\\
\times\frac{\vec{k}_{j\perp}^{2}+m_{j}^{2}-\frac{1}{4}q_{\perp}^{2}}{x_{j}}\psi_{n}(\\{x_{i},\vec{k}_{i,j\perp}^{-},s_{i}\\}),$
(66)

where,

$\displaystyle\vec{k}^{+}_{i,j\perp}=\,$
$\displaystyle\begin{cases}\vec{k}_{i\perp}+\frac{1}{2}x_{i}\vec{q}_{\perp},&\text{spectator:
}i\neq j\\\ \vec{k}_{i\perp}-\frac{1}{2}(1-x_{i})\vec{q}_{\perp},&\text{struck
parton: }i=j\\\ \end{cases}$ $\displaystyle\vec{k}^{-}_{i,j\perp}=\,$
$\displaystyle\begin{cases}\vec{k}_{i\perp}-\frac{1}{2}x_{i}\vec{q}_{\perp},&\text{spectator:
}i\neq j\\\ \vec{k}_{i\perp}+\frac{1}{2}(1-x_{i})\vec{q}_{\perp},&\text{struck
parton: }i=j\\\ \end{cases}$

The integrand is the $n$-body kinetic energy with an extra term
$-\frac{1}{4}q^{2}_{\perp}$. This term accounts for the recoil effect, which
has the same origin as the $D$-term as one can see from the free boson case
[49].

Using the coordinate space LFWFs, $M^{2}_{0}(q^{2}_{\perp})$ is dramatically
simplified as,

$M_{0}^{2}(q^{2}_{\perp})=\sum_{n}\sum_{\\{s_{i}\\}}\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}r_{i\perp}\big{]}_{n}\widetilde{\psi}_{n}^{*}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})\\\
\times\sum_{j}e^{i\vec{r}_{j\perp}\cdot\vec{q}_{\perp}}\frac{-\nabla_{j\perp}^{2}+m_{j}^{2}-\frac{1}{4}q_{\perp}^{2}}{x_{j}}\widetilde{\psi}_{n}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})\,.$
(67)

Its Fourier transform gives the OBD of the free mass squared,

$\mathcal{M}_{0}^{2}(r_{\perp})=\sum_{n}\sum_{\\{s_{i}\\}}\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}r_{i\perp}\big{]}_{n}\sum_{j}\delta^{2}(\vec{r}_{j\perp}-\vec{r}_{\perp})\\\
\times\bigg{\\{}\widetilde{\psi}_{n}^{*}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})\frac{-\nabla_{j\perp}^{2}+m_{j}^{2}}{x_{j}}\widetilde{\psi}_{n}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})\\\
+\frac{1}{4x_{j}}\nabla^{2}_{j\perp}\big{|}\widetilde{\psi}_{n}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})\big{|}^{2}\bigg{\\}}\,.$
(68)

For charmonium,

$M^{2}_{0}(q_{\perp}^{2})=\sum_{s,\bar{s}}\int\frac{{\mathrm{d}}x}{2x(1-x)}\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{3}}\Big{\\{}\\\
\psi^{*}_{s\bar{s}}(x,\vec{k}_{\perp}-{\frac{1}{2}}(1-x)\vec{q}_{\perp})\frac{k_{\perp}^{2}+m_{q}^{2}-\frac{1}{4}\vec{q}_{\perp}^{2}}{x}\\\
\times\psi_{s\bar{s}}(x,\vec{k}_{\perp}+{\frac{1}{2}}(1-x)\vec{q}_{\perp})+\psi^{*}_{s\bar{s}}(x,\vec{k}_{\perp}+{\frac{1}{2}}x\vec{q}_{\perp})\\\
\times\frac{k_{\perp}^{2}+m_{\bar{q}}^{2}-\frac{1}{4}\vec{q}_{\perp}^{2}}{1-x}\psi_{s\bar{s}}(x,\vec{k}_{\perp}-{\frac{1}{2}}x\vec{q}_{\perp})\Big{\\}}\,.$
(69)

In basis representation,

$M^{2}_{0}(q_{\perp}^{2})=\sum_{s,\bar{s},n,l,n^{\prime},l^{\prime}}\psi^{*}_{s\bar{s}}(n^{\prime},m^{\prime},l^{\prime})\psi_{s\bar{s}}(n,m,l)\int_{0}^{1}\frac{{\mathrm{d}}x}{4\pi}\\\
\times\chi_{l^{\prime}}(x)\chi_{l}(x)\Big{\\{}\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{2}}\phi^{*}_{n^{\prime}m^{\prime}}\big{(}\vec{k}_{\perp}-\frac{1}{2}\sqrt{\frac{1-x}{x}}\vec{q}_{\perp}\big{)}\\\
\times\big{[}(1-x)k_{\perp}^{2}+\frac{m_{q}^{2}-\frac{1}{4}q_{\perp}^{2}}{x}\big{]}\phi_{nm}\big{(}\vec{k}_{\perp}+\frac{1}{2}\sqrt{\frac{1-x}{x}}\vec{q}_{\perp}\big{)}\\\
+\phi^{*}_{n^{\prime}m^{\prime}}\big{(}\vec{k}_{\perp}+{\frac{1}{2}}\sqrt{\frac{x}{1-x}}\vec{q}_{\perp}\big{)}\big{[}xk_{\perp}^{2}+\frac{m_{\bar{q}}^{2}-\frac{1}{4}q_{\perp}^{2}}{1-x}\big{]}\\\
\times\phi_{nm}\big{(}\vec{k}_{\perp}-{\frac{1}{2}}\sqrt{\frac{x}{1-x}}\vec{q}_{\perp}\big{)}\Big{\\}}\,.$
(70)

###  IV.3 Gravitational form factor $D$: interaction part

The interaction part of the EMT in coordinate space can be formally
represented in Fock space using LFWFs. However, the expression is not
particularly useful in quantum field theories, since it involves particle
number changing. Furthermore, the interaction $v_{nm}$ also has to be
renormalized. The effective interactions, such as the one-gluon-exchange we
employed, are non-local in general. It may seem that a local interaction
density has to be obtained from an underlying field theory. On the other hand,
everything gravitates. As such, the EMT should exist for any system. As the
dynamics is fully determined by the Hamiltonian, the EMT including the
Hamiltonian density should be derivable from the Hamiltonian. In non-
relativistic quantum many-body theory, the Hamiltonian density can be obtained
by localizing the Hamiltonian operator, viz. inserting $\delta^{3}(r-r_{i})$
for the $i$-th constituent. For example, the charge density can be obtained by
localizing the charge,

$Q=\sum_{i}e_{i}\quad\to\quad\rho(r)=\sum_{i}e_{i}\delta^{3}(r-r_{i})\,.$ (71)

In relativistic quantum theory, particles cannot be consistently
localized222For mildly relativistic systems, Newton-Wigner operator can be
used to localize particles [168]. The application of this operator to QCD can
be found in Ref. [169].. Fortunately, there is one way to circumvent this
problem using light-front dynamics. In light-front dynamics, particles can be
localized on the transverse plane using the Dirac-$\delta$ prescription
$\delta^{2}(r_{\perp}-r_{i\perp})$, which is sufficient to extract the
transverse densities. For example, the transverse charge density can be
obtained by localizing the charge on the transverse plane,

$Q=\sum_{i}e_{i}\quad\to\quad\rho(r_{\perp})=\sum_{i}e_{i}\delta^{2}(r_{\perp}-r_{i\perp})\,.$
(72)

One can show that the charge density obtained in this way reproduces the
Drell-Yan-West formula. Similarly, the matter density can also be obtained
using particle localization on the light front.

The effective interaction we employed is a two-body potential,

$V=\frac{1}{2!}\sum_{i,j}v(x,r_{i\perp}-r_{j\perp}),$ (73)

where, $\vec{r}_{i\perp}$ is the transverse coordinate of the $i$-th parton.
Based on the localization prescription, we can write down the corresponding
energy density as,

$\mathcal{V}(r_{\perp})=\frac{1}{2!}\sum_{i,j}v(x,r_{i\perp}-r_{j\perp})\\\
\times\frac{1}{2}\big{\\{}\delta^{2}(r_{\perp}-r_{i\perp})+\delta^{2}(r_{\perp}-r_{j\perp})\big{\\}}\,.$
(74)

Here, the inter-particle energy density is evenly distributed between the two
interacting constituents. The above many-body expression can be converted into
the corresponding operator following the standard second quantization rule.

To obtain the inter-particle energy $v(x,r_{i\perp}-r_{j\perp})$, we first
observe that the interaction energy in the forward limit admits a diagonal
representation,

$\displaystyle M^{2}_{\text{int}}(0)=\,$
$\displaystyle\sum_{n}\sum_{\\{s_{i}\\}}\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}k_{i\perp}\big{]}_{n}\psi_{n}^{*}(\\{x_{i},\vec{k}_{i\perp},s_{i}\\})$
(75)
$\displaystyle\times\psi_{n}(\\{x_{i},\vec{k}_{i\perp},s_{i}\\})v_{n}(\\{x_{i},\vec{k}_{i\perp}\\}),$
(76) $\displaystyle=\,$
$\displaystyle\sum_{n}\sum_{\\{s_{i}\\}}\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}r_{i\perp}\big{]}_{n}\widetilde{\psi}_{n}^{*}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})$
(77) $\displaystyle\times
v_{n}(\\{x_{i},-i\vec{\nabla}_{i\perp}\\})\widetilde{\psi}_{n}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\}),$
(78)

where, the $n$-body potential energy is expressed in terms of the mass
eigenvalue and the kinetic energy

$v_{n}(\\{x_{i},\vec{k}_{i\perp}\\})=M^{2}-\sum_{j=1}^{n}\frac{\vec{k}_{j\perp}^{2}+m_{j}^{2}}{x_{j}}.$
(79)

This representation is exact and can be obtained directly from the Schrödinger
equation (2).

To generalize the above expression to the off-forward case, we need to specify
the locations of the interactions. We adopt the impulse ansatz that all
interactions happen at the same instant in light-front time. This ansatz is
certainly held for local interactions. In Ref. [49], we rigorously showed that
this ansatz is valid in a pion cloud model where all scattering appears at the
location of the nucleon. For small-size systems, such as charmonium, this
ansatz is expected to be a good approximation. For effective interactions,
this ansatz requires neglecting the propagation of the intermediate particles,
which is widely adopted in non-relativistic systems. Since we do not
incorporate dynamical gluons and sea quarks in this work, this ansatz is
consistent with our approximation.

The interaction form factor in this ansatz becomes,

$M^{2}_{\text{int}}(q^{2}_{\perp})=\sum_{n}\frac{1}{n}\sum_{\\{s_{i}\\}}\int\big{[}{\mathrm{d}}x_{i}{\mathrm{d}}^{2}r_{i\perp}\big{]}_{n}\widetilde{\psi}_{n}^{*}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\})\\\
\times\sum_{j}e^{i\vec{q}_{\perp}\cdot\vec{r}_{j\perp}}v_{n}(\\{x_{i},-i\vec{\nabla}_{i\perp}\\})\widetilde{\psi}_{n}(\\{x_{i},\vec{r}_{i\perp},s_{i}\\}),$
(80)

assuming the interaction is two-body. For charmonium, the above expression
becomes,

$M^{2}_{\text{int}}(q^{2}_{\perp})=\frac{1}{2}\sum_{s,\bar{s}}\int\frac{{\mathrm{d}}x}{4\pi
x(1-x)}\int{\mathrm{d}}^{2}r_{\perp}\widetilde{\psi}^{*}_{s\bar{s}}(x,\vec{r}_{\perp})\\\
\times\Big{[}e^{i\vec{q}_{\perp}\cdot\vec{r}_{1\perp}}+e^{i\vec{q}_{\perp}\cdot\vec{r}_{2\perp}}\Big{]}v(x,\vec{r}_{\perp},-i\nabla_{\perp})\widetilde{\psi}_{s\bar{s}}(x,\vec{r}_{\perp}),$
(81)

where,
$\vec{r}_{1\perp}=\vec{R}_{\perp}+(1-x)\vec{r}_{\perp}=(1-x)\vec{r}_{\perp}$,
$r_{2\perp}=\vec{R}_{\perp}-x_{1}\vec{r}_{\perp}=-x\vec{r}_{\perp}$, and

$v(x,\vec{r}_{\perp},-i\nabla_{\perp})=M^{2}-\frac{-\nabla_{\perp}^{2}+m_{q}^{2}}{x}-\frac{-\nabla_{\perp}^{2}+m^{2}_{\bar{q}}}{1-x}.$
(82)

In terms of the momentum-space LFWFs, the above expression can also be written
as,

$M^{2}_{\text{int}}(q^{2}_{\perp})=\frac{1}{2}\sum_{s,\bar{s}}\int\frac{{\mathrm{d}}x}{2x(1-x)}\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{3}}\\\
\times\Big{[}\psi^{*}_{s\bar{s}}\big{(}x,\vec{k}_{\perp}+(1-x)\vec{q}_{\perp}\big{)}+\psi^{*}_{s\bar{s}}(x,\vec{k}_{\perp}-x\vec{q}_{\perp})\Big{]}\\\
\times\Big{[}M^{2}-\frac{k_{\perp}^{2}+m_{q}^{2}}{x}-\frac{k_{\perp}^{2}+m^{2}_{\bar{q}}}{1-x}\Big{]}\psi_{s\bar{s}}(x,\vec{k}_{\perp})\,.$
(83)

In basis representation,

$M^{2}_{\text{int}}(q^{2}_{\perp})=\frac{1}{2}\sum_{s,\bar{s},n,l,n^{\prime},l^{\prime}}\psi^{*}_{s\bar{s}}(n^{\prime},m^{\prime},l^{\prime})\psi_{s\bar{s}}(n,m,l)\int_{0}^{1}\frac{{\mathrm{d}}x}{4\pi}\\\
\times\chi_{l^{\prime}}(x)\chi_{l}(x)\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{2}}\phi_{nm}\big{(}\vec{k}_{\perp}\big{)}\Big{(}M^{2}-k_{\perp}^{2}-\frac{m_{c}^{2}}{x(1-x)}\Big{)}\\\
\times\Big{[}\phi^{*}_{n^{\prime}m^{\prime}}\big{(}\vec{k}_{\perp}+\sqrt{\frac{1-x}{x}}\vec{q}_{\perp}\big{)}+\phi^{*}_{n^{\prime}m^{\prime}}\big{(}\vec{k}_{\perp}-\sqrt{\frac{x}{1-x}}\vec{q}_{\perp}\big{)}\Big{]}\,.$
(84)

In the forward limit ($q^{2}=0$), it is evident that

$M^{2}_{0}(0)+M^{2}_{\text{int}}(0)=M^{2},$ (85)

from the LFWF representation. From (65), this condition implies

$\lim_{q^{2}_{\perp}\to 0}q^{2}_{\perp}D(q^{2}_{\perp})=0.$ (86)

Therefore, the von Laue condition is satisfied.

##  V Numerical results

We substitute the charmonium LFWFs obtained from BLFQ [156]. As mentioned, the
form factors are evaluated with $N_{\text{max}}=L_{\text{max}}=8$. The
uncertainty associated with the basis resolution is quoted as the difference
between the $N_{\text{max}}=L_{\text{max}}=8$ results and the
$N_{\text{max}}=L_{\text{max}}=16$ results.

Figure 2: GFFs $A(Q^{2})$ (top) and $D(Q^{2})$ (bottom) of ground-state pseudoscalar charmonium $\eta_{c}$. The central value is based on $N_{\text{max}}=8$ wave functions and the difference between $N_{\text{max}}=8$ and $N_{\text{max}}=16$ results is shown as the bands to indicate the basis sensitivity. Table 1: BLFQ prediction of the $D$-term and the radii of charmonium $\eta_{c}$, $\eta_{c}^{\prime}$ and $\chi_{c0}$. We adopt the $N_{\text{max}}=8$ results as the central values, which corresponds to a UV resolution about the same as the hadron mass. In parentheses we quote the difference between $N_{\text{max}}=8$ and $N_{\text{max}}=16$ results as the basis sensitivity.  |  $D$ |  $r_{A}$ (fm) |  $r_{E}$ (fm) |  $r_{M^{2}}$ (fm) |  $r_{\theta}$ (fm) |  ${2\langle T\rangle}/{\langle V\rangle}$
---|---|---|---|---|---|---
$\eta_{c}$ |  $-5.05(8)$ |  $0.18(1)$ |  $0.26(3)$ |  0.33(3) |  0.39(1) |  $-1.7(6)$
$\eta_{c}^{\prime}$ |  $-7.1(3)$ |  $0.35(2)$ |  $0.40(4)$ |  0.452(1) |  0.50(2) |  $-14(16)$
$\chi_{c0}$ |  $-6.95(2)$ |  $0.24(2)$ |  $0.29(2)$ |  0.39(1) |  0.44(1) |  $-3.0(1.1)$

Figure 2 shows the GFFs $A(Q^{2})$, $D(Q^{2})$ for the ground-state
pseudoscalar charmonium $\eta_{c}$. From these results, we can extract $D$ and
radius $r_{A}$. The results are listed in Table 1. Also listed in Table 1 are
the radii, $r_{E}$, $r_{M^{2}}$, $r_{\theta}$, as combinations of these two
primary quantities. A somewhat surprising result is the large contribution
from terms proportional to $\lambda_{C}=1/M=0.066\,\mathrm{fm}$. Naïvely, for
quarkonium, this term is expected to be small. However, the matter radius is
roughly proportional to the inverse of the binding energy $r_{A}\sim
B^{-1}\sim\alpha^{-1}_{s}\lambda_{C}$, where the strong coupling constant
$\alpha_{s}\sim 0.3$ for charmonium. Therefore, $r_{A}$ and $\lambda_{C}$ are
actually of the same order of magnitude.

Figure 3 compares the GFFs for $\eta_{c}$ with those for its radial and
angular excitations, viz. $\eta^{\prime}_{c}$ and $\chi_{c0}$. Radial and
angular nodes appear for excited states as expected. The corresponding
$D$-term and the radii are extracted and collected in Table 1. From these
numbers, the matter radius $r_{A}$ of the $P$-wave charmonium $\chi_{c0}$ is
larger than that of $\eta_{c}$ by 30%. The radius $r_{A}$ of the
$\eta_{c}^{\prime}$, the radial excitation, almost doubles that of the ground
state.

Figure 3: Comparison of GFFs $A(Q^{2})$ (top) and $D(Q^{2})$ (bottom) of
charmonia $\eta_{c}$, $\eta_{c}^{\prime}$ and $\chi_{c0}$. The uncertainty
bands are determined as described in the caption to Fig. 2.

From these form factors, we can extract physical densities. One of the
advantages of the BLFQ approach is that the coordinate-space density can be
obtained in closed form. This avoids the problem of numeric fitting or
numerical Fourier transform, which introduces large and uncontrollable
numerical errors. The extracted pressure is shown in Fig. 4 for $\eta_{c}$,
$\eta^{\prime}_{c}$ [i.e. $\eta_{c}(2S)$] and $\chi_{c0}$. The total pressure
vanishes as a result of the von Laue condition. For all three states, the
pressure is positive (repulsive) at the center of the meson and negative
(attractive) near the edge, satisfying the stability conjecture [56]. The
pressure profile of the radial excited charmonium $\eta^{\prime}_{c}$ shows
one additional node. Note that the uncertainty of the pressure at small
$r_{\perp}$ is large. This is because our basis has a finite UV coverage. The
effective UV cutoff for $N_{\text{max}}=8$ is
$\Lambda_{\text{UV}}=\kappa\sqrt{N_{\text{max}}}\approx 2.8\,\mathrm{GeV}$,
corresponding to a spatial resolution of 0.07 fm. Smaller than this distance,
the pressure starts to show large uncertainty, i.e. stronger basis dependence.

Figure 4: Comparison of the pressure within $\eta_{c}$ (top),
$\eta_{c}^{\prime}$ (middle) and $\chi_{c0}$ (bottom). The uncertainty bands
are determined as described in the caption to Fig. 2.

We have mentioned the difference between the proper energy distribution
$\mathcal{E}(r_{\perp})$ defined from the instant form and the invariant mass
squared distribution $\mathcal{M}^{2}(r_{\perp})$ obtained on the light front.
The former is normalized to the total energy of a hadron at rest, i.e. its
rest mass $M$, whereas the latter is normalized to the invariant mass squared
$M^{2}$ of the hadron. Fig. 5 compares the normalized energy distribution of
$\eta_{c}$ and its normalized invariant mass squared distribution. The energy
distribution is positive. However, the $M^{2}$ distribution becomes negative
at short distance $r_{\perp}$. This tachyonic core within charmonium reveals a
tantalizing novel structure of the QCD bound state. Similar quantities can be
extracted and compared for $\eta^{\prime}_{c}$ and $\chi_{c0}$, as shown in
Fig. 6. Within basis uncertainty, the energy distributions
$\mathcal{E}(r_{\perp})$ are consistent with the weak energy condition that
the energy densities are positive.

Figure 5: (Top) Comparison of the normalized energy distribution and the
normalized distribution of the invariant mass squared within $\eta_{c}$. Note
that the invariant mass squared density is negative at the center of
$\eta_{c}$, shown in the bottom panel. This tachyonic core suggests a
tantalizing picture of the strong force. The uncertainty bands are determined
as described in the caption to Fig. 2.

Figure 6: Comparison of the energy distributions (top) and the distribution of
the invariant mass squared (bottom) of $\eta_{c}^{\prime}$ and $\chi_{c0}$.
The uncertainty bands are determined as described in the caption to Fig. 2.

We have so far introduced several densities: the matter density
$\mathcal{A}(r_{\perp})$, energy density $\mathcal{E}(r_{\perp})$, invariant
mass squared density $\mathcal{M}^{2}(r_{\perp})$ and scalar density
$\theta(r_{\perp})$. We have argued that a negative $D$ suggests a hierarchy
of the radii associated with these densities:

$r_{A}<r_{E}<r_{M^{2}}<r_{\theta}\,.$ (87)

Figure 7 compares these densities for $\eta_{c}$ after the normalization, i.e.
all densities are normalized to unity. Both the invariant mass squared density
$\mathcal{M}^{2}(r_{\perp})$ and the scalar density $\theta(r_{\perp})$ become
negative at the center of the meson. The emergence of a tachyonic core is
consistent with a relatively large radius.

Figure 7: Comparison of normalized densities of $\eta_{c}$. The densities are
computed with LFWFs with $N_{\text{max}}=8$ (see texts). The corresponding
radii are indicated by arrows.

Figure 8: (Top) Decomposition of the invariant mass squared density of
$\eta_{c}$ into the free part and the interacting part. (Middle) Comparison of
the free invariant mass squared densities for $\eta_{c}$, $\eta_{c}^{\prime}$
and $\chi_{c0}$. (Bottom) Comparison of the interacting invariant mass squared
densities for $\eta_{c}$, $\eta_{c}^{\prime}$ and $\chi_{c0}$. The uncertainty
bands are determined as described in the caption to Fig. 2.

As we mentioned, in the light-front Hamiltonian formalism, we are also able to
access the light-front kinetic energy and light-front potential energy
separately. The distributions of the free invariant mass squared as well as
the interacting invariant mass squared for $\eta_{c}$, $\eta_{c}^{\prime}$ and
$\chi_{c0}$ are shown in Fig. 8. Note that $\mathcal{M}_{0}^{2}(r_{\perp})$
does not remain positive for all distance for $\eta_{c}^{\prime}$ and
$\chi_{c0}$. This is caused by the recoil term $-{q^{2}_{\perp}}/{4x}$ in
(67). The interaction density $\mathcal{M}_{\text{int}}^{2}(r_{\perp})$ of
$\eta_{c}^{\prime}$ becomes positive at some large distance, consistent with
the positive confining potential $r^{2}_{\perp}$ dominant at large parton
separation.

Non-relativistically, for particles within a potential of the form
$V(r)=\alpha r^{n}$, the virial theorem333Lorcé et. al. generalized the virial
theorem to QFT and showed that it is equivalent to von Laue condition [165].
suggests, $n={2\langle T\rangle}/{\langle V\rangle}$. Inspired by this result,
we introduce a virial scaling index

$n\equiv\frac{2\langle T\rangle}{\langle V\rangle},$ (88)

where, we define the kinetic energy, i.e. the virial, as the free energy
$M^{2}_{0}\equiv M^{2}_{0}(Q^{2}=0)$ subtracting the square of the total quark
mass,

$\langle T\rangle\equiv M^{2}_{0}(Q^{2}=0)-(m_{q}+m_{\bar{q}})^{2}\,.$ (89)

The potential energy is simply the interacting energy,

$\langle V\rangle\equiv M^{2}_{\text{int}}(Q^{2}=0)\,.$ (90)

This index $n$ characterizes the shape of the inter-particle force under
spatial dilation. In non-relativistic limit, our effective interaction is a
harmonic oscillator $r^{2}$ at large distance plus a Coulomb $r^{-1}$ at short
distance. One would expect a virial index close to $-1$ for deep bound state,
e.g. $\eta_{c}$, and $+1$ for bound state with large size [170]. The obtained
indices are collected Table 1. Within the sizable basis uncertainty, these
numbers are in agreement with the prediction from the non-relativistic virial
theorem. Fig. 9 shows the decomposition of the invariant mass squared of
charmonia $\eta_{c}$, $\eta^{\prime}_{c}$ and $\chi_{c0}$ into quark mass
contribution, kinetic energy contribution and the potential energy
contribution.

Figure 9: Decomposition of charmonium invariant mass squared in terms of the
total quark mass squared $4m_{c}^{2}$, kinetic mass squared $\langle T\rangle$
and potential mass squared $\langle V\rangle$. See texts.

##  VI Summary and outlook

In this work, we investigate the gravitational form factors of the charmonium
system. Starting from the light-front wave functions obtained previously in
basis light-front quantization, we construct the wave function representation
of GFFs $A$ and $D$. The former is extracted from the operator $T^{++}$ while
the latter is from $T^{+-}$ with an impulse ansatz. The obtained GFFs satisfy
the known global constraints such as $A(0)=1$, and the von Laue condition.

From these two primary quantities, we also define physical densities on the
light front, including the pressure $\mathcal{P}(r_{\perp})$, the (proper)
energy density $\mathcal{E}(r_{\perp})$, invariant mass squared density
$\mathcal{M}^{2}(r_{\perp})$ and trace scalar density $\theta(r_{\perp})$. We
also identify the positive-defined density $\mathcal{A}(r_{\perp})$, the
Fourier transform of GFF $A(Q^{2})$, as the (convective) matter density. The
energy density $\mathcal{E}(r_{\perp})$ is also positive and is consistent
with the weak energy condition. These densities provide rich information of
the system. For example, we find that there is a hierarchy for the sizes of
the system: $r_{A}<r_{E}<r_{M^{2}}<r_{\theta}$, implying an onion-like
structure of hadrons. We also show that there is a tachyonic core inside
hadrons. Finally, we investigate the virial scaling index $2\langle
T\rangle/\langle V\rangle$, which provides a quantitative measure of the
strong force within charmonium. The obtained values are consistent with the
non-relativistic picture.

We note that the method we proposed in this work is general enough and is
applicable to states with a different spin such as $J/\psi$, and to other
systems, such as the nucleon. The same method was applied to the pion in the
context of holographic light-front QCD [52]. It is interesting to observe the
cancellation between the scalar and tensor glueballs and the emergence of
scalar meson dominance within the pion. It is interesting to note that the
scalar meson coupling to $D$ comes from the recoil term $-q^{2}_{\perp}/4x$. A
closely related set of observables are the gravitational transition form
factors [171, 172, 173, 174, 175, 176, 177, 178, 179], which describe hadron
production in gravity. These observables can also be accessed using the
methods proposed in this work.

One of the interacting questions is the comparison between GFF $D$ extracted
from $T^{+-}$ and from $T^{12}$ and from $T^{11}+T^{22}$. Operators
$T^{11},T^{22}$ also contain interactions and need to be renormalized.
However, for practical calculations on the light front, including perturbative
or non-perturbative, counterterms computed from the light-front Hamiltonian
are not sufficient to renormalize these operators. This is in sharp contrast
to $T^{+-}$, which shares the same counterterms as its conserved charge, the
light-front Hamiltonian. $T^{12}$ is interaction free. It is tempting to
extract GFF $D$ from $T^{12}$ [55]. However, covariant light-front dynamics
analysis shows that $T^{11},T^{22},T^{12}$ are all associated with a spurious
GFF – a Lorentz structure that breaks the Poincaré symmetry. The GFF $D$
extracted from $T^{12}$ is likely to differ from what we obtained from
$T^{+-}$, which does not suffer from the contamination of the spurious GFFs.

Our model does not contain dynamical gluons or sea quarks. Therefore, we are
not able to perform decomposition in terms of quarks and gluons. Incorporating
dynamical gluons and sea quarks are necessary next steps in basis light-front
quantization [138, 139, 140, 141, 142]. How to tame the computational
complexity as well as the non-perturbative renormalization arising in these
problems are critical challenges to be tackled.

## Acknowledgements

The authors acknowledge fruitful discussions with C. Mondal, S. Nair, K.
Serafin and J. More.

This work was supported in part by the National Natural Science Foundation of
China (NSFC) under Grant No. 12375081, by the Chinese Academy of Sciences
under Grant No. YSBR-101, and by the US Department of Energy (DOE) under Grant
No. DE-SC0023692. Y.L. is supported by the New faculty start-up fund of the
University of Science and Technology of China. 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. ZDBS-LY-7020, by the Natural Science Foundation of Gansu
Province, China, Grant No. 20JR10RA067, by the Foundation for Key Talents of
Gansu Province, by the Central Funds Guiding the Local Science and Technology
Development of Gansu Province, Grant No. 22ZY1QA006, by Gansu International
Collaboration and Talents Recruitment Base of Particle Physics (2023-2027), by
International Partnership Program of the Chinese Academy of Sciences, Grant
No. 016GJHZ2022103FN, by National Natural Science Foundation of China, Grant
No. 12375143, by National Key R&D Program of China, Grant No. 2023YFA1606903
and by the Strategic Priority Research Program of the Chinese Academy of
Sciences, Grant No. XDB34000000.

## References

  * [1] F. Gross, E. Klempt, S. J. Brodsky, A. J. Buras, V. D. Burkert, G. Heinrich, K. Jakobs, C. A. Meyer, K. Orginos and M. Strickland, et al. Eur. Phys. J. C 83, 1125 (2023) doi:10.1140/epjc/s10052-023-11949-2 [arXiv:2212.11107 [hep-ph]].
  * [2] H. Pagels, Phys. Rev. 144, 1250-1260 (1966) doi:10.1103/PhysRev.144.1250
  * [3] I. Y. Kobzarev and L. B. Okun, Zh. Eksp. Teor. Fiz. 43, 1904-1909 (1962)
  * [4] J. Y. Kim, Phys. Lett. B 834, 137442 (2022) doi:10.1016/j.physletb.2022.137442 [arXiv:2206.10202 [hep-ph]].
  * [5] S. Kumano, Q. T. Song and O. V. Teryaev, Phys. Rev. D 97, no.1, 014020 (2018) doi:10.1103/PhysRevD.97.014020 [arXiv:1711.08088 [hep-ph]].
  * [6] V. D. Burkert, L. Elouadrhiri and F. X. Girod, Nature 557, no.7705, 396-399 (2018) doi:10.1038/s41586-018-0060-z
  * [7] V. D. Burkert, L. Elouadrhiri and F. X. Girod, [arXiv:2104.02031 [nucl-ex]].
  * [8] K. Kumerički, Nature 570, no.7759, E1-E2 (2019) doi:10.1038/s41586-019-1211-6
  * [9] H. Dutrieux, C. Lorcé, H. Moutarde, P. Sznajder, A. Trawiński and J. Wagner, Eur. Phys. J. C 81, no.4, 300 (2021) doi:10.1140/epjc/s10052-021-09069-w [arXiv:2101.03855 [hep-ph]].
  * [10] V. D. Burkert, L. Elouadrhiri and F. X. Girod, [arXiv:2310.11568 [hep-ph]].
  * [11] B. Duran, Z. E. Meziani, S. Joosten, M. K. Jones, S. Prasad, C. Peng, W. Armstrong, H. Atac, E. Chudakov and H. Bhatt, et al. Nature 615, no.7954, 813-816 (2023) doi:10.1038/s41586-023-05730-4 [arXiv:2207.05212 [nucl-ex]].
  * [12] A. Accardi, J. L. Albacete, M. Anselmino, N. Armesto, E. C. Aschenauer, A. Bacchetta, D. Boer, W. K. Brooks, T. Burton and N. B. Chang, et al. Eur. Phys. J. A 52, no.9, 268 (2016) doi:10.1140/epja/i2016-16268-9 [arXiv:1212.1701 [nucl-ex]].
  * [13] D. P. Anderle, V. Bertone, X. Cao, L. Chang, N. Chang, G. Chen, X. Chen, Z. Chen, Z. Cui and L. Dai, et al. Front. Phys. (Beijing) 16, no.6, 64701 (2021) doi:10.1007/s11467-021-1062-0 [arXiv:2102.09222 [nucl-ex]].
  * [14] A. Accardi, P. Achenbach, D. Adhikari, A. Afanasev, C. S. Akondi, N. Akopov, M. Albaladejo, H. Albataineh, M. Albrecht and B. Almeida-Zamora, et al. [arXiv:2306.09360 [nucl-ex]].
  * [15] D. C. Hackett, D. A. Pefkou and P. E. Shanahan, [arXiv:2310.08484 [hep-lat]].
  * [16] D. C. Hackett, P. R. Oare, D. A. Pefkou and P. E. Shanahan, Phys. Rev. D 108, no.11, 114504 (2023) doi:10.1103/PhysRevD.108.114504 [arXiv:2307.11707 [hep-lat]].
  * [17] Y. B. Yang, J. Liang, Y. J. Bi, Y. Chen, T. Draper, K. F. Liu and Z. Liu, Phys. Rev. Lett. 121, no.21, 212001 (2018) doi:10.1103/PhysRevLett.121.212001 [arXiv:1808.08677 [hep-lat]].
  * [18] P. E. Shanahan and W. Detmold, Phys. Rev. D 99, no.1, 014511 (2019) doi:10.1103/PhysRevD.99.014511 [arXiv:1810.04626 [hep-lat]].
  * [19] P. E. Shanahan and W. Detmold, Phys. Rev. Lett. 122, no.7, 072003 (2019) doi:10.1103/PhysRevLett.122.072003 [arXiv:1810.07589 [nucl-th]].
  * [20] D. A. Pefkou, D. C. Hackett and P. E. Shanahan, Phys. Rev. D 105, no.5, 054509 (2022) doi:10.1103/PhysRevD.105.054509 [arXiv:2107.10368 [hep-lat]].
  * [21] M. V. Polyakov, Phys. Lett. B 555, 57-62 (2003) doi:10.1016/S0370-2693(03)00036-4 [arXiv:hep-ph/0210165 [hep-ph]].
  * [22] C. Cebulla, K. Goeke, J. Ossmann and P. Schweitzer, Nucl. Phys. A 794, 87-114 (2007) doi:10.1016/j.nuclphysa.2007.08.004 [arXiv:hep-ph/0703025 [hep-ph]].
  * [23] M. Mai and P. Schweitzer, Phys. Rev. D 86, 096002 (2012) doi:10.1103/PhysRevD.86.096002 [arXiv:1206.2930 [hep-ph]].
  * [24] M. Mai and P. Schweitzer, Phys. Rev. D 86, 076001 (2012) doi:10.1103/PhysRevD.86.076001 [arXiv:1206.2632 [hep-ph]].
  * [25] J. H. Jung, U. Yakhshiev and H. C. Kim, J. Phys. G 41, 055107 (2014) doi:10.1088/0954-3899/41/5/055107 [arXiv:1310.8064 [hep-ph]].
  * [26] M. Cantara, M. Mai and P. Schweitzer, Nucl. Phys. A 953, 1-20 (2016) doi:10.1016/j.nuclphysa.2016.04.032 [arXiv:1510.08015 [hep-ph]].
  * [27] J. Hudson and P. Schweitzer, Phys. Rev. D 97, no.5, 056003 (2018) doi:10.1103/PhysRevD.97.056003 [arXiv:1712.05317 [hep-ph]].
  * [28] M. V. Polyakov and H. D. Son, JHEP 09, 156 (2018) doi:10.1007/JHEP09(2018)156 [arXiv:1808.00155 [hep-ph]].
  * [29] I. V. Anikin, Phys. Rev. D 99, no.9, 094026 (2019) doi:10.1103/PhysRevD.99.094026 [arXiv:1902.00094 [hep-ph]].
  * [30] U. Özdem and K. Azizi, Phys. Rev. D 101, no.5, 054031 (2020) doi:10.1103/PhysRevD.101.054031 [arXiv:1912.06375 [hep-ph]].
  * [31] M. Varma and P. Schweitzer, Phys. Rev. D 102, no.1, 014047 (2020) doi:10.1103/PhysRevD.102.014047 [arXiv:2006.06602 [hep-ph]].
  * [32] D. Chakrabarti, C. Mondal, A. Mukherjee, S. Nair and X. Zhao, Phys. Rev. D 102, 113011 (2020) doi:10.1103/PhysRevD.102.113011 [arXiv:2010.04215 [hep-ph]].
  * [33] A. Metz, B. Pasquini and S. Rodini, Phys. Lett. B 820, 136501 (2021) doi:10.1016/j.physletb.2021.136501 [arXiv:2104.04207 [hep-ph]].
  * [34] J. More, A. Mukherjee, S. Nair and S. Saha, Phys. Rev. D 105, no.5, 056017 (2022) doi:10.1103/PhysRevD.105.056017 [arXiv:2112.06550 [hep-ph]].
  * [35] X. B. Tong, J. P. Ma and F. Yuan, JHEP 10, 046 (2022) doi:10.1007/JHEP10(2022)046 [arXiv:2203.13493 [hep-ph]].
  * [36] K. A. Mamo and I. Zahed, Phys. Rev. D 106, no.8, 086004 (2022) doi:10.1103/PhysRevD.106.086004 [arXiv:2204.08857 [hep-ph]].
  * [37] M. Fujita, Y. Hatta, S. Sugimoto and T. Ueda, PTEP 2022, no.9, 093B06 (2022) doi:10.1093/ptep/ptac110 [arXiv:2206.06578 [hep-th]].
  * [38] Z. Xing, M. Ding and L. Chang, Phys. Rev. D 107, no.3, L031502 (2023) doi:10.1103/PhysRevD.107.L031502 [arXiv:2211.06635 [hep-ph]].
  * [39] U. Özdem and K. Azizi, JHEP 03, 048 (2023) doi:10.1007/JHEP03(2023)048 [arXiv:2212.07290 [hep-ph]].
  * [40] A. Freese, A. Metz, B. Pasquini and S. Rodini, Phys. Lett. B 839, 137768 (2023) doi:10.1016/j.physletb.2023.137768 [arXiv:2212.12197 [hep-ph]].
  * [41] J. More, A. Mukherjee, S. Nair and S. Saha, Phys. Rev. D 107, no.11, 116005 (2023) doi:10.1103/PhysRevD.107.116005 [arXiv:2302.11906 [hep-ph]].
  * [42] B. Pasquini et al. [MAP (Multi-dimensional Analyses of Partonic distributions)], Phys. Rev. D 107, no.11, 114023 (2023) doi:10.1103/PhysRevD.107.114023 [arXiv:2303.01789 [hep-ph]].
  * [43] A. Garcia Martin-Caro, M. Huidobro and Y. Hatta, Phys. Rev. D 108, no.3, 034014 (2023) doi:10.1103/PhysRevD.108.034014 [arXiv:2304.05994 [nucl-th]].
  * [44] D. Chakrabarti, P. Choudhary, B. Gurjar, R. Kishore, T. Maji, C. Mondal and A. Mukherjee, Phys. Rev. D 108, no.1, 014009 (2023) doi:10.1103/PhysRevD.108.014009 [arXiv:2304.09908 [hep-ph]].
  * [45] A. Amor-Quiroz, W. Focillon, C. Lorcé and S. Rodini, Eur. Phys. J. C 83, no.11, 1012 (2023) doi:10.1140/epjc/s10052-023-12190-7 [arXiv:2304.10339 [hep-ph]].
  * [46] Y. Guo, X. Ji, Y. Liu and J. Yang, Phys. Rev. D 108, no.3, 034003 (2023) doi:10.1103/PhysRevD.108.034003 [arXiv:2305.06992 [hep-ph]].
  * [47] W. Kou and X. Chen, Phys. Lett. B 846, 138199 (2023) doi:10.1016/j.physletb.2023.138199 [arXiv:2305.19078 [hep-ph]].
  * [48] Z. Dehghan, K. Azizi and U. Özdem, Phys. Rev. D 108, no.9, 094037 (2023) doi:10.1103/PhysRevD.108.094037 [arXiv:2307.14880 [hep-ph]].
  * [49] X. Cao, Y. Li and J. P. Vary, Phys. Rev. D 108, no.5, 056026 (2023) doi:10.1103/PhysRevD.108.056026 [arXiv:2308.06812 [hep-ph]].
  * [50] H. Y. Won, H. C. Kim and J. Y. Kim, [arXiv:2310.04670 [hep-ph]].
  * [51] A. F. Krutov and V. E. Troitsky, Phys. Rev. D 108, no.9, 094043 (2023) doi:10.1103/PhysRevD.108.094043 [arXiv:2310.14287 [hep-ph]].
  * [52] Y. Li and J. P. Vary, Phys. Rev. D 109, no.5, L051501 (2024) doi:10.1103/PhysRevD.109.L051501 [arXiv:2312.02543 [hep-th]].
  * [53] H. Alharazin, Phys. Rev. D 109, no.1, 016009 (2024) doi:10.1103/PhysRevD.109.016009 [arXiv:2312.09675 [hep-ph]].
  * [54] Y. Hagiwara, X. B. Tong and B. W. Xiao, [arXiv:2401.12840 [hep-ph]].
  * [55] S. Nair, C. Mondal, S. Xu, X. Zhao, A. Mukherjee and J. P. Vary, [arXiv:2403.11702 [hep-ph]].
  * [56] M. V. Polyakov and P. Schweitzer, Int. J. Mod. Phys. A 33, no.26, 1830025 (2018) doi:10.1142/S0217751X18300259 [arXiv:1805.06596 [hep-ph]].
  * [57] V. D. Burkert, L. Elouadrhiri, F. X. Girod, C. Lorcé, P. Schweitzer and P. E. Shanahan, Rev. Mod. Phys. 95, no.4, 041002 (2023) doi:10.1103/RevModPhys.95.041002 [arXiv:2303.08347 [hep-ph]].
  * [58] N. Brambilla, S. Eidelman, B. K. Heltsley, R. Vogt, G. T. Bodwin, E. Eichten, A. D. Frawley, A. B. Meyer, R. E. Mitchell and V. Papadimitriou, et al. Eur. Phys. J. C 71, 1534 (2011) doi:10.1140/epjc/s10052-010-1534-9 [arXiv:1010.5827 [hep-ph]].
  * [59] 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, Phys. Rev. C 81, 035205 (2010) doi:10.1103/PhysRevC.81.035205 [arXiv:0905.1411 [nucl-th]].
  * [60] J. M. Namyslowski, Prog. Part. Nucl. Phys. 14, 49-104 (1985) doi:10.1016/0146-6410(85)90051-1
  * [61] W. M. Zhang, Chin. J. Phys. 32, 717-808 (1994) [arXiv:hep-ph/9412244 [hep-ph]].
  * [62] M. Burkardt, Adv. Nucl. Phys. 23, 1-74 (1996) doi:10.1007/0-306-47067-5_1 [arXiv:hep-ph/9505259 [hep-ph]].
  * [63] S. J. Brodsky, H. C. Pauli and S. S. Pinsky, Phys. Rept. 301, 299-486 (1998) doi:10.1016/S0370-1573(97)00089-6 [arXiv:hep-ph/9705477 [hep-ph]].
  * [64] J. Carbonell, B. Desplanques, V. A. Karmanov and J. F. Mathiot, Phys. Rept. 300, 215-347 (1998) doi:10.1016/S0370-1573(97)00090-2 [arXiv:nucl-th/9804029 [nucl-th]].
  * [65] T. Heinzl, Lect. Notes Phys. 572, 55-142 (2001) doi:10.1007/3-540-45114-5_2 [arXiv:hep-th/0008096 [hep-th]].
  * [66] G. A. Miller, Prog. Part. Nucl. Phys. 45, 83-155 (2000) doi:10.1016/S0146-6410(00)00103-4 [arXiv:nucl-th/0002059 [nucl-th]].
  * [67] B. L. G. Bakker, A. Bassetto, S. J. Brodsky, W. Broniowski, S. Dalley, T. Frederico, S. D. Glazek, J. R. Hiller, C. R. Ji and V. Karmanov, et al. Nucl. Phys. B Proc. Suppl. 251-252, 165-174 (2014) doi:10.1016/j.nuclphysbps.2014.05.004 [arXiv:1309.6333 [hep-ph]].
  * [68] J. R. Hiller, Prog. Part. Nucl. Phys. 90, 75-124 (2016) doi:10.1016/j.ppnp.2016.06.002 [arXiv:1606.08348 [hep-ph]].
  * [69] S. J. Brodsky, A. Deur and C. D. Roberts, Nature Rev. Phys. 4, no.7, 489-495 (2022) doi:10.1038/s42254-022-00453-3 [arXiv:2202.06051 [hep-ph]].
  * [70] P. A. M. Dirac, Rev. Mod. Phys. 21, 392-399 (1949) doi:10.1103/RevModPhys.21.392
  * [71] S. Weinberg, Phys. Rev. 150, 1313-1318 (1966) doi:10.1103/PhysRev.150.1313
  * [72] S. J. Chang and S. K. Ma, Phys. Rev. 180, 1506-1513 (1969) doi:10.1103/PhysRev.180.1506
  * [73] K. Bardakci and M. B. Halpern, Phys. Rev. 176, 1686-1699 (1968) doi:10.1103/PhysRev.176.1686
  * [74] J. B. Kogut and D. E. Soper, Phys. Rev. D 1, 2901-2913 (1970) doi:10.1103/PhysRevD.1.2901
  * [75] G. Tiktopoulos, Phys. Rev. 180, 1414-1416 (1969) doi:10.1103/PhysRev.180.1414
  * [76] S. J. Chang, R. G. Root and T. M. Yan, Phys. Rev. D 7, 1133-1148 (1973) doi:10.1103/PhysRevD.7.1133
  * [77] T. M. Yan, Phys. Rev. D 7, 1780-1800 (1973) doi:10.1103/PhysRevD.7.1780
  * [78] E. Tomboulis, Phys. Rev. D 8, 2736-2740 (1973) doi:10.1103/PhysRevD.8.2736
  * [79] F. Rohrlich and J. H. Ten Eyck, Phys. Lett. B 46, 102-104 (1973) doi:10.1016/0370-2693(73)90488-7
  * [80] A. Casher, Phys. Rev. D 14, 452 (1976) doi:10.1103/PhysRevD.14.452
  * [81] N. Nakanishi and K. Yamawaki, Nucl. Phys. B 122, 15-28 (1977) doi:10.1016/0550-3213(77)90424-2
  * [82] H. Leutwyler and J. Stern, Annals Phys. 112, 94 (1978) doi:10.1016/0003-4916(78)90082-9
  * [83] C. B. Thorn, Phys. Rev. D 19, 639 (1979) doi:10.1103/PhysRevD.19.639
  * [84] C. B. Thorn, Phys. Rev. D 20, 1934 (1979) doi:10.1103/PhysRevD.20.1934
  * [85] V. A. Franke, Y. V. Novozhilov and E. V. Prokhvatilov, Lett. Math. Phys. 5, 239-245 (1981) doi:10.1007/BF00420704
  * [86] M. Gell-Mann and M. Levy, Nuovo Cim. 16, 705 (1960) doi:10.1007/BF02859738
  * [87] M. Gell-Mann, Phys. Rev. 125, 1067-1084 (1962) doi:10.1103/PhysRev.125.1067
  * [88] M. Gell-Mann and Y. Ne’eman, Annals Phys. 30, 360-369 (1964) doi:10.1016/0003-4916(64)90122-8
  * [89] S. Fubini and G. Furlan, Physics Physique Fizika 1, no.4, 229-247 (1965) doi:10.1103/PhysicsPhysiqueFizika.1.229
  * [90] K. G. Wilson, Phys. Rev. 179, 1499-1512 (1969) doi:10.1103/PhysRev.179.1499
  * [91] H. Fritzsch and M. Gell-Mann, [arXiv:hep-ph/0301127 [hep-ph]].
  * [92] R. P. Feynman, Phys. Rev. Lett. 23, 1415-1417 (1969) doi:10.1103/PhysRevLett.23.1415
  * [93] J. D. Bjorken and E. A. Paschos, Phys. Rev. 185, 1975-1982 (1969) doi:10.1103/PhysRev.185.1975
  * [94] C. H. Llewellyn Smith, Nucl. Phys. B 17, 277-288 (1970) doi:10.1016/0550-3213(70)90165-3
  * [95] R. P. Feynman, M. Kislinger and F. Ravndal, Phys. Rev. D 3, 2706-2732 (1971) doi:10.1103/PhysRevD.3.2706
  * [96] R. P. Feynman,
  * [97] J. B. Kogut and L. Susskind, Phys. Rept. 8, 75-172 (1973) doi:10.1016/0370-1573(73)90009-4
  * [98] G. P. Lepage and S. J. Brodsky, Phys. Rev. D 22, 2157 (1980) doi:10.1103/PhysRevD.22.2157
  * [99] H. C. Pauli and S. J. Brodsky, Phys. Rev. D 32, 2001 (1985) doi:10.1103/PhysRevD.32.2001
  * [100] K. G. Wilson, T. S. Walhout, A. Harindranath, W. M. Zhang, R. J. Perry and S. D. Glazek, Phys. Rev. D 49, 6720-6766 (1994) doi:10.1103/PhysRevD.49.6720 [arXiv:hep-th/9401153 [hep-th]].
  * [101] H. Honkanen, P. Maris, J. P. Vary and S. J. Brodsky, Phys. Rev. Lett. 106, 061603 (2011) doi:10.1103/PhysRevLett.106.061603 [arXiv:1008.0068 [hep-ph]].
  * [102] X. Zhao, H. Honkanen, P. Maris, J. P. Vary and S. J. Brodsky, Phys. Lett. B 737, 65-69 (2014) doi:10.1016/j.physletb.2014.08.020 [arXiv:1402.4195 [nucl-th]].
  * [103] P. Wiecki, Y. Li, X. Zhao, P. Maris and J. P. Vary, Phys. Rev. D 91, no.10, 105009 (2015) doi:10.1103/PhysRevD.91.105009 [arXiv:1404.6234 [nucl-th]].
  * [104] Z. Hu et al. [BLFQ], Phys. Rev. D 103, no.3, 036005 (2021) doi:10.1103/PhysRevD.103.036005 [arXiv:2010.12498 [hep-ph]].
  * [105] S. Nair et al. [BLFQ], Phys. Lett. B 827, 137005 (2022) doi:10.1016/j.physletb.2022.137005 [arXiv:2201.12770 [hep-ph]].
  * [106] S. Nair et al. [BLFQ], Phys. Rev. D 108, no.11, 116005 (2023) doi:10.1103/PhysRevD.108.116005 [arXiv:2302.13645 [hep-ph]].
  * [107] Y. Li and J. P. Vary, Phys. Lett. B 825, 136860 (2022) doi:10.1016/j.physletb.2021.136860 [arXiv:2103.09993 [hep-ph]].
  * [108] Y. Li and J. P. Vary, Phys. Rev. D 105, no.11, 114006 (2022) doi:10.1103/PhysRevD.105.114006 [arXiv:2202.05581 [hep-ph]].
  * [109] Y. Li, P. Maris and J. P. Vary, Phys. Lett. B 836, 137598 (2023) doi:10.1016/j.physletb.2022.137598 [arXiv:2203.14447 [hep-th]].
  * [110] Y. Li, P. Maris, X. Zhao and J. P. Vary, Phys. Lett. B 758, 118-124 (2016) doi:10.1016/j.physletb.2016.04.065 [arXiv:1509.07212 [hep-ph]].
  * [111] Y. Li, P. Maris and J. P. Vary, Phys. Rev. D 96, 016022 (2017) doi:10.1103/PhysRevD.96.016022 [arXiv:1704.06968 [hep-ph]].
  * [112] S. Leitão, Y. Li, P. Maris, M. T. Peña, A. Stadler, J. P. Vary and E. P. Biernat, Eur. Phys. J. C 77, no.10, 696 (2017) doi:10.1140/epjc/s10052-017-5248-0 [arXiv:1705.06178 [hep-ph]].
  * [113] M. Li, Y. Li, P. Maris and J. P. Vary, Phys. Rev. D 98, no.3, 034024 (2018) doi:10.1103/PhysRevD.98.034024 [arXiv:1803.11519 [hep-ph]].
  * [114] L. Adhikari, Y. Li, M. Li and J. P. Vary, Phys. Rev. C 99, no.3, 035208 (2019) doi:10.1103/PhysRevC.99.035208 [arXiv:1809.06475 [hep-ph]].
  * [115] S. Tang, Y. Li, P. Maris and J. P. Vary, Phys. Rev. D 98, no.11, 114038 (2018) doi:10.1103/PhysRevD.98.114038 [arXiv:1810.05971 [nucl-th]].
  * [116] J. Lan, C. Mondal, M. Li, Y. Li, S. Tang, X. Zhao and J. P. Vary, Phys. Rev. D 102, no.1, 014020 (2020) doi:10.1103/PhysRevD.102.014020 [arXiv:1911.11676 [nucl-th]].
  * [117] S. Tang, Y. Li, P. Maris and J. P. Vary, Eur. Phys. J. C 80, no.6, 522 (2020) doi:10.1140/epjc/s10052-020-8081-9 [arXiv:1912.02088 [nucl-th]].
  * [118] S. Tang, S. Jia, P. Maris and J. P. Vary, Phys. Rev. D 104, no.1, 016002 (2021) doi:10.1103/PhysRevD.104.016002 [arXiv:2011.05454 [hep-ph]].
  * [119] Y. Li, M. Li and J. P. Vary, Phys. Rev. D 105, no.7, L071901 (2022) doi:10.1103/PhysRevD.105.L071901 [arXiv:2111.14178 [hep-ph]].
  * [120] M. Li, Y. Li, G. Chen, T. Lappi and J. P. Vary, Eur. Phys. J. C 82, no.11, 1045 (2022) doi:10.1140/epjc/s10052-022-10988-5 [arXiv:2111.07087 [hep-ph]].
  * [121] Z. Wang, M. Li, Y. Li and J. P. Vary, Phys. Rev. D 109, no.3, 3 (2024) doi:10.1103/PhysRevD.109.L031902 [arXiv:2312.02604 [hep-ph]].
  * [122] S. Jia and J. P. Vary, Phys. Rev. C 99, no.3, 035206 (2019) doi:10.1103/PhysRevC.99.035206 [arXiv:1811.08512 [nucl-th]].
  * [123] J. Lan, C. Mondal, S. Jia, X. Zhao and J. P. Vary, Phys. Rev. Lett. 122, no.17, 172001 (2019) doi:10.1103/PhysRevLett.122.172001 [arXiv:1901.11430 [nucl-th]].
  * [124] J. Lan, C. Mondal, S. Jia, X. Zhao and J. P. Vary, Phys. Rev. D 101, no.3, 034024 (2020) doi:10.1103/PhysRevD.101.034024 [arXiv:1907.01509 [nucl-th]].
  * [125] W. Qian, S. Jia, Y. Li and J. P. Vary, Phys. Rev. C 102, no.5, 055207 (2020) doi:10.1103/PhysRevC.102.055207 [arXiv:2005.13806 [nucl-th]].
  * [126] C. Mondal et al. [BLFQ], Phys. Rev. D 104, no.9, 094034 (2021) doi:10.1103/PhysRevD.104.094034 [arXiv:2109.02279 [hep-ph]].
  * [127] L. Adhikari et al. [BLFQ], Phys. Rev. D 104, no.11, 114019 (2021) doi:10.1103/PhysRevD.104.114019 [arXiv:2110.05048 [hep-ph]].
  * [128] Y. Li, P. Maris and J. P. Vary, Phys. Lett. B 836, 137598 (2023) doi:10.1016/j.physletb.2022.137598 [arXiv:2203.14447 [hep-th]].
  * [129] Z. Zhu et al. [BLFQ], Phys. Lett. B 839, 137808 (2023) doi:10.1016/j.physletb.2023.137808 [arXiv:2301.12994 [hep-ph]].
  * [130] C. Mondal, S. Xu, J. Lan, X. Zhao, Y. Li, D. Chakrabarti and J. P. Vary, Phys. Rev. D 102, no.1, 016008 (2020) doi:10.1103/PhysRevD.102.016008 [arXiv:1911.10913 [hep-ph]].
  * [131] S. Xu et al. [BLFQ], Phys. Rev. D 104, no.9, 094036 (2021) doi:10.1103/PhysRevD.104.094036 [arXiv:2108.03909 [hep-ph]].
  * [132] Y. Liu et al. [BLFQ], Phys. Rev. D 105, no.9, 094018 (2022) doi:10.1103/PhysRevD.105.094018 [arXiv:2202.00985 [hep-ph]].
  * [133] Z. Hu et al. [BLFQ], Phys. Lett. B 833, 137360 (2022) doi:10.1016/j.physletb.2022.137360 [arXiv:2205.04714 [hep-ph]].
  * [134] T. Peng et al. [BLFQ], Phys. Rev. D 106, no.11, 114040 (2022) doi:10.1103/PhysRevD.106.114040 [arXiv:2208.00355 [hep-ph]].
  * [135] Z. Zhang et al. [BLFQ], Phys. Rev. D 109, no.3, 034031 (2024) doi:10.1103/PhysRevD.109.034031 [arXiv:2312.00667 [hep-th]].
  * [136] Z. Zhu et al. [BLFQ], Phys. Rev. D 108, no.3, 036009 (2023) doi:10.1103/PhysRevD.108.036009 [arXiv:2304.05058 [hep-ph]].
  * [137] Z. Kuang et al. [BLFQ], Phys. Rev. D 105, no.9, 094028 (2022) doi:10.1103/PhysRevD.105.094028 [arXiv:2201.06428 [hep-ph]].
  * [138] J. Lan et al. [BLFQ], Phys. Lett. B 825, 136890 (2022) doi:10.1016/j.physletb.2022.136890 [arXiv:2106.04954 [hep-ph]].
  * [139] S. Xu, C. Mondal, X. Zhao, Y. Li and J. P. Vary, [arXiv:2209.08584 [hep-ph]].
  * [140] S. Xu et al. [BLFQ], Phys. Rev. D 108, no.9, 094002 (2023) doi:10.1103/PhysRevD.108.094002
  * [141] B. Lin et al. [BLFQ], Phys. Lett. B 847, 138305 (2023) doi:10.1016/j.physletb.2023.138305 [arXiv:2308.08275 [hep-ph]].
  * [142] S. Kaur et al. [BLFQ], Phys. Lett. B 851, 138563 (2024) doi:10.1016/j.physletb.2024.138563 [arXiv:2401.03480 [hep-ph]].
  * [143] S. J. Brodsky, G. F. de Teramond, H. G. Dosch and J. Erlich, Phys. Rept. 584, 1-105 (2015) doi:10.1016/j.physrep.2015.05.001 [arXiv:1407.8131 [hep-ph]].
  * [144] S. Weinberg, Cambridge University Press, 2005, ISBN 978-0-521-67053-1, 978-0-511-25204-4 doi:10.1017/CBO9781139644167
  * [145] C. Lorcé and P. Lowdon, Eur. Phys. J. C 80, no.3, 207 (2020) doi:10.1140/epjc/s10052-020-7779-z [arXiv:1908.02567 [hep-th]].
  * [146] B. R. Barrett, P. Navratil and J. P. Vary, Prog. Part. Nucl. Phys. 69, 131-181 (2013) doi:10.1016/j.ppnp.2012.10.003
  * [147] K. Hornbostel, S. J. Brodsky and H. C. Pauli, Phys. Rev. D 41, 3814 (1990) doi:10.1103/PhysRevD.41.3814
  * [148] N. Anand, A. L. Fitzpatrick, E. Katz, Z. U. Khandker, M. T. Walters and Y. Xin, [arXiv:2005.13544 [hep-th]].
  * [149] S. S. Chabysheva and J. R. Hiller, Phys. Lett. B 711, 417-422 (2012) doi:10.1016/j.physletb.2012.04.032 [arXiv:1103.0037 [hep-ph]].
  * [150] S. D. Glazek and K. G. Wilson, Phys. Rev. D 48, 5863-5872 (1993) doi:10.1103/PhysRevD.48.5863
  * [151] M. M. Brisudova, R. J. Perry and K. G. Wilson, Phys. Rev. Lett. 78, 1227-1230 (1997) doi:10.1103/PhysRevLett.78.1227 [arXiv:hep-ph/9607280 [hep-ph]].
  * [152] M. M. Brisudova, S. Szpigel and R. J. Perry, Phys. Lett. B 421, 334-340 (1998) doi:10.1016/S0370-2693(97)01556-6 [arXiv:hep-ph/9709479 [hep-ph]].
  * [153] S. D. Głazek, M. Gómez-Rocha, J. More and K. Serafin, Phys. Lett. B 773, 172-178 (2017) doi:10.1016/j.physletb.2017.08.018 [arXiv:1705.07629 [hep-ph]].
  * [154] S. S. Chabysheva and J. R. Hiller, Annals Phys. 337, 143-152 (2013) doi:10.1016/j.aop.2013.06.016 [arXiv:1207.7128 [hep-ph]].
  * [155] G. F. de Teramond and S. J. Brodsky, Phys. Rev. D 104, no.11, 116009 (2021) doi:10.1103/PhysRevD.104.116009 [arXiv:2103.10950 [hep-ph]].
  * [156] Li, Yang (2019), “Heavy quarkonium light front wave functions from basis light-front quantization with a running coupling”, Mendeley Data, V2, doi: 10.17632/cjs4ykv8cv.2
  * [157] A. Freese and G. A. Miller, Phys. Rev. D 103, 094023 (2021) doi:10.1103/PhysRevD.103.094023 [arXiv:2102.01683 [hep-ph]].
  * [158] A. Freese and G. A. Miller, Phys. Rev. D 104, no.1, 014024 (2021) doi:10.1103/PhysRevD.104.014024 [arXiv:2104.03213 [hep-ph]].
  * [159] A. Freese and G. A. Miller, Phys. Rev. D 105, no.1, 014003 (2022) doi:10.1103/PhysRevD.105.014003 [arXiv:2108.03301 [hep-ph]].
  * [160] Y. Li, W. b. Dong, Y. l. Yin, Q. Wang and J. P. Vary, Phys. Lett. B 838, 137676 (2023) doi:10.1016/j.physletb.2023.137676 [arXiv:2206.12903 [hep-ph]].
  * [161] A. Freese and G. A. Miller, Phys. Rev. D 108, no.3, 034008 (2023) doi:10.1103/PhysRevD.108.034008 [arXiv:2210.03807 [hep-ph]].
  * [162] A. Freese and G. A. Miller, Phys. Rev. D 108, no.9, 094026 (2023) doi:10.1103/PhysRevD.108.094026 [arXiv:2307.11165 [hep-ph]].
  * [163] O. V. Teryaev, Front. Phys. (Beijing) 11, no.5, 111207 (2016) doi:10.1007/s11467-016-0573-6
  * [164] L. Rezzolla and O. Zanotti, Oxford University Press, 2013, ISBN 978-0-19-174650-5, 978-0-19-852890-6 doi:10.1093/acprof:oso/9780198528906.001.0001
  * [165] C. Lorcé, Eur. Phys. J. C 78, no.2, 120 (2018) doi:10.1140/epjc/s10052-018-5561-2 [arXiv:1706.05853 [hep-ph]].
  * [166] S. J. Brodsky, D. S. Hwang, B. Q. Ma and I. Schmidt, Nucl. Phys. B 593, 311-335 (2001) doi:10.1016/S0550-3213(00)00626-X [arXiv:hep-th/0003082 [hep-th]].
  * [167] X. Cao, Y. Li and J.P. Vary, in prepration
  * [168] M. H. L. Pryce, Proc. Roy. Soc. Lond. A 195, 62-81 (1948) doi:10.1098/rspa.1948.0103
  * [169] C. Lorcé, Eur. Phys. J. C 78, no.9, 785 (2018) doi:10.1140/epjc/s10052-018-6249-3 [arXiv:1805.05284 [hep-ph]].
  * [170] A. P. Trawiński, S. D. Głazek, S. J. Brodsky, G. F. de Téramond and H. G. Dosch, Phys. Rev. D 90, no.7, 074017 (2014) doi:10.1103/PhysRevD.90.074017 [arXiv:1403.5651 [hep-ph]].
  * [171] K. L. Ng, Phys. Rev. D 47, 5187-5190 (1993) doi:10.1103/PhysRevD.47.5187 [arXiv:gr-qc/9305002 [gr-qc]].
  * [172] U. Özdem and K. Azizi, Phys. Rev. D 101, no.5, 054031 (2020) doi:10.1103/PhysRevD.101.054031 [arXiv:1912.06375 [hep-ph]].
  * [173] M. V. Polyakov and A. Tandogan, Phys. Rev. D 101, no.11, 118501 (2020) doi:10.1103/PhysRevD.101.118501
  * [174] H. Alharazin, D. Djukanovic, J. Gegelia and M. V. Polyakov, Phys. Rev. D 102, no.7, 076023 (2020) doi:10.1103/PhysRevD.102.076023 [arXiv:2006.05890 [hep-ph]].
  * [175] J. Y. Kim, Phys. Lett. B 834, 137442 (2022) doi:10.1016/j.physletb.2022.137442 [arXiv:2206.10202 [hep-ph]].
  * [176] U. Özdem and K. Azizi, JHEP 03, 048 (2023) doi:10.1007/JHEP03(2023)048 [arXiv:2212.07290 [hep-ph]].
  * [177] J. Y. Kim, H. Y. Won, J. L. Goity and C. Weiss, Phys. Lett. B 844, 138083 (2023) doi:10.1016/j.physletb.2023.138083 [arXiv:2304.08575 [hep-ph]].
  * [178] H. Alharazin, B. D. Sun, E. Epelbaum, J. Gegelia and U. G. Meißner, JHEP 03, 007 (2024) doi:10.1007/JHEP03(2024)007 [arXiv:2312.05193 [hep-ph]].
  * [179] H. Alharazin, Phys. Rev. D 109, no.1, 016009 (2024) doi:10.1103/PhysRevD.109.016009 [arXiv:2312.09675 [hep-ph]].