October 17, 2025
Using \((10.087\pm0.044)\times10^9\) \(J/\psi\) events collected with the BESIII detector at the \(e^+e^-\) BEPCII collider, we present the first amplitude analysis of \(J/\psi\to\gamma p\bar{p}\) with the \(p\bar p\) invariant mass in the \(\eta_c\) mass region \([2.70,3.05]\) GeV/\(c^2\). The product branching fraction \(\mathcal{B}(J/\psi\to\gamma\eta_c)\times\mathcal{B}(\eta_c\to p\bar{p})\) is determined to be \((2.11\pm0.02_{\rm stat}\pm0.07_{\rm syst})\times10^{-5}\) with precision improved by one order of magnitude. Combining with the product branching fractions \(\mathcal{B}(\eta_c\to p\bar{p})\times\mathcal{B}(\eta_c\to \gamma\gamma)\) and \(\mathcal{B}(J/\psi\to\gamma\eta_c)\times\mathcal{B}(\eta_c\to \gamma\gamma)\), the branching fractions of \(\mathcal{B}(J/\psi\to\gamma\eta_c)\) and \(\mathcal{B}(\eta_c\to\gamma\gamma)\) are calculated to be \((2.29\pm0.01_{\rm stat}\pm0.04_{\rm syst}\pm0.18_{\rm opbf})\%\) and \((2.28\pm0.01_{\rm stat}\pm0.04_{\rm syst}\pm0.18_{\rm opbf})\times10^{-4}\), respectively, which are consistent with the latest lattice quantum chromodynamics calculations. Here, opbf is the uncertainty from the other product branching fractions used in the calculation.
The transition between heavy quarkonium systems presents an ideal laboratory to investigate the theory of the strong interaction, quantum chromodynamics (QCD), in both the perturbative and non-perturbative regions. The magnetic dipole (M1) transition between the two lowest-lying charmonium states, \(J/\psi\to\gamma\eta_c\), is of great interest. The predicted transition width in the non-relativistic limit [1] is found to be significantly larger than the experimental results [2] by a factor of between 2 and 3. Several theoretical studies have attempted to resolve this long-standing puzzle, including dispersion sum rules [3], QCD sum rules [4], relativistic quark models [5], non-relativistic potential models [6], [7], effective field theories [1], [8], [9], light-cone sum rules [10], and lattice QCD (LQCD) [11]–[18]. However, a significant discrepancy remains between experimental measurements and theoretical predictions. In particular, LQCD calculations are systematically larger than the Particle Data Group (PDG) average [2] by approximately a factor of two.
Experimental measurements of \(\mathcal{B}(J/\psi\to\gamma\eta_c)\) were early reported by CLEO-c [19], KEDR [20], [21], and Crystal Ball [22] via inclusive hadronic decays of \(\eta_c\). Although these works successfully revealed the asymmetric lineshape of \(\eta_c\) and the necessity of damping factors, the potential interference between the \(\eta_c\) and non-resonant (NR) amplitudes was either ignored or insufficient considered. This neglect could potentially result in a bias of up to dozens of percent. Later, BESIII presented some measurements of the product branching fraction (BF) \(\mathcal{B}(J/\psi\to\gamma\eta_c)\times\mathcal{B}(\eta_c\to f)\) in various exclusive final states \(f\) [23]–[25], but still with the interference being ignored. However, as demonstrated in Ref. [26], even if taking into account this interference in one-dimensional fit to the \(\eta_c\) mass spectrum, a considerable uncertainty up to dozens of percent is inevitable due to the unknown NR contributions other than \(J^{PC}=0^{-+}\). Besides, two solutions with different interference patterns are found to be indistinguishable in Ref. [26], which lead to further ambiguities. Therefore, an amplitude analysis, such as that performed in Refs. [27], [28], is highly desired to incorporate all available information into the fit and provide a more reliable description of the interference, hence benefits precise determinations of the BF and resonance parameters.
The BESIII collaboration recently reported a BF measurement of \(\eta_c\to\gamma\gamma\) in \(J/\psi\to\gamma\eta_c\) [25]. The product \(\mathcal{B}(J/\psi\to\gamma\eta_c)\times\mathcal{B}(\eta_c\to \gamma\gamma)\) is in good agreement with the latest LQCD calculations [17], [29], while the BF of \(\eta_c\to \gamma\gamma\), determined using the BF of \(J/\psi\to\gamma\eta_c\) from the PDG [2], significantly differs from these LQCD calculations and the PDG global fit value [2] by more than 3\(\sigma\). Therefore, independent and precise measurements of the BF of \(J/\psi\to\gamma\eta_c\) are crucial for clarifying these discrepancies and testing the theoretical models.
Additionally, the hyperfine mass splitting between the \(J/\psi~(1^3 S_1)\) and \(\eta_c~(1^1 S_0)\) \(c\bar{c}\) systems is of importance for our understanding of interquark potential for charmonium system. Various Lattice-QCD calculations [30]–[33] have been reported, and most of which are in good agreement with each other within small theoretical uncertainties. Hence, precise determination of \(\eta_c\) resonance parameters in experiments is highly desired for the validation of relevant calculations.
In this Letter, by analyzing \((10.087\pm0.044)\times10^9\) \(J/\psi\) events [34] collected with the BESIII detector at the symmetric \(e^+e^-\) collider BEPCII, we present the first amplitude analysis of \(J/\psi\to\gamma p\bar{p}\) with the \(p\bar p\) invariant mass \(M_{p\bar{p}}\) in the \(\eta_c\) mass region \([2.70,3.05]\) GeV/\(c^2\), based on which precise measurements on \(\eta_c\) contributions and resonance parameters are performed.
Details about the design and performance of the BESIII detector are provided in Refs. [35]–[38]. The inclusive Monte Carlo (MC) sample, which includes both the production of the \(J/\psi\) resonance and the continuum processes incorporated in the kkmc [39] generator, are employed to study potential background contributions. All particle decays are modeled with the evtgen tool [40] using branching fractions either taken from the PDG [2], when available, or otherwise estimated with the lundcharm model [41]. Final state radiation (FSR) from charged final state particles is incorporated using photos [42]. The simulations of exclusive MC samples are described below.
Candidates for \(J/\psi\to\gamma p\bar{p}\) must have two charged tracks with zero net charge. Charged tracks detected in the main drift chamber (MDC) are required to be within a polar angle (\(\theta\)) range of \(|\rm{cos\theta}|<0.93\), where \(\theta\) is defined with respect to the \(z\)-axis, which is the symmetry axis of the MDC. Their distance of closest approach to the interaction point must be less than 10 cm along the \(z\)-axis, and less than 1 cm in the transverse plane. Particle identification (PID) is performed using the specific ionization energy loss and time of flight information, and the resultant likelihood for protons is required to be greater than those for pions and kaons. Photon candidates are chosen from isolated clusters in the electromagnetic calorimeter (EMC). Their energies are required to be greater than 25 MeV in the barrel (\(\vert\!\cos\theta\vert<0.8\)) region and 50 MeV in the end-cap (\(0.86<\vert\!\cos\theta\vert<0.92\)) region. Reconstructed clusters due to electronic noise or beam backgrounds are suppressed by requiring the EMC timing to be within [0, 700] ns after the event start time. To suppress background photons produced by hadronic interactions in the EMC, and secondary photons from bremsstrahlung radiation, clusters within cone angles of \(20^\circ\) and \(30^\circ\) around the extrapolated positions in the EMC of protons and anti-protons, respectively, are rejected [43]. At least one photon candidate is required for further analysis.
A four constraint (4C) kinematic fit, requiring energy and momentum conservation between the initial and final states, is imposed under the hypothesis of \(e^+e^-\to J/\psi\to\gamma p\bar{p}\). If there is more than one combination due to multiple photon candidates, the combination with the minimum \(\chi^2_{\rm 4C}\) is selected. The \(\chi^2_{\rm 4C}\) is required to be less than 23 based on the optimization of the Figure-of-Merit defined as \(\frac{S}{\sqrt{S+B}}\) [44], where \(S\) in the numerator is the signal yield from MC simulation and \(S+B\) in the denominator is the number of events from the data sample. Only candidates within the \(\eta_c\) mass region are kept for the amplitude analysis.
After applying all selection criteria mentioned above, 479,652 candidate events for \(J/\psi\to\gamma p\bar{p}\) survive in the data sample. The \(J/\psi\) inclusive MC simulation contains two dominant background components: \(J/\psi\to p\bar{p}\gamma^{\rm F}\) and \(J/\psi\to p\bar{p}\pi^0\), where \(\gamma^{\rm F}\) is a FSR photon. The process \(J/\psi\to p\bar{p}\gamma^{\rm F}\) is well simulated in the inclusive MC sample with photos [42], which contributes to the background at a level of 22.8% of the total events. The consistency between data and MC simulation is checked with the control sample \(J/\psi\to p\bar{p}\), where the energy spectra of FSR photons show good agreement. Another exclusive MC sample is simulated for the \(J/\psi\to p\bar{p}\pi^0\) decay based on a preliminary amplitude analysis result, and the corresponding background level is estimated to be 4.2%. The contribution from other \(J/\psi\) background processes is predicted to be less than 0.2% and modeled with inclusive MC simulation in the amplitude fit. The non-\(J/\psi\) background is studied with the data sample taken at \(\sqrt{s}=3.080\) GeV, corresponding to an integrated luminosity of \((167.4\pm0.1)\) pb\(^{-1}\) [34]. Its fraction, after taking into account the difference in luminosities [34], is estimated to be less than 0.7% and ignored in further analysis.
Figure 1 shows the distributions of the \(p\bar{p}\), \(p\gamma\), and \(\bar{p}\gamma\) invariant masses, as well as those of \(\cos(\theta_\gamma)\), \(\cos(\theta_p)\) and \(\phi_p\). Here, \(\theta_{\gamma}\) is the polar angle of the photon in the \(J/\psi\) rest frame with the \(z\)-axis defined as the direction of the \(e^+\) beam and \((\theta_{p},\phi_{p})\) are the polar and azimuthal angles, respectively, in the \(p\bar{p}\) helicity frame. No significant structures in the \(M_{p\bar{p}}\) or \(M_{p\gamma}\) two-body invariant mass spectra are seen, other than the \(\eta_c\) meson. As a check, potential contributions from \(N^*(1440/1520/1535/1650)\to\gamma p\) are estimated with N_J/[J/N^*|p, N^*p^0+c.c.], where \(\mathcal{B}[J/\psi\to N^*\bar{p},~N^*\to p\pi^0+c.c.]\) is taken from Ref. [45]; \(\mathcal{B}(N^*\to p\gamma/\pi^0)\) is from the PDG [2]; and \(\varepsilon\) is the detection efficiency of \(J/\psi\to N^*\bar{p}~(N^*\to p\gamma)+c.c.\) determined using MC simulation. The sum of background fractions of \(N^*\) baryons is less than 1.4%, and the \(N^{*}\) background is ignored in the amplitude analysis.
The covariant tensor amplitude constructed in Ref. [46] is applied in the amplitude analysis, which is expressed as
&A^(s)=_(p,m_J/)e^*_(q,m_)
&__s(p_p,S_p;p_|p,S_|p)_i_i U_i^_s.
Here, \(\psi_\mu(p,m_{J/\psi})\) is the polarization four-vector of the \(J/\psi\) with a spin projection \(m_{J/\psi}\) and four momentum \(p\); \(e_\nu(q,m_\gamma)\) is the polarization four-vector of the photon with spin projections \(m_{\gamma}\) and four momentum \(q\); \(\psi_{\lambda_s}(p_p,S_p;p_{\bar{p}},S_{\bar{p}})\) is the spin wave function of the proton and anti-proton system with polarizations \(S_{p,\bar{p}}\) and momenta \(p_{p,\bar{p}}\), where the index \(s\) is the total spin of the \(p\bar{p}\) system; \(U_i^{\mu\nu\lambda_s}\) is the \(i\)-th partial wave amplitude with a coupling strength determined by a complex parameter \(\Lambda_i\). The form of the \(\psi\) four-vector is detailed in Ref. [46].
Summing over the polarizations, the squared amplitude is given as
&||^2 = _S_p,S_|p=_m_J=_m_=|A^(s)|^2
&= -_i,j_i^*_j^2_U_i^_s g^()_^U_j^*^^_s_S_p,S_|p^*__s_^_s,
with \(-g^{(\perp\perp)}_{\mu\nu}=\sum_{m_\gamma}{ e^*_{\mu}(q,m_{\gamma})e_\nu(q,m_{\gamma})}\) [46]. For \(J/\psi\to\gamma\)“\(0^{-+}\)"\(\to\gamma p\bar{p}\), \(U_{i}^{\mu\nu\lambda_s}=\epsilon^{\mu\nu\rho\sigma}p_{\mu}q_{\sigma}B_1(Q_b)R\), where \(\epsilon^{\mu\nu\rho\sigma}\) is the Levi-Civita tensor, \(B_1(Q_b)\) is the Blatt-Weisskopf barrier factor [47] with angular momentum \(L=1\) and \(Q_{b}\) is the momentum of \(X\) in \(J/\psi\to\gamma X\) with \(X=\eta_c\) or non-resonant, \(R\) describes the line shape, and \(i\) and \(j\) are iterated over all possible processes. Explicit expressions of \(U^{\mu\nu\lambda_s}\) for other spin-parity cases are available in Ref. [46].
The line shape \(R\) of \(\eta_c\) is described by a relativistic Breit-Wigner function \(\frac{1}{M_{\eta_c}^2-M^2_{p\bar{p}}-i M_{\eta_c}\Gamma_{\eta_c}}\), where the mass \(M_{\eta_c}\) and width \(\Gamma_{\eta_c}\) of \(\eta_c\) vary freely in the fit. Additionally, the Blatt-Weisskopf barrier factor in \(J/\psi\to\gamma\eta_c\) is replaced by the square root of the damping factor \(f_d\) [19], [20], which is widely used to suppress the divergent long tail of \(\eta_c\). Two well-known damping factors \(e^{-E^2_\gamma/{8\beta^2}}\) of CLEO-c [19] with \(\beta\) floating and \(E^2_{\gamma 0}/{\left(E_{\gamma_0}E_\gamma+(E_\gamma-E_{\gamma_0})^2\right)}\) of KEDR [20], are considered in this work. Here, \(E_{\gamma}=\frac{M^2_{J/\psi}-M^2_{p\bar{p}}}{2 M_{J/\psi}}\) is the energy of the radiative photon and \(E_{\gamma_0}\) is the photon energy under the assumption \(M_{p\bar{p}}=M_{\eta_c}\). To account for detector resolution effects, the product \(R_{i}\times R^*_{j}\) in \(|\mathcal{M}|^2\) is convolved with a Gaussian function \(G(\delta_M,\sigma_M)\). The mass shift \(\delta_M=(1.01\pm0.07)\) MeV/\(c^2\) and resolution \(\sigma_M=(3.93\pm0.05)\) MeV/\(c^2\) are determined by studying the control sample \(\psi(3686)\to\gamma\chi_{c1},\chi_{c1}\to p\bar{p}\). For the non-resonant contributions, \(R\) is modeled by a constant.
The complex coupling constants \(\Lambda_i\) and the resonance parameters of \(\eta_c\) are determined with a maximum likelihood fit. The log-likelihood function is constructed as = _dt - _bg, where \(\ln \mathcal{L}_{\rm dt(bg)}\) sums over all the data or simulated background events \(N_{\rm dt(bg)}\) and is defined as _dt(bg) = ^N_dt(bg)_k=1, where \(p\) is the momentum of the final state particles and \(\Phi_{3}\) is the phase space factor. Integral of \(\int\epsilon(p)|\mathcal{M}(p)|^2 \Phi_{3}(p)dp\) is calculated numerically using MC events as (p)|(p)|^2 _3(p)dp ^N_MC_k_MC=1|(p^k_MC)|^2. Here, \(N_{\rm MC}=2.3\times10^{6}\) is the number of phase space events that survive the data selection. The fit fractions with detection efficiency are estimated with f_i = / , where \(\mathcal{M}_{i}\) is the amplitude of contribution \(i\) alone.
First the CLEO-c and KEDR damping factors are tested with \(\eta_c\) and all the potential NR contributions included. Since the KEDR damping factor results in a better log-likelihood value with a statistical significance \(\sqrt{2\times\Delta\ln\mathcal{L}}=4.7\sigma\), as well as improved fit quality \((\Delta\chi^2/{\rm nbin}=10.4/100)\) in the \(M_{p\bar{p}}\) projection compared to the CLEO-c form, the KEDR form is used in this analysis. An additional test is performed by excluding the damping factor, which yields a much worse fit quality with \(\Delta\chi^2/{\rm nbin}=49.1/100\) in the \(M_{p\bar{p}}\) projection. All 13 potential NR contributions, each of which corresponds to different angular momentum and spin coupling in the transitions, are excluded from the solution one at a time, and the corresponding statistical significances are calculated based on the change of log-likelihood function \(\Delta\ln\mathcal{L}\) and the number of free parameters \(\Delta N_{\rm par}=2\). Five components with statistical significance greater than 3\(\sigma\), including three waves with \(J^{P}=\) \(0^{-}\), \(1^+\), and \(2^+\), as well as two waves with \(J^{P}=\) \(2^-\), are kept in the final solution. By performing a scan on the phase angle \(\phi\) between \(\eta_c\) and \(0^{-+}\) NR, two local minima are found. The best one exhibits a superior log-likelihood value with a statistical significance of \(7.7\sigma\), hence we only consider it in the further analysis.
Figure 1 shows the projections of the nominal amplitude analysis result. The mass and width of \(\eta_c\) are determined to be \(M_{\eta_c}=(2984.55\pm0.09_{\rm stat})\) MeV/\(c^2\) and \(\Gamma_{\eta_c}=(29.74\pm0.17_{\rm stat})\) MeV. The fit fraction of \(\eta_c\) is determined to be \(f_{\eta_c}=(30.94\pm0.24)\%\) and the product BF, \(\mathcal{B}(J/\psi\to\gamma\eta_c)\times\mathcal{B}(\eta_c\to p\bar{p})=\frac{(N_{\rm dt}-N_{\rm bg})\cdot f_{\eta_c}}{\varepsilon_{\eta_c}\cdot N_{J/\psi}}\), is calculated to be \((2.11\pm0.02_{\rm stat})\times10^{-5}\). Here, \(N_{\rm dt}-N_{\rm bg}\) is the net number of signal events, \(N_{J/\psi}=(10.087\pm0.044)\times10^{9}\) is the number of \(J/\psi\) events [34], and \(\varepsilon_{\eta_c}=50.55\%\) is the signal efficiency determined with MC simulation based on the amplitude analysis result.
As a comparison, we have also tried to extract the \(\eta_c\) signal yield via one-dimensional fit to the \(M_{p\bar{p}}\) spectrum following previous publications. Without the interference between \(\eta_c\) and NR contributions [19]–[22], the fit model can not provide an acceptable description of the data sample, and the fitted \(\eta_c\) signal yield deviates from our nominal result by more than 30%. With the interference considered following Ref. [26], two indistinguishable solutions with distinct interference patterns and \(\eta_c\) signal yields are observed as expected. Additionally, the fraction of NR components other than \(0^{-+}\) are found to be about 15% in the amplitude analysis, which can not be determined in the one-dimensional fit and thereby causes large uncertainties [26]. In contrast, the amplitude analysis successfully overcome these two issues.
| Source | \(\mathcal{B}\) (%) | \(M_{\eta_c}\) (MeV/\(c^2\)) | \(\Gamma_{\eta_c}\) (MeV) |
|---|---|---|---|
| \(N_{J/\psi}\) | 0.5 | — | — |
| Tracking | 0.2 | — | — |
| PID | 0.3 | — | — |
| Photon | 1.0 | — | — |
| 4C kinematic fit | 0.5 | — | — |
| NR line shape | 0.4 | 0.38 | 0.14 |
| Insig. NR waves | 1.8 | \(\star\) | 0.13 |
| Background | 0.8 | \(\star\) | 0.06 |
| Mass calibration | \(\star\) | 0.37 | 0.24 |
| Fit bias | 2.0 | 0.55 | 0.10 |
| Total | 3.1 | 0.77 | 0.33 |
The systematic uncertainties on the product \(\mathcal{B}(J/\psi\to\gamma\eta_c)\times\mathcal{B}(\eta_c\to p\bar{p})\) and the resonance parameters of \(\eta_c\) are summarized in Table 1. The uncertainty on the total number of \(J/\psi\) events is 0.5% [34]. The systematic uncertainties due to the tracking and PID of \(p(\bar{p})\) are studied using the control sample \(J/\psi\to p\bar{p}\pi^+\pi^-\), and are determined to be 0.2% and 0.3%, respectively. The systematic uncertainty for photon reconstruction is assigned to be 1.0% [48]. The systematic uncertainty associated with the 4C kinematic fit is estimated by performing corrections on the charged track helix parameters in the MC simulation. The difference between the detection efficiencies obtained with and without the helix parameter correction [49] is taken as the systematic uncertainty.
The systematic uncertainty due to the NR line shape is estimated by modeling the \({\rm NR}(0^{-+})\) component with an alternative line shape, the magnitude and phase of which are allowed to vary linearly as a function of \(M_{p\bar{p}}\). The systematic uncertainty from NR waves considered insignificant is estimated by including all 13 potential waves. The systematic uncertainty in the background estimation is studied by varying the contribution from \(J/\psi\to p\bar{p}\pi^0\) within the uncertainty of the quoted BF, and by including the simulated \(N^*\to p\gamma\) contribution as background. The systematic uncertainty due to the mass calibration parameters, mass shift \(\delta_M\) and resolution \(\sigma_M\), is estimated by varying them within their statistical uncertainties and by using the control sample \(\psi(3686)\to\gamma\chi_{c2},\chi_{c2}\to p\bar{p}\) as an alternative. For all four sources listed above in this paragraph, the amplitude fit is reperformed, and the largest difference for each source is assigned as its systematic uncertainty. The systematic uncertainty caused by a possible fit bias is studied by performing input and output checks with toy MC samples. The difference between the input and the averaged output values, predominantly attributed to detector effects and the statistical fluctuations of the phase space MC sample used in amplitude fit, is conservatively assigned as this uncertainty.
In summary, the first amplitude analysis of \(J/\psi\to\gamma p\bar{p}\) with \(M_{p\bar{p}}\) in the \(\eta_c\) mass region \([2.70,3.05]\) GeV/\(c^2\) is performed. Compared to the one-dimensional fits to the \(\eta_c\) spectrum [19]–[22], [26], our amplitude analysis approach provides more reliable estimations of the different non-resonant components, and their interference with \(\eta_c\), hence, avoids several 10%–level biases and uncertainties. The mass and width of \(\eta_c\) are measured to be \((2984.55\pm0.09_{\rm stat}\pm0.77_{\rm syst})\) MeV/\(c^2\) and \((29.74\pm0.17_{\rm stat}\pm0.33_{\rm syst})\) MeV, respectively, which are in good agreement with the PDG world averages [2]. The product BF \(\mathcal{B}(J/\psi\to\gamma\eta_c)\times\mathcal{B}(\eta_c\to p\bar{p})\) is determined to be \((2.11\pm0.02_{\rm stat}\pm0.07_{\rm syst})\times10^{-5}\), whose precision is improved by one order of magnitude compared to the previous measurements [2].
Given most measurements used in the PDG global fits ignored interference effects [2], the PDG fitted \(\mathcal{B}(\eta_c\to p\bar{p})\) is not quoted for the extraction of \(\mathcal{B}(J/\psi\to\gamma\eta_c)\). Instead, two unique processes \(J/\psi\to\gamma\eta_c,\eta_c\to\gamma\gamma\) and \(\gamma\gamma\leftrightarrow p\bar{p}\) around \(\eta_c\) peak, which benefit from limited interference effects [25], [50], are quoted. Combining this result with the products \(\mathcal{B}(\eta_c\to p\bar{p})\times\mathcal{B}(\eta_c\to \gamma\gamma)=(2.1\pm0.3)\times10^{-7}\) averaged based on Refs. [50]–[52] and the recently reported \(\mathcal{B}(J/\psi\to\gamma\eta_c)\times\mathcal{B}(\eta_c\to \gamma\gamma)=(5.23\pm0.40)\times10^{-6}\) [25], we obtain
(J/_c)&==(2.29)%,
(_c)&==(2.28)^-4,
(_cp|p)&==(0.92)^-3,
where the uncertainties are statistical, systematic, and those from the other product BFs used in the calculation. Some systematic uncertainties are correlated between our work and Ref. [25], which are dominated by the photon detection and the form of the damping factor. Figures 2 and 3 show comparisons of \(\mathcal{B}(J/\psi\to\gamma\eta_c)\) and \(\mathcal{B}(\eta_{c}\to\gamma\gamma)\) determined in this study, with various theoretical calculations and other measurements. Our results deviate from the PDG global fit values [2] by \(3\sigma\), but are in good agreement with the latest LQCD calculations [17], [18], [29]. This helps resolve a long-standing puzzle. Furthermore, the obtained BF \(\mathcal{B}(\eta_c\to p\bar{p})\) deviates from the PDG global fit values by 3\(\sigma\). Additionally, using \(M_{J/\psi}=3096.9\) MeV/\(c^2\) [2], the mass splitting between the \(J/\psi\) and \(\eta_c\) \(c\bar{c}\) systems is determined to be \((112.35\pm0.77)\) MeV/\(c^2\), which is consistent with the results calculated in Refs. [30]–[33].
The BESIII Collaboration thanks the staff of BEPCII (https://cstr.cn/31109.02.BEPC) and the IHEP computing center for their strong support. This work is supported in part by National Key R&D Program of China under Contracts Nos. 2025YFA1613900, 2023YFA1606000, 2023YFA1606704; National Natural Science Foundation of China (NSFC) under Contracts Nos. 11635010, 11935015, 11935016, 11935018, 12025502, 12035009, 12035013, 12061131003, 12192260, 12192261, 12192262, 12192263, 12192264, 12192265, 12221005, 12225509, 12235017, 12361141819; the Chinese Academy of Sciences (CAS) Large-Scale Scientific Facility Program; CAS under Contract No. YSBR-101; 100 Talents Program of CAS; The Institute of Nuclear and Particle Physics (INPAC) and Shanghai Key Laboratory for Particle Physics and Cosmology; German Research Foundation DFG under Contract No. FOR5327; Istituto Nazionale di Fisica Nucleare, Italy; Knut and Alice Wallenberg Foundation under Contracts Nos. 2021.0174, 2021.0299; Ministry of Development of Turkey under Contract No. DPT2006K-120470; National Research Foundation of Korea under Contract No. NRF-2022R1A2C1092335; National Science and Technology fund of Mongolia; Polish National Science Centre under Contract No. 2024/53/B/ST2/00975; Swedish Research Council under Contract No. 2019.04595; U. S. Department of Energy under Contract No. DE-FG02-05ER41374