Measurement of the Born cross section for $e^{+}e^{-}\to\eta h_{c}$ at center-of-mass energies between 4.1 and 4.6 GeV

In recent years, several new vector resonances, including the $Y(4230)$, $Y(4360)$ and $Y(4660)$, have been observed in the charmonium region at B-factories and $\tau$-charm factories [1, 2, 3]. While the potential models [4] can accommodate the $\psi(4040)$, $\psi(4160)$, and $\psi(4415)$, the new states appear to be supernumerary. Several models have been proposed to explain them as exotic non-$c\bar{c}$ mesons [5, 6, 7, 8, 9]. The nature of the first observed vector charmonium-like state, the $Y(4230)$ (also known as $\psi(4230)$) is still mysterious. It is regarded as a good candidate for a hybrid state because its mass is close to the vector hybrid state predicted by lattice QCD [10] and also because of its small electronic width [11, 12] and decay pattern [13]. But it is also interpreted as a conventional charmonium [14, 15, 16], hadronic molecule [17, 18], non-resonant enhancement [19], etc. More experimental measurements are valuable to elucidate the nature of the $Y(4230)$.

Recently, BESIII performed a precise measurement of the cross section of $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ at center-of-mass (c.m.) energies from 3.774 to 4.600 GeV, in which two resonances were observed, namely the $Y(4230)$ and the $Y(4320)$ [20]. This observation has challenged our initial understanding of the $Y(4230)$. A similar double-resonance structure is also observed in other hadronic modes, such as the $e^{+}e^{-}\to\pi^{+}\pi^{-}h_{c}$ [21] and $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ [22] processes. Studying the transitions from vector states to the $h_{c}$ meson is particularly interesting because strong coupling to $h_{c}$ is indicative of a hybrid state with a $c\bar{c}$ pair in a spin-singlet configuration [23].

The process $e^{+}e^{-}\to\eta h_{c}$ was previously observed for the first time at $\sqrt{s}=$ 4.226 GeV by BESIII using data collected in 2012-2014 [24]. A hint of a resonance around 4.200 GeV was observed in the c.m. energy-dependent cross section. Recently, the BESIII experiment collected more data at c.m. energies from 4.129 to 4.600 GeV. In this Letter we use the new data sets to update the study of $e^{+}e^{-}\to\eta h_{c}$ with $h_{c}\to\gamma\eta_{c}$ and $\eta\to\gamma\gamma$ to substantially improve our measurement of the cross section. The total integrated luminosity is measured to be $\rm\sim 15\,fb^{-1}$ using large-angle Bhabha scattering events with an uncertainty of $1.0\%$ [25], and the c.m. energies are measured using the di-muon process [26].

The BESIII detector is described in detail elsewhere [27]. The determination of the detection efficiency and estimation of physics backgrounds are carried out with Monte Carlo (MC) samples. Geant4-based [28, 29] detector simulation software is used to model the detector response. The signal MC events are simulated for each decay mode of the $\eta_{c}$ meson with kkmc [30] and besevtgen [31], in which the line shape of $\eta_{c}$ is a Breit-Wigner (BW) function. The $e^{+}e^{-}\to\eta h_{c}$ process is assumed to be dominated by the S wave, and the E1 transition $h_{c}\to\gamma\eta_{c}$ is simulated using the dedicated helicity formalism [32]. In order to study potential backgrounds, inclusive MC samples are simulated at each c.m. energy with kkmc. These MC samples consist of charmed meson production, initial state radiation (ISR) production of the low-mass vector charmonium states, QED events and continuum processes. The known decay modes of the resonances are simulated with besevtgen with branching fractions set to the world average values [33], and the remaining events associated with charmonium decays are simulated with lundcharm [34]. Other hadronic events are simulated with pythia [35].

In this measurement, we first reconstruct the E1 photon and bachelor $\eta$, and then the $\eta_{c}$ is reconstructed with sixteen hadronic final states with a total branching fraction of about 40%: $p\bar{p}$, $2(\pi^{+}\pi^{-})$, $2(K^{+}K^{-})$, $K^{+}K^{-}\pi^{+}\pi^{-}$, $p\bar{p}\pi^{+}\pi^{-}$, $3(\pi^{+}\pi^{-})$, $K^{+}K^{-}2(\pi^{+}\pi^{-})$, $K^{+}K^{-}\pi^{0}$, $p\bar{p}\pi^{0}$, $K^{0}_{S}K^{\pm}\pi^{\mp}$, $K^{0}_{S}K^{\pm}\pi^{\mp}\pi^{\pm}\pi^{\mp}$, $\pi^{+}\pi^{-}\eta$, $K^{+}K^{-}\eta$, $2(\pi^{+}\pi^{-})\eta$, $\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, and $2(\pi^{+}\pi^{-})\pi^{0}\pi^{0}$. The $\eta/\pi^{0}$ candidates are reconstructed using two photons and the $K_{S}^{0}$ candidates are reconstructed via the $\pi^{+}\pi^{-}$ decay channel. For the selected candidates, we apply a fit to the distribution of the $\eta$ recoil mass to obtain the $\eta h_{c}$ signal yield.

The selection criteria for charged tracks and photon candidates as well as the $\pi^{0}$, $\eta$ and $K_{S}^{0}$ reconstruction are described in Ref. [24]. After this selection, a four-constraint (4C) kinematic fit is performed for each event imposing overall energy-momentum conservation, and the $\chi^{2}_{\rm 4C}$ is required to be less than 25 to suppress background events. The best candidates of $\pi^{0}$, $\eta$ and $K_{S}^{0}$ as well as the particle identification (PID) assignments of charged tracks in an event are determined by minimizing $\chi^{2}\equiv\chi^{2}_{\rm 4C}+\chi^{2}_{\rm 1C}+\chi^{2}_{\rm PID}+\chi^{2}_{\rm vertex}$, where $\chi^{2}_{\rm 1C}$ is the overall $\chi^{2}$ of the 1C fit for all $\pi^{0}$ and $\eta$ candidates, $\chi^{2}_{\rm PID}$ is the sum over all charged tracks of the $\chi^{2}$ of the PID hypotheses, and $\chi^{2}_{\rm vertex}$ is the $\chi^{2}$ of the $K_{S}^{0}$ secondary-vertex fit. If there is no $\pi^{0}$ ($K_{S}^{0}$) in an event, the corresponding $\chi^{2}_{\rm 1C}$ ($\chi^{2}_{\rm vertex}$) is set to zero. If more than one $\eta$ candidate with recoil mass in the $h_{c}$ pre-selection signal region [3.480, 3.600] GeV$/c^{2}$ is found, the one with mass of the $\eta_{c}$ candidate closest to $\eta_{c}$ known mass is selected.

The requirements on $\chi_{\rm 4C}^{2}$ and the mass windows for $\eta$ and $\eta_{c}$ selection are determined by maximizing the figure-of-merit, which is defined as $S/\sqrt{S+B}$. Here, $S$ and $B$ are the signal and background yields, respectively. The optimization was performed using the combined statistics of four samples with high integrated luminosity taken at $\sqrt{s}$ = 4.179, 4.189, 4.199 and 4.209 GeV.

After applying all the above criteria and using a $h_{c}$ mass window of [3.510, 3.540] ${\rm GeV}/c^{2}$, a clear $\eta_{c}$ signal is observed in the invariant mass spectrum of hadrons in the data sample taken at $\sqrt{s}=$ 4.179 GeV, which is shown in Fig. 1. A clear $h_{c}$ signal is also seen in the $\eta$ recoil mass distribution when applying a mass window of [2.944, 3.024] ${\rm GeV}/c^{2}$ to the invariant mass distribution of hadrons from $\eta_{c}$ decay, as shown in Fig. 2.

In the $\eta$ recoil mass spectra, as depicted in Fig. 2, no peaking structures are observed within the $\eta_{c}$ sidebands, which are delineated as [2.870, 2.910] ${\rm GeV}/c^{2}$ and [3.050, 3.090] ${\rm GeV}/c^{2}$, as demonstrated in Fig. 1. Likewise, the hadronic invariant mass spectrum resulting from $\eta_{c}$ decay within the $h_{c}$ sidebands, specified as [3.490, 3.505] ${\rm GeV}/c^{2}$ and [3.545, 3.560] ${\rm GeV}/c^{2}$, exhibits no peak formations, as illustrated in Fig. 1. In addition, inclusive MC samples simulated at $\sqrt{s}=4.179\,\rm{GeV}$ are analyzed with the same event selection as applied to data, and the dominant background is found to be the continuum processes. Other sources just contribute negligible background. The comparison between data and the inclusive MC sample is shown in Fig. 2.

To obtain the $h_{c}$ yield, the sixteen $\eta$ recoil mass distributions are fitted simultaneously with an unbinned maximum likelihood method. In the fit, the signal shape is determined by MC simulation, and the background is described by an ARGUS function [36]. The truncation point is set to be the same for all channels and fixed according to simulation at each c.m. energy. The total signal yield, $N_{sig}$, is the sum of the yields from all the sixteen $\eta_{c}$ decay channels. The signal yield of the $i$-th channel is set to be $N_{\rm sig}\cdot f_{i}$, where $f_{i}$ is the weight factor $f_{i}=\epsilon_{i}\mathcal{B}_{i}/\sum\epsilon_{i}\mathcal{B}_{i}$, $\mathcal{B}_{i}$ denotes the branching fraction of $\eta_{c}$ decays to the $i$-th final state and $\epsilon_{i}$ is the corresponding efficiency. The fraction of the background in each mode is free. The sum of the fit results at $\sqrt{s}=4.179$ GeV is shown in Fig. 2. The corresponding $\chi^{2}$ per degree of freedom (dof) for this fit is $\chi^{2}/\rm{dof}=43.8/46=0.95$. The total signal yield is $104\pm 16$ with a statistical significance of $9.6\sigma$.

Using the same method, we also study the data samples taken at other c.m. energies. Significant signals (more than $5\sigma$) are observed at $\sqrt{s}=$ 4.179, 4.189 and 4.226 GeV, and evidence (between $3\sigma$ and $5\sigma$) is found at $\sqrt{s}=$ 4.209, 4.358 and 4.436 GeV.

The Born and “Dressed” cross sections are calculated by the following formula:

where $\mathcal{L}$ is the integrated luminosity of the data sample taken at each c.m. energy. The radiative correction factor $(1+\delta)$ at each c.m. energy is calculated iteratively [24]. The term $|1+\Pi|^{2}$ is the vacuum-polarization (VP) correction factor and is calculated according to Ref. [37]. In addition, $\mathcal{B}(\eta\to\gamma\gamma)$ and $\mathcal{B}(h_{c}\to\gamma\eta_{c})$ are the branching fractions of the $\eta$ and $h_{c}$ decays, respectively, and $\epsilon_{i}$ is the efficiency determined with MC simulation. The measured Born cross sections and the quantities that enter the calculation at all c.m energies are listed in the Supplemental Material [38].

The cross section of $e^{+}e^{-}\to\eta h_{c}$ from $\sqrt{s}=$ 4.129 to 4.600 GeV is parameterized as:

The relativistic BW amplitude for a resonance $Y\to\eta h_{c}$ in the fit is written as:

where $M$, $\Gamma_{\rm tot}$, $\Gamma_{ee}$, and ${\cal B}(Y\to\eta h_{c})$ are the mass, full width, electronic partial width , and branching fraction of the corresponding resonance, respectively, and $\sqrt{{\rm PS}(\sqrt{s})/{\rm PS}(M)}$ is the two-body phase space factor [33]. It is important to note that the definition of $\Gamma_{ee}$ includes vacuum polarization effects. Consequently, the dressed cross sections are fitted rather than the Born cross sections. A maximum likelihood fit is used to obtain the parameters of these three resonances. In the fit, the parameters of the second BW function are fixed to that of the $Y(4360)$ due to the large uncertainty of the cross section in this region, while the other two BW functions are free. Our analysis is primarily concerned with the measurement of ${\rm BW}_{1}$. Consequently, the interference effects between ${\rm BW}_{1}$ and ${\rm BW}_{3}$, as well as between ${\rm BW}_{2}$ and ${\rm BW}_{3}$, have been neglected in the nominal fit because of statistics for ${\rm BW}_{2}$ and ${\rm BW}_{3}$. There could be multiple solutions in the fit due to interference, but we vary the initial values of $\phi$ and only find one solution.

The fit results are shown in Fig. 3(a) and the fitted parameters of the BW functions are listed in Table 1, where the uncertainties are statistical only. The statistical significance of the first resonance is calculated to be $8\sigma$, which is obtained by comparing the change of the log-likelihood value $\Delta(-{\rm ln}L)=41$ and degrees of freedom $\Delta{\rm dof}=4$ with and without this resonance in the fit.

To check the fit stability of the first resonance, we also use many different parameterizations, including changing the second resonant parameters to different hypotheses: $Y(4320)$ [20], $Y(4380)$ [22], $Y(4390)$ [21] and removing the second resonance in the model. In addition, we also use the sum of a BW function and phase space shape to fit the cross section line shape and use a model that takes the interference between all three resonances into account. The comparison of the fitted line shapes of c.m. energy-dependent cross sections is shown in Fig. 3(b). The choice of the fit model leads to the dominant systematic uncertainties on the mass and width parameters. The significance of the first resonance remains above 7$\sigma$ in all the alternative fits.

The systematic uncertainties for the measured Born cross sections are determined as follows. The integrated luminosity is measured using Bhabha events, with an uncertainty of $1.0\%$ [25]. To estimate the uncertainty due to the data/MC mass resolution difference in the fit to the recoil mass of $\eta$, the simulated signal shape is shifted and convolved with a Gaussian function to match the shape in data. The parameters of the Gaussian function are obtained by a control sample of $e^{+}e^{-}\to\eta J/\psi$. The average change of the cross section with and without this correction among all the data sets is taken as the common systematic uncertainty. To estimate the uncertainty due to the background shape, a second order Chebyshev function instead of the ARGUS function is used as an alternative model. The average fit result difference between these two background shapes in all data sets is adopted as the systematic uncertainty. The systematic uncertainty from the fit range is determined by changing the fit range randomly and redoing the fit, and the average difference from the nominal result among all data sets is taken as the systematic uncertainty. The branching fraction of $h_{c}\to\gamma\eta_{c}$ is taken from Ref. [39], and its uncertainty, which is 15.7%, will propagate to the cross section measurement. The ISR correction factor at a given c.m. energy is determined using the cross section line shape from the threshold to the c.m. energy of interest. The nominal energy-dependent cross section is parameterized with the sum of three BW functions as shown in Fig. 3(a). The uncertainty of the line shape is estimated by using different models as discussed in the fit to the c.m. energy-dependent cross section (shown in Fig. 3(b)). The line shape of the cross section also affects the efficiency. To consider this effect, we use the method introduced in Ref. [24]. To investigate the uncertainty due to the vacuum polarization factor, we use two available VP parameterization [37, 40]. The difference between them is 0.3% and is taken as the systematic uncertainty. In the simultaneous fit to the $\eta$ recoil mass distributions, $\epsilon_{i}\mathcal{B}_{i}$ is used to constrain the weights between different $\eta_{c}$ decay modes, so the uncertainty from $\epsilon_{i}\mathcal{B}_{i}$ will affect the signal yield. To consider this issue, we use the refitting method introduced in Ref. [24] to estimate the uncertainty due to $\epsilon_{i}\mathcal{B}_{i}$. The systematic uncertainties related to efficiency including charged track, photon, $K^{0}_{S}$, $\pi^{0}$ and $\eta$ reconstruction, PID, kinematic fit, $\eta_{c}$ tag and cross feed are estimated with the same method as described in Ref. [24]. The angular momentum between the $\eta$ and $h_{c}$ mesons is investigated from data and is found to be between S and D waves. The uncertainty is estimated by comparing the efficiencies from the pure S wave MC and mixture of S and D wave MC. The uncertainties related to efficiency at $\sqrt{s}=4.179$ GeV are given in the Supplemental Material [38].

The systematic uncertainties from different sources are listed in Table 2. All sources are treated as uncorrelated, so the total systematic uncertainty is obtained by summing them in quadrature. For the data sets without significant $\eta h_{c}$ signals, an upper limit of the cross section at the 90% confidence level is obtained using a Bayesian method. The systematic uncertainties are taken into account by convolving the probability density function of the measured cross section with a Gaussian function [41].

The systematic uncertainties for the fit parameters to the cross section line shape are described as follows. To estimate the uncertainty for the fit model, we use all the variations described in the fit to the c.m. energy-dependent cross section as shown in Fig. 3(b). The largest deviation is taken as the systematic uncertainty. The c.m. energies of all data sets are measured using di-muon events with uncertainty $\pm 1$ MeV. The uncertainty of the c.m. energy measurement will propagate to the mass of the resonances directly. The uncertainty due to the beam energy spread is found to be negligible. The systematic uncertainty for the cross section measurements can be classified to two categories: the uncertainties due to the fit, which are data-set independent and the remaining uncertainties, which are treated as correlated among all data-sets. The first kind are studied by changing the configuration in the fit to the $\eta$ recoil mass. We vary the signal shape, background shape and fit range for each data set as discussed before and then redo the fit to these alternatively measured cross sections, and the largest deviations of the fitted parameters of ${\rm BW}_{1}(s)$ are taken as the systematic uncertainties. The second kind of uncertainties would not affect mass and width of each resonance, but will propagate to the $\Gamma_{ee}^{Y}\mathcal{B}$ by the same amount directly. The average of the correlated systematic uncertainty for $\Gamma_{ee}^{Y}\mathcal{B}$ is 19%. Table 3 summarizes the uncertainties of parameters for the resonance around $\sqrt{s}=$ 4.200 GeV from the c.m. energy-dependent cross sections.

In summary, the $e^{+}e^{-}\to\eta h_{c}$ process is studied with the data samples taken at c.m. energies from 4.129 to 4.600 GeV. The corresponding Born cross sections or upper limits of Born cross sections at each c.m. energy are obtained. In the cross section lineshape, a resonant structure near 4.200 GeV is observed with a statistical significance of 7$\sigma$. The parameters of this resonance are measured to be: $M=4188.8\pm 4.7\pm 8.0\,{\rm MeV/}c^{2}$ and $\Gamma=49\pm 16\pm 19\,{\rm MeV}$. Here, the first uncertainties are statistical and second systematic. This result is consistent with the parameters of $\psi(4160)$ but also not far away from the $\psi(4230)$ observed from the $\pi^{+}\pi^{-}J/\psi$ process. The $1^{--}$ hybrid charmonium state predicted by the BOEFT model [42] has a mass of $4.15\pm 0.15$ GeV, therefore, our measurement is also consistent with this prediction.