Measurement of $e^+ e^- \rightarrow \phi \eta^{\prime}$ cross sections at center-of-mass energies between 3.508 and 4.600 GeV

We present a measurement of the dressed cross sections for $e^+ e^- \rightarrow \phi \eta^{\prime}$ at different center-of-mass energies between 3.508 and 4.600 GeV based on 15.1 fb$^{-1}$ of $e^+ e^-$ annihilation data collected with the BESIII detector operating at the BEPCII collider. In addition, a search for the decay $Y(4230) \to \phi \eta^{\prime}$ is performed. No clear signal is observed and the corresponding upper limit is provided.

We present a measurement of the Born cross sections for e + e − → φη ′ at different center-ofmass energies between 3.508 and 4.600 GeV based on 15.1 fb −1 of e + e − annihilation data collected with the BESIII detector operating at the BEPCII collider. In addition, a search for the decay Y (4230) → φη ′ is performed. No clear signal is observed and the corresponding upper limit is provided.

I. INTRODUCTION
With the observation of a series of new charmoniumlike states, the so-called XY Z states, charmonium physics has seen a resurgence of interest from both theory and experiment. These charmonium-like states do not fit in the conventional cc charmonium spectroscopy and could be exotic states that lie outside the naive quark model [1]. A better understanding of these states would shed light on the non-perturbative regime of the strong interaction. The first observed Y state, the Y (4260), was found by the BaBar collaboration in the initial state radiation (ISR) process e + e − → γ ISR J/ψπ + π − [2]; it was confirmed by CLEO-c [3], Belle [4] and in an updated analysis by the BaBar collaboration [5]. In later experiments, the Y (4260) was also observed in a se-ries of processes measured by the BESIII collaboration, such as e + e − → π 0 π 0 J/ψ [6], ωχ c0 [7], π + π − h c [8], π + π − J/ψ [9], π + π − ψ(3686) [10], and π + D 0 D * − [11]. Furthermore, evidence for transitions from the Y (4260) to other charmonium-like states, such as the X(3872) and Z c (3900), have been reported [12,13]. These new measurements at the BESIII experiment also led to a downward shift in the mass of the Y (4260), so it has been renamed the Y (4230) [14]. In the remainder of this paper, we will use Y (4230) to represent this state * .
In this paper, we utilize data samples collected by the BESIII detector to search for a new Y (4230) decay mode by measuring the Born cross sections for e + e − → φη ′ at center-of-mass energies between 3.508 and 4.600 GeV, as summarized in Table I. This measurement is also an extension of a previous measurement of the same process in a lower energy region, performed in the vicinity of the φ(2170) by the BESIII collaboration [32]. The φ meson is reconstructed through its K + K − decay mode, and the η ′ through both its γπ + π − decay (mode I) and its decay to ηπ + π − with η → γγ (mode II). The sum of data or Monte Carlo (MC) simulated samples at all 20 energy points are used hereafter unless explicitly stated.

II. BESIII DETECTOR AND DATA SAMPLES
The BESIII detector is a magnetic spectrometer [33] located at the Beijing Electron Positron Collider (BEPCII) [34]. The cylindrical core of the BESIII detector consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over 4π solid angle. The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the resolution of specific ionization measured in the MDC (dE/dx) is 6% for the electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region. The time resolution of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps. The end cap TOF system was upgraded in 2015 with multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [35].
MC simulated data samples are used to determine the detection efficiency and to estimate the backgrounds. They are produced with a geant4-based [36] software package that includes the geometric description of the BESIII detector and the detector response. To study the efficiency of each final state, a sample of 1 × 10 5 events is generated at each energy point, and the simulation includes the beam energy spread and ISR in the e + e − annihilations modeled with kkmc [39]. For the signal process, e + e − → φη ′ is generated using the vvs pwave decay model [37]; η ′ → γπ + π − and φ → K + K − are simulated with the ρ 0 − ω − box anomaly model considered [38] and the vss decay model [37], respectively; other decay modes are generated with phase space (PHSP) distributions. For the background study, three generic MC samples at the energies of 3.773, 4.178 and 4.226 GeV are generated. The known decay modes of the charmonium states are produced with evtgen [40,41] according to the world average branching fraction (BF) values [14], while the unknown decay modes are generated with lundcharm [42,43]. Final state radiation (FSR) from charged final state particles is incorporated using the photos package [44]. Continuum hadronic events are generated with kkmc [39] and QED processes such as Bhabha scattering, µ + µ − , τ + τ − , and γγ events are generated with kkmc [39] and babayaga [45].

A. Event selections and background analysis
Charged tracks are required to have a polar angle θ with respect to the detector axis within the MDC acceptance | cos θ| < 0.93, and a distance of closest approach to the interaction point within 10 cm along the beam direction and 1 cm in the plane perpendicular to the beam direction. The particle type for each charged track is determined by selecting the hypothesis with the highest probability, which is calculated with the combination of time information from the TOF and dE/dx for different particle hypotheses.
Photon candidates are reconstructed from isolated electromagnetic showers in the EMC. The energy of a photon candidate is required to be larger than 25 MeV (50 MeV) in the barrel (end cap) region, corresponding to an angular coverage of | cos θ| < 0.80 (0.86 < | cos θ| < 0.92). The electromagnetic shower time from the EMC has to be within 700 ns of the event start time to suppress electronic noise and energy deposition unrelated to the event of interest. To eliminate the showers associated with charged particles, the opening angle between a photon candidate and the extrapolated position of the closest charged track must be larger than 10 degrees.
The signal candidates for the e + e − → φη ′ process are selected by requiring four charged tracks with net charge zero and identified as K + , K − , π + and π − , as well as at least one (two) photon(s) for mode I (II). To improve the resolution and suppress backgrounds, a four-constraint (4C) kinematic fit is performed for the decay mode I, constraining the total four-momentum of the final-state particles to the total initial four-momentum of the colliding beams. For the mode II, a five-constraint (5C) kinematic fit is performed with an additional constraint of the invariant mass of the two photons to the world average η mass [14]. If there is more than one combination in an event, the one with the smallest kinematic fit χ 2 is selected. The χ 2 of the candidate events are required to be less than 50 for both modes.
The signal yields N are determined from the distributions of M K + K − versus M γπ + π − for mode I (a) and and N E represent the numbers of events observed in the equal-area regions X, A, B, C, D, and E, as shown in Fig. 2. The boundaries of these regions correspond to the previous definitions of the signal and sidebands. The backgrounds due to misreconstruction of the η ′ are assumed to be linear in the M γπ + π − and M ηπ + π − distributions, and are estimated using the number of events in the regions A and B. The background with correctly reconstructed η ′ but no φ is estimated using the M K + K − sideband region D. The regions C and E represent the non-resonant background without a φ or an η ′ .
The ratio of non-φ backgrounds under the M K + K − peak over that in the sideband region of M K + K − is defined as r, which is evaluated to be 0.66 and 0.39 for mode I and mode II, respectively. To obtain the value of r for the two modes, maximum likelihood fits are performed on the M K + K − distributions for the two modes, as shown in Fig. 3. The probability density function of the M K + K − spectra for the φ is obtained from a P -wave Breit-Wigner function convolved with a Gaussian function that accounts for the detector resolution. The Pwave Breit-Wigner function is defined as FIG. 1. Invariant mass spectra of the φ candidates after selecting the η ′ signal regions for mode I (a) and for mode II (b), and the invariant mass spectra for the η ′ candidates after selecting the φ signal regions for mode I (c) and for mode II (d). The black dots with error bars are experimental data, the red line is the signal MC, the slate blue arrows mark the signal regions, and the lime green arrows denote the sidebands. The MC simulation is arbitrarily normalized to data.
where m 0 is the nominal φ mass as specified in the Particle Data Group [14], p is the momentum of the kaon in the rest frame of the K + K − system, p ′ is the momentum of the kaon at the nominal mass of the φ, and Γ 0 is the width of the φ. The angular momentum (ℓ) is assumed to be equal to one, which is the lowest allowed given the parent and daughter spins, B(p) is the Blatt-Weisskopf form factor, and R is the radius of the centrifugal barrier, whose value is taken to be 3 GeV/c −1 [46]. The background shape is described by an ARGUS function [47]. The parameters of the Gaussian function and the ARGUS function are left free in the fit. The data is described by an incoherent sum of signal and background contributions.
To validate our analysis method, we also perform the same treatment for generic MC samples to validate the reliability of the method at the energies 3.773 and 4.178 GeV. The signal obtained from Eq. (1) is consistent with zero, as expected for the background-only generic MC sample. 1.08 , where the red rectangle shows the signal region, and the slate blue and lime green rectangles show the 2D sidebands.
To determine the statistical uncertainties of the number of signal events, we construct a likelihood function L(s) as Eq.  signal events N sig , and the asymmetric statistical uncertainties σ + sig and σ − sig can be obtained from the following formulas: The Born cross section is calculated by where L is the integrated luminosity; ǫ is the detection efficiency; B is the product of branching fractions, i.e. B(φ → K + K − ) · B(η ′ → γπ + π − ) for mode I, and 1 |1−Π(s)| 2 is the vacuum polarization correction factor, which is taken from Ref. [48]; (1 + δ ISR ) is the ISR correction factor, which is obtained from a QED calculation taking the line shape of e + e − → φη ′ cross sections at the 20 energy points shown in Table I as input, and is calculated in an iterative procedure until the variation of the correction factor is below 1% compared to the previous iteration. The Born cross sections are measured with asymmetric statistical uncertainties separately for the two decay modes η ′ → γπ + π − and η ′ → ηπ + π − . The cross sections are taken as the parameters of a variable Gaussian as defined in Ref. [49] to construct the likelihood function to give a combined result. The variable Gaussian form can be expressed as where σ represents the combined cross section, and σ i , σ + i and σ − i are the cross sections and their asymmetric statistical uncertainties for the decay mode i. The asymmetric statistical uncertainties of combined cross section is obtained according to ∆ ln L = −1/2. The results of the Born cross section measurements are summarized in Table I.
Since there is no obvious structure in the line shape of the Born cross sections of e + e − → φη ′ , as shown in Fig. 4, the upper limit of Y (4230) → φη ′ is determined by fitting the line shape using a coherent sum of the continuum and the resonance Y (4230) amplitude with Eq. (11) where f con and n are the fit parameters for the continuum process; φ is the relative phase between the continuum and resonant amplitudes; Γ and Γ ee are the total width and partial width to e + e − , respectively; B φη ′ is the branching fraction for the resonance decay into φη ′ , and M is the mass of the resonance. The mass and total width of the Y (4230) are set to the world average values 4222.7 ± 2.6 MeV and 49 ± 8 MeV [14]. We vary the product Γ ee · B φη ′ with fixed step size. For each value, the correlations among different data points are considered in the fit with a fitting estimator Q 2 constructed as Eq. (12), which is minimized by minuit [50].
Here σ Bi and σ Bi fit are the measured and fitted Born cross section of the i th energy point, respectively; δ i is the energy-dependent part of the total uncertainty at each energy point, which includes the statistical uncertainty and the energy-dependent contribution to the systematic uncertainty; δ c is the relative systematic uncertainty corresponding to the energy-independent part, and h is a free parameter introduced to take into account the correlations [51]. Then we construct the likelihood by L = e −0.5Q 2 , whose normalized distribution is used to get the upper limit of Γ ee ·B φη ′ at the 90% confidence level (C.L.). Two solutions with the same minimum value of Q 2 are found with different interference between the two amplitudes, where the second solution can also be derived from the first solution by a numerical method [52]. The fit results are shown in Fig. 4 (the line shapes of the two solutions are identical) and summarized in Table II. Since we cannot determine which solutions is real, we only set the larger one as the upper limit of Γ ee · B φη ′ to be 0.50 eV conservatively.

IV. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties on the cross sections mainly come from the center-of-mass energy, luminosity, tracking efficiency, PID, photon and η reconstruction, ISR correction factor, vacuum polarization factor, quoted BFs and kinematic fit, which are energy independent. The uncertainty of the center-of-mass energy is less than 0.01%, measured by analyzing the di-muon process e + e − → γ ISR/FSR µ + µ − , and is negligible [53]. The uncertainty on the integrated luminosity is estimated to be 1.0% using events from large-angle Bhabha scattering [54]. The tracking and PID differences in the efficiencies between data and MC simulation are studied using the control samples J/ψ → K 0 S K ± π ∓ , J/ψ → π + π − π 0 and J/ψ → K + K − π 0 , and the systematic uncertainties of tracking and PID are both determined to be 4% (1% per track) [55]. The systematic uncertainty from the photon detection efficiency is determined to be 1% per photon by utilizing a control sample of J/ψ → ρ 0 π 0 with ρ 0 → π + π − and π 0 → γγ [56]. The uncertainty from the η selection is 1% per η, which is determined from a control sample of J/ψ → ηpp [56]. For the ISR correction factor, we use a power function to parameterize the line shape of the cross sections, then change the line shape by using or not using the low energy points from 2.900 to 3.080 GeV [32], and take the difference in the cross sections as the systematic uncertainty. The uncertainty associated with the vacuum polarization factor is found to be negligible [57]. The uncertainties of the quoted BFs are taken from [14]. The systematic uncertainty due to the kinematic fit is estimated by correcting the track helix parameters of charged tracks and the corresponding covariance matrix for the signal MC sample to improve the agreement between data and MC simulation. The detailed method can be found in Ref. [58]. The resulting change of the detection efficiency with respect to the one obtained without the corrections is taken as -0.5 0.2 MC statistics 0.7 0.8 0.5 Total 6.8 6.6 6.5 the systematic uncertainty. The systematic uncertainties from the mass interval for the signal and sideband regions are estimated for all energy points with the largest data sample at 4.178 GeV. The systematic uncertainties associated with the signal regions of the φ and η ′ are estimated by changing the φ region from (1.010, 1.034) to (1.0076, 1.0364) and the η ′ region from (0.940, 0.975) to (0.93825, 0.97675) GeV/c 2 . The uncertainties due to the sideband regions are determined by changing the mass windows to M K + K − ∈ (1.0588, 1.0852) GeV/c 2 , M γπ + π − and M ηπ + π − ∈ (0.8832, 0.9218) ∪ (0.9933, 1.0318) GeV/c 2 . The differences in the efficiencies are taken as the corresponding systematic uncertainties. The uncertainty in MC statistics is obtained by 1 where the ǫ is detection efficiency and N is the total number of the generated MC events.
Due to correlations between the two decay modes, the combined systematic uncertainties ζ i are calculated with ζ = (w 1 ζ 1 ) 2 + (w 2 ζ 2 ) 2 + 2w 1 w 2 ρ 12 ζ 1 ζ 2 where w 1 , w 2 represent B(η ′ → γπ + π − ) · ǫ 2 and B(η ′ → ηπ + π − ) · B(η → γγ) · ǫ 1 respectively. The corresponding uncertainties for mode I and mode II are denoted ζ 1 and ζ 2 , respectively. The correlation coefficients ρ 12 are taken as 0 for items that are uncorrelated between the two modes, such as η reconstruction, B(η ′ → ηπ + π − ), B(η ′ → γπ + π − ), B(η → γγ) and MC statistics. Other systematic effects are correlated between the two modes and ρ 12 is taken as 1. Table III summarizes all the systematic uncertainties related to the cross section measurement of the e + e − → φη ′ process for the individual decay modes and the combined one, respectively. The overall systematic uncertainties are obtained by adding all the sources of systematic uncertainties in quadrature, assuming they are all uncorrelated.

V. SUMMARY
The Born cross sections of e + e − → φη ′ are measured with data samples collected with the BESIII detector operating with the BEPCII collider at 20 energy points between √ s = 3.508 and 4.600 GeV. The line shape of the Born cross sections is consistent with the continuum process. A fit with an additional resonance is performed to search for the decay Y (4230) → φη ′ . No clear resonant structure is observed, and an upper limit on the product Γ ee × B(Y (4230) → φη ′ ) at the 90% C.L. is determined to be less than 0.50 eV.