Stress out of charmonia

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 light-front energy distributions of the system. We find that hadrons are multi-layer systems.


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 µν [2].These operators also dictate how the system gravitates, and conversely, how it responds to the gravitational field [3].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 γ * γ → π 0 π 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/ψ [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 π = 170 MeV [15,16] (cf.[17][18][19][20]).The hadronic EMT were also investigated in various theories and phenomenological models, e.g., .See Refs.[56,57] for recent reviews.
One of the challenges to access the hadronic EMT is that T µν 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.Since m c ≫ Λ QCD and the strong coupling α s (m c ) ≪ 1, phenomena occurring at the scale ∼ m c may be treated perturbatively.On the other hand, with the binding energy α s m c ≳ Λ 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 nonperturbative physics of QCD.It is sometimes dubbed as the "hydrogen atom of QCD" [58].
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 [52,.For charmonium, the effective Hamiltonian combines confinement from holographic lightfront QCD [142] for the long-distance physics as well as one-gluon exchange interaction from light-front QCD [102] for the short-distance physics [109].The computed arXiv:2404.06259v3[hep-ph] 19 Aug 2024 mass spectrum is shown to be within 40 MeV deviation with the PDG compilation [110,114,116,117].For the ground state η c , the mass deviation is about ∼ 15% of the total binding energy.The obtained lightfront wave functions (LFWFs) were used to investigate a variety of observables, including the decay constants [109,110], electromagnetic form factors [113], parton distribution functions [115], and GPDs [113].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 [118][119][120].
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 |ψ h ⟩ can be obtained by solving the light-front Schrödinger equation [63], Here, we adopt the light-front coordinates, x ± = x 0 ± x 3 , ⃗ x ⊥ = (x 1 , x 2 ).The hadronic state vector |ψ h (p, j, m j )⟩ can be further classified by its momentum p µ , total spin j and spin magnetic projection m j , viz.[63,143], where, P µ is the 4-momentum operator, ⃗ S 2 = W 2 /P 2 is the total spin operator, constructed from the Pauli-Lubanski vector W µ = 1 2 ε µνρσ P ν J ρσ , where J αβ 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 ⊥ + 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], Here, ω µ = (ω + , ω − , ⃗ ω ⊥ ) = (0, 2, ⃗ 0 ⊥ ) 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., where, T µν is the EMT.In particular, the density of the light-front interaction energy H int is related to the +− component of the EMT, H int = 1 2 T +− int .Momentum conservation (2) imposes constraints on the forward HMEs of the EMT.It is straightforward to show, Similar constraints can be obtained from the conservation of the angular momentum (3-4) as shown by Lorcé et.al. [144].
In QCD, the force is strong and we need a nonperturbative approach to solve the eigenvalue problems (2)(3)(4).Inspired by the success of ab initio nuclear structure theory, especially the no-core shell model (NCSM) [145], 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, j where H ij = ⟨i|H|j⟩ and ψ i = ⟨i|ψ⟩.Here |i⟩ is the basis state and |ψ⟩ is the state vector of a hadron, and ⊥ is equivalent to the light-front Hamiltonian operator P − since P + and ⃗ P 2 ⊥ 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 [146], conformal basis functions employed in Hamiltonian conformal truncation method [147], and coherent basis in light-front coupled cluster method [148].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 [149][150][151][152][153]. In the bottom-up approach, one starts with a phenomenological effective Hamiltonian that embodies the key physics of the system [109,110].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 [142].
In this work, we adopt the phenomenological effective Hamiltonian with a confining interaction from lightfront holographic QCD.The form is remarkably simple, V q q ∝ x(1 − x)r 2 ⊥ , where, x = p + q /P + h is the longitudinal momentum fraction of the quark, and ⃗ r ⊥ 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 [106,107,154,155].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 [102].
The basis functions are the eigen-functions of the confining interactions.Their analytic forms are, where, N l = 4π(2l The basis parameters α = 2m q (m q + m q )/κ 2 , β = 2m q (m q + m q )/κ 2 , where κ is the strength of the confining interaction and m q = m q = m c is the effective quark mass.Then effective Hamiltonian for the cc system is, Here, I, J, I ′ , J ′ each represent a collection of singleparticle 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 2 c + 2κ 2 (2n + |m| + l + 3 2 ) + (κ 4 /4m 2 c )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 cc Fock sector.
BLFQ adopts the N max basis truncation scheme to render the basis space finite.Specifically for charmonium, we impose the conditions, N max defines the UV and IR regulators of the basis space as Λ uv = √ N max κ, and λ ir = κ/ √ N max .L max defines the resolution of the longitudinal momentum fraction δx = L −1 max .For simplicity, we have tied the basis scale parameter to the confining strength κ.
We diagonalized the Hamiltonian with different N max , L max to obtain the mass squared eigenvalues and the LFWFs.The model parameters κ, and m c were fit to the charmonium mass spectrum below the DD threshold, which is a good approximation for states with narrow resonance.Since radiative corrections are neglected, we adopt N max = L max = 8 for observables that are sensitive to the short-distance quantum fluctuations, which corresponds to a UV resolution Λ uv ≈ M cc .To estimate the uncertainty associated with the basis resolution, we quote the difference between the N max = L max = 8 results and the N max = L max = 16 results.The details of the model can be found in Ref. [110].
From the basis functions, we can reconstruct the valence LFWFs, where, ψ (mj ) ss/h (n, m, l) are the basis coefficients obtained from the diagonalization.The state vector can be represented by the LFWFs as, Here, P µ = (P − , P + , ⃗ P ⊥ ) is the 4-momentum of the bound state h.⃗ k ⊥ = ⃗ p ⊥ − x ⃗ P ⊥ 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 η c and χ c0 , can be parametrized in terms of covariant tensors [56], 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 [49,[157][158][159][160][161][162][163], where, t αβ (⃗ q ⊥ ; P ) is the normalized hadronic matrix elements, The factor 1/(2P + ) is due to the normalization of the state vector: Note that t αβ (hence T αβ ) depends on the c.m. momentum P µ of the hadron [159].The density of the lightfront longitudinal momentum P + is [161,164], Similarly, the density of the light-front transverse momentum ⃗ P ⊥ is, From these expressions, the Fourier transform of GFF A(q 2 ⊥ ), can be interpreted as the matter density.In Refs.[161,162,164], this density is called the internal light-front longitudinal momentum density.We emphasize that A(r ⊥ ) describes the convective distribution of all three momenta.Note that for spin-1 2 particles, such as the nucleon, there is an extra contribution from spin, Here, ⃗ S(r ⊥ ) = 2J (r ⊥ )⃗ s is the spin density, ⃗ s = sẑ, and J (r ⊥ ) is the 2D Fourier transform of the GFF J(Q 2 ), which exists for spin- Momentum conservation (8) implies the conservation of matter, The 3D r.m.s radius of the matter distribution is, The distribution of the light-front energy P − = (P 2 ⊥ + M 2 )/P + in a hadron state with mass M reads [49], where, The first term A(r ⊥ )P 2 ⊥ /P + represents the light-front energy due to the motion of the particle.Here, we have introduced a new form factor M 2 (q 2 ⊥ ) and its 2D lightfront density M 2 (r ⊥ ).In Refs.[161,162], this density [Eq.( 27)] is termed the internal light-front energy density.We find it is better to interpret M 2 (r ⊥ ) as the distribution of the invariant mass squared.Hereafter, we will use these two terms interchangeably.For spin-1 2 particles, the spin term also contributes to the light-front energy density T +− (r ⊥ ), The last term vanishes upon integration.
Light-front energy conservation (8) implies [49], From (28), the above condition is equivalent to lim As we will see later, this forward limit constraint, also known as the von Laue condition, is a necessary condition for the force equilibrium of a hadron [27,56].It is useful to define the r.m.s.radius of the invariant mass squared distribution as, Here, λ 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 E(x), pressure P(x) and (traceless) shear Π αβ (x) as [49,166], Here, u α ≡ P α / √ P 2 is the 4-velocity vector.The energy density can be extracted as, where, the energy form factor E(q 2 ⊥ ) is defined as, For spin-1 2 particles, the energy form factor should include a spin contribution, This is in contrast to M 2 (q 2 ⊥ ), the light-front energy form factor, in which, the spin contribution is absent.The energy-momentum conservation ( 24) & (30) [or (31)] implies i.e. energy density is normalized to the total energy in local rest frame, i.e. the mass.The 3D r.m.s.radius of the energy distribution is, The pressure can be extracted as, The form factor q 2 ⊥ D(q 2 ⊥ ) is also known as the hadronic cosmological constant [165], in analogy to the quantum expectation value of the EMT with respect to the vacuum, ⟨0|T µν |0⟩ = g µν Λ.One can readily identify Λ = 1 2 q 2 D(q 2 ).The von Laue condition (31) implies that hadrons are in mechanical equilibrium, The shear is also related to the GFF D(q 2 ⊥ ), where, ri ⊥ = r i ⊥ /r ⊥ .The shear is "traceless", ∆ αβ Π αβ = 0, where ∆ αβ = u α u β − g αβ is the spatial metric tensor.In fact, pressure and shear are the trace and traceless parts of the Cauchy stress tensor C αβ = Π αβ − P∆ αβ .
Note that the energy density E(r ⊥ ) is related to the invariant mass squared density M 2 (r ⊥ ) through, Physically, the proper energy density E is generally assumed positive, a constraint known as the weak energy condition, which provides a constraint on D: These conditions are trivially satisfied for negative D, which is conjectured for stable systems such as the proton [167,168].
A related quantity of interest is the scalar density, θ = T α α , which obeys the classical relation [169], The scalar radius is ).The D-term is related to the second moment of the pressure, It is conjectured that D < 0 for a mechanically stable system [167,168].Intuitively, a negative D represents a repulsive core along with an attractive periphery.A negative D-term suggests a chain of inequalities, For charmonium, as we will show below, D ∼ −5.Then, r E > r A for this system.Then, we obtain a layered picture of charmonium as shown in Fig. 1.
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 A(r ⊥ ).Therefore, at the core of a charmonium are the valence quarks.For the invariant mass squared density M 2 (r ⊥ ), small-x partons, known as "wee partons", also have a significant contribution since the light-front energy p − ∝ (p 2 ⊥ + m 2 )/xP + .The fact that r M 2 > r A implies that the distribution of the wee partons is broader FIG. 1.A 3D picture of the normalized matter density A(r ⊥ ), energy density E(r ⊥ ), invariant mass squared density M 2 (r ⊥ ), aka.internal light-front energy density, scalar density θ(r ⊥ ) and pressure P(r ⊥ ) on the light front.All of these densities, except the pressure, are normalized to unity upon integration over the transverse plane.For the pressure, we plot P(r ⊥ )/M .
than the valence quarks.The scalar trace T α α 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 [169,170].In principle, the light-front energy 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 total internal lightfront energy M 2 is a Lorentz scalar.In the previous work [49], some of us proposed to decompose the light-front energy 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.
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, where, recall, M 2 (q 2 ⊥ ) = (M 2 + 1 4 q 2 ⊥ )A(q 2 ⊥ ) + 1 2 q 2 ⊥ D(q 2 ⊥ ).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 [171], where, And the n-body integration measure is defined as, where, S n is the symmetry factor.Using the coordinate space representation developed in our previous work [49], A(q 2 ⊥ ) becomes, x j e i⃗ r j⊥ •⃗ q ⊥ , ( where, Then, the matter distribution becomes a one-body density (OBD), And the corresponding r.m.s.matter radius is, 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, Similarly, the matter density reads, In basis representation, Note that m = m ′ = m j − s − 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 A(⃗ r ⊥ ) reads, B. Gravitational form factor D: kinetic part The GFF D(q 2 ⊥ ) and the light-front energy form factor M 2 (q 2 ⊥ ) can be extracted from either the light-front energy density T +− , or the transverse stress T 11 + T 22 and T 12 .In practical model calculations, such as this work, truncations and approximations are introduced which result in the loss of some Poincaré symmetries.One then deals with a system with a reduced number of symmetries, and GFFs extracted from different current components are no longer equivalent.One possible way to extract the form factors is to parametrize the HMEs in terms of the reduced symmetry.Effectively, one introduces additional Lorentz tensor structures.Form factors associated with non-covariant structures are called spurious ones.See Ref. [49] and the references therein for a discussion on this approach to the EMT.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 condition1 .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: In the Breit frame P ⊥ = 0, only the light-front energy form factor M2 (q 2 ⊥ ) survives.Form factor M 2 (q 2 ⊥ ) is related to GFFs A and D as given in (28).
In Ref. [49], we propose to split T +− into two parts T +− = T +− 0 + T +− int : a free part T +− 0 and an interaction part T +− int .We adopt the light-cone gauge A + = 0 and in this gauge, the kinetic part T +− 0 is simple.The corresponding invariant mass squared form factor is also split into two pieces, M 2 (q 2 ⊥ ) = M 2 0 (q 2 ⊥ ) + M 2 int (q 2 ⊥ ).The free energy density T +− 0 is diagonal in Fock space, and the corresponding invariant mass squared form factor admits an exact LFWF representation, where, The integrand is the n-body kinetic energy with an extra term − 1 4 q 2 ⊥ .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 [2] (cf.[49]).
Using the coordinate space LFWFs, M 2 0 (q 2 ⊥ ) is dramatically simplified as, Its Fourier transform gives the OBD of the free light-front energy, For charmonium, In basis representation, C. 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, in non-relativistic quantum many-body theory, the Hamiltonian density can also be obtained by localizing the Hamiltonian operator, viz.inserting δ 3 (r − r i ) for the i-th constituent.For example, the charge density can be obtained by localizing the charge, The above method can be generalized to the relativistic quantum many-body systems.In relativistic quantum theory, particles cannot be consistently localized 2 .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-δ prescription δ 2 (r ⊥ − r i⊥ ), 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, One can show that the charge density obtained in this way reproduces the Drell-Yan-West formula [173,174].
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, where, ⃗ r i⊥ is the transverse coordinate of the i-th parton.
Based on the localization prescription, we can write down the corresponding energy density as, 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⊥ − r j⊥ ), we first observe that the interaction energy in the forward limit admits a diagonal representation, = where, the n-body potential energy is expressed in terms of the mass eigenvalue and the kinetic energy 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, assuming the interaction is two-body.For charmonium, the above expression becomes, where In terms of the momentum-space LFWFs, the above expression can also be written as, In basis representation, In the forward limit (q 2 = 0), it is evident that from the LFWF representation.From (28), this condition implies lim Therefore, the forward limit constraint is satisfied.

V. NUMERICAL RESULTS
We substitute the charmonium LFWFs obtained from BLFQ [156].As mentioned, the form factors are evaluated with N max = L max = 8.The uncertainty associated with the basis resolution is quoted as the difference between the N max = L max = 8 results and the N max = L max = 16 results.
Figure 2 shows the GFFs A(Q 2 ), D(Q 2 ) for the ground-state pseudoscalar charmonium η c .From these results, we can extract D and radius r A .The results are listed in Table I.Also listed in Table I are the radii, r E , r M 2 , r θ , as combinations of these two primary quantities.A somewhat surprising result is the large contribution from terms proportional to λ C = 1/M = 0.066 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 ∼ B −1 ∼ α −1 s λ C , where the strong coupling constant α s ∼ 0.3 for charmonium.Therefore, r A and λ C are actually of the same order of magnitude.
c and χc0.The uncertainty bands are determined as described in the caption to Fig. 2.
Figure 3 compares the GFFs for η c with those for its radial and angular excitations, viz.η ′ c and χ c0 .Radial and angular nodes appear for excited states as expected.The corresponding D-term and the radii are extracted and collected in Table I.From these numbers, the matter radius r A of the P -wave charmonium χ c0 is larger than that of η c by 30%.The radius r A of the η ′ c , the radial excitation, almost doubles that of the ground state.
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 η c , η ′ c [i.e.η c (2S)] and χ 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 [167,168].The pressure profile of the radial excited charmonium η ′ c shows two additional nodes [23].Note that the uncertainty of the pressure at small r ⊥ is large.This is because our basis has a finite UV coverage.The effective UV cutoff for N max = 8 is 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.
We have mentioned the difference between the proper energy distribution E(r ⊥ ) defined from the instant form and the invariant mass squared distribution M 2 (r ⊥ ) 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 η c and its normalized invariant mass squared distribution.The energy distribution is positive, consistent with the weak energy condition.However, the invariant mass squared distribution becomes negative (i.e.tachyonic) at short distance r ⊥ , although at the short distance, the results suffer from slow basis convergence as indicated by the wide bands.Similar quantities can be extracted and compared for η ′ c and χ c0 , as shown in Fig. 6.Within basis uncertainty, the energy distributions E(r ⊥ ) are consistent with the weak energy condition that the energy densities are positive.
We have so far introduced several densities: the matter density A(r ⊥ ), energy density E(r ⊥ ), invariant mass squared density M 2 (r ⊥ ) and (trace) scalar density θ(r ⊥ ).Both E(r ⊥ ) and A(r ⊥ ) are strictly positive as required by the weak and null energy conditions.The energy conditions impose constraints on D (46), which are satisfied by a negative D. We have argued that a negative D suggests a hierarchy of the radii associated with these densities: Figure 7 compares these densities for η c after the normalization, i.e. all densities are normalized to unity.Both the invariant mass squared density M 2 (r ⊥ ) and the scalar density θ(r ⊥ ) become negative at the center of the meson.
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 internal light-front energy (i.e.free invariant mass squared) as well as the interacting light-front energy for η c , η ′ c and χ c0 are shown in Fig. 8.Note that M 2 0 (r ⊥ ) does not remain positive for all distance for η ′ c and χ c0 .This is caused by the recoil term −q 2 ⊥ /4x in (66).The interaction density M 2 int (r ⊥ ) of η ′ c becomes positive at some large distance, consistent with the positive confining potential r 2 ⊥ dominant at large parton separation.
Non-relativistically, for particles within a potential of the form V (r) = αr n , the virial theorem3 suggests, n = 2⟨T ⟩/⟨V ⟩.Inspired by this result, we introduce a virial scaling index where, we define the kinetic energy, i.e. the virial, as the free energy M 2 0 ≡ M 2 0 (Q 2 = 0) subtracting the square of the total quark mass, The potential energy is simply the interacting energy, 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.η c , and +1 for bound state with large size [178].The obtained indices are collected Table I.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 M 2 of charmonia η c , η ′ c and χ c0 into quark mass contribution, kinetic energy contribution and the potential energy contribution.

VI. SUMMARY AND OUTLOOK
In this work, we investigate the gravitational form factors of the charmonium system.Starting from the lightfront wave functions obtained previously in basis lightfront 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 P(r ⊥ ), the (proper) energy density E(r ⊥ ), invariant mass squared density (aka.internal light-front energy density) M 2 (r ⊥ ) and trace scalar density θ(r ⊥ ).We also   identify the positive-defined density A(r ⊥ ), the Fourier transform of GFF A(Q 2 ), as the (convective) matter density.The energy density E(r ⊥ ) 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 θ , implying an onion-like structure of hadrons.Finally, we investigate the virial scaling index 2⟨T ⟩/⟨V ⟩, 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/ψ, 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 ⊥ /4x.A closely related set of observables are the gravitational transition form factors [30,39,53,[179][180][181][182][183][184], which describe hadron production in gravity.These observables can also be accessed using the methods proposed in this work.
One of the interesting 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 lightfront 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 lightfront 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   all associated with a spurious GFF -a 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.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 [137][138][139][140][141]. How to tame the computational complexity as well as the nonperturbative renormalization arising in these problems are critical challenges to be tackled.

Our model does not contain dynamical gluons or sea
FIG. 2. GFFs A(Q 2 ) (top) and D(Q 2 ) (bottom) of groundstate pseudoscalar charmonium ηc.The central value is based on Nmax = 8 wave functions and the difference between Nmax = 8 and Nmax = 16 results is shown as the bands to indicate the basis sensitivity.

- 1 ]FIG. 5 .
FIG. 5. Comparison of the normalized energy distribution and the normalized distribution of the invariant mass squared within ηc.The uncertainty bands are determined as described in the caption to Fig. 2.

2 [fm - 1 ]FIG. 6 .
FIG. 6.Comparison of the energy distributions (top) and the distribution of the invariant mass squared (bottom) of η ′ c and χc0.The uncertainty bands are determined as described in the caption to Fig. 2.

2 [fm - 1 ]FIG. 8 .
FIG. 8. (Top) Decomposition of the internal light-front energy density, i.e. invariant mass squared density, of ηc into the free part and the interacting part.(Middle) Comparison of the free invariant mass squared densities for ηc, η ′ c and χc0.(Bottom) Comparison of the interacting invariant mass squared densities for ηc, η ′ c and χc0.The uncertainty bands are determined as described in the caption to Fig. 2.

FIG. 9 .
FIG. 9.Decomposition of charmonium invariant mass squared M 2 in terms of the total quark mass squared 4m 2 c , kinetic light-front energy ⟨T ⟩ and potential light-front energy ⟨V ⟩.See texts.

TABLE I .
BLFQ prediction of the D-term and the radii of charmonium ηc, η ′ c and χc0.We adopt the Nmax = 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 Nmax = 8 and Nmax = 16 results as the basis sensitivity.