Introduction. Charm family hadrons, according to the color confinement principle, are believed to be capable of containing one, two, three, four or even more charm quarks. Although the \(J/\psi\) meson composed of a charm-anticharm pair was discovered fifty-one years ago [1], [2] and the D meson composed of a charm-antilight quark two years later [3], it was not until 2020 that the LHCb experiment found the first evidence of a pure charm tetraquark state, the \(X(6900)\), in the \(J/\psi\)-pair invariant mass spectrum—using a data sample corresponding to an integrated luminosity of \(9fb^{-1}\)—thereby providing a signal for a hadron with multiple charm quarks [4].

Using a data sample fifteen times larger but with a focus on the small rapidity region, both the ATLAS and CMS collaborations have confirmed the narrow \(X(6900)\) state and observed another narrow structure around \(7100MeV\) and a broad structure above twice the \(J/\psi\) mass [5], [6].

Several particularly interesting aspects of these exotic states in the \(J/\psi\)-pair spectrum include: (1)the spin-parity quantum numbers for all observed data distributions are found to be compatible with the \(J^{PC}=2^{++}\) hypothesis from the very recent CMS analysis [7], which constrains the possible internal quark configuration for these exotic structures; (2)a Regge trajectory is predicted and observed from their mass spectrum [8], [9], which gives hints of a radial excitation tetraquark family; (3)examination of data from three experiments reveals that the production rate of the fully charm tetraquark state is unaffected by rapidity, yet its signal is strongly enhanced at high transverse momentum, which can be inferred that the dependence on rapidity across low to medium values should be mild, and that the transverse momentum distribution curves for \(J/\psi\)-pair and for the tetraquark state themselves should be distinct.

Addressing and clarifying the above three points requires both a consistent theoretical framework and precise calculations. From the perspective of group theory, the fully charm tetraquark states can always be represented by a complete set of basis states. In terms of color space, one can choose either the color-symmetry-antisymmetry \(\mathbf{6}\otimes \bar{\mathbf{6}}\)/\(\bar{\mathbf{3}}\otimes \mathbf{3}\) basis or the color-singlet–octet \(\mathbf{1}\otimes \mathbf{1}\)/\(\mathbf{8}\otimes \mathbf{8}\) basis. In color-symmetry-antisymmetry basis, the two charm quarks can form the “good” axial-vector diquark with antisymmetric color \(\bar{3}\) and symmetric spin-flavor configuration or “bad” scalar diquark with symmetric color \(6\) and symmetric spin-flavor configuration. The two charm antiquarks do similarly. Both charm diquark and antidiquark carry colour charge and form a colour-charge-neutral and tightly-bound hadron from a QCD confining atrractive force. Interestingly, charm diquark here is distinctly different from light-diquark system, due to the additional constraint of isospin symmetry present in the latter. In color-singlet–octet basis, the two pair of charm quark-antiquark can form a color-singlet charmonium or a color-octet cluster separately, and then the two charmonia or coloured clusters form a colour-charge-neutral and loosely-bound molecule hadron from either a residual strong interaction or a repulsive force within QCD. The two basis expansions are equivalent in principle.

A true physical system may reside in one of the basis states mentioned above, or in a superposition of these basis states. For example, we can expand the fully heavy tetraquark state within the color-symmetry-antisymmetry basis as \[\begin{align} |T_{4Q}\rangle =& \sum_{L=0}C^s_L |(Q_1Q_2)_0;(\bar{Q}_3\bar{Q}_4)_0;L\rangle \nonumber\\&+C^a_L |(Q_1Q_2)_1;(\bar{Q}_3\bar{Q}_4)_1;L\rangle, \end{align}\] where \(C^s_L\) and \(C^a_L\) are the symmetric and antisymmetric configuration coefficients with \(\sum_{L=0}|C^s_L|^2+|C^a_L|^2=1\), respectively. Orbital excitation effects are marked with \(L\). The higher Fock states with gluon content are omitted here but can be included in a similar manner.

In this Letter, we calculate the differential and total cross sections for the hadronic production of S-wave fully charm tetraquark states based on the color-symmetry-antisymmetry basis and nonrelativisitic QCD (NRQCD) effective theory. Our calculations include contributions from all possible configurations, such as \(J^{PC}=0^{++}\) or \(2^{++}\), as well as mixing effects arising from identical quantum numbers. On hadron colliders, gluon radiation effects are non-negligible; therefore, we provide a complete treatment of QCD corrections up to the next-to-leading order (NLO). When the transverse momentum of the fully charm tetraquark state is much lower than its typical momentum scale, we perform resummation at the next-to-leading logarithmic (NLL) accuracy to mitigate large logarithmic terms induced by soft gluon emissions in the low transverse momentum region. These physical results will aid in identifying the quantum numbers and microscopic structures of fully charm tetraquark states, while also facilitating the understanding of the physical mechanisms behind the observed signal enhancement at high transverse momenta.

QCD factorization formulae. The typical Feynman diagrams for fully charm tetraquark hadronic production are plotted in Fig. 1. The hadronic production of fully charm tetraquark can be factorized as \[\begin{align} &\frac{d\sigma}{ dy dP_\perp^2}(p+p\to T_{4c}+X) \nonumber\\&~~= \sum_i \int_{x_A}^1 d\xi_1 f_{i/A} (\xi_1, \mu_F)\sum_j \int_{x_B}^1 d\xi_2 f_{j/B} (\xi_2, \mu_F) \nonumber\\ &~~ \times \int d^2b e^{iP_\perp \cdot b} W(b, \xi_1,\xi_2, M)\hat{\sigma}^{(0)}_{ij} C_{iA} \left( \frac{x_A}{\xi_1},\alpha_s(\mu_F) \right) \nonumber\\ & ~~\times C_{jB} \left( \frac{x_B}{\xi_2}, \alpha_s(\mu_F) \right) + Y(P_\perp, M, x_A, x_B), \end{align}\] where \(y\) and \(P_\perp\) are the rapidity and transverse momentum for fully charm tetraquark. \(f_{i/A}\) and \(f_{j/B}\) are parton distribution functions of parton in the proton. \(W(b, \xi_1,\xi_2, M)\) is the Fourier transformation of the differential cross section at low \(P_\perp\) into the impact parameter b-space considering the resummation contribution. \(C_{iA}\) and \(C_{jB}\) are the matching coefficients. \(Y(P_\perp, M, x_A, x_B)\) is the regular term contributing the medium and large \(P_\perp\) region.

In the factorization formuale, \(\hat{\sigma}^{(0)}_{ij}\) is related to the leading-order (LO) partonic cross section \(\hat{\sigma}^{(LO)}(i+j\to T_{4c})=\hat{\sigma}^{(0)}_{ij}\delta(1-z)\) with \(z=M^2/\hat{s}\) where \(M\) is the fully charm tetraquark mass and \(\hat{s}\) is the center-of-mass energy in partonic frame. We have \[\begin{align} \hat{\sigma}^{(0)}_{ij} = \frac{\pi}{M^4} \frac{1}{c_i c_j s_i s_j} \left| \mathcal{M}^{(0)}_{ij} \right|^2, \end{align}\] where the color and spin average factor are \(c_g=N_c^2-1\), \(c_q=f_{\bar{q}}=N_c\), \(s_g=2-2\epsilon\), and \(s_q=s_{\bar{q}}=2\). The Feynman amplitude squared from certain operator contribution can be written as \[\begin{align} \left| \mathcal{M}^{(n)}_{ij} \right|_{O^J_{g_1,g_2}}^2 = \frac{M\pi ^4 \alpha_s^4}{128 m_c^8} c^{(n)}_{ij,O^J_{g_1,g_2}}\langle0|O^{J}_{g_1,g_2}|0\rangle. \end{align}\] In this expression, \(m_c\) is the charm quark mass. \(O^{J}_{g_1,g_2}\) are the NRQCD long-distance operators for tetraquarks with spin quantum number \(J\) \[\begin{align} O_{g_1,g_2}^{J}=\mathcal{O}_{g_1}^{(J)} \sum_X\left|T_{4 c}^J+X\right\rangle\left\langle T_{4 c}^J+X\right| \mathcal{O}_{g_2}^{(J) \dagger}. \end{align}\] Therein \(g_i\) can be either \(\mathbf{6}\otimes \overline{\mathbf{6}}\) or \(\overline{\mathbf{3}} \otimes \mathbf{3}\) in color symmetry-antisymmetry basis. Then we have [10], [11] \[\begin{align} \mathcal{O}_{\overline{\mathbf{3}} \otimes \mathbf{3}}^{(0)} & =-\frac{1}{\sqrt{3}}\left[\psi_a^T\left(i \sigma^2\right) \sigma^i \psi_b\right]\left[\chi_c^{\dagger} \sigma^i\left(i \sigma^2\right) \chi_d^*\right] \mathcal{C}_{\overline{\mathbf{3}} \otimes \mathbf{3}}^{a b ; c d},\nonumber \\ \mathcal{O}_{\overline{\mathbf{3}} \otimes \mathbf{3}}^{(2;i j )} & =\left[\psi_a^T\left(i \sigma^2\right) \sigma^m \psi_b\right]\left[\chi_c^{\dagger} \sigma^n\left(i \sigma^2\right) \chi_d^*\right] \Gamma^{i j ; m n} \mathcal{C}_{\overline{\mathbf{3}} \otimes \mathbf{3}}^{a b ; c d},\nonumber \\ \mathcal{O}_{\mathbf{6} \otimes \overline{\mathbf{6}}}^{(0)} & =\left[\psi_a^T\left(i \sigma^2\right) \psi_b\right]\left[\chi_c^{\dagger}\left(i \sigma^2\right) \chi_d^*\right] \mathcal{C}_{\mathbf{6} \otimes \overline{\mathbf{6}}}^{a b ; c d}, \end{align}\] where the rank-4 color tensor is defined as \(\mathcal{C}_{\mathbf{6} \otimes \overline{\mathbf{6}}}^{cd;ef} = \frac{1}{2\sqrt{6}}\left(\delta_{ce}\delta_{df}+\delta_{cf}\delta_{de}\right)\) and \(\mathcal{C}_{\overline{\mathbf{3}} \otimes \mathbf{3}}^{cd;ef} = \frac{1}{2\sqrt{3}}\left(\delta_{ce}\delta_{df}-\delta_{cf}\delta_{de}\right)\). The rank-4 Lorentz tensor is defined as \(\Gamma^{\alpha\beta ; \mu\nu}=\frac{1}{2}g^{\mu \alpha}g^{\nu \beta}+\frac{1}{2}g^{\mu \beta}g^{\nu \alpha}-\frac{1}{3}g^{\mu \nu}g^{\alpha\beta}\).

At LO, the short-distance coefficients for gluon-gluon fusion are \[\begin{align} c^{(0)}_{gg,O^J_{g_1,g_2}}=\{\frac{15488}{27},2304,\frac{1408\sqrt{6}}{3}, \frac{200704}{27}\}, \end{align}\] for \(O^J_{g_1,g_2}=O^0_{\mathbf{6}\otimes \overline{\mathbf{6}},\mathbf{6}\otimes \overline{\mathbf{6}}},O^0_{\overline{\mathbf{3}} \otimes \mathbf{3},\overline{\mathbf{3}} \otimes \mathbf{3}},O^0_{\mathbf{6}\otimes \overline{\mathbf{6}},\overline{\mathbf{3}} \otimes \mathbf{3}},O^2_{\overline{\mathbf{3}} \otimes \mathbf{3},\overline{\mathbf{3}} \otimes \mathbf{3}}\), respectively. The result for \(O^0_{\overline{\mathbf{3}} \otimes \mathbf{3},\mathbf{6}\otimes \overline{\mathbf{6}}}\) is the same as for \(O^0_{\mathbf{6}\otimes \overline{\mathbf{6}},\overline{\mathbf{3}} \otimes \mathbf{3}}\). We observe that the non-valishing contribution in quark-antiquark annihilation at LO is for spin-2 tetaquark with \(c^{(0)}_{q\bar{q},O^2_{\overline{\mathbf{3}} \otimes \mathbf{3},\overline{\mathbf{3}} \otimes \mathbf{3}}}=\frac{32768}{9}\).

At LO, there are 62 and 4 tree-level non-vanishing Feynman diagrams contributing to \(g+g\rightarrow T_{4c}\) and \(q+\bar{q}\rightarrow T_{4c}\). Since it is a 2-to-4-body process at the partonic level, this process can generate quite a few topological diagrams. At NLO, there are 2008 and 170 one-loop non-vanishing Feynman diagrams contributing to \(g+g\rightarrow T_{4c}\) and \(q+\bar{q}\rightarrow T_{4c}\). Also we need to consider the real corrections and find 618, 98 and 98 tree-level non-vanishing Feynman diagrams contributing to \(g+g\rightarrow T_{4c}+g\), \(q+\bar{q}\rightarrow T_{4c}+g\) and \(q+g\rightarrow T_{4c}+q\), respectively. Dimension regularization with \(D=4-2\epsilon\) is employed for both UV and IR divergences.

We have generated the tree-level and one-loop amplitude by using the package FeynArts [12]. After substituting the fermion chains, the package FeynCalc [13], [14] is used to simplify the Dirac matrices. Then the amplitudes can be expressed by the several form factor, \[\begin{align} \mathcal{M}_{ij}^{(1)} = \sum_{k=1}^{10} d_kF_k. \end{align}\] Then the corresponding coefficients \(\{d_k\}\) are the linear combinations of one-loop scalar integrals, which can be reduced to a set of basis integrals called master integrals (MIs) due to the identities from integration by parts (IBP) [15], [16] with the package Kira [17]. The analytical results of one-loop MIs can be easily obtained by package FeynHelpers [18] based on Package-X [19]. In addition, the ward identity is used to check our result.

The NLO virtual correction to the partonic cross-section can be expressed by \[\begin{align} 2{|\mathcal{M}_{ij}^{(0)}\mathcal{M}_{ij}^{(1)\dagger}|} &= |\mathcal{M}_{ij}^{(0)}|^2 \left(\frac{\alpha_s}{\pi}\right) \Gamma(1+\epsilon)(4\pi)^\epsilon K^\text{V}_{ij}. \end{align}\] As a concrete demonstration, the analytic expression for the \(2^{++}\) tetraquark from gluon-gluon fusion is given here

where the \(n_l (n_h)\) are the number of the light (heavy) quarks and \(\mu_R\) is the renormalization scale.

The NLO real correction including the gluon radiation to the partonic cross-section can be written as \[\begin{align} &\hat{\sigma}^{\mathrm{NLO,R}}\left(i+j\rightarrow T_{4c}+k\right) \nonumber\\&~~= \frac{1}{2\hat{s}} \frac{4^\epsilon}{\Gamma(1-\epsilon)}\left(\frac{4 \pi \mu_R^2}{\hat{s}}\right)^\epsilon \frac{1}{16 \pi}(1-z)^{1-2 \epsilon} \nonumber\\&~~~~\times \int_{-1}^1 d y_\theta\left(1-y_\theta^2\right)^{-\epsilon} \frac{1}{c_i c_j s_i s_j} \left|\mathcal{M}_{ijk}^{ (1)}\right|^2\left(z, y_\theta\right), \end{align}\] where \(y_\theta=\cos \theta\) with the angle \(\theta\) between the momenta of initial parton \(i\) and the tetraquark in lab frame. The explicit expressions for \(\left|\mathcal{M}_{ijk}^{ (1)}\right|^2\) and the virtual correction \(K^\text{V}_{ij}\) are listed in the Supplemental material.

Constraints on long-distance operator matrix elements for fully charm tetraquark. The aforementioned factorization formula is generally regarded as a model-independent effective theory derived from QCD, which factorizes the short-distance hard scattering processes from the non-perturbative interactions desribed by the long-distance matrix elements (LDMEs). The LDMEs are process-independent input parameters. There are four LDMEs for fully charm tetraquark in this work. We rewrite them as \[\begin{align} \langle\mathcal{O}_{6\bar{6}}^{T4c}\rangle\left(^1S_0\right)&=\langle0|O^0_{\mathbf{6}\otimes \overline{\mathbf{6}},\mathbf{6}\otimes \overline{\mathbf{6}}}|0\rangle, \nonumber\\ \langle\mathcal{O}_{3\bar{3}}^{T4c}\rangle\left(^1S_0\right)&=\langle0|O^0_{\overline{\mathbf{3}} \otimes \mathbf{3},\overline{\mathbf{3}} \otimes \mathbf{3}}|0\rangle,\nonumber\\ \langle\mathcal{O}_{mix}^{T4c}\rangle\left(^1S_0\right)&=\langle0|O^0_{\mathbf{6}\otimes \overline{\mathbf{6}},\overline{\mathbf{3}} \otimes \mathbf{3}}|0\rangle, \end{align}\] for \(0^{++}\) fully charm tetraquark, and \[\begin{align} \langle\mathcal{O}_{3\bar{3}}^{T4c}\rangle\left(^5S_2\right)&=\langle0|O^2_{\overline{\mathbf{3}} \otimes \mathbf{3},\overline{\mathbf{3}} \otimes \mathbf{3}}|0\rangle, \end{align}\] for \(2^{++}\) fully charm tetraquark.

In LHCb experiment, the production cross-section of the \(X(6900)\) structure relative to that of \(J/\psi\)-pair, times the branching fraction \({\cal B}(X(6900)\to J/\psi+J/\psi)\) is measured to be \({\cal R}=[1.1\pm 0.4(stat)\pm0.3(syst)]\%\) in the proton-proton collision data at \(\sqrt{S}=13TeV\) [4]. The prompt cross-section of \(J/\psi\)-pair at the same center-of-mass energy has also measured to be \(16.36\pm0.28(stat)\pm0.88(syst)nb\) in the transverse momentum range \(0<P_\perp<14GeV\) and the rapidity range \(2.0<y<4.5\) [20]. Considering the cross-section of \(J/\psi\)-pair in all \(P_\perp\) and \(y\) range, there is a factor around 4 to 4.8 from the theoretical calculation [21]. Thus the total cross-section of the \(X(6900)\) at \(\sqrt{S}=13TeV\) can be estimated taking the biggest uncertainty into account

\[\begin{align} \sigma(X(6900)){\cal B}(X(6900)\to J/\psi+J/\psi)&=0.72^{+0.79}_{-0.48}nb. \end{align}\] Assuming the \(X(6900)\) is a \(2^{++}\) tetraquark from the observation from CMS experiment [7] and its \(J/\psi\)-pair decay branching ratio around 1, we plot the NLO cross-section results for the fully charm tetraquark for various center-of-mass energy. In numerical calculation, the PDF4LHC21 combination of global PDF fits for the LHC Run III is employed [22], [23]. The mass of \(X(6900)\) is setted the newly value \(6847MeV\). The \(2^{++}\) LDMEs is then determined as \(\langle\mathcal{O}_{3\bar{3}}^{T4c}\rangle\left(^5S_2\right)=6.16\times 10^{-5}GeV^9\). The ratio among various LDMEs are also considered, i.e. \(\langle\mathcal{O}_{6\bar{6}}^{T4c}\rangle\left(^1S_0\right):\langle\mathcal{O}_{3\bar{3}}^{T4c}\rangle\left(^1S_0\right):\langle\mathcal{O}_{mix}^{T4c}\rangle\left(^1S_0\right):\langle\mathcal{O}_{3\bar{3}}^{T4c}\rangle\left(^5S_2\right)\), is equal to [1, 2.7, 1.6, 5.6] (Set-I) and [1, 1.3, -1.2, 4.5](Set-II) [11], [24]–[26]. By calculation, we find that the NLO K factor for the cross section at 13TeV is 2.31 and 1.09 for the \(0^{++}\) and \(2^{++}\) configuration, respectively. Both the NLO K factors increase when the center-of-mass energy decreases.

Discussions and conclusion

We also plot the rapidity and transverse momentum distribution in Figs 3 and 4 for different quark configurations. The differential cross section exhibits a relatively flat profile in rapidity, particularly across the low and intermediate rapidity ranges. The contribution from the \(2^{++}\) configuration is significantly enhanced—for instance, it is nearly two orders of magnitude higher than that of the \(6\bar{6}(0^{++})\) configuration. This enhancement stems not only from differences in LDMEs, but also from the difference in short-distance coefficients, a finding which is also consistent with previous theoretical and experimental works [7], [8], [27]–[30].

In NLO calculations, UV and IR divergences arise from loop diagrams and real corrections. The UV divergences can be rigorously canceled through the renormalization procedure, while the soft IR divergences cancel between the loop and real correction terms, leaving only the collinear divergences. The divergence-free pattern is also first observed in the quadripole picture of fully charm tetraquark, then we do not need introduce any NRQCD operator renormalization here. The final collinear divergences are obsorbed into the parton distribution functions.

The effects of soft and collinear gluon radiation induce singular structures, specifically terms proportional to \(1/P_\perp\) and \(\ln(P^2_\perp/M^2)/P_\perp\), in the transverse momentum distribution. This undermines the reliability of the theoretical prediction at low \(P_\perp\). By the standard CSS resummation procedure, these singularity at the Fourier transformation impact parameter b-space will be canceled between the loop and real correction terms while the large logrithems can be resummed at NLL accuracy [31]. We then also plot the resummation results for the transverse momentum distribution in Fig 4. We observe that the transverse momentum distribution is steeper than the rapidity distribution. However, it is flatter than the \(J/\psi\)-pair transverse momentum distribution, particularly at \(0<P_\perp < 20 GeV\). This relative flatness can explain the LHCb data with \(X(6900)\) signal enhancement at high transverse momenta.

We would like to conclude this work with the following remarks. First, our theoretical framework is built upon the fundamental principles of QCD and systematic symmetry analysis, ensuring the rigor of both the theoretical framework and the analytical results. Second, by utilizing LHCb and CMS experimental data, we have for the first time constrained and extracted the non-perturbative LDMEs for X(6900). These results can be widely applied to other reaction processes involving X(6900). Third, while the current discussion is limited to S-wave states without orbital excitation, even if contributions from P-wave [32] or higher partial waves need to be considered in the future, the S-wave results presented in this work should remain essential and critical for understanding these exotic hadrons. This is because we believe the basis expansion method based on symmetry is effective and covers these exotic states with corresponding coefficients. With the increasing data samples collected by the LHC, we also hope that experimental colleagues will perform more precise measurements of (differential) cross sections, rapidity and transverse momentum distributions for exotic hadrons in \(J/\psi\)-pair mass spectrum. Such efforts will provide further insights into clarifying the structural properties of these states, particularly regarding the existence of compact clustering mechanisms.

Acknowledgments. Thanks to Profs. Liu-Pan An and Zhi-Guo He for useful discussions. This work is supported by the National Natural Science Foundation of China Grants No.12322503, No.12405117 and No.12075124.

References↩︎