Spin-orbit correlation of quarks within quarkonium

Understanding the distribution of spin and orbital angular momentum (OAM) within hadrons remains a central challenge in quantum chromodynamics (QCD) [1]. Motivated largely by the longstanding proton spin puzzle [2], decades of experimental campaigns have established that quark and gluon OAM are essential, yet elusive, components of the hadronic spin sum rule [3−9]. As the field advances toward precision tomographic mapping at the future Electron-Ion Collider (EIC) through generalized parton distributions (GPDs) and transverse-momentum-dependent distributions (TMDs) [10, 11], it is equally important to develop clear theoretical observables that isolate specific spin-orbit dynamics. Extending these theoretical investigations beyond the proton to systems such as heavy quarkonia provides a particularly clean environment for investigating the underlying quantum many-body dynamics and fundamental spin-orbit correlations.

A critical step toward understanding hadronic spin lies in the energy-momentum tensor (EMT) $ T^{\mu\nu} $, which rigorously encodes the spatial distribution of energy, momentum, and stress within the hadron [12, 13]. Although the standard parity-even EMT is foundational for extracting total angular momentum through Ji's sum rule [14], decomposing this total $ J^{q,g} $ into distinct spin and OAM contributions remains theoretically subtle because of scheme and gauge dependencies [15, 16]. To circumvent these ambiguities and directly probe the interplay between spin and spatial dynamics, one can investigate the parity-odd (P-odd) EMT. Unlike its P-even counterpart, the P-odd EMT defines a robust spin-orbit correlation (SOC) observable that cleanly isolates the quantum entanglement between spin and orbital motion [17]. Evaluating this correlation naturally connects formal field-theoretical operators with nonperturbative tools, such as light-front wave function decompositions [18] and Wigner distributions [19].

At a fundamental level, SOCs $ \vec L \cdot \vec S $ are ubiquitous across physical systems, from atomic fine structure to nuclear shell models, and play a similarly important role in QCD. Remarkably, even in spin-0 hadrons such as the pion or $ B_c $, a nontrivial spin-orbit structure can arise through the correlation $ \langle L_z S_z\rangle $. Despite the vanishing total angular momentum of these states, this correlation encodes rich internal partonic dynamics [19−27]. This concept was first articulated by Lorcé for the nucleon [28, 29] and subsequently extended to spin-0 systems [30], motivating explicit calculations for both mesons [31−33] and the proton [23]. Formally, this SOC is intimately connected to the aforementioned P-odd EMT, defined as $ \widetilde{T}^{\mu\nu} = T_R^{\mu\nu} - T_L^{\mu\nu} $, where $ T_{L,R}^{\mu\nu} $ represent the EMTs associated with left- and right-handed quark fields, respectively. Arising from the chiral separation of the quark EMT, the matrix elements of $ \widetilde{T}^{\mu\nu} $ between hadronic states directly map the spatial distribution of the spin-orbit coupling, thereby linking it to the anomalous tensor charge and chiral-odd gravitational form factors (FFs) [28].

At the parton level, the SOC is encoded in twist-2 Wandzura-Wilczek contributions together with genuine twist-3 structures [34]. These correlations are central to novel momentum sum rules, analogous to the Jaffe-Manohar spin sum rule, that rigorously decompose hadronic energy into kinetic contributions and potential energies driven by color Lorentz forces [34]. Furthermore, in the small-x limit, the SOC exhibits a universal anti-alignment behavior: $ C_q(x) \approx -\frac{1}{2}q(x) $ and $ C_g(x) \approx -g(x) $. This striking asymptotic relation highlights a state of maximal quantum entanglement between partonic spin and spatial degrees of freedom in this regime [24, 25, 35].

In this paper, we investigate spin-orbit coupling in heavy quarkonia, particularly charmonium ($ c\bar c $) and $ B_c $ ($ b\bar c $) states, as ideal theoretical laboratories for probing hadronic spin structure [1]. Quarkonia provide a controlled setting: at leading order, sea-quark and gluon contributions are suppressed, enabling a clean separation of relativistic and nonrelativistic dynamics. Moreover, their mass scale permits complementary treatments using lattice QCD [36], effective field theories such as potential nonrelativistic QCD (pNRQCD) [37], and light-front Hamiltonian approaches [38−42].

This work has two primary objectives. First, we derive a quantum many-body expression for the quark SOC density, $ {\cal{C}}_z^q(r_\perp) $, within the framework of light-front dynamics. We further decompose the SOC along the longitudinal direction to establish a rigorous connection to its partonic interpretation, $ C_z^q(x) $. Second, we explicitly compute $ {\cal{C}}_z^q(r_\perp) $ for charmonia and $ B_c $ states using light-front wave functions obtained via basis light-front quantization (BLFQ, see Ref. [43] for a recent review). These evaluations yield quantitative predictions for future experiments and provide theoretical insights relevant to the ongoing debate surrounding spin-orbit correlations in spin-0 hadrons.

The remainder of the article is organized as follows. Section II connects the field-theoretical SOC operator $ \hat C_z^q $, defined via the P-odd EMT, to the canonical quantum many-body correlation $ \sum\limits_i S_i^z L_i^z $. In Sec. III, we analyze the covariant structure of the P-odd EMT hadronic matrix elements and demonstrate the proper extraction of both the SOC and the relevant physical FFs. Section IV establishes the partonic interpretation of the SOC. We then present our explicit numerical calculations for heavy quarkonia in Sec. V, before summarizing our conclusions in Sec. VI.

Within quantum many-body theory, the SOC is defined as

where $ \vec S $ and $ \vec L $ denote the intrinsic spin and OAM, respectively. The quantum many-body state $ |\Psi \rangle $ can be represented explicitly by the wave function $ \Psi(\vec r_1, \vec r_2, \cdots, \vec r_n) $. For example, in the hydrogen atom, $ C_{LS} \approx 0 $ for the 1S ground state, whereas $ C_{LS} \approx -1 $ for the $ {\rm{P}}_{1/2} $ state. The local spatial density of the SOC is formally obtained by inserting a Dirac delta function:

Hadrons are inherently relativistic QCD bound states. To extend the concept of SOC to QCD, we employ the light-front Hamiltonian formalism, in which the complete quantum information of the system is encoded in the light-front wave functions (LFWFs). In a relativistic system, the spin-orbit operator $ \vec L\cdot\vec S $ is generally not conserved. However, in the infinite-momentum frame, the transverse components of $ \vec L $ and $ \vec S $ are suppressed, while the longitudinal components remain invariant under longitudinal boosts. Consequently, it is sufficient to evaluate the SOC solely along the z-direction on the light front:

Here, $ \eta_j $ represents the intrinsic parity of the j-th parton, with $ \eta_j = \pm 1 $ for quarks and antiquarks, respectively; $ s_j $ is the spin projection, and $ L_{j}^z = -\hat z \cdot (\vec r_{j\perp} \times i\overleftrightarrow \nabla_{j\perp})/2 $ is the longitudinal OAM operator, with $ \overleftrightarrow \nabla = \vec\nabla - \overleftarrow \nabla $. The n-body LFWF, denoted by $ \widetilde\psi_n(\{x_i,\vec{r}_{i\perp}, s_i\}) $, depends on the transverse coordinates $ \vec r_{i\perp} $ and the longitudinal momentum fractions $ x_i = p^+_i/P^+ $, where $ p_i^+ $ and $ P^+ $ are the longitudinal momenta of the i-th parton and the hadron, respectively. The light-front coordinates are defined as $ x^\pm = x^0 \pm x^3 $ and $ \vec x_\perp = (x^1, x^2) $. The integration measure $ [{{\mathrm{d}}} x_i{{\mathrm{d}}}^2 r_{i\perp}]_n $ covers the appropriate light-front phase space, as detailed in Appendix A. Thus, Eq. (4) provides a direct LFWF representation of the SOC.

As in the nonrelativistic limit, transverse boosts are Galilean and strictly kinematical on the light front. This crucial property allows us to localize the SOC precisely in the transverse plane:

Integrating this density over the transverse coordinates yields the total SOC.

This integrated formulation, however, does not yet determine the distribution of $ C_{LS}^z $ along the longitudinal axis.

In connecting this framework to field theory, the quark contribution to the SOC is closely related to the P-odd EMT, as discussed in Ref. [28]:

where the parity-odd kinematic EMT $ \widetilde T^{\mu\nu}_q(x) $ is defined as

Here, $ \overleftrightarrow D^\mu = \vec\partial^\mu - {{{\overleftarrow{{{\partial^\mu}}}}}} - 2ig_sA^\mu $ is the covariant derivative, and $ T^{\mu\nu}_{R,L} =(1/4)\overline q(x) \gamma^{\mu} (1 \pm \gamma_5) i \overleftrightarrow {D}^{\nu} q(x) $ are the corresponding chiral projections, i.e., $ T^{\mu\nu} = T^{\mu\nu}_R + T^{\mu\nu}_L $. In the present calculation, we restrict ourselves to the valence $q\bar q$ Fock sector and do not include explicit gluon degrees of freedom. It follows directly that the hadronic expectation value $ C_z^q \equiv \langle \Psi | \hat C_z^q |\Psi\rangle $ obeys the same LFWF representation derived for $ 2C_{LS}^z $ in Eq. (6) when the quark and antiquark contributions are summed, namely:

The relative sign encoded by $ \eta_j $ follows from the charge-conjugation property of $ \hat C_z^q $, $ {\cal{C}} \hat C_z^q {\cal{C}}^{-1}=-\hat C_z^q $, where $ {\cal{C}} $ is the charge-conjugation operator. Consequently, the full flavor-current matrix element vanishes for a hadron with charge symmetry.

To evaluate the spatial distribution of this correlation, we define the quark SOC density analogously to the integral operator in Eq. (8):

where the light-front transverse density of the P-odd EMT in the Drell–Yan frame ($ p'^+ = p^+ $) is given by

Physically, the Drell-Yan kinematic condition $ \Delta^+ \equiv p'^+ - p^+ = 0 $ is equivalent to integrating over the longitudinal coordinate $ x^- $. The factor of $ 2p^+ $ arises naturally from the normalization of the hadronic states, $ \langle p' | p\rangle = 2p^+(2\pi)^3 \delta^3(p-p') $. By construction, $ {\cal{C}}_z^q(r_\perp) $ satisfies the LFWF representation established for $ 2{\cal{C}}_{LS}^z(r_\perp) $.

Thus, this transverse density is the two-dimensional Fourier transform of the hadronic matrix element $ C_z^q(\Delta_\perp^2) $ evaluated in the Drell-Yan frame:

The underlying hadronic matrix element is extracted using

with the average momentum $ P = (p'+p)/2 $ and momentum transfer $ \Delta = p'-p $. In the Drell–Yan frame, $ \Delta_\perp^2 \equiv \vec \Delta_\perp^2 = -\Delta^2 $. The total quark SOC is then obtained from the forward limit, $ C_z^q = C_z^q(0) $.

Following Eq. (8), the SOC is directly determined by the hadronic matrix element of the P-odd EMT, $ \langle p' | \widetilde T^{\mu\nu}_q |p \rangle $. By Lorentz invariance, this matrix element can be parameterized in terms of fundamental P-odd covariant structures. For spin-0 hadrons, such as the scalar quarkonia studied in this work, the Lorentz decomposition takes the form [30],

Here, $ P = (p'+p)/2 $, $ \Delta = p'-p $, $ t = -\Delta^2 $, $ \varepsilon^{\mu\nu\Delta P} = \varepsilon^{\mu\nu\rho\sigma} \Delta_\rho P_\sigma $, and $ \varepsilon^{\mu\nu\rho\sigma} $ is the Levi-Civita tensor. We use the convention $ \varepsilon^{0123} = - \varepsilon_{0123} = +1 $, which corresponds to $ \varepsilon^{-+12} = +2 $ in light-front coordinates. Note that $ \varepsilon^{\mu\nu\Delta P} $ is antisymmetric.

In their study of pionic SOC [31], Tan and Lyu augmented the standard Lorentz tensor $ i \varepsilon^{\mu\nu\Delta P} $ with an additional P-odd tensor structure, $ {P^{[\mu}\varepsilon^{\nu]+\Delta P}} $, a form originally proposed by Lorcé for the nucleon [28]. Because this term explicitly depends on the light-front “$ + $” direction, it violates manifest Lorentz covariance for spin-0 hadrons, as recently emphasized by Lorcé and Song [30]. Nevertheless, such symmetry-breaking terms often arise in practical calculations employing LFWFs based on finite Hamiltonian truncations, including Fock-space truncations ¹. Under such truncations, the full Lorentz group is explicitly reduced to the kinematic subgroup of the light front. To rigorously analyze and isolate these truncation-induced structures, one can employ covariant light-front dynamics (CLFD) [44, 45]. The CLFD framework explicitly tracks this frame dependence by introducing the light-front null vector $ \omega^\mu = (1, 0, 0, -1) $, which specifies the orientation of light-front time ($ \omega_\mu \omega^\mu = 0 $, $ \omega_\mu a^\mu = a^+ $). It is particularly convenient to construct the scaled null vector $ N^\mu = M^2\omega^\mu/(\omega \cdot P) $, where $ M^2 = p^2 = p'^2 $ denotes the squared hadron mass. This vector has the appropriate mass dimension and remains manifestly invariant under arbitrary longitudinal rescaling, $ \omega^\mu \to \lambda \omega^\mu $.

Incorporating this explicit ω-dependence, the generalized covariant decomposition takes the form:

where $ a^{[\mu}b^{\nu]} = (a^\mu b^\nu - a^\nu b^\mu)/2 $ and $ a^{\{\mu}b^{\nu\}} = (a^\mu b^\nu + a^\nu b^\mu)/2 $. Consequently, spurious FFs arise from the explicit dependence on the light-front orientation $ N^\mu $, which violates manifest Lorentz covariance. Specifically, the FFs $ \widetilde A_{1,2}(t) $ in the antisymmetric part of $ \widetilde T^{\mu\nu}_q $, together with $ \widetilde S_{1,2,3}(t) $ in its symmetric part, are entirely spurious. These artifacts must vanish when full Lorentz symmetry is restored. Furthermore, the alternative P-odd structure proposed in Ref. [31] is mathematically related to the basis structures identified above through the Schouten identities.

The physical FF $ \widetilde F_q(t) $ can be extracted from the antisymmetric part of the transverse components of $ \widetilde T^{+i} $. Evaluated in the symmetric Drell-Yan frame ($ \Delta^+=0, \vec P_\perp = 0 $), this expression is

Here, $ \epsilon^{ij}_T $ denotes the antisymmetric tensor in the transverse plane, with $ \epsilon^{12}_T = -\epsilon^{21}_T = +1 $. In contrast, the symmetric part of $ \widetilde T^{+i}_q $ need not vanish identically for model wave functions and is directly related to the spurious FF $ \widetilde S_1 $:

The complete set of hadronic matrix elements is collected in Appendix B. Returning to the SOC, the FF $ C_z^q(\Delta^2_\perp) $ is defined directly in terms of the hadronic matrix element of the full transverse component $ \widetilde T^{+i}_q = \widetilde T^{[+i]}_q + \widetilde T^{\{+i\}}_q $ and therefore includes the spurious contribution:

Lorcé and Song showed that the FF $ \widetilde F_q(t) $ is related to two additional FFs [30],

a relation that follows directly from the QCD operator identity [28]

Here, the FFs $ F^q(t) $ and $ H^q(t) $ are defined in terms of local hadronic matrix elements,

Since $ \widetilde F_q $ is extracted solely from the antisymmetric part of the hadronic matrix element, as in Eq. (20), this underlying relation is naturally preserved even in the presence of spurious contributions. In our model calculation, the matrix-element relation implied by Eq. (24) is therefore maintained. When explicit gauge-field contributions are omitted, it reduces at the one-body quark level to the Dirac spinor identity

This identity follows directly from the on-shell Dirac equations for the constituent spinors.

Having established a rigorous procedure for isolating the physical FF $ \widetilde F_q $ from the macroscopic EMT using CLFD, we now turn to the underlying microscopic picture. In the next section, we examine the partonic interpretation of the SOC, linking this formal matrix element to the internal dynamics of the bound state through Wigner distributions and GPDs.

The SOC has a natural partonic interpretation in terms of polarized Wigner distributions [28],

where $ \widetilde f_q $ denotes the polarized quark Wigner distribution:

In this expression, $ W_\pm $ denotes a Wilson line connecting the two spacetime points $ \pm z/2 $. The choice of the Wilson-line path determines which definition of the SOC is obtained. A straight Wilson line corresponds to the kinetic SOC in Ji's spin decomposition, whereas a staple-shaped Wilson line corresponds to the canonical SOC in the Jaffe-Manohar decomposition [28]. In the present calculation, the Wilson line reduces to the identity. Consequently, the path dependence disappears, and the two definitions coincide.

Retaining the explicit dependence on the impact parameter $ \vec b_\perp $ in Eq. (28) yields the SOC transverse density $ {\cal{C}}_z^q(b_\perp) $ introduced in Eq. (12). Similarly, integrating over the transverse phase space while leaving the longitudinal momentum fraction x unintegrated yields the parton distribution of the SOC in the longitudinal direction:

Retaining the impact parameter $ \vec b_\perp $ during this integration naturally leads to the definition of the generalized parton distribution:

Although the above definition is intuitive, it requires evaluating a gauge link in the transverse direction. Alternatively, the purely longitudinal distribution $ C_z^q(x) $ can be determined within the GPD formalism in the valence Fock sector [30, 34]. From Eqs. (20–23), $ C_z^q(x) $ is directly related to the leading-twist GPDs $ {\cal{H}}^q(x, \xi, t) $ and $ {\cal{H}}^q_1(x, \xi, t) $, as well as to the twist-3 GPD $ {\cal{G}}_2^q(x, \xi, t) $:

where the GPDs are defined in terms of bilocal light-cone operators:

Using these definitions, the LFWF representation of $ C_z^q(x) $ can be obtained straightforwardly. Moreover, when the spurious contributions vanish, the twist-3 GPD $ {\cal{G}}_2^q $ is rigorously related to the leading-twist GPDs $ {\cal{H}}^q $ and $ {\cal{H}}_1^q $ through the QCD equations of motion. Under this condition, $ C_z^q(x) $ depends exclusively on $ {\cal{G}}_2^q $, as first demonstrated in Ref. [30]. The integral of the distribution $ C_z^q(x) $ gives the hadron's SOC.

In this section, we evaluate the SOC for charmonia and $ B_c $ ($ b\bar c $) states. The $ B_c(nS) $ and $ B_c(nP) $ states considered here correspond to spin-zero states with $ J^P=0^- $ and $ J^P=0^+ $, respectively. We evaluate the SOC for both spin-zero and spin-one mesons using the many-body expression in Eq. (4). Because the valence Fock sector dominates the wave functions of heavy quarkonia [46, 47], we evaluate the observables within this leading sector. Truncating Eq. (4) accordingly, the valence contribution to the SOC takes the explicit form:

where x and $ \vec r_\perp $ denote the longitudinal momentum fraction and the quark's relative transverse coordinate, respectively. Because the relevant operator $ \hat C_z^q $ is odd under charge conjugation, its full matrix element, including both quark and antiquark contributions, vanishes identically for charge-symmetric states such as charmonia. Consequently, for these states, our analysis focuses exclusively on the isolated quark (or antiquark) contributions. By contrast, for asymmetric systems such as the $ B_c $ meson, the total SOC is generally nonzero. In the limit $ x\ll 1 $, the valence expression in our model asymptotically reduces to:

We adopt heavy quarkonium wave functions obtained from the BLFQ approach [39, 40]. In this framework, the wave functions are determined by diagonalizing an effective light-front Hamiltonian in a holographic basis. Practical calculations require a basis truncation governed by the cutoff parameter $ N_\text{max} $. This truncation determines the effective ultraviolet (UV) resolution through $ \Lambda_\text{UV}\sim\kappa\sqrt{N_\text{max}} $, where $ \kappa\approx 1\ \text{GeV} $ is the basis scale parameter. Guided by our previous studies [39, 40, 48−50], we use $ N_\text{max}=8 $ to compute the central values for charmonium, which corresponds to a UV resolution comparable to the charmonium mass ($ \Lambda_\text{UV}\approx M $). The sensitivity to this basis truncation is estimated from the difference between the $ N_\text{max}=8 $ and $ N_\text{max}=16 $ results. For the $ B_c $ system, we similarly use the $ N_\text{max}=32 $ results as our central values, while the basis sensitivity is estimated as twice the difference between the $ N_\text{max}=32 $ and $ N_\text{max}=24 $ results. The specific values of the BLFQ model parameters used for each hadronic system are given in Refs. [39, 40].

Table I summarizes the charm-quark SOC for several low-lying charmonia, together with the total SOC for the $ B_c $ states. A simple estimate based on the non-relativistic quark model (NRQM) is also presented as a baseline. As shown in the table, the SOC of the quark is consistently negative, indicating that the quark spin and OAM are anti-aligned within the quarkonium system. In contrast, the SOC of the antiquark is positive owing to the explicit negative sign associated with the antiquark in the definition. Consequently, charge conjugation symmetry causes the total SOC for charmonia to vanish exactly.

Furthermore, Table I shows that the BLFQ results for S-wave states are not exactly zero, contrary to expectations from the NRQM. This deviation is a distinct relativistic effect. Specifically, it originates from dynamically generated P-wave mixing, which is forbidden by non-relativistic symmetry rules but naturally accommodated in relativistic dynamics, including light-front dynamics [39, 52]. For P-wave and D-wave mesons, our BLFQ results are of the same order of magnitude as the NRQM estimates, consistent with the predominantly non-relativistic nature of these heavy bound states.

The spatial distribution of the SOC in the transverse plane, $ {\cal{C}}_z^q(r_\perp) $, is shown in Figs. 1–2. Figure 1 compares these transverse densities for selected ground and orbitally excited states, namely $ \eta_c(1S) $, $ \chi_{c0}(1P) $, $ B_c(1S) $, and $ B_c(1P) $. Consistent with physical intuition, the magnitude of the SOC is substantially enhanced in the P-wave systems relative to their S-wave counterparts. Extending this analysis to radial excitations, Fig. 2 shows the corresponding transverse densities for higher radially excited quarkonia. In these states, a clear nodal structure appears in the density, reflecting the sign-changing structure of the radially excited wave functions, together with the competition among different helicity components and, where applicable, the quark and antiquark contributions.

Figure 1.  (color online) Spin–orbit correlation distributions for the ground and excited states $ \eta_c $, $ \chi_{c0} $, $ B_c(1S) $, and $ B_c(1P) $. Panel (a) shows the antiquark contribution alone, whereas panel (b) shows the corresponding quark contribution. The shaded bands indicate the theoretical sensitivity to the basis truncation (see text for details).

Figure 2.  (color online) Transverse distributions of the spin–orbit correlation, $ 2\pi r_\perp {\cal{C}}^q_z(r_\perp) $, for (a) S-wave charmonia ($ \eta_c, \eta_c' $, and $ \eta_c'' $); (b) P-wave charmonia ($ \chi_{c0} $ and $ \chi_{c0}' $); (c) S-wave $ B_c $ mesons ($ B_c(1S), B_c(2S) $); and (d) P-wave scalar $ B_c $ mesons ($ B_c(1P), B_c(2P) $). The shaded bands represent the theoretical uncertainty arising from basis sensitivity (see text for details).

The parton distribution of the SOC along the longitudinal direction, $ C_z^q(x) $, is shown in Fig. 3. For the charge-symmetric charmonium system, the quark and antiquark distributions both peak around $ x=1/2 $. In contrast, the unequal constituent masses in the $ B_c $ meson produce a relative shift in the peaks of the b- and c-quark distributions, resulting in a distinctive zigzag profile. To provide a complete spatial picture, the fully unintegrated three-dimensional mixed distribution of the SOC on the light front is visualized in Fig. 4.

Figure 3.  (color online) Parton distributions of the spin–orbit correlation, $ C^q_z(x) $, for (a) $ \eta_c $; (b) $ \eta'_c $; (c) $ \chi_{c0} $; (d) $ \chi_{c0}' $; (e) $ B_c(1S) $; (f) $ B_c(2S) $; (g) $ B_c(1P) $; and (h) $ B_c(2P) $. The shaded bands indicate the theoretical uncertainty arising from basis sensitivity (see text for details).

Figure 4.  (color online) Mixed distributions of the spin-orbit correlation $ C_z^q(x,r_\perp) $ as functions of x and the impact parameter $ r_\perp $ are shown for (a) $ \eta_c $, (b) $ \eta'_c $, (c) $ \chi_{c0} $, and (d) $ \chi_{c0}' $, and for (e) $ B_c(1S) $, (f) $ B_c(2S) $, (g) $ B_c(1P) $, and (h) $ B_c(2P) $. For the charmonium states in (a–d), only the quark contribution is presented. For the $ B_c $ states in (e–h), the total contribution is shown.

The results for the FF $ \widetilde F_q(Q^2) $ extracted from the antisymmetric tensor in Eq. (20) are shown in Figs. 5–6. The FF $ \widetilde F_q(Q^2) $ for the 1S charmonium state $ \eta_c $ is presented in Fig. 5(a), while the corresponding results for the $ B_c(1\text{S}) $ meson are shown in Fig. 5(b). Figures 5(c) and 5(d) show the results for the P-wave states, $ \chi_{c0} $ and $ B_c(1\text{P}) $, respectively, which are significantly less sensitive to the basis truncation.

Figure 5.  (color online) Calculated form factors (FFs) $ \widetilde F_q(Q^2) $ for selected S-wave and P-wave systems. The panels show results for (a) the S-wave charmonium state $ \eta_c $, (b) the S-wave $ B_c $ meson, (c) the P-wave charmonium state $ \chi_{c0} $, and (d) the P-wave $ B_c $ meson. The shaded bands indicate the theoretical sensitivity to basis truncation (see text).

Figure 6.  (color online) FF $ \widetilde F_q(Q^2) $ evaluated for selected radial excitations. The panels are organized by system and orbital angular momentum: (a) S-wave charmonia ($ \eta_c $, $ \eta'_c $, and $ \eta''_c $), (b) P-wave charmonia ($ \chi_{c0} $ and $ \chi_{c0}' $), (c) S-wave $ B_c $ mesons ($ B_c(1S) $ and $ B_c(2S) $), and (d) P-wave $ B_c $ mesons ($ B_c(1P) $ and $ B_c(2P) $). The shaded bands indicate the theoretical uncertainty associated with sensitivity to basis truncation (see text).

In the covariant limit, one expects $ \widetilde F_q(0) = C_z^q $. Although the P-wave states approximately satisfy this relation, the extracted $ \widetilde F_q(0) $ for the S-wave states is an order of magnitude larger than $ C_z^q $, as shown in Table II. This discrepancy indicates a substantial spurious contribution in the S-wave states, suggesting that the small nonzero $ C_z^q $ is strongly affected by truncation effects. Consequently, the magnitude of $ C_z^q $ may carry a sizable systematic uncertainty as an estimate of the physical SOC.

Figure 6 compares the quark FFs $ \widetilde F_q(Q^2) $ across different radial excitations. Specifically, Fig. 6(a) shows $ \widetilde F_q(Q^2) $ for the S-wave charmonia ($ \eta_c $, $ \eta_c' $, and $ \eta_c'' $). Figure 6(b) presents the corresponding quark FFs for the P-wave charmonia ($ \chi_{c0} $ and $ \chi_{c0}' $). Finally, the total FF $ \widetilde F_q(Q^2) $ for the various radially excited $ B_c $ states is summarized in Figs. 6(c)–6(d).

In this work, we systematically investigated the spin-orbit correlation (SOC) of quarks using a relativistic quantum many-body approach. By deriving the many-body representation of $ C_z^q $ within the valence Fock sector, we demonstrated that it is equal to twice the longitudinal SOC on the light front, i.e., $ \hat C_z^q = 2 \sum\nolimits_i \eta_i L_i^z S_i^z $. This physical interpretation extends naturally to coordinate space, where the density $ {\cal{C}}_z^q(r_\perp) $ is identified as twice the SOC one-body density. Finally, we provided a rigorous partonic interpretation of the SOC by linking it directly to the formalism of generalized parton distributions (GPDs).

Prior work by Lorcé et al. established a fundamental link between the SOC $ C_z^q $ and the hadronic matrix element of the P-odd energy-momentum tensor (EMT) $ \widetilde T^{\mu\nu}_q $. Nevertheless, tensions persist in the literature regarding the correct procedure for extracting the associated FF $ \widetilde F_q $ in practical calculations. In this work, we resolved this ambiguity by performing a covariant light-front dynamics (CLFD) analysis, which yielded the complete parametrization of the matrix element $ \langle p' | \widetilde T^{\mu\nu}_q(0)| p\rangle $ for hadronic states described by partonic degrees of freedom. We rigorously demonstrated that the genuine physical FF $ \widetilde F_q $ is accessed exclusively through the antisymmetric P-odd EMT component, $ \widetilde T^{[+i]}_q $. Our framework not only generalizes the decomposition introduced by Tan and Lyu but also smoothly recovers the standard covariant decomposition in the appropriate limit. Consequently, this extraction method provides a robust theoretical tool for practical light-front calculations in which manifest covariance is broken.

As a practical application of this formalism, we computed the SOC for charmonia and $ B_c $ mesons. As anticipated, the numerical results are consistent with estimates based on the non-relativistic quark model (NRQM). We also compared the transverse and longitudinal SOC distributions in detail across various S-wave, P-wave, and D-wave quarkonia, including radially excited states, which exhibit distinct and intriguing spatial topologies.

Experimental access to the SOC presents a unique phenomenological challenge. Formally, the most direct pathway to the SOC is through the framework of GPDs. As we have demonstrated, the physical P-odd EMT FF is rigorously connected to specific twist-3 GPDs within the valence Fock sector, providing a profound 3D tomographic map of the spin-orbit correlation in impact-parameter space. For stable hadrons such as the nucleon, such GPDs are systematically probed via deeply virtual Compton scattering (DVCS). However, the unstable nature of heavy quarkonia precludes their use as fixed targets, rendering standard DVCS measurements unfeasible. Nevertheless, this exact mapping to GPDs remains highly advantageous because it tightly links the macroscopic SOC to the underlying microscopic light-front wave functions (LFWFs). Consequently, precision measurements of hard exclusive processes, such as two-photon fusion ($ \gamma^* \gamma^* \to \eta_c, \chi_{c0} $) or exclusive double charmonium production at $ e^+ e^- $ colliders (e.g., BES III, Belle II, and the proposed Super Tau-Charm Factory), can constrain these phenomenological LFWFs and provide a data-driven method for reconstructing quarkonium GPDs and extracting the internal SOC.

As a complement to this spatial GPD picture, the partonic SOC can be comprehensively probed in full phase space through generalized transverse momentum distributions (GTMDs). As the exact quantum-mechanical phase-space distributions of the bound state, GTMDs contain the most direct macroscopic representation of the SOC. Recently, as proposed by Bhattacharya et al., highly promising avenues for directly accessing the proton's SOC via GTMDs at the EIC have been established using exclusive dijet production ($ \gamma^* p \to p' + \text{jet} + \text{jet} $) and semi-inclusive diffractive deep inelastic scattering (SIDDIS) [51].

Our current numerical evaluation is restricted to the valence quark sector, which provides a robust approximation for heavy quarkonia. However, the SOC observable is fundamentally capable of probing the full many-body content of QCD bound states. Consequently, investigating the SOC contributions from sea quarks and gluons remains a primary goal for future study. For scalar mesons, a particularly intriguing paradox arises regarding these gluon contributions. On the one hand, the many-body representation can be directly extended to incorporate gluons, namely, $ \hat C_z^{q,g} = 2\sum\nolimits_{i\in q, g} L_i^z S_i^z $. On the other hand, Lorcé and Song demonstrated the absence of an appropriate gluonic operator for the P-odd EMT, implying that the macroscopic gluonic contribution to scalar mesons must vanish identically. Reconciling this apparent tension between the microscopic many-body picture and the macroscopic EMT formalism will be a critical step toward unraveling the complex gluonic structure of QCD.

This work was supported by the National Natural Science Foundation of China (NSFC) under Grant No. 12375081 and by the Chinese Academy of Sciences under Grant No. YSBR-101. A portion of the computational resources was also provided by the Dongjiang Yuan Intelligent Computing Center.

For an n-body light-front-quantized Fock state

we introduce the relative momenta as,

where $ P^+ = \sum\nolimits_i p_i^+ $ and $ \vec P_\perp = \sum\nolimits_i \vec p_{i\perp} $. The relative momenta satisfy

The light-front wave functions depend only on these relative variables.

The n-body measure is defined as,

where $ S_n $ is the symmetry factor.

The transverse coordinate-space wave function is defined as

Here, we have denoted,

Finally, the n-body measure over the transverse coordinates is defined as follows:

The covariant decomposition of the parity-odd tensor matrix element for spin-0 hadrons is given in Eqs. (17–18):

In the symmetric Drell–Yan frame ($\Delta^+=\Delta^-=0$, $\vec{P}_\perp=0$), the light-front components are:

The transverse Levi-Civita symbol is defined by $ \varepsilon^{12}_T = +1 $.