Observation of the $Y(4220)$ and $Y(4360)$ in the process $e^{+}e^{-} \to \eta J/\psi$

The cross sections of the process $e^{+}e^{-} \to \eta J/\psi$ at center-of-mass energies ($\sqrt{s}$) between 3.81 and 4.60 GeV are measured with high precision by using data samples collected with the BESIII detector operating at the BEPCII storage ring. Three structures are observed by analyzing the lineshape of the measured cross sections, and a maximum-likelihood fit including three resonances is performed by assuming the lowest lying structure is the $\psi(4040)$. For the other resonances, we obtain masses of $(4218.7 \pm 4.0 \pm 2.5)$ and $(4380.4 \pm 14.2 \pm 1.8)$ MeV/c$^{2}$ with corresponding widths of $(82.5 \pm 5.9 \pm 0.5)$ and $(147.0 \pm 63.0 \pm 25.8)$ MeV, respectively, where the first uncertainties are statistical and the second ones systematic. The measured resonant parameters are consistent with those of the $Y(4220)$ and $Y(4360)$ from pr evious measurements of different final states. For the first time, we observe the decays of the $Y(4220)$ and $Y(4360)$ into $\eta J/\psi$ final states.

The cross sections of the process e + e − → ηJ/ψ at center-of-mass energies ( √ s) between 3.81 and 4.60 GeV are measured with high precision by using data samples collected with the BESIII detector operating at the BEPCII storage ring. Three structures are observed by analyzing the lineshape of the measured cross sections, and a maximum-likelihood fit including three resonances is performed by assuming the lowest lying structure is the ψ(4040). For the other resonances, we obtain masses of (4218.7 ± 4.0 ± 2.5) and (4380.4 ± 14.2 ± 1.8) MeV/c 2 with corresponding widths of (82.5 ± 5.9 ± 0.5) and (147.0 ± 63.0 ± 25.8) MeV, respectively, where the first uncertainties are statistical and the second ones systematic. The measured resonant parameters are consistent with those of the Y (4220) and Y (4360) from previous measurements of different final states. For the first time, we observe the decays of the Y (4220) and Y (4360) into ηJ/ψ final states.
The process e + e − → ηJ/ψ was studied using ISR by Belle [19]. Unlike the process e + e − → π + π − J/ψ [12], two resonant structures at 4040 and 4160 MeV/c 2 , regarded as the ψ(4040) and ψ(4160), respectively, were observed by studying the cross section dependence on the center-of-mass (c.m.) energy. Furthermore, a hint of an enhancement at 4360 MeV/c 2 was reported. Using data samples at 17 c.m. energies from 3.81 to 4.60 GeV, BESIII reported more accurate measurements of cross sections of the e + e − → ηJ/ψ process [20]. The BESIII data agree well with that of Belle. However, due to limited statistics it was not possible to establish any potential involved Y states.
In this Letter, we present an updated analysis of e + e − → ηJ/ψ at c.m. energies between 3.81 and 4.60 GeV, where J/ψ is reconstructed with ℓ + ℓ − (ℓ = e/µ) final states, and η is reconstructed via its γγ (mode I) and π + π − π 0 (mode II) decay modes. The samples used in this analysis include a set of high luminosity data samples with more than 50 pb −1 at each c.m. energy adding up to a total integrated luminosity of 13.1 fb −1 (referred to as "XY Z data") [21]. Compared to the previous analysis [20], 10 high luminosity data sets [21], which are of luminosity greater than 500 pb −1 individually and with c.m. energy around 4200 MeV/c 2 , are added. A set of data samples of about 7-9 pb −1 at each c.m. energy with a total integrated luminosity of 0.8 fb −1 (named as "scan data") [21] was used in this study, in addition, which is not available in the earlier study of Ref. [20]. The additional data obtained from mode II, which was absent from the previous measurement [20], are about one quarter of the total statistics collected with mode I.
Details on the features and capabilities of the BEPCII collider and the BESIII detector can be found in Ref. [22]. The geant4-based [23] Monte Carlo (MC) simulation software package boost [24], which includes the geometric description of the BESIII detector and the detector response, is used to optimize event selection criteria, determine the detection efficiencies, and estimate the background events. Signal MC samples of e + e − → ηJ/ψ with the corresponding J/ψ and η decay modes are generated using helamp [25] and evtgen [26] at each c.m. energies. The ISR is simulated with kkmc [27] by requiring a maximum energy of the ISR photon corresponding to the ηJ/ψ mass threshold. Final-state radiation (FSR) is simulated with photos [28]. Possible background contributions are studied with the inclusive MC samples generated by kkmc with comparable luminosity to the XY Z data, where the known decay modes are simulated by evtgen [26] with branching fractions taken from the PDG [1], and the remaining unknown are simulated with the lundcharm model [29].
Candidate events are required to have two (with zero net charge) charged tracks for mode I and four (with zero net charge) charged tracks for mode II. The charged tracks are reconstructed with the hit information in the main drift chamber (MDC). Each charged track is required to be within the polar angle (θ) ranging | cos θ| < 0.93 and to have a point of closest approach to the interaction point (IP) within ±10 cm along the beam direction and 1 cm in the plane perpendicular to the beam. Since pions and leptons have distinct momenta for the signal processes, we assigned charged tracks with momenta (p) larger than 1.0 GeV/c to be leptons, otherwise as pions. The separation of electrons from muons is realized by taking into account the deposited energy (E) in the electromagnetic calorimeter (EMC), i.e., a muon is required to have E < 0.5 GeV, while electron to be E/p > 0.8. The signal candidates are required to have a pair of lepton with same flavor but opposite charge. In mode II, two additional pions with opposite charge are further required. Photon candidates are reconstructed by isolated clusters in the EMC, which are at least 10 • away from the nearest charged track. The photon energies are required to be at least 25 MeV in barrel region (| cos θ| < 0.8), or 50 MeV in end-cap regions (0.86 < | cos θ| < 0.92). A requirement on the EMC timing (0 < t < 700 ns) is implemented to suppress electronic noise and energy deposition unrelated to the event. Candidate events with at least two photons are taken for further analysis.
To improve the kinematic resolution and to suppress the background events, a four-constraint (4C) kinematic fit imposing energy-momentum conservation with the hypothesis of e + e − → γγℓ + ℓ − is applied for the candidates of mode I, while a five-constraint (5C) kinematic fit is performed under the hypothesis of e + e − → γγπ + π − ℓ + ℓ − with additional π 0 mass constraint for the photon pair of the mode II. For events with more than two photon candidates, all photon pairs are tested in the kinematic fit and the combination with the smallest χ 2 4C/5C is retained. The surviving events are further required to satisfy χ 2 4C < 40 or χ 2 5C < 80. To further suppress the background events from the radiative Bhabha and dimuon events associated with a random photon candidate for events of the mode I, the energy of each of two selected photons is required to be larger than 80 MeV. Figure 1 presents the distributions of the invariant mass of the ℓ + ℓ − pair (M (ℓ + ℓ − )) versus that of the γγ pair (M (γγ)) or π + π − π 0 combination (M (π + π − π 0 )) for the surviving events at √ s = 4.1780 GeV after applying the previously described selection criteria. Clear accumulations of candidate events of the signal process e + e − → ηJ/ψ are observed around the intersections of the J/ψ and η mass regions. Signal candidates are required to be within the J/ψ mass region, defined as [3.067, 3.127] GeV/c 2 , which is approximately 3σ of the resolution of the invariant mass distributions of lepton pairs. The events in the J/ψ mass sideband regions, defined as [3.027, 3.057] and [3.137, 3.167] GeV/c 2 , are used to estimate the non-J/ψ background, and non-peaking background are observed in the M (γγ) and M (π + π − π 0 ) distributions. A significantly larger non-J/ψ background is observed in the e + e − mode than in the µ + µ − mode in mode I, which is due to the large Bhabha cross section.
The Born cross section is obtained from where N sig is the signal yield, which will be explained below, L int is the integrated luminosity, (1 + δ) r is the ISR correction factor, (1 + δ) v is the vacuum polarization factor taken from a QED calculation [30], Br is the product of the branching fractions of the subsequent decays of intermediate states quoted from the PDG [1], and ǫ is the detection efficiency obtained from a MC simulation. The ISR correction factor is obtained by using iteratively the QED calculation as described in Ref. [31], where the last measured cross section is taken as the input lineshape. For the XY Z data, an unbinned maximum-likelihood fit is performed on the distributions of M (γγ) and M (π + π − π 0 ) to extract the signal yields, where the signal is described by a MC-simulated shape convolved with a Gaussian function, representing the resolution difference between data and MC simulation, and the back- (Color online) Scatter plots of M (ℓ + ℓ − ) versus M (γγ/π + π − π 0 ) (a, b, e, f) and the spectra of the M (γγ/π + π − π 0 ) distribution (c, d, g, h) in the J/ψ signal region for data at √ s = 4.1780 GeV. The upper panels correspond to the mode I and those at bottom for the mode II.
In the scatter plots, the solid (red) lines denote the signal region, and the dashed (blue) lines for the sideband region. For the mass spectra plots, the dots with error bars represent data. The solid curves correspond to the fit results whereby the long dashed (red) curves for the signal and the short dashed (blue) curves for background.
ground is described by a linear function. A simultaneous fit is performed by considering the four processes, i.e. two observation variables M (γγ) and M (π + π − π 0 ), as well as two J/ψ decay modes e + e − and µ + µ − for the 27 data samples at different c.m. energies. In the fit, the different processes are constrained by the same Born cross section σ B , and the expected signal yields are N sig = σ B · L int · (1 + δ) r · (1 + δ) v · Br · ǫ. Among the different data sets we used common fit parameters for the mean and width of the Gaussian function representing differences between data and MC. For center of mass energies where the signal is not significant, we compute the upper limits on the cross sections at 90% C.L. using a Bayesian method with a flat prior. The optimized likelihoods L are presented as a function of the cross section. The upper limits on the cross section σ UP at 90% C.L. are the values that yield 90% of the likelihood integral over σ from zero to infinity: The systematic uncertainty is taken into account by smearing the posterior distribution.
For the scan data sets, the signal yields are determined by counting the number of events in the signal region after subtracting the background estimated by the normalized number of events in the J/ψ mass sideband region. Only the mode I is considered to extract the Born cross sections by using Eq. (1).
The measured Born cross sections at the different c.m. energies for both XY Z and scan data are shown in the top and bottom panels of Fig. 2, respectively. Clear structures are observed. The numbers used in the calculation of the Born cross section (upper limit at 90% C.L.) are summarized in Tables I and II in the supplemental material [21]. The following sources of systematic uncertainty are considered in the cross section measurements. The uncertainty of the integrated luminosity is 1% measured by analyzing events of the Bhabha scattering process [32]. The uncertainty related with the efficiencies of leptons, pions and photons is 1% for each particle [33,34]. The uncertainties related to the J/ψ mass window requirement and kinematic fit are estimated by tuning the MC sample for the J/ψ mass resolution and the helix parameters of charged tracks [35] according to data, and taking the resulting changes in efficiency as the uncertainties. The uncertainty associated with ISR correction factor is taken to be the difference of (1 + δ) r · ǫ between the last two iterations in the cross section measurement. The uncertainties of the branching fractions of the intermediate states are taken from the PDG [1]. As described above, the signal yields are extracted by performing a simultaneous fit, thus, those uncertainties, which are correlated (i.e. luminosity, lepton and photon efficiencies), are directly propagated to the measured cross sections. Otherwise, we repeated the simultaneous fits by changing the corresponding value by ±1σ, individually, and the largest changes in the results are taken as the uncertainties. To extract the uncertainties associated with the fit procedure, we perform alternative fits by replacing the linear function with a second-order polynomial function for the background, fixing the width of the Gaussian function for the signal to be its nominal value and in addition changing its uncertainty and varying the fit range, individually, and the relative changes in the results are taken as the uncertainties. The efficiencies for the other selection criteria, the trigger simulation, the event start time determination, and the FSR simulation, are quite high (>99%), and their systematic errors are estimated to be less than 1%. Assuming all sources of uncertainties are independent, the total uncertainties in the ηJ/ψ cross section measurement are determined to be 3.5 -13.7% depending on the c.m. energy. In general, the systematical errors are much smaller than the statistical ones. For details, we refer to Table III of the supplemental material [21].
To extract the resonant parameters of the structures observed in the measured cross sections, a simultaneous maximum-likelihood fit is performed to the results extracted from the XY Z and scan data. The fit function is a coherent sum of a P-wave phase space component (P-PHSP) (Φ( √ s)) of the process e + e − → ηJ/ψ and three Breit-Wigner functions (B i=1,2,3 ) for the structures observed around 4040, 4230 and 4360 MeV/c 2 , respectively: where φ i is the relative phase of a resonance (i) to the P-PHSP component, C 0 is a free parameter. The parameterizations of the P-PHSP and Breit-Wigner components are formualted as where q, M i , Γ i , Γ e + e − i and Br i are the daughter momentum in the rest frame of its parent, the mass, width, partial width of the decay in e + e − , and the branching fraction to ηJ/ψ mode for the resonance i, respectively.
The fit is carried out by incorporating the statistics uncertainties only, where the number of events for the scan data are assumed to be Poisson distribution, and those for the XYZ data are Gaussian distribution. Additionally, the beam energy spread of BEPCII (1.6 MeV) is considered by convolving with a Gaussian function whose width is 1.6 MeV [14,36]. The structure around 4040 MeV/c 2 is assumed to be the ψ(4040), and its mass and width are fixed to those given in the PDG [1], due to a lack of data sets at this energy region. Three solutions are found with equal fit quality and with identical masses and widths for the structures around 4220 and 4360 MeV/c 2 . The fit quality is χ 2 /n.d.f. = 107.7/120, estimated by a χ 2 -test approach, where n.d.f. is the number of degree of freedom. The fit results are summarized in Table I and the fit curves (Solution 1) are exhibited in Fig. 2.  The systematic uncertainties of the resonant parameters of the structures at 4220 and 4360 MeV/c 2 and of the product of Br i and Γ e + e − i are discussed as follows. The uncertainties associated with the systematic of measured cross section are estimated by incorporating the correlated and uncorrelated systematic uncertainties of measured Born cross section in the fit. The uncertainty associated with the c.m. energy (0.8 MeV) [12] is common for all data samples and propagates directly to the mass measurement. The uncertainty associated with the fit range is investigated by excluding the last energy point √ s = 4.60 GeV in the fit. The uncertainty from the ψ(4040) resonant parameters is studied by varying the parameters within its uncertainties. We performed the alternative fits with above scenarios, individually, the resultant difference are taken as the systematic uncertainties, and are summarized in Table II.
The structure at 4360 MeV/c 2 is firstly observed in the process e + e − → ηJ/ψ, the corresponding significance is studied by performing an alternative fit without this structure included. The significance is 6.0 σ, calculated with the change of likelihood values and of n.d.f. relative to the nominal fit and incorporated all the uncertainties discussed above.
In summary, we measured the Born cross sections It should be noted that we found a resonant structure with a mass around 4220 MeV/c 2 that is significantly higher than the one (4160 MeV/c 2 ) observed by the Belle experiment [19]. A comparison of masses versus widths for the structures in this measurement as well as those obtained from the processes e + e − → π + π − J/ψ [12], π + π − ψ(3686) [13], π + π − h c [14], ωχ c0 [15] and π + D 0 D * − [16] by the BESIII experiment are presented in Fig. 3. The measured resonant parameters of the two observed structures are consistent with or close to those of previous measurements, however, the intrinsic scenario for the difference on width is still unknown. Assuming that the two observed structures are the Y (4220) and Y (4360), our result would be the first mass, width, and branching fraction measurements of the two Y states decaying into the ηJ/ψ final state. The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center and the supercomputing center of USTC for their strong support.