Observation of a new structure near 10.75 GeV in the energy dependence of the e+e−→ ϒ (nS)π+π− (n = 1, 2, 3) cross sections

We report a new measurement of the e+e−→ ϒ(nS)π+π− (n = 1, 2, 3) cross sections at energies from 10.52 to 11.02 GeV using data collected with the Belle detector at the KEKB asymmetric-energy e+e− collider. We observe a new structure in the energy dependence of the cross sections; if described by a Breit-Wigner function its mass and width are found to be M=10752.7±5.9−1.1+0.7MeV/c2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ M=\left(10752.7\pm {5.9}_{-1.1}^{+0.7}\right)\mathrm{MeV}/{c}^2 $$\end{document} and Γ=35.5−11.3−3.3+17.6+3.9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Gamma =\left({35.5}_{-11.3\kern0.5em -3.3}^{+17.6+3.9}\right) $$\end{document}MeV, where the first error is statistical and the second is systematic. The global significance of the new structure including systematic uncertainty is 5.2 standard deviations. We also find evidence for the e+e−→ ϒ (1S)π+π− process at the energy 10.52 GeV, which is below the BB¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{B} $$\end{document} threshold.


Introduction
The observed vector bottomonium states above the B B threshold, Υ(4S), Υ(10860), and Υ(11020), have properties that are unexpected for pure b b bound states [1].Their transitions to lower bottomonia with the emission of light hadrons have much higher rates compared to expectations for ordinary bottomonium, and some of these transitions strongly violate Heavy Quark Spin Symmetry.A possible explanation of these unusual properties is a contribution of hadron loops or, equivalently, the presence of a B ( * ) (s)
In the region of the Υ(4S, 5S, 6S) states, quark models also predict the Υ(3D, 4D) levels [5,6].The electron widths of the D-wave states arise from mixing with the S-wave states, which is expected to be quite small for bottomonia below the B B threshold [7].However, it can be significantly enhanced for the states above open-flavor thresholds because of B-meson loops [8].Thus, the "dressed" Υ(3D, 4D) states might be formed abundantly in e + e − annihilations.
The unexpected properties of the Υ(10860) and Υ(11020) could also be due to the presence of other exotic states, e.g., compact tetraquarks [9] or hadrobottomonia [10].To understand the nature of the already known Υ states above the B B threshold and to search for additional states expected in this energy region, it is of interest to study the energy dependence of the e + e − → (b b...) cross sections, where (b b...) denotes exclusive final states containing the b and b quarks, both with open and hidden flavor.

B( * )
s [14] in the energy region from 10.63 to 11.02 GeV.The shapes of the Υ(nS)π + π − and h b (nP )π + π − cross sections show prominent Υ(10860) and Υ(11020) signals with no significant non-resonant contributions.Production of the h b (nP )π + π − is found to proceed entirely via the exotic charged bottomonium-like states Z b (10610) and Z b (10650) [15].The B * s B * s cross section shows a prominent Υ(10860) signal, while the B s Bs and B s B * s cross sections are relatively small and do not show any significant structures.The uncertainties in the χ bJ (1P ) π + π − π 0 cross section at various energies are too large to draw conclusions about its shape.No evidence is found for new structures in any of these cross sections except possibly near 10.77 GeV in the Υ(nS)π + π − final states.It is of interest to study more channels and to improve the accuracy of the previously measured cross sections.
In this paper, we report an updated measurement of the σ[e + e − → Υ(nS)π + π − ] (n = 1, 2, 3) energy dependence.We improve the accuracy by reconstructing Υ(nS) → e + e − in addition to µ + µ − , and by extracting signal yields via fits to the π + π − recoil mass distributions, instead of counting events with inverse efficiency weights in the signal and sideband regions, as was done previously [11].We also use the initial-state-radiation (ISR) process in the high statistics Υ(10860) on-resonance data to obtain additional information about the cross section shapes.As a result of these improvements, we observe a new structure in the energy dependence of the e + e − → Υ(nS)π + π − (n = 1, 2, 3) cross sections with a mass near 10.75 GeV.We also find evidence for the e + e − → Υ(1S)π + π − process below the B B threshold at 10.52 GeV.This implies that the e + e − → Υ(nS)π + π − cross sections have non-resonant contributions.
The paper is organized as follows.In Section 2 we present the data sets used in the analysis, and briefly describe the Belle detector and Monte-Carlo (MC) simulation.In Section 3 we list selection requirements for signal events.Section 4 is devoted to the center-of-mass (c.m.) energy calibration using the e + e − → µ + µ − process.In Section 5 we present an analytical calculation of the e + e − → Υ(nS)π + π − signal shapes in the M recoil (π + π − ) distributions and calibration of momentum resolution using the Υ(2S, 3S) → Υ(1S, 2S)π + π − transitions in the Υ(2S) and Υ(3S) data.Fits to the M recoil (π + π − ) distributions are described in Section 6, while the results of the measurement of the Born cross sections and c.m. energy calibration using combined e + e − → µ + µ − and e + e − → Υ(nS)π + π − processes are presented in Section 7. The fit to the cross section energy dependence, determination of the parameters of the intermediate states, and of the significance of the new structure are presented in Section 9.In Section 10 we give conclusions and mention possible interpretations of the new structure.

Data sets, Belle detector and simulation
The analysis is based on data collected by the Belle detector at the KEKB asymmetricenergy e + e − collider [16,17].We use energy scan data with approximately 1 fb −1 per point collected in the energy range from 10.63 GeV to 11.02 GeV.These data were collected during two running periods, one with six energy points in 2007 and the other with sixteen energy points in 2010.We also use the Υ(10860) on-resonance data sample, with a total luminosity of 121 fb −1 .These data were collected in five running periods with slightly different c.m. energies between 10.864 GeV and 10.868 GeV.Finally, we use 60 fb −1 of the continuum data sample collected at 10.52 GeV.Thus, in total, there are 28 energy points where we calibrate c.m. energies and measure cross sections.Energies and luminosities of various data samples are presented in Table 1.
The Belle detector is a large-solid-angle magnetic spectrometer that consists of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), a barrel-like arrangement of time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter (ECL) comprised of CsI(Tl) crystals located inside a superconducting solenoid that provides a 1.5 T magnetic field.An iron flux return located outside the coil is instrumented to detect K 0 L mesons and to identify muons (KLM).A more detailed description of the detector can be found in Ref. [16].
The integrated luminosity is measured with barrel Bhabha events.The systematic error in the luminosity measurement is about 1.4% and is dominated by the theoretical uncertainty in the Bhabha cross section; the statistical error is usually small compared to the systematic error.

Event selection
We select events of the type e + e − → Υ(nS)π + π − (n = 1, 2, 3) with Υ(nS) → µ + µ − or e + e − .A preselected sample containing µ + µ − or e + e − pairs with invariant masses greater than 8 GeV/c 2 is used.Muons are identified by their range and transverse scattering in the KLM [25].Electrons are identified by the presence of an ECL cluster matching a track in position and energy and having a transverse energy profile consistent with an electromagnetic shower; the ionization loss measurement in the CDC and the response of the ACC are also used [26].We require the presence of two additional charged tracks that are identified as pions.Identification is based on the information from the TOF and ACC, combined with the ionization loss measurement in the CDC.We also apply an electron veto.The total efficiency of the identification requirements is at the level of 99% per pion.All tracks are required to originate from the interaction point (IP) region; this requirement helps eliminate poorly reconstructed tracks.Multiple candidates occur in about 0.5% of the events.We select the best candidate based on the smallest distance to IP in the plane transverse to the beam direction.
To suppress background from converted photons, we require cos θ π + π − < 0.95, where θ π + π − is the pion opening angle in the laboratory frame.In the e + e − π + π − final state we apply additional requirements, M recoil (e + e − ) > 350 MeV/c 2 and cos θ e − < 0.82, where θ e − is the angle between the e − momentum in the c.m. frame and the electron beam.The e + e − recoil mass is defined as where E c.m. is the c.m. energy, E e + e − and p e + e − are the e + e − energy and momentum measured in the c.m. frame.Figure 1 (left) shows the M (µ + µ − ) vs. M recoil (π + π − ) distribution for a high-statistics data sample (No. 9 in Table 1).The clusters along the diagonal are due to e + e − →  Υ(nS)π + π − (n = 1, 2, 3).Below the diagonal there are events in which some final-state particles are not reconstructed.The fully reconstructed events are selected with a requirement: where ℓ = µ or e.Since M recoil (π 2) corresponds to an energy balance requirement.The M recoil (π + π − ) distributions are presented in Fig. 1 (right) for both µ + µ − π + π − and e + e − π + π − events.
4 Calibration of the center-of-mass energies with e + e − → µ + µ − In this section we describe the c.m. energy calibration with e + e − → µ + µ − .Subsequently, the E c.m. calibration is improved using the e + e − → Υ(nS)π + π − (n = 1, 2, 3) processes, as described in Section 6.For the Υ(10860) on-resonance data we use only e + e − → Υ(nS)π + π − .We select e + e − → µ + µ − events with the same requirements as described in Section 3 except requiring the presence of the π + π − pair.We fit the M (µ + µ − ) − E c.m. distributions for all data samples to a Gaussian in the range ±70 MeV/c 2 from the peak position, which corresponds to about ±1 σ.Here and throughout the paper, we use a binned maximumlikelihood fit unless stated otherwise.An example of a fitted M (µ + µ − ) distribution in the scan data sample is shown in Fig. 2. The statistical uncertainty in the peak position is There is a difference between the M (µ + µ − ) peak position and E c.m. due to radiative effects.This difference is determined to be −8.3 ± 0.6 MeV using the e + e − → Υ(nS)π + π − processes in the Υ(10860) on-resonance data as described in Section 6.The error assigned is based on the scatter of the measurements in different on-resonance data samples; it corresponds to the uncertainty due to long-term stability of the detector performance.The dependence of the shift on E c.m. is determined from MC simulation.The shift increases in absolute value by 0.3 MeV over the energy scan range 10.52 − 11.02 GeV.
To estimate the systematic uncertainty in E c.m. , we vary the fit interval from ±50 MeV/c 2 to ±80 MeV/c 2 and take the largest deviation as an error; it is in the range 0.1 − 0.7 MeV.
We investigate the stability of the K 0 S → π + π − invariant mass over the data-taking period and assign a 1 MeV uncertainty due to a possible variation of the Belle solenoid magnetic field during energy scan.The overall uncertainty is determined as a quadrature sum of all contributions; the values are in the range 1.2 − 1.4 MeV.
5 Signal density functions in the M recoil (π + π − ) distributions The c.m. energy calibration using the e + e − → Υ(nS)π + π − (n = 1, 2, 3) processes and the measurement of the corresponding yields is done by fitting the π + π − recoil mass distributions.In this Section we describe the determination of signal density functions in the M recoil (π + π − ) distributions, which is a key tool of this analysis.
The density functions for the Υ(nS)π + π − signals in the M recoil (π + π − ) distributions are calculated as a sequence of convolutions that take into account the π + π − momentum resolution, ISR, and the beam energy spread.
The momentum resolution function is a sum of a symmetric part and a tail on the high M recoil (π + π − ) side due to FSR, decays in flight, and secondary interactions.The symmetric part is described by a sum of five Gaussians.Such a parameterization has enough flexibility to describe non-Gaussian tails and is fast to compute, which is crucial for performing convolutions.The sum of the first three Gaussians contains about 96% of the signal, its standard deviation is between 2 MeV/c 2 and 4 MeV/c 2 for various c.m. energies and channels.The FSR contribution is modeled with a photon energy threshold of 0.1 MeV.
To calibrate the momentum resolution and to verify the simulation of its tails we use the high-statistics Υ(2S) → Υ(1S)π + π − signal in the 25 fb −1 data sample collected by Belle at the Υ(2S) peak.We use the same selection requirements as described in Section 3, with the upper boundary in Eq. (3.2) released to 250 MeV.The M recoil (π + π − ) − m(Υ(1S)) distribution is presented in Fig. 3. Due to the small total width of the Υ(2S) state, the contributions from the beam energy spread and ISR are negligibly small.Therefore we describe the signal by the momentum resolution function; its floated parameters are the normalization, the overall shift in M recoil (π + π − ), and a scale factor f for the width of the symmetric component (we multiply the σ parameters of all Gaussians by f ).The non-peaking background is parameterized with the threshold function where x = M recoil (π + π − ); A, x 0 , p, and c are parameters that are floated.We also consider the peaking backgrounds due to the Υ(2S) → Υ(1S)η transitions with η → π + π − π 0 and π + π − γ.The shapes of these peaking backgrounds are determined from MC simulations and their normalizations are fixed relative to the signal.Fit results are presented in Fig. 3, which shows that the fit function describes the data well.The value of f is found to be We have used Υ(3S) → Υ(1S, 2S)π + π − transitions in the 3 fb −1 data sample collected at the Υ(3S), and find no indication of a dependence of f on π + π − energy.The shifts in For all transitions, we find that the shifts are consistent with zero; thus, no corrections to the recoil mass scale are applied.
For the ISR probability, we use a calculation to second order in α by Kuraev and Fadin [28].We multiply it by the relative change of the cross section with ∆E c.m. , the shift in c.m. energy due to the emission of a photon, and by the relative change in the reconstruction efficiency.The efficiency slowly decreases with ∆E c.m. due to the π + π − becoming softer.To convert from ∆E c.m. to ∆M recoil (π + π − ) we use the approximate relation ∆E c.m. = −∆M recoil (π + π − ).This approximation does not produce any visible effects in our analysis.
The contribution of the c.m. energy spread is modeled by a Gaussian multiplied by the cross section energy dependence.The σ parameter of the Gaussian is found from a fit to the e + e − → Υ(nS)π + π − (n = 1, 2, 3) signals in the combined on-resonance data sample; this fit is described below.The measured value is at the Υ(10860) on-resonance energy of E c.m. = 10.866GeV.To find the spread at other c.m. energies, we assume that it is proportional to E c.m. .We note that if the cross section changes rapidly with E c.m. , then, because of the energy spread, the average energy of the produced Υ(nS)π + π − combination is different from the average energy of the colliding beams.This can result in a shift of the visible peak position that is as large as a few MeV/c 2 .
After the momentum resolution, ISR and energy spread convolutions are performed, we multiply the resulting signal density function by the ∆M recoil (π + π − ) efficiency dependence of the energy balance requirement, which is a step function smeared by the M (µ + µ − ) resolution, which is typically σ = 60 MeV/c 2 .The smearing is described with the error function.
The integrals of the momentum resolution and energy spread functions are normalized to unity; therefore the integral of the signal density function corresponds to that of the ISR function and the measured signal yield already includes the ISR correction, 1 + δ ISR .Such a normalization of the signal function was used in Ref. [12].
The calculation of the signal density function is performed inside the fit function that allows to float various parameters that influence the signal, such as the peak position, the c.m. energy spread, etc.The energy dependence of the cross section is required for the calculation of the signal density function; therefore the analysis is performed iteratively: we compute the signal density functions using results of the fit to the cross section energy dependence from the previous iteration.
6 Fits to M recoil (π + π − ) distributions To determine the e + e − → Υ(nS)π + π − (n = 1, 2, 3) signal yields and to calibrate the c.m. energies, we fit the M recoil (π + π − ) distributions in different data samples.The µ + µ − π + π − and e + e − π + π − candidates are fitted simultaneously, while the fits in different data samples are independent.The fit function is the sum of the signal components, plus non-peaking and peaking backgrounds.
Signal density functions are calculated as described in Section 5.The yields in the µ + µ − π + π − final state are floated, while the ratio of the yields in the e + e − π + π − and µ + µ − π + π − is fixed from MC simulation for all data samples except for the Υ(10860) onresonance data.This latter sample is used to tune the MC simulation.The masses of the Υ(1S), Υ(2S) and Υ(3S) that enter the calculation of the signal density functions are floated with the constraints of their world average values [27].We also introduce a common shift for all peaks, that represents a possible c.m. energy miscalibration.This common shift is floated without constraints for the on-resonance data and with the constraints from the e + e − → µ + µ − calibration for the scan and continuum data.
The µ + µ − π + π − final-state background originates from the QED production of four tracks, e.g.two-photon e + e − annihilations in which one virtual photon produces a highmass µ + µ − pair, and the other produces a π + π − , µ + µ − or e + e − pair, with the latter two misidentified as π + π − .There is a small peaking component in this type of processes due to the high-mass lepton pair being produced via the intermediate Υ(nS) (n = 1, 2, 3).In case of the e + e − π + π − final state there are also photon exchange diagrams and, thus, the background is higher.
The non-peaking component is described by an empirical function where A, x 0 and p are parameters, P 3 (x) is a 3rd-order Chebyshev polynomial with a constant term set to unity.For high-statistics data samples all background parameters are floated.For low-statistics samples, only the normalization A is floated, while the other parameters are fixed using the results of the fit to the combined on-resonance data sample; the threshold x 0 is recalculated for each of these data samples based on the energy of each sample.The peaking background is determined from MC simulation [23] separately for each data sample.Its influence on the measured cross sections is found to be small.In case of the lowest energy, 10.52 GeV, reflections from the Υ(3S) → Υ(1S)π + π − and Υ(3S) → Υ(2S)π + π − transitions leak into the signal band in the M (µ + µ − ) vs. M recoil (π + π − ) plane due to the finite resolution of the µ + µ − mass.The reflections are described in the fit by Gaussians with all parameters floated.Examples of the fits for the Υ(10860) on-resonance and scan data samples are shown in Figs. 1 and 4. The   1. Top and bottom panels correspond to µ + µ − π + π − and e + e − π + π − final states, respectively.Results of the fit (red solid histogram) and their background components (black dotted histogram) are also shown.
M recoil (π + π − ) spectra are fitted with 1 MeV/c 2 bins, but we present them with larger bin sizes for clarity.
7 Results for the center-of-mass energies and the Υ(nS)π + π − cross sections The calibrated c.m. energies are obtained from the fit results for the common shifts of the Υ(nS)π + π − signals.We do not find significant deviations from the constraints of the e + e − → µ + µ − calibration; the e + e − → µ + µ − and e + e − → Υ(nS)π + π − energy calibrations agree well.The calibrated E c.m. values are presented in Table 1.
Based on the measured signal yields we calculate Born cross sections: No.  where N is the signal yield (since the signal density function has been appropriately normalized, N includes the ISR correction, as discussed in Section 5), |1 − Π| 2 is a vacuum polarization correction taken from Ref. [29], ε is the reconstruction efficiency (Fig. 5), L is the integrated luminosity of each data sample (Table 1), and B is the branching fraction for Υ(nS) → µ + µ − [27].Measured e + e − → Υ(nS)π + π − cross sections are presented in Table 1.
To study systematic uncertainties in E c.m. and the cross sections, we vary the E c.m. spread, the scale factor f , and the e + e − π + π − to µ + µ − π + π − efficiency ratio by ±1 σ.We also increase the polynomial order in the background parameterization for high-luminosity data samples, while for low luminosity we release the coefficient of the linear term of the Chebyshev polynomial.Uncertainties associated with the cross section energy dependence are estimated using MC pseudoexperiments that are generated according to the fit results described in Section 9.The signal density functions are computed based on the fit results of each pseudoexperiment.
The systematic uncertainties are estimated based on the deviations of the measured quantities from their nominal values under the above variations of the analysis.In order to avoid overestimation of the relative uncertainties among the Υ(10860) on-resonance points, we separate uncertainties into correlated and uncorrelated parts.For each variation, we first find the average deviation over all the on-resonance points; this is used to estimate the correlated uncertainty.Then, we subtract the average deviation from deviations at each energy point.These relative deviations are used to estimate the uncorrelated uncertainties.In the case of the cross section energy dependence, we take the root mean squares of the deviations (for both correlated and uncorrelated parts).In all other cases, we take the maximal deviation as the uncertainty.
Uncertainties in efficiency associated with the DP model are treated as uncorrelated.
The total uncorrelated systematic uncertainty is the sum in quadrature of all the contributions; dominant among them is the cross section energy dependence.The systematic uncertainty is small compared to the statistical uncertainty, as can be seen in Fig. 6.We add the two contributions in quadrature to find the total uncorrelated uncertainty.The c.m. energies and Born cross sections with total uncorrelated uncertainties are presented in Table 1.The correlated systematic errors include the uncertainties in the efficiency due to possible differences between data and MC in track reconstruction (0.35% per high and 1% per low momentum tracks, which are muons and pions, respectively) and muon identification, and the uncertainties in the luminosity and the branching fractions for Υ(nS) → µ + µ − decays.The summary of the correlated errors is presented in Table 2.
The energy-dependent cross sections (Fig. 6) show clear Υ(10860) and Υ(11020) peaks that were seen in previous publications [11,30].Due to the improved precision an additional structure around E c.m. = 10.75 GeV, hinted at by the measurements in Ref. [11], is visible.
The obtained expectations for the Υ(2S) → Υ(1S)π + π − and Υ(3S) → Υ(2S)π + π − tails are quite large because matrix elements of these three-body decays are dominated by the terms proportional to M 2 (π + π − ) (see, e.g., [22]).Corresponding partial widths are calculated as integrals of the matrix elements over phase space; they increase rapidly with the c.m. energy as high values of M (π + π − ) become kinematically allowed.The dependence of the cross sections for various tails on E c.m. is shown in Fig. 7.The cross sections start to increase above a certain energy.Uncertainties of the prediction of the tail contributions are discussed in the next Section.To verify the hypothesis of the Υ(2S, 3S) tails, we note that the M (π + π − ) distribution for the QED background is quite different from the expectations for the tails.The background is studied using the µ + µ − π + π − events in the Υ(1S) sideband region 8.7 < M recoil (π + π − ) < 9.4 GeV/c 2 .The background M (π + π − ) distribution (Fig. 8) is enhanced at low values and in the region of the ρ meson; the MC simulation describes the shape and the normalization of the sideband data quite well.In contrast, the Υ(2S, 3S) tail Events / 20 MeV/c events are expected to concentrate near the upper kinematic boundary.To suppress the QED background, we apply a requirement M (π + π − ) > 0.85 GeV/c 2 , which is optimized using the figure of merit S 3/2+ √ B [32], where S and B are the number of signal and background events in the simulation, respectively.The π + π − recoil mass distribution with this additional requirement for the µ + µ − π + π − final state is shown in Fig. 9.A clear signal for Events / 10 MeV/c the Υ(1S)π + π − process is evident; its significance is estimated using Wilks' theorem [33] to be 3.6 σ.The signal is stable against variations of the fit interval and the order of the polynomial used for the background parameterization.We conclude that the significance including systematic uncertainties is larger than 3.5 σ.The expected M (π + π − ) requirement efficiency is 72% and is consistent with the reduction of the signal yield in the data.
Based on the signal yield measured with the M (π + π − ) requirement, we determine the e + e − → Υ(nS)π + π − cross section at E c.m. = 10.52 GeV to be 0.042 +0.017 −0.015 pb, where the uncertainties are statistical only.This result is model-dependent because of the unknown M (π + π − ) distribution.In the fit to the cross section energy dependence, we use the value measured without the M (π + π − ) requirement.

Fit to the energy dependence of the cross sections
Since the numbers of signal events in some scan data samples are small, the errors in the measured cross sections might be non-Gaussian.To address this issue, we scan −2 ln L of the M recoil (π + π − ) fits in the signal yields.The −2 ln L dependence on the signal yield is then converted into the −2 ln L dependence on the cross section and is parameterized using an empirical function where p 1 , p 2 and p 3 are fit parameters and P 6 is a 6th-order Chebyshev polynomial, with the zeroth-order parameter set to unity and all other parameters floated.To account for the systematic uncertainty, we convolve the distributions in L with Gaussian functions that represent uncorrelated systematic uncertainties.We find that there is good agreement between the scan results and the asymmetric Gaussian errors in the ±1 σ region; however, there are noticeable discrepancies for larger deviations.Therefore Gaussian errors should give a good approximation for a default fit that describes data well; while the −2 ln L scan results are important for significance calculation, in which case some points deviate from the alternative fit function by several standard deviations.In the fits to the cross section energy dependence, we use the −2 ln L scan results for all scan data points outside the Υ(10860) peak (points 2 to 8 and 20 to 28 in Table 1).
The ISR tails of the Υ(nS)π + π − signals in the Υ(10860) on-resonance data are sensitive to the cross section shapes in the region of the new structure at 10.75 GeV.Therefore we include the M recoil (π + π − ) distribution for the µ + µ − π + π − final state, where background is lower, into the fit of the cross section energy dependence.Thus we perform a simultaneous fit to the Υ(nS)π + π − (n = 1, 2, 3) cross sections and the M recoil (π + π − ) distribution.
To parameterize the cross sections energy dependence, we consider contributions from the Υ(10860) and Υ(11020) resonances, the new structure, and the Υ(2S, 3S) tails.Each contribution is represented by a Breit-Wigner amplitude: where M and Γ are the mass and total width, Γ 0 ee is the "bare" electron width related to the physical width by Γ 0 ee = Γ ee |1 − Π(s)| 2 , and B f and Γ f (s) are branching fraction and energy-dependent partial width of the decay to the Υ(nS)π + π − final state.The values of Γ f (s) at various energies are computed numerically by integrating the three-body decay matrix element over the DPs.
The new structure might have either resonant or non-resonant origin.In some cases, non-resonant effects produce peaks and have phase motion similar to resonances, as discussed, e.g., in Ref. [34].Thus, we use the Breit-Wigner parameterization for the new structure amplitude.
In the default model the DP distributions of the Υ(10860), Υ(11020), and new structure decays are assumed to be the same, thus their amplitudes are added coherently.For the Υ(2S) → Υ(1S)π + π − and Υ(3S) → Υ(2S)π + π − tails we use the three-body matrix elements measured by CLEO [22].The interference terms between the tails and the rest of the amplitude are multiplied by "decoherence factors" that are calculated as overlap integrals of the DP matrix elements [11,12], and can take values that range from 0 (incoherence) to 1 (full coherence).The values of the decoherence factors are found to be typically 0.7-0.8.
Complex phases of the Υ(11020), the new structure, and the Υ(2S, 3S) amplitudes relative to the Υ(10860) amplitude are all floated in the fit individually for the three channels.We also float the M , Γ and Γ 0 ee × B f parameters of the Breit-Wigner amplitudes except for M and Γ of the Υ(2S) and Υ(3S).The fit function to the Υ(nS)π + π − cross section at each c.m. energy contains a convolution with a Gaussian to account for the energy spread.The fit function to the M recoil (π + π − ) distribution is the same as described in Section 6.Energy dependence of the e + e − → Υ(nS)π + π − cross sections (n = 1, 2, 3 from top to bottom).The points with error bars are data with combined statistical and uncorrelated systematic uncertainties; the curves are the results of the simultaneous fit to these distributions and the M recoil (π + π − ) distribution in the on-resonance data (Fig. 11).Shown are the results of the default fit (red solid curve) and the fit with the new structure excluded in all channels (blue dotted curve).The hatched histograms and the points with error bars at their maxima show the Υ(4S) line shapes and the cross sections measured at the Υ(4S) peak in Ref. [31].These points are not used in the fit to the cross section energy dependence.Insets show a zoom of the low energy region, while right side panels show a zoom of the Υ(10860) on-resonance region.
The fit results are presented in Figs. 10 and 11.For illustration, we show in Fig. 10 the cross sections measured at the Υ(4S) peak and the expected Υ(4S) line shape; these measurements are not used in the fit.The fit results for masses and widths of the Υ(10860), Υ(11020), and new structure are presented in Table 3.The results are slightly shifted from the previous Belle measurements based on the same data sample [11] because we updated The M recoil (π + π − ) distribution for the µ + µ − π + π − events in the Υ(10860) onresonance data; the right panel is a zoom of the left panel.Histograms show the results of the simultaneous fit to this distribution and the energy dependence of the e + e − → Υ(nS)π + π − (n = 1, 2, 3) cross sections (Fig. 10).The red solid histogram corresponds to the default fit, while the blue dashed histogram in the right panel corresponds to the fit with the new structure excluded.The black dotted histogram shows the background component of the default fit, the hatched onepeaking background from η decays.Table 3. Measured masses and widths of the Υ(10860), Υ(11020) and the new structure.The first uncertainty is statistical, the second is systematic.the cross section values at each energy, we fit now Born cross sections instead of visible ones, the new structure and the non-resonant contributions, i.e. the Υ(2S, 3S) tails, are included in the fit model, and the energy dependence of the Γ f is taken into account.
A sum of several Breit-Wigner amplitudes produces multiple solutions that have the same values of −2 ln L, the same masses and widths, but different normalizations and relative phases.To search for multiple solutions, we create points with fitted values of the cross sections and small uncertainties.We then fit these points, one Υ(nS)π + π − channel at a time, repeating the fit many times with randomly generated initial values of the fit parameters.We find four or eight solutions in various channels, as expected for the sum of three or four Breit-Wigner amplitudes, respectively.The ±1 σ intervals of various solutions typically overlap, therefore we present only ranges in Γ ee ×B f from the lowest to the highest solution (Table 4).
To estimate the significance of the new structure in a single Υ(nS)π + π − channel, we repeat the fit with the new structure excluded in that channel.The significance is found based on the change in the −2 ln L of the fit using Wilks' theorem.The values are 2.7 and 5.4 standard deviations (σ) in the channels Υ(1S)π + π − and Υ(2S)π + π − , respectively.Using large sample of MC pseudoexperiments, we find that Wilks' theorem provides slightly conservative estimation in our case.This is related to the fact that the number of experimental points in the region of the new structure is quite limited.We also estimate the significance for all three Υ(nS)π + π − channels combined by performing the fit with the new structure excluded in all channels simultaneously.The −2 ln L difference between the default fit and the null hypothesis fit is 65.8, with the cross section scan results contributing 51.7 and recoil mass distribution contributing 14.1.This −2 ln L difference corresponds to a local significance of 7.0 σ, estimated using Wilks' theorem.The global significance is estimated using the method described in Ref. [35].We find that the p-value of the fluctuation increases by a factor 4.5 due to the "look-elsewhere effect", and the resulting global significance is 6.8 σ.
To estimate systematic uncertainties in the Υ(10860), Υ(11020), and new structure parameters, we vary the fit procedure.As an alternative parameterization of the new structure amplitude we use the threshold function, Eq. (5.1).For its parameters we find x 0 ≈ 10.73 GeV/c 2 , p ≈ 0.17 and c ≈ −10 (GeV/c 2 ) −1 .The quality of this fit is comparable to the default fit, with −2 ln L worse by 3.4.Such a parameterization represents a thresholdlike contribution without variation of a complex phase with energy.It has no clear physical meaning and we use it to conservatively estimate systematic uncertainty in the Υ(10860) and Υ(11020) parameters.
We study the influence of thresholds on the line shape of the Υ(10860) resonance.We consider the thresholds in the region of the new structure that show strong coupling to the Υ(10860).The factors 2  3 and 1 3 correspond approximately to the ratio of the B B * π and B * B * π cross sections [36].We set the weights x and y to various values of 0, 0.2, 0.4, and 0.6 with the restriction x + y ≤ 0.8.We find that by introducing the effect of the Z b π thresholds we always increase the significance of the new structure.The reason for this is that the thresholds suppress the lower mass tails of the Υ(10860) resonance, and this leads to worse description of data under the null hypothesis.
In the default model the tails of the Υ(2S) → Υ(1S)π + π − and Υ(3S) → Υ(2S)π + π − transitions are calculated under an assumption that matrix elements of the three-body decays grow with M (π + π − ) as M 2 (π + π − ).This growth can be damped due to the presence of some form factor related to details of the hadronization process [37].We consider a damping factor of the decay amplitude in the form with the parameter Λ taking the values 1, 2 and 4 GeV.The Γ f (s) dependence for the Υ(2S, 3S) tails is recalculated for each value of Λ.In all cases we find that the significance of the new structure increases.Thus, the default model corresponds to the fastest possible growth of the Υ(2S, 3S) tails with energy and this gives a conservative estimate of the significance.We note that another reason why the growth of the non-resonant contribution with the energy could be damped is the crossing of various open-flavor thresholds.The tail of the Υ(3S) resonance can also contribute to the Υ(3S)π + π − final state.Thus, we include it into the fit using the same parametrization as for the Υ(2S)π + π − channel.To calculate its energy dependence and decoherence factors we assume that corresponding distribution over DP is uniform.
The significance of the new structure in the Υ(3S)π + π − channel in the default model is 3.9 σ.The reason why it is so high is evident from Fig. 10 (lowest panel): in the absence of the new structure the description of the Υ(3S)π + π − data in the region of the Υ(10860) peak becomes very poor.Thus, the new structure increases flexibility of the fit function in the region of the Υ(10860) peak due to interference.Obviously, a similar effect of interference can be achieved with a non-resonant contribution instead of the new structure.Indeed, the significance of the new structure in the Υ(3S)π + π − channel drops to 0.6 σ once the non-resonant contribution is added.The production threshold in the Υ(3S)π + π − channel is above the continuum energy of 10.52 GeV, therefore the information about the non-resonant contribution is limited, and available data do not allow to study the interplay between the new structure and the non-resonant contribution.
We consider the maximum deviation of the result as a systematic uncertainty associated with a given source.The different contributions and the total uncertainty obtained as their quadrature sum are presented in Table 5.We find that the significance of the new structure in the Υ(2S)π + π − channel and global significance combined over all channels remain above 5.1 σ and 5.2 σ, respectively, for all the variations that were introduced to study systematic uncertainties.The lowest significances are reached when we include the non-resonant contribution into the Υ(3S)π + π − channel.
To visualize the ISR contribution to the measurement of the cross section energy dependence we estimate the e + e − → Υ(nS)π + π − cross sections based on the ISR tails.For this, we divide the background-subtracted M recoil (π + π − ) distributions in the ISR tail regions by the ISR luminosity function and efficiency, and apply other corrections from Eq. (7.1).Technically, to obtain the correction function we compute the Υ(nS)π + π − signal shape under an assumption that the cross section is constant with the c.m. energy and divide the obtained shape by the Υ(10860) on-resonance cross section.The results are presented in Fig. 12.The cross sections estimated via ISR are compatible with the scan results and provide support for the new structure.However, they are not intended to be used in the fit, because the uncertainties are statistical only and do not include contributions due to the background subtraction.Also, the ISR luminosity changes rapidly in the studied energy region, which complicates accurate description of the resolution effects.It is only in the fit to the ISR tails in the M recoil (π + π − ) distribution that all these issues are addressed rigorously.
The new structure could have a resonant origin and correspond to a signal for the not yet observed Υ(3D) state provided S − D mixing is enhanced [8], or an exotic state, e.g. a compact tetraquark [9] or hadrobottomonium [10].It could also be a non-resonant effect due to some complicated rescattering.Information on the cross section energy dependence for more channels, with both hidden and open b flavor, is needed to clarify the nature of the new structure.nology of Taiwan; and the United States Department of Energy and the National Science Foundation.

Figure 3 .
Figure 3.The M recoil (π + π − ) − m(Υ(1S)) distribution for the data collected at the Υ(2S) resonance (points with error bars) and fit results (black solid histogram).Also shown are combinatorial background (cross-hatched histogram), background from η decays (hatched histogram) and various contributions to the signal function: FSR with E γ > 0.1 MeV (blue dashed), decays in flight (magenta dotted) and secondary interactions (red dash-dotted).Both panels show the same distribution, but differ in bin sizes, linear or logarithmic scales of the vertical axes and ranges of the horizontal axes.

Figure 4 .
Figure 4. M recoil (π + π − ) distributions in the data samples No. 4 (left) and 5 (right); the numbering is given in Table1.Top and bottom panels correspond to µ + µ − π + π − and e + e − π + π − final states, respectively.Results of the fit (red solid histogram) and their background components (black dotted histogram) are also shown.
Figure 10.Energy dependence of the e + e − → Υ(nS)π + π − cross sections (n = 1, 2, 3 from top to bottom).The points with error bars are data with combined statistical and uncorrelated systematic uncertainties; the curves are the results of the simultaneous fit to these distributions and the M recoil (π + π − ) distribution in the on-resonance data (Fig.11).Shown are the results of the default fit (red solid curve) and the fit with the new structure excluded in all channels (blue dotted curve).The hatched histograms and the points with error bars at their maxima show the Υ(4S) line shapes and the cross sections measured at the Υ(4S) peak in Ref.[31].These points are not used in the fit to the cross section energy dependence.Insets show a zoom of the low energy region, while right side panels show a zoom of the Υ(10860) on-resonance region.
Figure 11.The M recoil (π + π − ) distribution for the µ + µ − π + π − events in the Υ(10860) onresonance data; the right panel is a zoom of the left panel.Histograms show the results of the simultaneous fit to this distribution and the energy dependence of the e + e − → Υ(nS)π + π − (n = 1, 2, 3) cross sections (Fig.10).The red solid histogram corresponds to the default fit, while the blue dashed histogram in the right panel corresponds to the fit with the new structure excluded.The black dotted histogram shows the background component of the default fit, the hatched onepeaking background from η decays.
These are B B * π, B * B * π and B * s B * s at 10.75 GeV, 10.79 GeV and 10.83 GeV, respectively.The B * s B * s channel is the only among the B a prominent Υ(10860) signal [14].Production of the B B * π and B * B * π channels proceeds entirely via intermediate Z b (10610)π and Z b (10650)π states, respectively [36], while the Z b π channels show prominent Υ(10860) signals in the h b (nP )π + π − final state [12].Thus we multiply the constant width Γ by an energy-dependent factor 1 − x − y + x (p 1 /p p i (i = 1, 2, 3) are momenta in the B * s B * s , Z b (10610)π and Z b (10650)π pairs, respectively, and superscript (0) denotes momentum calculated for the nominal Υ(10860) mass.

Table 4 .
The ranges of the Γ ee × B values from the multiple solutions (in eV).