Evidence for $e^+e^-\to\gamma\chi_{c1, 2}$ at center-of-mass energies from 4.009 to 4.360 GeV

Using data samples collected at center-of-mass energies of $\sqrt{s}$ = 4.009, 4.230, 4.260, and 4.360 GeV with the BESIII detector operating at the BEPCII collider, we perform a search for the process $e^+e^-\to\gamma\chi_{cJ}$ $(J = 0, 1, 2)$ and find evidence for $e^+e^-\to\gamma\chi_{c1}$ and $e^+e^-\to\gamma\chi_{c2}$ with statistical significances of 3.0$\sigma$ and 3.4$\sigma$, respectively. The Born cross sections $\sigma^{B}(e^+e^-\to\gamma\chi_{cJ})$, as well as their upper limits at the 90% confidence level are determined at each center-of-mass energy.

Johannes Gutenberg University of Mainz, Johann-Joachim-Becher-Weg 45, D-55099 Mainz, Germany 22  the BESIII detector operating at the BEPCII collider, we perform a search for the process e + e − → γχcJ (J = 0, 1, 2) and find evidence for e + e − → γχc1 and e + e − → γχc2 with statistical significances of 3.0σ and 3.4σ, respectively. The Born cross sections σ B (e + e − → γχcJ ), as well as their upper limits at the 90% confidence level are determined at each center-of-mass energy.

I. INTRODUCTION
The charmonium-like state Y (4260) was first observed in the initial state radiation (ISR) process e + e − → γ ISR π + π − J/ψ by BaBar [1], and later confirmed by the CLEO [2] and Belle [3] experiments. Recently, both BaBar and Belle updated their results with full data sets, respectively, and further confirmed the existence of the Y (4260) [4,5].
Since it is produced through ISR in e + e − annihilation, the Y (4260) has the quantum numbers J P C = 1 −− . However, there seems to be no cc slot available for the Y (4260) in the conventional charmonium family [6]. In addition, a number of unusual features, such as a strong coupling to hidden-charm final states, suggest that the Y (4260) is a non-conventional cc meson. Possible interpretations of this state can be found in Refs. [7][8][9][10][11], but all need further experimental input.
Most previous studies of the Y (4260) have utilized hadronic transitions. Besides the clear signal observed in the π + π − J/ψ decay mode, the Belle experiment failed to find evidence of the Y (4260) via the e + e − → γ ISR ηJ/ψ process [12]. Based on 13.2 pb −1 of e + e − data collected at √ s = 4.260 GeV, the CLEO experiment investigated fourteen hadronic decay channels, but only few decay modes had a significance more than 3σ [13]. The BESIII Collaboration first observed the process e + e − → γX(3872) using data samples taken between √ s = 4.009 and 4.420 GeV [14], which strongly supports the existence of the radiative transition decays of the Y (4260). To further understand the nature of the Y (4260) state, an investigation into the radiative transitions between the Y (4260) and other lower mass charmonium states, like the χ cJ (J = 0, 1, 2), is important [15,16]. The cross sections of e + e − → γχ cJ have been evaluated theoretically within the framework of NRQCD [16]. Experimentally, the only existing investigation comes from the CLEO experiment [13], which did not observe a signal. The large data sample collected with the BESIII detector provides a good opportunity to deeply investigate these decay modes, which may shed more light on the properties of the Y (4260).
In this paper, we report on a search for e + e − → γχ cJ (J = 0, 1, 2) based on the large e + e − annihilation data samples collected with the BESIII detector at center-of-mass energies (CME) √ s = 4.009, 4.230, 4.260, and 4.360 GeV, where the χ cJ is reconstructed by its γJ/ψ decay mode, and the J/ψ is by its decay to µ + µ − . The decay J/ψ → e + e − is not considered in this analysis due to the high background of Bhabha events. The corresponding luminosities of the data samples at different CME used in this analysis are listed in Table I.

II. BESIII DETECTOR AND MONTE CARLO
The BESIII detector at the BEPCII collider [17] is a large solid-angle magnetic spectrometer with a geometrical acceptance of 93% of 4π solid angle consisting of four main components. The innermost is a small-cell, helium-based (40% He, 60% C 3 H 8 ) main drift chamber (MDC) with 43 layers providing an average single-hit resolution of 135 µm. The resulting charged-particle momentum resolution for a 1 T magnetic field setting is 0.5% at 1.0 GeV/c, and the resolution on the ionization energy loss information (dE/dx) is better than 6%. The next detector, moving radially outwards, is a time-of-flight (TOF) system constructed of 5 cm thick plastic scintillators, with 176 detectors of 2.4 m length in two layers in the barrel and 96 fan-shaped detectors in the end-caps. The barrel (end-cap) time resolution of 80 ps (110 ps) provides a 2σ K/π separation for momenta up to 1.0 GeV. Continuing outward, we have an electromagnetic calorimeter (EMC) consisting of 6240 CsI(Tl) crystals in a cylindrical barrel structure and two end-caps. The energy resolution at 1.0 GeV is 2.5% (5%) and the position resolution is 6 mm (9 mm) in the barrel (end-caps). Finally, the muon counter (MUC) consists of 1000 m 2 of Resistive Plate Chambers (RPCs) in nine barrel and eight end-cap layers, which provides a 2 cm position resolution.
A GEANT4 [18] based Monte Carlo (MC) simulation software, which includes the geometric description of the detector and the detector response, is used to optimize the event selection criteria, determine the detection efficiency, and estimate the potential backgrounds. Signal MC samples of e + e − → γχ cJ are generated for each CME according to the electric-diplole (E1) transition assumption [19]. Effects of ISR are simulated with KKMC [20] by assuming that γχ cJ is produced via Y (4260) decays, where the Y (4260) is described by a Breit-Wigner function with resonance parameters from the Particle Data Group (PDG) [21]. For the background studies, an 'inclusive' Y (4260) MC sample equivalent to 500 pb −1 integrated luminosity is generated which includes the Y (4260) resonance, ISR production of the known vector charmonium states, and events driven by QED processes. The known decay modes are generated with EvtGen [19] with branching fractions set to their world average values in the PDG [21], and the remaining events are generated with Lundcharm [22] or PYTHIA [23].

III. EVENT SELECTION
Charged tracks are reconstructed in the MDC. For each good charged track, the polar angle must satisfy | cos θ| < 0.93, and the point of closest approach to the interaction point must be within ±10 cm in the beam direction and within ±1 cm in the plane perpendicular to the beam direction. The number of good charged tracks is required to be two with a zero net charge. Charged tracks are identified as muons if they have E/p < 0.35 and p > 1.0 GeV/c, where E is the energy deposited in the EMC and p is the momentum measured by the MDC.
Photons are reconstructed from isolated showers in the EMC that are at least 20 degrees away from any of the charged tracks. To improve the reconstruction efficiency and the energy resolution, the energy deposited in the nearby TOF counters is included. Photon candidates are required to have energy greater than 25 MeV in the EMC barrel region (| cos θ| < 0.8), and 50 MeV in the end-cap region (0.86 < | cos θ| < 0.92).
In order to suppress electronic noise and energy deposits that are unrelated to the event, the EMC time t of the photon candidates must be in coincidence with collision events within the range 0 ≤ t ≤ 700 ns. At least two photon candidates in the final state are required.
To improve the momentum resolution and to reduce backgrounds, a kinematic fit with five constraints (5C-fit) is performed under the e + e − → γγµ + µ − hypothesis, imposing overall energy and momentum conservation and constraining the invariant mass of µ + µ − to the nominal J/ψ mass. Candidates with a χ 2 5C < 40 are selected for further analysis. If more than one candidate occurs in an event, the one with the smallest χ 2 5C is selected. Due to the kinematics of the reaction, the first radiative photon from e + e − → γχ cJ has a high energy while the second radiative photon from χ cJ → γJ/ψ has a lower energy at √ s = 4.230, 4.260, and 4.360 GeV. The invariant mass of the low energy photon and J/ψ, M γJ/ψ , is used to search for χ cJ signals. However, for the data sample taken at √ s =4.009 GeV, there is an overlap of the energy distributions of the photons from e + e − → γχ c1,2 and from χ c1,2 decays, as shown in Fig. 1. To separate the overlapping photon spectra, the energy of photons from χ c1,2 decays is further required to be less than 0.403 GeV at √ s =4.009 GeV.

IV. BACKGROUND STUDY
The potential backgrounds from e + e − → P + J/ψ, P → γγ (P = π 0 , η, or η ′ ) can be rejected by requiring where M γγ is the invariant mass of two selected photons. The background from e + e − → γ ISR ψ(3686), ψ(3686) → γχ cJ is rejected by applying the 5C kinematic fit. After imposing all the selection criteria above, the remaining dominant background is from radiative dimuon events, which is not expected to peak in the M γJ/ψ distribution. This has been validated by a dedicated simulation study. For other remaining backgrounds, such as e + e − → π 0 π 0 J/ψ, only few events (normalized to data luminosity) survive and can be neglected.

V. FIT TO THE MASS SPECTRUM
The resulting M γJ/ψ distributions, after applying the above selection criteria, at √ s = 4.009, 4.230, 4.260 and 4.360 GeV are shown in Figure 2. An unbinned maximum likelihood fit of the M γJ/ψ distribution is performed to extract the numbers of χ cJ signal events. In the fit, the shapes of the χ cJ signals are described by double Gaussian functions, where the means and the standard deviations of the double Gaussian functions are determined from a fit to the corresponding signal MC samples at √ s = 4.260 GeV. These shapes are also used for the other three CME points, as the resolution varies only mildly between √ s = 4.009 − 4.360 GeV. This has been validated by MC simulation. Since the dominant background comes from radiative dimuon events, the corresponding MC simulation is used to represent the background shape. To reduce the effect of statistical fluctuations, the dimuon MC shape is smoothed before it is taken as the background function. Figure 2 also shows the fitted results for the M γJ/ψ distribution at different CME. The number of fitted χ cJ signal events, as well as the corresponding statistical significances (calculated by comparing the fit log likelihood values with and without the χ cJ signal) at the four CME points are listed in Table II. The same fit is applied to the sum of M γJ/ψ distributions of the four CME points. The statistical significances for χ c0 , χ c1 and χ c2 are found to be 1.2σ, 3.0σ and 3.4σ, respectively. As a test, we perform similar analyses to control samples from the J/ψ sideband regions, 2.917 < M µ + µ − < 3.057 GeV/c 2 and 3.137 < M µ + µ − < 3.277 GeV/c 2 , by constraining the invariant mass of µ + µ − to 3.047 or 3.147 GeV/c 2 in 5C-fit, and find no obvious χ cJ signals.

VI. RESULTS
The Born cross section at different CME is calculated with where N obs is the number of observed events obtained from the fit, L is the integrated luminosity, 1 + δ r is the radiative correction factor for χ cJ with the assumption that the e + e − → γχ cJ cross section follows the Y (4260) Breit-Wigner line shape [24], 1 + δ v is the vacuum polarization factor [25], B is the combined branching ratio of χ cJ → γJ/ψ and J/ψ → µ + µ − , and ǫ is the detection efficiency. The detection efficiencies, radiative correction factors as well as the calculated Born cross sections at different CME are shown in Table II. The much lower efficiencies for χ c1,2 at √ s = 4.009 GeV are due to the requirement on the photon energy used to separate the overlapping photon spectra as described in Section III.
Since the χ cJ signals are not statistically significant at the individual CME points, we also give in Table II the by assuming the non-existence of signals, the upper limits on the Born cross sections at the 90% confidence level (C. L.) under the assumption that no signals are present. The upper limits are derived using a Bayesian method [21], where the efficiencies are lowered by a factor of (1 − σ sys ) to take systematic uncertainties into account.
We also perform a simultaneous fit to the M γJ/ψ distribution at √ s = 4.009, 4.230, 4.260, and 4.360 GeV, assuming the production cross section of e + e − → γχ cJ at different CME point follows the line shape of the Y(4260) state. In the fit, the line shapes of the χ cJ signals and the background are as same as those in previous fits, and the number of χ cJ events at each CME point is expressed as a function of ǫ c.m. L c.m. R c.m. (1 + δ r ), where ǫ c.m. and L c.m. are the detection efficiency and luminosity, respectively, and R c.m. is the ratio of the cross section calculated with the Y (4260) line shape (a Breit-Wigner function with parameters fixed to the PDG values) at different CME points to that at √ s = 4.260 GeV. The corresponding fit result is shown in Fig. 3, and the statistical significances for χ c0 , χ c1 and χ c2 signals are 0.0σ, 2.4σ and 4.0σ, respectively.  The solid lines show the total fit results. The χcJ signals are shown as dashed lines, dotted lines, and dash-dotted lines, respectively, and the backgrounds are indicated by red dashed lines.

VII. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties in the cross section measurements of e + e − → γχ cJ caused by various sources are par-tially in common for all channels. The common sources of systematics include the luminosity measurement, reconstruction efficiencies for charged tracks and photons, the vacuum polarization factor, kinematic fit and branching fractions of TABLE II. The results on e + e − → γχcJ Born cross section measurement. Shown in the table are the significance σ, detection efficiency ǫ, number of signal events from the fits N obs , radiative correction factor (1 + δ r ), vacuum polarization factor (1 + δ v ), upper limit (at the 90% C.L.) on the number of signal events N UP , Born cross section σ B and upper limit (at the 90% C.L.) on the Born cross section σ UP at different CME point. The first uncertainty of the Born cross section is statistical, and the second systematic.  The systematic uncertainty due to the luminosity measurement is estimated to be 1.0% using Bhabha events. The uncertainty related to the track reconstruction efficiency of high-momentum muons is 1.0% per track [26]. The systematic uncertainty related to the photon detection is estimated to be 1.0% per photon [14]. The systematic uncertainty due to 5C-fit is 0.6%, obtained by studying a control sample of ψ ′ → ηJ/ψ decays. The uncertainty related to the branching fractions of χ cJ and J/ψ decays are taken from the PDG [21]. The uncertainty for the vacuum polarization factor is 0.5% [25].
The other systematic uncertainties arising from the χ cJ mass resolution, the shift of the χ cJ reconstructed mass, the MC model, the shape of background, the radiative correction factor and the fit range at different CME points are discussed below.
The ψ ′ → γχ cJ channel is employed as a control sample to extract the differences on the mass resolution of the χ cJ signal by fitting the M γJ/ψ spectrum. The differences in the mass resolutions between data and MC are found to be 0.02%, 0.01%, 0.2% for χ cJ (J = 0, 1, 2). A similar fit is performed, in which the signal shapes are smeared to compensate for the mass resolution difference, and the differences on the yields of χ cJ signal are taken as the systematic uncertainties due to An alternative fit is performed shifting the mean of χ cJ signal shapes by one standard deviation of the PDG values, and the deviations of the signal yields to the nominal values are taken as the systematic uncertainties due to the uncertainties of the signal line shapes.
The detection efficiency is evaluated using MC samples based on the E1 transition assumption [19] for Y (4260) → γχ cJ . Another set of MC samples is generated where the Y (4260) → γχ cJ decay is modeled using a phase space distribution, and the differences of the detector efficiencies between the two sets of MC samples are treated as systematic uncertainties from the MC model.
To estimate the systematic uncertainty related to the background shape, a control sample is selected from the data by requiring a µ + µ − pair and at least one photon. An alternative background shape is then extracted by re-weighting the γµ + µ − invariant mass spectrum of the control sample, where the weights are the efficiency ratio of e + e − → (nγ)µ + µ − MC simulated events surviving the signal selection criteria to the same selection criteria for the control sample. A fit with the alternative background shape is performed, and the differences between the yields of χ cJ signal to the nominal ones are taken as the systematic uncertainties due to the shape of background.
The possible distortions of the Y (4260) line shape due to interference effects with nearby resonances could introduce uncertainties in the radiative correction factor ǫ × (1 + δ r ). To estimate the related systematic uncertainties, we instead assume that e + e − → γχ cJ are produced via ψ(4040) decays at √ s = 4.009 GeV, ψ(4160) decays at √ s = 4.229 and 4.260 GeV, and ψ(4415) decays at √ s = 4.360 GeV. The variations in the factor ǫ × (1 + δ r ) are taken as the systematic uncer-tainties due to the radiative correction factor.
A series of similar fits are performed in different ranges of the M γJ/ψ distribution, and the largest differences on the signal yields to the nominal values are taken as systematic uncertainties.
All the systematic uncertainties from the different sources are summarized in Table III. The total systematic uncertainties are calculated as the quadratic sum of all individual terms.

VIII. SUMMARY
Using data samples collected at CME of √ s = 4.009, 4.230, 4.260, and 4.360 GeV with the BESIII detector, we perform a search for e + e − → γχ cJ (J = 0, 1, 2) with the subsequent decay χ cJ → γJ/ψ and J/ψ → µ + µ − . We find evidence for the processes e + e − → γχ c1 and e + e − → γχ c2 with statistical significances of 3.0σ and 3.4σ, respectively. No evidence of e + e − → γχ c0 is observed. The corresponding Born cross sections of e + e − → γχ cJ at different CME are calculated and listed in Table II. Under the assumption of the absence of χ cJ signals, the upper limits on the Born cross sections at the 90% C.L. are calculated and listed in Table II, too. These upper limits on the Born cross section of e + e − → γχ cJ are compatible with the theoretical prediction from an NRQCD calculation [16].