Observation of $e^+e^- \to \gamma \chi_{c1}$ and search for $e^+e^- \to \gamma \chi_{c0}, \gamma \chi_{c2},$ and $\gamma\eta_c$ at $\sqrt{s}$ near 10.6 GeV at Belle

Using data samples of 89.5~fb$^{-1}$, 711.0~fb$^{-1}$, and 121.4~fb$^{-1}$ collected with the Belle detector at the KEKB asymmetric-energy $e^+e^-$ collider at center-of-mass energies 10.52 GeV, 10.58 GeV, and 10.867 GeV, respectively, we study the exclusive reactions $e^+e^-\to\gamma\chi_{cJ}$ ($J = 0,~1,~2$) and $e^+e^-\to\gamma\eta_c$. A significant $\gamma \chi_{c1}$ signal is observed for the first time at $\sqrt{s}=10.58$ GeV with a significance of $5.1\sigma$ including systematic uncertainties. No significant excesses for $\gamma \chi_{c0}$, $\gamma \chi_{c2}$, and $\gamma \eta_c$ final states are found, and we set 90\% credibility level upper limits on the Born cross sections ($\sigma_{\rm B}$) at 10.52 GeV, 10.58 GeV, and 10.867~GeV. Together with cross sections measured by BESIII at lower center-of-mass energies, the energy dependency of $\sigma_{\rm B}(e^+e^-\to\gamma\chi_{c1})$ is obtained.

Using data samples of 89.5 fb −1 , 711.0 fb −1 , and 121.4 fb −1 collected with the Belle detector at the KEKB asymmetric-energy e + e − collider at center-of-mass energies 10.52 GeV, 10.58 GeV, and 10.867 GeV, respectively, we study the exclusive reactions e + e − → γχcJ (J = 0, 1, 2) and e + e − → γηc. A significant γχc1 signal is observed for the first time at √ s = 10.58 GeV with a significance of 5.1σ including systematic uncertainties. No significant excesses for γχc0, γχc2, and γηc final states are found, and we set 90% credibility level upper limits on the Born cross sections (σB) at 10.52 GeV, 10.58 GeV, and 10.867 GeV. Together with cross sections measured by BESIII at lower center-of-mass energies, the energy dependency of σB(e + e − → γχc1) is obtained. The production of heavy quark pairs in high-energy lepton collisions is described well by perturbative Quantum Chromodynamics (pQCD). Yet a description of these pairs forming quarkonium-charmonium or bottomonium-is theoretically challenging. Quarkonium formation is governed by nonperturbative longdistance effects [1]. Nonrelativistic Quantum Chromodynamics (NRQCD) factorization was used to compute the cross section for several processes, including the doublecharmonium production cross section [2,3], e + e − → γχ cJ (J = 0, 1, 2) [4][5][6][7][8] and e + e − → γη c [5,6,9,10] at B factories with relativistic and higher-order corrections included.
Electromagnetic quarkonium production is relatively simpler than other production mechanisms, and therefore it serves as a good testing ground for such NRQCD predictions. The BESIII experiment measured e + e − → γχ cJ cross sections at √ s = 4.01 GeV, 4.23 GeV, 4.26 GeV, and 4.36 GeV and e + e − → γη c cross section at the same energies and additionally at 4.42 GeV and 4.60 GeV [11,12]. At none of the individual energy points does the statistical significance for production of χ cJ or η c exceed 3σ, and when the data from all energy points are combined, the statistical significances for χ c1 , χ c2 , and η c are 3.0σ, 3.4σ, and greater than 3.6σ, respectively.
In this paper, we report cross-section measurements for the exclusive reactions e + e − → γχ cJ and γη c with data recorded at √ s ∼ 10.6 GeV by the Belle experiment at the KEKB asymmetric-energy e + e − collider [17,18]. The data used in this analysis corresponds to 89.5 fb −1 of integrated luminosity at 10.52 GeV, referred to as the continuum sample; 711 fb −1 at 10.58 GeV, referred to as the Υ(4S) sample; and 121.4 fb −1 at 10.867 GeV, referred to as the Υ(5S) sample.
The Belle detector [19,20] is a large solid-angle magnetic spectrometer that consists of a silicon vertex detector, a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), a barrellike arrangement of time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter comprised of CsI(Tl) crystals (ECL) located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux-return yoke instrumented with resistive plate chambers (KLM) located outside the coil is used to detect K 0 L mesons and to identify muons. The Belle detector is described in detail elsewhere [19,20].
We determined event-selection criteria using a large sample of Monte Carlo (MC) signal events (100k) for e + e − → γχ cJ and γη c at √ s = 10.52, 10.58, and 10.867 GeV generated with EvtGen [21]. In the generator, the polar angle of the transition photon in the e + e − C.M. system (θ γ ) is distributed according to (1 + cos 2 θ γ ) for γχ c0 and γη c production, and (1 + 0.63 cos 2 θ γ ) for γχ c1 production [22]. No definite model exists for the distribution of θ γ in γχ c2 production because the combination of tensor-meson production and γ emission is theoretically complicated and requires experimental input. So we model the production of this channel as evenly distributed in phase space and account for differences from (1 ± cos 2 θ γ ) distributions as systematic uncertainties.
Corrections due to initial-state radiation (ISR) are taken into account in all studied channels, where we assume σ(e + e − → γχ cJ /η c ) ∼ 1/s n in the calculation of the radiative-correction factor. The values of n, determined from Refs. [8,10] s , and e + e − → qq (q = u, d, s, c) are checked with a MC sample four times larger than the data sample, and are also gen-erated with EvtGen [21]. GEANT3 [23] is used to simulate the detector response to all MC events. The χ cJ candidates are reconstructed from their decays to γJ/ψ with J/ψ → µ + µ − , and the η c candidates are reconstructed from five hadronic decays into K 0 S K + π − , π + π − K + K − , 2(π + π − ), 2(K + K − ), and 3(π + π − ) [24].
We define a well-reconstructed charged track as having impact parameters with respect to the nominal interaction point of less than 0.5 cm and 4 cm perpendicular to and along the beam direction, respectively. For a e + e − → γχ cJ candidate event, we require the number of well-reconstructed charged tracks, N trk , to be two, and the net charge be zero. For e + e − → γη c , we require N trk = 6 for the 3(π + π − ) final state and N trk = 4 for the other final states, also with a zero net charge. For the particle identification (PID) of a well-reconstructed charged track, information from different detector subsystems, including specific ionization in the CDC, time measurement in the TOF, and the response of the ACC, is combined to form a likelihood L i [25] for particle species i. Tracks with R K = L K /(L K + L π ) < 0.4 are identified as pions with an efficiency of 96%, while 9% of kaons are misidentified as pions; tracks with R K > 0.6 are identified as kaons with an efficiency of 98%, while 8% of pions are misidentified as kaons.
For muons from J/ψ → µ + µ − , we require at least one of the paired tracks to have R µ = L µ /(L µ + L K + L π ) > 0.95; if one track has R µ < 0.95, it must have associated hits in the KLM agreeing with the extrapolated trajectory provided by the CDC [26]. The efficiency of muonpair identification is 94%.
Using a multivariate analysis with a neural network [27] based on two sets of input variables [28], a K 0 S candidate is reconstructed from a pair of oppositelycharged tracks that are treated as pions. An ECL cluster with energy higher than 50 MeV is treated as a photon candidate if it does not match the extrapolation of any charged track. The photon with the maximum energy in the e + e − C.M. system is taken as the transition photon. Since there are two photons in the e + e − → γχ cJ channel, the transition photon is denoted as γ h , and the one with the second highest energy is denoted as γ l and is taken as the photon from the χ cJ → γJ/ψ decay. We require E(γ l ) > 300 MeV to suppress the backgrounds from fake photons.
A four-constraint (4C) kinematic fit constraining the four-momenta of the final-state particles to the initial e + e − collision system is performed. In e + e − → γχ cJ , an additional constraint is used, constraining the mass of the µ + µ − pair to the J/ψ nominal mass, giving a five-constraint (5C) kinematic fit. Kinematic fits with χ 2 5C < 25 for e + e − → γχ cJ and χ 2 4C < 30 for e + e − → γη c are required to improve the resolutions of the momenta of charged tracks and the energies of photons, and to suppress backgrounds with more than two photons, such as ISR processes.
The invariant mass distribution of the µ + µ − pair from the continuum, Υ(4S), and Υ(5S) data samples, prior to the application of the 5C fit, is shown in Fig. 1, together with the result of fitting the data to the sum of a Gaussian function for the J/ψ and a first order polynomial for the background. In the plot, a clear J/ψ signal is observed. We define the J/ψ signal region as |M µ + µ − − m J/ψ | < 48 MeV/c 2 corresponding to three times the detector resolution, where m J/ψ is the J/ψ mass [29]. After all of the above requirements, some non-peaking background events are observed in the processes e + e − → γχ cJ and γη c at the studied C.M. energy points. Figure 2 shows the γ l J/ψ invariant mass distributions for the data. A clear χ c1 signal is observed in the Υ(4S) data sample, but is not evident in the other two data samples. Unbinned extended maximum-likelihood fits to the M γ l J/ψ distributions are performed to extract the χ cJ signal yields. The χ cJ signal shapes in the fits are a Breit-Wigner (BW) function convolved with a Log-normal [30] function with all the values of the χ cJ resonance parameters fixed from the fits to MC simulations; second-order polynomial functions are used to describe the background distributions. The MC-simulated χ cJ signals have mass resolutions around 6 MeV/c 2 with small low-mass tails due to the measurement of E(γ l ). The results from the fits are listed in Table I. The statistical significances of the χ cJ signals are calculated using −2 ln(L 0 /L max ), where L 0 and L max are the maximized likelihoods of the fits without and with the χ cJ signal, respectively. The statistical significance of the χ c1 signal in the Υ(4S) sample is 5.2σ. The signal significance remains at 5.1σ when convolving the likelihood profile with a Gaussian function of width equal to the total systematic uncertainty discussed below. The χ cJ signals in the continuum sample and the Υ(5S) sample are not significant, as indicated in Table I.  Figure 3 shows the η c invariant mass distributions for the five hadronic final states combined. Clear signals resulting from the production of J/ψ by ISR are present, while no significant η c signal is evident.
We perform a simultaneous fit to the five η c final states, in which the ratio of the yields in each channel is fixed to the ratio of B i ε i , where i is the η c decay-mode index, B i is the branching fraction taken from the Particle Data Group (PDG) [29], and ε i is the reconstruction efficiency determined from MC simulation. In the fit, we use a BW function convolved with a Gaussian resolution function to describe the η c signal; the values of all parameters are fixed from the fits to MC simulations. A Gaussian function with free parameters is used to describe the J/ψ signal, and a second-order Chebyshev polynomial function is used for the backgrounds. The fit results are shown in Fig. 3 and summarized in Table I. The Born cross section for e + e − → γX is given by the formula where N obs is the number of signal events obtained from the fit, L is the integrated luminosity of the data sample, B i and ε i are the branching fraction and the detection efficiency of the i-th X decay mode (χ cJ is reconstructed in one decay mode and η c in five decay modes). (1+δ) ISR is the radiative-correction factor, calculated using the formula given in Ref. [31], and |1 − | 2 is the vacuum polarization factor, calculated according to Ref. [32]. The obtained Born cross sections for e + e − → γχ cJ and γη c are listed in Table I together with all the parameter results needed for the cross section calculation. For all processes but e + e − → γχ c1 at √ s = 10.58 GeV, upper limits at 90% credibility level (C.L.) [33] on the numbers of signal events (N UL ) and the Born cross sections (σ UL B ) are determined by solving the equation where x is the assumed signal yield or Born cross section, and F likelihood (x) is the corresponding maximized likelihood from a fit to the data. To take into account the systematic uncertainties discussed below, the likelihood is convolved with a Gaussian function whose width equals the corresponding systematic uncertainty.
Combining the measurement of σ B (e + e − → γχ c1 ) from BESIII [11] and this analysis, we show the cross section as a function of √ s in Fig. 4. We fit these data points with a function proportional to 1/s n assuming that the reaction e + e − → γχ c1 proceeds through the continuum process only: from a fit to the seven points for e + e − → γχ c1 , we find n = 2.1 +0. 3 −0.4 . The significance of the fitted n is 2.2σ, calculated using χ 2 0 − χ 2 min = 2.2, where χ 2 0 is the χ 2 with n fixed at 0, and χ 2 min is the mini-  mum χ 2 with the value of n free, respectively. Adding an additional possible resonance, such as ψ(4040), ψ(4160), Y (4260), or Υ(4S), the largest change in the fitted value of n is 0.3. The result is consistent with the prediction by NRQCD with all leading relativistic corrections included in Ref. [8]. Due to the large uncertainties, we do not fit the √ s dependence of e + e − → γχ c0 , e + e − → γχ c2 , or e + e − → γη c . There are several sources of systematic uncertainty in the Born cross section measurements, including detection efficiency, the statistical error of the MC efficiency, trigger simulation, intermediate state branching fractions, resonance parameters, the distribution of θ γ for e + e − → γχ c2 , fit uncertainty, the s dependence of the cross sections, and the integrated luminosity. The systematic uncertainty for detection efficiency is a finalstate-dependent combined uncertainty for all the different types of particles detected, including tracking efficiency, PID, K 0 S selection, and photon reconstruction. Based on a study of D * + → D 0 (→ K 0 S π + π − )π + , the uncertainty in tracking efficiency is taken to be 0.35% per track. The uncertainties in PID are studied via γγ → ℓ + ℓ − for leptons and a low-background sample of D * decay for charged kaons and pions. The studies show uncertainties of 2.2% for each muon, 1.0% for each charged kaon, and 1.2% for each charged pion.
Comparison of the K 0 S selection efficiencies determined from data and MC results in 1 − ε data εMC = (1.4 ± 0.3)%; 1.7% is taken as a conservative systematic uncertainty. The uncertainty in the photon reconstruction is 2.0% per photon, according to a study of radiative Bhabha events. For each final state, the final detection efficiency uncertainty is obtained by adding all sources in quadrature.
The statistical uncertainty in the determination of efficiency from MC is less than 1.0%. We include uncertainties of 4.8% and 0.6% from trigger simulations for e + e − → γχ cJ and γη c , respectively. The uncertainties from the intermediate decay branching fractions are taken from Ref. [29]. For e + e − → γχ cJ , the total uncertainties from the branching fractions are obtained by adding all relative uncertainties in quadrature. For e + e − → γη c , the total uncertainty from the branching fraction is obtained by summing in quadrature over the five decay modes with weight factors equal to the corresponding efficiency. The uncertainties from the resonance parameters are estimated by changing the values of mass and width of a resonance by 1σ in the fits [29]. Additionally, for the mode e + e − → γχ c2 , the uncertainty from simulating the θ γ dependence is estimated to be 8.2% by comparing the difference between a phase space distribution and the angular distributions of (1±cos 2 θ γ ).
In determining the number of signal events from the fits to data, the fit range and the choice of the function to describe the backgrounds are the main sources of systematic uncertainty. For the latter, the background shapes are replaced by an exponential form or a higherorder Chebyshev polynomial, and the largest difference compared to the nominal fit result is taken as the related systematic uncertainty. Changing the s dependence of the cross sections from fitted values of n to a large number, e.g., n = 4, gives very small differences in the radiative-correction factor (< 1%). The total luminosity is determined to 1.4% precision using wide-angle Bhabha scattering events. All the uncertainties are summarized in Table II and, assuming all the sources are independent, summed in quadrature for the total systematic uncertainties.
In summary, we perform measurements of e + e − → γχ cJ (J = 0, 1, 2) and γη c at √ s = 10.52 GeV, 10.58 GeV, and 10.867 GeV using a 921.9 fb −1 data sample taken by the Belle detector. A clear χ c1 signal is observed at 10.58 GeV with a statistical significance of 5.2σ, and the Born cross section is measured to be (17.3 +4.2 −3.9 (stat.)±1.7(syst.)) fb. For the cases where a χ cJ or η c signal is not evident, upper limits on the Born cross sections are determined at 90% C.L. Using the cross sections measured at three different √ s in this analysis and from BESIII at much lower √ s and assuming the reaction e + e − → γχ c1 proceeds through the continuum process only, we determine the cross section s-dependence to be 1/s 2.1 +0.3 −0.4 ±0.3 for e + e − → γχ c1 . We thank the KEKB group for the excellent operation of the accelerator; the KEK cryogenics group for the efficient operation of the solenoid; and the KEK computer group, the National Institute of Informatics, and the Pacific Northwest National Laboratory (PNNL) Environmental Molecular Sciences Laboratory (EMSL) computing group for valuable computing and Science Information NETwork 5 (SINET5) network support. We acknowledge support from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, the Japan Society for the Promotion of Science (JSPS), and the Tau