Form Factors and Generalized Parton Distributions of Heavy Quarkonia in Basis Light Front Quantization

We calculate the electromagnetic (charge, magnetic and quadrupole) form factors and the associated static moments of heavy quarkonia (charmonia and bottomonia) using the Basis Light Front Quantization (BLFQ) approach. For this work, we adopt light front wavefunctions (LFWFs) generated by a holographic QCD confining potential and a one-gluon exchange interaction with fixed coupling. We compare our BLFQ results with the limiting case of a single BLFQ basis state description of heavy quarkonia and with other available results. These comparisons provide insights into relativistic effects. Using the same LFWFs generated in the BLFQ approach, we also present the generalized parton distributions (GPDs) for selected mesons including those for radially excited mesons such as $\psi^\prime$ and $\Upsilon^\prime$. Our GPD results establish the foundation within BLFQ for further investigating hadronic structure such as probing the spin structure of spin-one hadrons in the off-forward limit.

therefore natural to ask what one can learn about the spin-one hadronic structure from an investigation using the recently developed non-perturbative Basis Light Front Quantization (BLFQ) approach. In this paper, we present the EM FFs and the associated static moments, the charge radii, magnetic moments and the quadrupole moments, for a selection of heavy quarkonia. We also calculate corresponding quantities in a simplified basis, which is referred to as the Single Basis Limit (SBL) approach in this paper, and interpret the differences with the BLFQ results to uncover dynamical effects arising from the different interactions. We also compare our results with corresponding results from the Contact Interaction (CI) [16][17][18], Dyson-Schwinger Equation (DSE) [4,5] and Lattice [3] approaches wherever available.
The generalized parton distribution (GPD) has emerged as a powerful tool to describe hadrons in terms of quark and gluon degrees of freedom. Several reviews cover the GPDs and their connections to experiments [25][26][27][28][29][30][31][32][33][34][35][36][37]. In particular, there are several investigations on the spin-zero and spin-one GPDs [26,[30][31][32][33][34][35][36]. For example, in Refs. [33,36], the pion GPDs have been calculated in the LF phenomenological models with both the valence and non-valence contributions. Similarly, in Ref. [26], the angular momentum sum rule for the spin-one system within the gaugeinvariant decomposition framework has been investigated. That study has also discussed the connections between the deeply virtual Compton scattering (DVCS) amplitudes and the total quark angular momentum for the ground-state vector meson GPDs. In Ref. [35], The deuteron GPDs have been investigated within the impulse approximation in a framework with non-zero longitudinal momentum transfer. The recent study of the deuteron in Ref. [31] has also presented the GPDs in the framework of holographic QCD.
These studies, however, do not address the heavy mesons (for example, cc and bb). Hence, such basic properties as the momentum-transfer dependence of the valence quarks for the radially excited states ψ and Υ are not previously available. Furthermore, we are motivated by the feasibility of experiments to investigate the hadronic structure of the (pseudo) scalar and vector meson GPDs in the forward limit (zero momentum transfer limit). Such measurements provide connections with unpolarized parton distributions [30].
In this work, we calculate the EM FFs and GPDs through the corresponding matrix elements which are defined by the overlap integrals of light front wavefunctions (LFWFs) in the Drell-Yan frame. The non-perturbative solutions for the LFWFs are provided by a recent BLFQ study [20] of heavy quarkonia. This work implements a transverse confining potential from light front holography and a longitudinal confining interaction which has a similar shape in the non-relativistic limit. It also includes the one-gluon exchange interaction with a fixed coupling to generate the spin structure of the charmonium and bottomonium systems. Note a recent investigation [38] following the quarkonia study [20] with a running coupling [21] provides comparisons of the mass spectrum and decay constant between the results obtained from the BLFQ and that from the covariant spectator theory (CST). The CST treatment [38] is an independent, fully relativistic approach and the comparison of the results of CST and BLFQ showed favorable correspondence.
For this work, we adopt the BLFQ approach which is developed for solving bound-state problems in quantum field theory [20,39,40]. This approach not only provides easy conversion between the transverse coordinates and momentum space [21,41], but also connects mass spectroscopy with other observables [20,42]. The BLFQ is a Hamiltonian-based formalism that uses the advantages of LF dynamics [43] with advances in solving many body bound-state problems [44]. It has been successfully applied to the single electron problem in quantum electrodynamics (QED) [45], the strong coupling bound-state positronium problem [40,41] and the running coupling quarkonium problem [21] . Furthermore, the BLFQ approach has been extended to time-dependent strong external field problems such as non-linear Compton scattering [46]. The reviews related to BLFQ and its application are available in Refs. [20-23, 38-42, 44-50].
We organize this paper as follows: In Sec. II, FFs and GPDs are defined through the local and non-local matrix elements of the plus component of the current operator, respectively. Then, the FFs and GPDs are expressed in terms of the overlap integrals of LFWFs. In Sec. III, we briefly introduce BLFQ along with SBL, the simplified case of a single BLFQ basis state, and in Sec. IV, we present our results. Finally, we present the summary of this work and outlook for further research in Sec. V.

ON THE LIGHT FRONT
The Lorentz-invariant, elastic FFs F i (t) for spin-one hadrons are defined by the local matrix elements of the current operator J µ [ ψ (0)γ µ ψ(0)] as [1,12,51] where ψ(ψ) is the quark (anti-quark) field operator, p (p ) is the momentum of the initial (final) state of the hadron, J is total angular momentum for the hadron, m J (m J ) is the total angular-momentum projection in the initial (final) state of the hadron, t ≡ (p − p) 2 , M is the mass of the hadron, P = p + p, = (p, m J ) and = (p , m J ) are the polarization vectors of the hadron in the initial and final helicity states, respectively, satisfying · p = · p = 0, and We adopt the following conventions [12] to calculate helicity amplitudes I m J ,m J (t) in the Drell-Yan equivalent frame. where is the light front variables in this paper, q = √ −t and τ ≡ −q 2 /(4M 2 ).
There is only one helicity amplitude I 0,0 (t) for J = m J = 0 that can be computed from the plus component of the current defined in Eq. (1), and the charge form factor for spin-zero hadron is therefore defined by G C (t) ≡ I 0,0 (t).
But, in the case of J = 1 with m J (and m J )= +, 0, −, there are nine helicity amplitudes I m J ,m J (t) that can be computed for the same current. One can reduce them to four amplitudes I +,− (t), I +,+ (t), I +,0 (t) and I 0,0 (t) using the light front parity and the charge conjugation symmetries in LF dynamics.
In the case of the spin-one hadrons, there are three Lorentz-invariant, elastic FFs F i (t), and hence three EM FFs, but there are four helicity amplitudes. Studies presented in Refs. [1,6,9,10,12,13] have claimed that computing the EM FFs is more feasible than the elastic FFs F i (t). The four helicity amplitudes in LF dynamics and the three EM FFs create an ambiguity on how to compute them in the case where the current conservation is not preserved, and therefore the relations that define the EM FFs from the helicity amplitudes are not unique. There are several choices in which the four helicity amplitudes can be combined to extract the EM FFs. One can find the most popular choices in Refs. [1,[8][9][10]13]. The studies available in Refs. [1,6,12] suggest to adopt the prescription defined by Grach and Kondratyuk (GK) available in Refs. [7,13] because this prescription does not contain any contribution from the helicity amplitude I 0,0 (t) showing the prescription free from the zero-mode contributions. In this work, we therefore use the GK prescription to calculate EM FFs.
Following the GK prescription, one can define the three EM FFs, the charge FF G C (t), the magnetic FF G M (t) and the quadrupole FF G Q (t), in terms of the four helicity amplitudes as The charge root-mean-squared (r.m.s.) radius r 2 , magnetic moment µ and the quadrupole moment Q are defined by [9] with normalization G C (t = 0) = 1. Note that heavy quarkonium is charge symmetric, thus the total charge of the system is zero. We therefore calculate the form factors by considering only "the quark" contribution. Although this case is fictitious, it is well-defined and can be compared with related spin-one-hadron theoretical work. It is noted that, in LFWF representation [52], the r.m.s. radii can also be related to the impact parameter b ⊥ ≡ (1 − x)r ⊥ [29] by r 2 = (3/2) b 2 ⊥ [21,22]. One can define a total of nine real GPDs for the spin-one hadrons through the non-local matrix elements of the (axial) vector current on the LF. Five of them are computed from the non-local matrix elements of the same current operator (plus component) which is used as a local operator in Eq. (1) whereas the remaining four are computed from that of the axial current [30,35]. Although there are nine non-local matrix elements that can be computed from the plus component of the current, only five of them are linearly independent because of the constraints from parity invariance. Thus, there are five real GPDs that can be calculated from the five linearly independent non-local matrix elements. In this paper, we only present the (pseudo) scalar and vector meson GPDs that are computed from the current with no quark helicity flip because the meson GPDs with no helicity flip are the ones most readily compared with phenomenological applications [30]. It is however straightforward to calculate helicity-flip GPDs using the same method that is used for helicity-non-flip GPDs.
The five vector meson GPDs for the spin-one hadron are defined through the non-local matrix elements of the vector current on the LF as [30,35] Here, n = (1, 0, 0, 1) is a null vector perpendicular to the light front direction. We choose Ji's convention [53] to define arguments x, ξ and t of the GPDs H i , i = 1, 2, . . . , 5, where x is the momentum fraction carried by the quark in the longitudinal direction and ξ is the skewness parameter. In this work, we choose the Drell-Yan frame ∆ + = 0, or It is straightforward to extract the five GPDs in terms of the five linearly independent non-local matrix elements using Eq. (2) in Eq. (10). Due to the time reversal symmetry on matrix elements V m J ,m J (x, 0, t), one can write [30,35], and we therefore choose those independent non-local matrix elements to be Thus, the expressions for the GPDs read Note our expressions are consistent with those from Ref. [35] in the limit ξ = 0. It is interesting to observe that the integrations of H 4 (x, 0, t) and H 5 (x, 0, t) over x do not correspond to F i (t) of the local current [see Eq. (1)] and therefore vanish. This arises from the time reversal constraints in the case of H 4 (x, 0, t). In the case of H 5 (x, 0, t), this arises because of the term n µ n ν /(P · n) 2 whose analog is absent in the decomposition of the local current [Eq. (1)] as a consequence of Lorentz-invariance [30,35]. Here, we point out that the right-hand side of Eq. (15), after integrating over x, is widely known and cited in the spin-one FF calculations as the angular condition [12,13]. Thus, the first moments of the GPDs can be related to F i (t) for the spin-one hadrons by the first set of sum rules on the LF as [30,35] stress tensor decomposition) as defined in Refs. [26,54].
In the Drell-Yan frame, within the impulse approximation, the helicity amplitudes I m J ,m J (t) and the non-local matrix elements V m J ,m J (x, 0, t) in the region 0 ≤ x ≤ 1 can be written as overlap integrals between LFWFs. The expression for V m J ,m J (x, 0, t) reads [28,33] and that for I m J ,m J (t) reads [11,23,55] where k ⊥ and k ⊥ = k ⊥ +(1−x)∆ ⊥ are the respective relative transverse momenta of the quark before and after being struck by the virtual photon, λ q (λq) is the helicity of the quark (anti-quark). Note that integrating V m J ,m J (x, 0, t) over x yields the local matrix elements (helicity amplitudes) I m J ,m J (t). Here, the LFWFs are truncated to only the valence Fock sector. The valence sector LFWF is normalized according to [20] λq,λq Note that this convention for normalization is introduced in Ref. [20] and the non-perturbative solutions of LFWFs are generated accordingly.

A. Basis Light Front Quantization (BLFQ)
A recent study of heavy quarkonia [20], in a LF Hamiltonian approach [39], presents the effective Hamiltonian based, in part, on the LF holographic QCD [56] as where m q is the mass of the quark. V T is the "soft-wall" light front holography in the transverse direction and is defined as where ζ is holographic variable [56], κ is the confining strength, and r ⊥ is the transverse separation between the quark and the anti-quark. The longitudinal confining potential reads V g is the one-gluon exchange term and in the momentum space, it reads [20,40] where C F = 4/3 is the color factor for the color singlet state, α s is the fixed coupling constant, and 2 is the average momentum squared carried by the exchanged gluon.
In the BLFQ approach, if quarkonium is described by state vectors |ψ J m J , the eigenvalue equations can be defined by and solved non-perturbatively to obtain eigenfunctions that represent the LFWFs ψ J m J (k ⊥ , x, λ q , λq) for heavy quarkonium. To solve Eq. (25), two functions φ nm and χ l are adopted to form the basis in which to evaluate the Hamiltonian matrix. In the transverse direction, 2-dimensional (2D) harmonic oscillator (HO) functions are adopted and are de- where v = |v ⊥ |, θ = arg v ⊥ , n and m are the radial and angular quantum numbers, L |m| n (z) is the associated Laguerre polynomial and b is the HO basis scale with dimension of mass. In the longitudinal direction, the basis functions are defined by where P (α,β) l (z) is the Jacobi polynomial, α = β = 4m 2 q /κ 2 are dimensionless basis parameters, and we drop α and β from the arguments of χ l hereafter.
Using Eqs. (26) and (27) as basis functions, the expansion of momentum-space LFWFs reads [20,41] where n, m, l, λ q , λq|ψ J m J are the LFWFs in the BLFQ basis, obtained by diagonalizing the truncated Hamiltonian matrix [20]. The following truncation is applied to restrict the quantum numbers.
It is clear from the truncation that L max controls the basis resolution in the longitudinal direction whereas N max controls the transverse momentum covered by 2D-HO functions. In the BLFQ approach, the total angular momentum J is only an approximate quantum number due to the breaking of the rotational symmetry by the Fock sector truncation and the basis truncation. However, the total angular momentum projection (m J ) for the system is conserved.
Inserting Eq. (28) in Eqs. (19) and (18) yields the integral over the product of the two 2D-HO functions with different arguments, and that is simplified using the TM coefficients [61] to reduce it to an integral over one 2D-HO function [41]. Then the integral is calculated numerically. Readers are referred to Refs. [20,41] for further details of the BLFQ approach. we surmise that the gluon-exchange dynamics plays a significant role.

IV. RESULTS AND DISCUSSION
In this section, we present and discuss our results for FFs, associated static moments and GPDs. The details of the Hamiltonian's parameters used in calculations are summarized in Table I. The fixed gluon mass µ g = 0.02 GeV is introduced to regularize the singularity present in Eq. (25) [20]. The convergence study of mass eigenvalues with different µ g keeping N max = L max fixed in Ref. [20] suggested that the mass eigenvalues are well converged with respect to µ g . Therefore, the gluon mass is kept fixed in these calculations. Similarly, the HO basis scale b is chosen to be equal to the confining strength κ at the given N max = L max value and at the fixed gluon mass µ g . Fixed, but flavor-dependent, coupling constants α s are used to produce results presented in this work.
Our masses are obtained from the mass eigenvalue equations for total angular momentum projection m J = 0 at the given N max = L max truncation. An important issue that arises in a LF Hamiltonian approach, such as BLFQ, concerns the relative sign between different eigenstates. In particular, since the relative sign between two states with different m J is not fixed by the diagonalization (though the signs of all basis states are fixed by our basis state conventions), we control the overall sign of each eigenfunction to have positive derivative at the origin in coordinate space.

A. The EM FFs and the associated static moments
In this subsection, we present results for the EM FFs and the associated static moments. We start by presenting the charge FFs G C (t) for (pseudo) scalar and (axial) vector mesons in Fig. 1
-t (GeV 2 )    of the left panel in Fig. 1, η c , χ c0 , η c , η b , χ b0 and η b , Eq. (19) directly produces the charge FFs as G C (t) ≡ I 0,0 (t), whereas for the vector mesons of the right panel in Fig. 1, J/ψ, χ c1 , ψ , Υ, χ b1 and Υ , the GK prescription [Eq. (4)] is used to calculate the charge FFs. The FFs for the radially excited charmonia, η c and ψ , exhibit a tendency to develop a node while the corresponding states in bottomonium show this tendency only at larger values of −t (not shown). Nodes in FFs are common features for excited states in non-relativistic systems.
We present the charge FF results for four selected mesons, η c , η b , J/ψ and Υ at a sequence of N max = L max = 8, 16, and 24 values to gain a perspective on their convergence. On the left panel of Fig. 2, we present the convergence of −t G C (t) for the pseudo scalar mesons, and in the right panel, we present the same observable for the vector mesons.
The results show a good convergence trend over this range of −t as evident by finding that the N max = L max = 24 and N max = L max = 16 results are nearly coincident with each other in contrast with the N max = L max = 8 results presented in Fig. 2. This observed convergence in the FFs is reassuring since the charmonia and bottomonia spectroscopy are also reasonably well converged at N max = L max = 24 [20]. Therefore, we only present our FF and  GPD results calculated with N max = L max = 24. The difference between the N max = L max = 24 and 8 values are presented as our uncertainty estimate.
We now turn our attention to the charge mean squared radii of charmonia and bottomonia calculated both in the BLFQ and the SBL approaches. We note again that the charge mean squared radius is an artificial quantity defined with the neglect of the contribution of the anti-quark to the form factor.  bottomonia. This relative relationship is found for both BLFQ and SBL results as well as for the available CI results.
This observation about the relative radii can be understood simply from the tendency towards the non-relativistic limit with increasing quark mass. It is also noted from Table II that the charge radius of η c is smaller than that of χ c0 both in the BLFQ and SBL approaches, and this relationship is consistent with the Lattice results [3].
We note Tables II and III show significant differences among the results calculated in different formalisms which is reasonable considering the major distinctions among the formalisms. For example, in Ref. [4], the DSE results were calculated describing J/ψ by the solutions of the homogeneous Bethe-Salpeter equations (BSE) in rainbowladder truncation. The DSE results also reflect the adoption of an effective running coupling via one-gluon exchange.
Among the many differences with our BLFQ results we note our use of a fixed coupling. Furthermore, in Refs. [16][17][18], the CI results were calculated using contact interactions within the framework of the DSE and BSE. Despite several differences between the CI and BLFQ approaches, there is however a reasonable agreement among the resulting charge radii for the mesons η c and η b . We also observe that for each meson the radius calculated in the SBL approach is larger than the radius calculated in the BLFQ approach. This observation can be understood from the fact that SBL results are produced by only taking the leading basis function into account, which means that the radius is controlled by the dominant mode and by the confining length scale, while the BLFQ includes the gluon exchange, an attractive interaction.
Next, we present the magnetic FFs G M (t) and the quadrupole FFs G Q (t) of vector mesons calculated with N max = L max = 24 in the BLFQ approach. Figure 3   from either of these amplitudes to the magnetic and quadrupole moments, and that is because only the leading basis function contribution(s) of the LFWFs is (are) taken into account. The BLFQ magnetic moments for the mesons J/ψ, Υ and Υ are below 2.0 while the results from the cited literature are above 2.0 where available. This led us to make additional checks of our calculations to confirm the accuracy of our results. It is interesting to note that theoretical results for the rho meson are often below 2.0 as well [57][58][59][60]. For example, in Ref. [59], the investigation in the framework of a covariant extension of the LF formalism has found the rho meson magnetic moment to be 1.83.
Another investigation in the LF quark model [60] has found it to be 1.92 and an investigation in the framework of QCD sum rules [58] has found it to be 1.5±0.3.   1.
x H(x, 0, t) Now, before we close this subsection we present comparisons of selected EM FFs calculated in the SBL and BLFQ approaches. However, we first note that there is likely to be a dominant effect from the difference in the r.m.s. radii between these two approaches. To reduce the impact of this simple difference, we can scale the momentum transfer variable by the appropriate ratio of the charge radii. For this comparison, we select the EM FFs for J/ψ and Υ.
Following the logic for scaling the momentum transfer for the SBL results, the −t values of the SBL charge FF of a vector meson has been scaled so that its slope equals to that of corresponding quantity calculated in the BLFQ approach while keeping both quantities at t = 0 fixed. Then, the −t values of the SBL magnetic and quadrupole FFs are multiplied by the same factor that sets the slopes of the SBL and BLFQ charge FFs for the given meson. The scale factor applied to the SBL results for −t is found to be 1.71 for the case of J/ψ and 1.84 for the case of Υ. Figure 6 presents the resulting comparisons of the J/ψ EM FFs (left panel) and the Υ EM FFs (right panel). The BLFQ magnetic and quadrupole FFs in Fig. 6 are very similar to the corresponding scaled SBL quantities with deviations becoming somewhat visible above approximately −t = 2.4 GeV 2 in the case of J/ψ and above approximately −t = 6 GeV 2 in the case of Υ. This suggests that the dominant role of gluon exchange dynamics for these form factors is a re-scaling of the size of the system from the size dictated by the confinement scale.

B. Generalized parton distributions
In this subsection, we present GPDs for a selection of heavy quarkonia starting with the (pseudo) scalar GPDs. For the (pseudo) scalar mesons such as η c , χ c0 , η c , η b , χ b0 , and η b , Eq. (18) directly produces the GPDs, H(x, ξ = 0, t) ≡ V 0,0 (x, 0, t) . In the previous paper [41], we have presented 3D plots of (pseudo) scalar GPDs of positronium with the one photon exchange (where the longitudinal confining term in the Hamiltonian is absent, of course). We first present the (pseudo) scalar GPDs at fixed |t| to observe the x-dependence of the GPDs in the non-zero momentum transfer limit. Figure.  Let us turn our attention to vector meson GPDs for the ground state identified as 1 3 S 1 of heavy quarkonia 1 in the BLFQ approach with N max = L max = 24. Each vector meson has five GPDs, but we have seen from Eq. (14) that only four of them are non-zero on the LF. Figures 8 and 9 show the non-zero GPDs for J/ψ and Υ, respectively.
It is interesting to note that for both J/ψ and Υ, the integration of GPD H 5 (x, 0, t) over x does not vanish except at |t| = 0 which contradicts the consequence of Lorentz-invariance [see Eq. (17) and associated text]. As we have pointed out before, the x-integration of GPD H 5 (x, 0, t) is the angular condition, the same term that is widely used in spin-one FF on the LF. It is not surprising that due to Lorentz symmetry breaking, the angular condition is not satisfied since we obtain a non-vanishing result for the x-integration of GPD H 5 (x, 0, t) except at the forward limit    (15)] for ψ (2 3 S1) in the BLFQ approach.
(|t| = 0). Repairing this deficiency requires the inclusion of higher Fock sectors in BLFQ, which is a subject for future research. The x-dependence of the ground-state vector meson GPDs is comparable with the corresponding quantities presented in Ref. [32]. The work in Ref. [32] presents the GPDs for the charged ρ meson in a light front constituent quark model. Although the ρ meson is light compared to the heavy quarkonia, the peak somewhere between x = 0.4 to x = 0.6 is the feature that both results have in common.
Next, we present GPDs for the radially excited meson state identified as 2 3 S 1 (ψ and Υ ). Figures 10 and 11 show the non-zero vector meson GPDs for ψ and Υ , respectively. We again comment that, for both ψ and Υ , the integration of GPD H 5 (x, 0, t) over x does not vanish except at |t| = 0 similar to what we found for the case of ground-state vector meson. We note that the decaying trend of the vector meson GPDs (Figs. 10, 11) is more rapid with increasing |t| for the radially excited state compared to the corresponding GPDs for the ground state (Figs. 8,9). This trend is consistent with the fact that ψ (Υ ) has narrower radial extension in momentum space compared to that of J/ψ (Υ) . Similarly, these differences correlate with the relative sizes of these mesons as seen in the results of Table III. That is, larger charge r.m.s. radii correlate with smaller spread in momentum space, as expected. Furthermore, in the forward limit t = 0, the x-dependence of the vector meson GPDs changes character significantly for the radially excited states compared to that of the corresponding ground states. This observation is useful as the x-dependence of the GPDs in the forward limit is directly connected to the partonic interpretation of the hadronic spin [26,30,63].
From the GPDs presented in this work it can be observed the decaying trend of the vector meson GPDs with x is rapid for bottomonia compared to the GPDs for their counterparts in charmonia, and this trend can be understood from considering the relative proximity to the non-relativistic limit where we expect that increasing quark mass leads to a sharper peak in x. Similarly, the rapid fall-off trend of the GPDs with x for heavy quarkonia reflects the consequence of the impulse approximation. This follows the notion that the single quark cannot account for very large longitudinal momentum fraction for equal mass quark constituents. This observation is consistent with properties of the deuteron vector GPDs available in Ref. [35].
Note the vector meson GPDs, investigated in this work, play important roles in various applications. The second moment of GPD H 2 (x, 0, 0) gives the spin-one angular momentum via a sum rule [26,54], and GPD H 5 (x, 0, 0) is equal to b 1 (x), the Deep Inelastic Scattering (DIS) structure function, for the spin-one target such as the deuteron [26,30,34,64]. There is growing interest in these quantities since the announcements of the experimental measurements on b 1 (x) from HERMES [37,63]. Our GPD results, presented in the off-forward limit in this work provide insight to further investigate angular momentum and the structure functions for the spin-one target within BLFQ.

V. SUMMARY AND OUTLOOK
We have calculated the EM FFs for a selection of heavy quarkonia. We have compared the charge radii, the magnetic moments and quadrupole moments calculated in both BLFQ and SBL approaches with the results from other approaches available in the literature. The differences between the BLFQ and SBL results for selected mesons highlight the dynamics of the internal structure of heavy quarkonia. We have also studied the convergence of BLFQ