Observation of $e^+e^-\to\pi^+\pi^-\pi^0\chi_{b1,2}(1P)$ and search for $e^+e^-\to\phi\chi_{b1,2}(1P)$ at $\sqrt{s}=$10.96---11.05 GeV

We report searches for the processes $e^+e^-\to\pi^+\pi^-\pi^0\chi_{bJ}$ and $e^+e^-\to\phi\chi_{bJ}$ (J=1,2) in the vicinity of the $\Upsilon(11020)$ resonance using energy scan data collected by the Belle experiment at the KEKB collider. We observe $e^+e^-\to \pi^+\pi^-\pi^0\chi_{b1}$ and find evidence for $e^+e^-\to \omega\chi_{b2}$ processes for data with center-of-mass energies from 10.96 to 11.05~GeV for the first time. The statistical significance for $\pi^+\pi^-\pi^0\chi_{b1}$ and $\omega\chi_{b2}$ are $6.1\sigma$ and $4.0\sigma$, respectively. We investigate energy dependence of the $e^+e^-\to\pi^+\pi^-\pi^0\chi_{bJ}$ cross section and find that it is consistent with production of $\Upsilon(10860)$ and $\Upsilon(11020)$ resonances without significant non-resonant contribution. Assuming $e^+e^-\to \pi^+\pi^-\pi^0\chi_{bJ}$ proceeds via $\Upsilon(10860)$ and $\Upsilon(11020)$, the branching fraction $\mathcal{B}(\Upsilon(11020)\to\pi^+\pi^-\pi^0\chi_{bJ})=(8.6\pm4.1\pm6.1^{+4.5}_{-2.5})\times10^{-3}$ is measured for the first time, where the first error is statistical, the second is systematic, and the third from the branching fraction of $\Upsilon(11020)\to e^+e^-$. The signals for $e^+e^-\to \phi\chi_{bJ}$ are not significant, and the upper limits of the Born cross sections at the 90\% confidence level are $0.6$ and $1.0$~pb for $e^+e^-\to \phi\chi_{b1}$ and $\phi\chi_{b2}$, respectively, for center-of-mass energies from 10.96 to 11.05~GeV.

PACS numbers: 13.25.Gv, 12.38.-t Hadronic transitions among heavy quarkonium states serve as a key source of information for better understanding of the strong interaction between a quark-antiquark pair, and thus quantum chromodynamics (QCD). The heavy quarkonium systems are in general non-relativistic and the hadronic transitions to lower lying states have long been described using the QCD multipole expansion [1]. However, the existence of anomalously large hadronic transition rates from the Υ(10860) as reported by the Belle experiment [2][3][4][5][6][7][8][9] challenge the theoretical calculations as well as the pure bottomonium nature of the Υ(10860) and Υ(11020) [10][11][12].
The processes e + e − → ωχ bJ were observed recently [4] using data samples taken at energies near the Υ(10860) peak, but the dependence of the e + e − → ωχ bJ cross section versus energy was not measured. Therefore, it is unclear whether this process occurs from resonant Υ(10860) or continuum process. Nevertheless, the result has been investigated extensively by theorists to understand the dynamics of these transitions, producing studies of Sand D-wave mixing for the observed heavy quark symmetry violation between ωχ b1 and ωχ b2 [13], the possible contribution of Υ(10860) → πZ b → πρΥ(1S) [14], a molecular component in the Υ(10860) wave function [14], and hadronic loop effects [15].
In this Letter, we report the results of a search for Υ(11020) → ωχ bJ and φχ bJ using the Υ(10860) and Υ(11020) energy scan data collected with the Belle detector. Hereinafter, the Υ(10860) and Υ(11020) are referred to, for brevity, as the Υ(5S) and Υ(6S) according to the potential model assignment. The data that we are using consist of 121.4 fb −1 from three energy points near the Υ(5S) peak, 19 data samples with integrated luminosity of about 1 fb −1 per point, and 18 data samples of about 50 pb −1 per point taken in 5 MeV steps between 10.96 and 11.05 GeV [17]. The corresponding luminosities of the Υ(5S) peak data and the 19 samples with high luminosities are listed in Table I as well as the center of mass (c.m.) energies. We use χ bJ → γΥ(1S), Υ(1S) → ℓ + ℓ − (l = e, µ), ω → π + π − π 0 to reconstruct the e + e − → ωχ b1,2 signal; while for e + e − → φχ bJ signal, we reconstruct φ with its decays to K + K − and determine transitions to χ bJ by studying the K + K − recoil mass.
The Belle detector, located at the KEKB asymmetricenergy e + e − collider [18] is described in Ref. [19]. The EVTGEN [20] generator is used to produce simulated events using Monte Carlo (MC) mothods. The nominal parameters of the states in the decay chains are quoted from Ref. [21]. To take the initial state radiation (ISR) into consideration, we use the VECTORISR model [20] in EVTGEN. Generic MC samples at Υ(5S) peak and continuum data, taken at 10.52 GeV, are also used to study the background shape. I. Energy-dependent Born cross sections for e + e − → π + π − π 0 χ bJ at different c.m. energy (in GeV) with corresponding integrated luminosity.
Ec.m. (GeV) L( fb −1 ) σ B (π + π − π 0 χ bJ ) (pb) 10 For charged tracks, the impact parameters perpendicular to and along the beam direction with respect to the interaction point are required to be less than 1.0 and 3.5 cm, respectively. The transverse momentum is restricted to be higher than 0.1 GeV/c. A likelihood L(X) for each charged track is obtained from different detector subsystems for a particle hypothesis X ∈ e, µ, π, K, p (PID). Tracks with a likelihood ratio R(K) = L(K)/(L(K) + L(π)) < 0.4 are identified as pions while those with R(K) > 0.6 are identified as kaons. Similarly, we define the likelihood ratios R(e) and R(µ) for identification of electrons and muons, respectively, with R(e) > 0.01 and R(µ) > 0.1. A neutral cluster in the electromagnetic calorimeter is reconstructed as a photon if it does not match the extrapolated position of any charged track and its energy is greater than 30 MeV.
To select e + e − → π + π − π 0 χ bJ candidates, we require that there be exactly four tracks, of which two are positively identified as pions and the other two as leptons. At least three photons are required in the event and a π 0 -list is created with invariant mass of the photon pairs satisfying M (γγ) ∈ [0.12, 0.15] GeV/c 2 , which covers nearly ±3σ around the π 0 peak. To improve momentum resolution, photon energy resolution, and reduce background, a five-constraint (5C) kinematic fit is performed, where the invariant mass of the π 0 candidate is constrained to be the nominal mass and the energy and momentum of the final state system are constrained to those of the initial e + e − system. The momenta after the 5C kinematic fit are kept for further analysis. The χ 2 5C /ndf is required to be less than 20 for both Υ(1S) → µ + µ − and e + e − modes. Here ndf = 5 is the number of degrees of freedom. If there are multiple π 0 combinations in an event, the one with the smallest χ 2 5C /ndf is kept. The invariant mass of ℓ + ℓ − is required to be in the region [9.42, 9.60] GeV/c 2 .
The χ bJ candidates are reconstructed with the selected Υ(1S) and the photon not in π 0 combination. The invariant mass of π + π − π 0 (M (π + π − π 0 )) versus that of γΥ(1S) Fig. 1 for the sum of the data samples in the Υ(6S) energy region, which is defined as E c.m. > 10.96 GeV. Clusters of events for the production of χ bJ can be seen both when M (π + π − π 0 ) is in ω mass region ([0.75, 0.81] GeV/c 2 ) and at higher masses (> 0.81 GeV/c 2 ). For events having M (π + π − π 0 ) in the ω mass region, the χ b2 signal is dominant while for signal events with higher π + π − π 0 masses, the χ b1 signal is dominant. Background mainly comes from fake π 0 's from photon candidate combinatoric.
An unbinned two-dimensional (2D) extended maximum likelihood fit to M (π + π − π 0 ) and M (γΥ(1S)) is applied to extract the numbers of ωχ bJ and π + π − π 0 χ bJ events. In the fit, the shapes of ωχ bJ and π + π − π 0 χ bJ obtained from MC simulation are used to describe the signals, and a 2D function f (x, y) = ax + by (x = M (γΥ(1S)) and y = M (π + π − π 0 )), is used to fit the background. Here the π + π − π 0 χ bJ MC sample is generated assuming π + π − π 0 follow a phase space (PHSP) distribution, and this process is denoted (π + π − π 0 ) non-ω χ bJ . The statistical significance for π + π − π 0 χ b1 and π + π − π 0 χ b2 , including both ω and non-ω contributions, are 6.4σ and 3.5σ, respectively. The significances are calculated based on the changed likelihood and number of degrees of freedom with or without signals throughout the paper. The signal yields for ωχ b1 and (π + π − π 0 ) non-ω χ b2 are consistent with zero in the fitting results, and thus are set to zero and the fit is repeated. The projections of the fit results for events in χ bJ signal region (M (γΥ(1S)) ∈ [9.87, 9.93] GeV/c 2 ), in ω signal region , and higher than ω mass region are also shown in Fig. 1. The signal yields (statistical significances) for (π + π − π 0 ) non-ω χ b1 and ωχ b2 are 19.6 ± 5.3 (6.1σ) and 7.8 ± 3.2 (4.0σ), respectively. We also use other forms as alternative background descriptions, such as the product of two independent second-order polynomial functions, the 2D function with additional correlative component cxy or constant component. Changes in the fit results are negligible.
In order to study the line shape of π + π − π 0 χ b1 and π + π − π 0 χ b2 events, we extract the signal yields N sig at each energy scan point. Because of the limited statistics for most energy points, we do not perform a 2D fit as for the summed sample, nor do we separate π + π − π 0 into ω and non-ω, nor γΥ(1S) into χ b1 and χ b2 . The number of χ bJ signal events is computed using the formula: N sig = N obs − N side , where N obs is the number of events in χ bJ signal region and N side is that in the sideband region. Here the signal region is defined as M (γΥ(1S)) ∈ [9.852, 9.952] GeV/c 2 , while the sideband region is [9.77, 9.82] and [9.98, 10.03] GeV/c 2 .
The energy-dependent cross sections for e + e − → π + π − π 0 χ bJ are listed in Table I and plotted in Fig. 2. Assuming the π + π − π 0 χ bJ signal comes from Υ(5S) and Υ(6S) decay, a maximum likelihood fit of the cross sections is performed. The likelihood for the three Υ(5S) peak data samples is calculated assuming the Gaussian distribution, while for the other samples, the likelihood is calculated assuming Poisson distribution scaled to the cross section values because of the limited statistics. The fit function is a coherent sum of two BW amplitudes, i.e., Υ(5S) and Υ(6S), and the masses and widths of the Υ(5S) and Υ(6S) are fixed to their world average values [21] while the corresponding branching fractions are left free. The fit results are shown in Fig. 2. Two solutions are found that differ in phase, but the resulting Γ ee · B are nearly the same. The product branching fractions are B(Υ(5S) → e + e − ) · B(Υ(5S) → π + π − π 0 χ bJ ) = (15.3 ± 3.7) × 10 −9 , B(Υ(6S) → e + e − ) · B(Υ(6S) → π + π − π 0 χ bJ ) = (18.3 ± 9.0) × 10 −9 , where the errors are statistical. We also try to introduce a continuum component into the fit, but its significance is only 1.4σ.
There are several sources of systematic error in the cross section measurements. Most of the uncertainties are similar FIG. 2. Fitting to the cross sections of e + e − → π + π − π 0 χ bJ . Red boxes with error bars are the cross sections of e + e − → π + π − π 0 χ bJ and solid blue curve is fitting curve.
to the previous work [4]. The uncertainty from tracking efficiency is 1.0% per pion and kaon track and 0.35% per lepton. The uncertainty from PID efficiency is 1.3% per pion and 1.6% per lepton. The uncertainty in the calibration of the photon energy resolution is less than 1.1%. The uncertainty in the selection of π 0 candidates is estimated by comparing control samples of η → π 0 π 0 π 0 and η → π + π − π 0 in data, and amounts to 2.2%. The uncertainty due to the 5C kinematic fit is 4.2% and a 3.0% uncertainty is assigned to the trigger simulation. All above systematic uncertainties are quoted from Ref. [4]. The uncertainty from luminosity is 1.5% [9]. Comparing the reconstruction efficiency with the ISR process in EVTGEN and with it removed in the generator, but still corrected for with the ISR correction factor, yields an uncertainty of 1.0% due to the reconstruction efficiency of ISR events. The corresponding uncertainty from the branching fraction is 8.2% quoted from the PDG [21]. The total systematic uncertainty, 11.9%, is obtained by adding all the above results in quadrature.
The systematic uncertainties in the lineshape fit mainly come from the parametrization of the BW function, PHSP factor, resonance parameters, and the fitting function used to describe the line shape. The first is estimated by replacing the constant width with an energy dependent width Γ tot = Γ 0 tot · Φ( √ s)/Φ(M ). The second source is estimated by replacing the approximate two-body PHSP factor of π + π − π 0 χ bJ with the two-body PHSP factor of ωχ bJ . The third source is estimated by varying the resonance parameters Υ(5S) and Υ(6S) within ±1σ. The final systematic uncertainty is estimated by adding a coherent amplitude (A con · e iφ ′ ) to the fit function. The changes of the resonance parameters and branching fractions are taken as the systematic uncertainty. The details are listed in Table II.   TABLE II. Summary of the absolute systematic uncertainties in branching fractions (in 10 −3 ) from fitting to the cross sections, where B(5, 6S) represent B(Υ(5, 6)S) → π + π − π 0 χ bJ . , we obtain B(Υ(5S) → π + π − π 0 χ bJ ) = (2.5 ± 0.6 ± 2.0 ± 0.7) × 10 −3 , B(Υ(6S) → π + π − π 0 χ bJ ) = (8.7 ± 4.3 ± 6.1 +4.5 −2.5 ) × 10 −3 , where the first errors are statistical, the second are systematic errors combined from the cross sections measurement and line shape fit, and the third result from the branching fractions of Υ(5S) and Υ(6S) → e + e − .
To reconstruct e + e − → φχ bJ , we require at least two kaons in one event. There is no requirement on the number of photons, but a list of photon candidates is created for the events with photons satisfying |M (γγ 2 ) − m π 0 | > 13 MeV/c 2 , where γ 2 is the second photon in the event with E γ2 > 0.1 GeV, and m π 0 is the nominal mass of π 0 . The data are divided into two categories, one where events with one of the photons in the above list satisfies M (γK + K − ) recoil ∈ [9.42, 9.50] GeV/c 2 , i.e., in the Υ(1S) mass region, used to tag χ bJ → γΥ(1S) events, the other including all other events, used to tag χ bJ →non-γΥ(1S) events.
We use the figure of merit, S/ √ S + B, to optimize the φ signal window requirement, where S is the reconstructed number of signal events obtained from MC simulation in the signal region, [9.88, 9.93] GeV/c 2 , normalized to the theoretical branching fraction of Υ(6S) → φχ bJ [16], and B is the number of background events in the signal region in the Υ(5S) generic MC sample with the c.m. energy shifted to 11.022 GeV. For the first category, we require M (K + K − ) within m φ ± 7.5 MeV/c 2 , and for category two, we require M (K + K − ) within m φ ± 7.0 MeV/c 2 , where m φ is the nominal mass of φ [21]. The φ mass sideband region is defined as M (K + K − ) ∈ [1.000, 1.005] or [1.035, 1.040] GeV/c 2 . There is no evidence for peaking background using the φ mass sideband events, nor in the generic MC sample mentioned above.
After applying all the selection criteria, the recoil mass spectra of φ from both data categories are shown in Fig. 3 for combined data in energy region √ s = 10.96-11.05 GeV. We perform a simultaneous unbinned maximum likelihood fit to the mass spectra with the signals described using MC simulated shapes and a background shape obtained by summing up the expected background shapes at each individual energy point normalized to the luminosity. The expected background shape is obtained from Υ(5S) onresonance data, where, in calculating the K + K − recoil mass, the c.m. energy is changed to that of each individual data point. The ratios of the numbers of χ bJ in the two categories are fixed according to the expected branching fractions [21] and efficiencies. The fit results, which yield χ 2 /ndf = 104.2/55 = 1.9, are shown in Fig. 3. According to the fit, (1.5 ± 0.5) × 10 3 χ b1 and (2.4 ± 0.5) × 10 3 χ b2 events are produced. The statistical significances are found to be 3.3σ and 4.8σ for χ b1 and χ b2 , respectively.
When we vary the background shape by multiplying the nominal background shape with a first, second, or third order polynomial, the smallest significances of the χ b1 and χ b2 signals are found to be 2.6σ and 2.1σ, respectively. The most conservative upper limits on the numbers of produced signal events in all the above tests are reported. After considering the systematic uncertainty which we discuss later, the upper limits for the produced numbers of φχ b1 and φχ b2 signal events are determined to be 2.0 × 10 3 and 3.1 × 10 3 at 90% confidence level (C.L.), respectively, corresponding to upper limits of the Born cross sections of e + e − → φχ b1 and φχ b2 as 0.6 and 1.0 pb, respectively, averaged over the Υ(6S) region, specifically √ s = 10.96-11.05 GeV.
The sources of systematic uncertainties in the φχ bJ cross section measurement are similar to the π + π − π 0 χ bJ modes, including the tracking efficiency, PID, photon detection, luminosity, trigger simulation, ISR correction, φ mass window, and corresponding branching fraction. Most of these have been discussed in the π + π − π 0 χ bJ analysis. The uncertainties from the φ mass window requirement is found to be negligible by studying the consistency of the K + K − invariant mass between data and MC simulation. The uncertainty from the branching fraction of φ → K + K − is 1.0% [21]. The total systematic uncertainty for the cross section measurement is thus, combining all uncertainties in quadrature, 5.5% for both modes.
The processes e + e − → φχ bJ are also searched for in data within √ s = 10.96-11.05 GeV, with no significant signals being observed. We report upper limits on the Born cross sections of e + e − → φχ b1 and φχ b2 are 0.6 and 1.0 pb at 90% C.L., respectively. Compared with the total cross section of e + e − → Υ(6S), these upper limits correspond to Υ(6S) decay branching fractions of order 10 −3 , well above the theoretical predictions of order 10 −6 [16].