Measurement of the energy dependence of the $e^+e^-\to B\bar{B}$, $B\bar{B}^*$ and $B^*\bar{B}^*$ exclusive cross sections

We report the first measurement of the exclusive cross sections $e^+e^-\to B\bar{B}$, $e^+e^-\to B\bar{B}^*$, and $e^+e^-\to B^*\bar{B}^*$ in the energy range from 10.63 GeV to 11.02 GeV. The $B$ mesons are fully reconstructed in a large number of hadronic final states and the three channels are identified using a beam-constrained-mass variable. The shapes of the exclusive cross sections show oscillatory behavior with several maxima and minima. The results are obtained using data collected by the Belle experiment at the KEKB asymmetric-energy $e^+e^-$ collider.

The total e + e − → b b cross section at various energies above the B B threshold is shown in Fig. 1 (left) [1].It exhibits peaks of Υ(4S), Υ(10860), and Υ(11020), possibly a dip in the region of Υ(10775), and also dips at the B B * and B * B * thresholds.The exclusive two-body cross sections e + e − → B B, e + e − → B B * , and e + e − → B * B * , that saturate the total cross section below the Υ(10860) peak and give a dominant contribution also at higher energy, are expected to show much more pronounced behaviour, as shown in Fig. 1 (right) [2].The expected oscillatory behavior of the exclusive cross sections might be due to the nodes of the Υ(4S), Υ(10860), and Υ(11020) wave functions [2].These cross sections provide important information about the interactions in this energy region and, in particular, about the structure of the Υ(4S), Υ(10860), and Υ(11020) states.This topic is of special interest since the above states show anomalies, primarily in the pattern of transitions to lower bottomonium states, that are currently not well understood (for a review see, e.g., Ref. [3]).
Here we report the first measurement of the energy dependence of the e + e − → B B, e + e − → B B * , and e + e − → B * B * exclusive cross sections. 1Our approach is to perform a full reconstruction of one B meson in hadronic channels, and then to identify the B B, B B * , and B * B * signals using the M bc distribution, M bc = (E cm /2) 2 − p 2 B , where E cm is the center-of-mass (c.m.) energy and p B is the B-candidate momentum measured in the c.m. frame.The M bc distribution for B B events peaks at the nominal B-meson mass, m B , while the distributions for B B * and B * B * events peak approximately at m B − ∆m B * 2 and m B − ∆m B * , respectively, where ∆m B * is the mass difference of the B * and B mesons [4].If the B meson originates from a B * → Bγ decay, there is an additional broadening of the signal due to the photon recoil momentum.
To reconstruct B mesons in a large number of hadronic final states we apply the Full Event Interpretation (FEI) package of the Belle II software that was developed primarily for tagging B mesons in the Υ(4S) → B B decays [5].This package uses multivariate analysis for the event selection and provides high flexibility in choosing the B decay channels and the input variables for the classifier.
Before going into details, we describe how the paper is organized and give an overview of the analysis.In Section 2 we briefly describe the Belle detector, the data samples and the simulation.Selection of events is described in Section 3.For the FEI classifier, we choose input variables that are not correlated with the B candidate momentum, which helps to avoid distortion of the background in the M bc distribution and to keep efficiency approximately independent of E cm .We do not include the energy of the B candidate, E B , into the FEI training and use sidebands in the M bc versus ∆E plane to study background, where ∆E = E B − E cm /2.We find that there is a peaking background in the M bc distribution which is primarily due to misreconstructed soft photons.We first apply FEI to the Υ(4S) data sample (Section 4).We construct the M bc fit function in which the signal shape is calculated based on the E cm spread, the cross section energy dependence, and the momentum resolution.We also calibrate simulation of the peaking background and determine the total B meson yield which is later used to determine the efficiency.We proceed with the analysis of the Υ(10860) data sample (Section 5) with the aim to verify the fit procedure and measure the signal yield.We also study the distribution of B mesons in the polar angle (Appendix B).The fits to the data samples at various energies are presented in Section 6.We measure the efficiency at the Υ(4S) and Υ(10860) energies (Section 7).To determine the total numbers of B mesons in the Υ(10860) data sample, we use five clean B + and B 0 decay channels reconstructed without FEI.Determination of the cross sections, parameterization of the cross section energy dependence, and estimation of the systematic uncertainties are presented in Section 8.The results are discussed in Section 9.As a byproduct, we measure f s , the fraction of the B ( * ) s B( * ) s events at Υ(10860) (Section 10).We conclude in Section 11.

Belle detector and data sets
The analysis is based on data collected by the Belle detector at the KEKB asymmetricenergy e + e − collider [6].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-offlight 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).Two different inner detector configurations were used.For the first sample of 156 fb −1 , a 2.0 cm-radius beam pipe, and a 3-layer silicon vertex detector were used; for the latter sample of 833 fb −1 , a 1.5 cm-radius beam pipe, a 4-layer silicon vertex detector (SVD2), and a small-cell inner drift chamber were used.This analysis is based only on data collected with the SVD2 configuration.Detailed description of the detector can be found in Ref. [7,8].
We use energy scan data with approximately 1 fb −1 per point, six points collected in 2007 and 16 points collected in 2010.We use also the Υ(5S) on-resonance data with a total luminosity of 121 fb −1 collected at five points with energies from 10.864 GeV to 10.868 GeV.The E cm calibration of these data is reported in Ref. [9].We combine the data samples with similar energies so that finally we obtain 23 energy points.The energies and integrated luminosities of these 23 data samples are presented in Table 5 below.We use also the SVD2 part of the Υ(4S) data sample with the primary goal to calibrate the reconstruction efficiency; its integrated luminosity is 571 fb −1 .

B( * )
s , e + e − → q q (q = u, d, s, c) events are generated using EvtGen [10].The detector response is simulated using GEANT [11].The Monte-Carlo (MC) simulation includes run-dependent variations in the detector performance and background conditions.

Event selection
The event selection is performed primarily by FEI.Our strategy is not to include the ∆E variable in the training of the FEI classifier so that the ∆E sidebands (or, more precisely, the ∆E sidebands with the ∆E variable defined as ∆E ≡ ∆E + M bc − m B ) are available.We use the sidebands to study the smooth background component.We perform the FEI training and then apply further channel-dependent selection criteria on the FEI output variable and the ∆E variable, as described below.
We reconstruct B + and B 0 in the decay channels , where D denotes D0 and D − .We do not use B-decay channels with π 0 as their energy resolution is rather poor.The D 0 , D + , and D + s mesons are reconstructed in the final states with K ± , K 0 S , π ± , up to one π 0 , and multiplicity up to five.The complete list of the B-and D-meson decay channels is given in Appendix A. We reconstruct D * in all possible decay modes: Dπ and Dγ.J/ψ are reconstructed in both µ + µ − and e + e − channels.To improve momentum resolution, we apply a mass-constrained fit to π 0 , J/ψ, and D * ; mass-vertex-constrained fit to D and D + s ; and vertex fit to K 0 S and B. In FEI, the Fast Boosted Decision Tree algorithm [12] is used to discriminate signal and background events.First, the final-state particles are classified, and then all the candidates for unstable particles are obtained as combinations.Thus, FEI training is performed in stages, and the results of the previous stage are used in the training of the current one.The training is performed using MC simulation at the Υ(5S) energy.The classifier output is the probability that a given candidate is the signal.
As training variables of charged pions, kaons, and leptons, we use particle identification information, momentum, and transverse momentum.Prior to training we select the candidates that originate from the interaction-point (IP) region: we require dr < 0.5 cm and dz < 3 cm, where dr and dz are cylindrical coordinates of the point of the closest approach of the track to the beam.For kaons we apply an additional requirement of L(K)/(L(K) + L(π)) > 0.1, where L(K) and L(π) are likelihoods of the kaon and pion hypotheses formed from the measurements in the ACC and TOF systems, as well as energy loss measurement in CDC.The efficiency of this requirement is 98% and the probability to misidentify a pion as a kaon is about 20%.
Photon candidates are clusters in ECL with the energies above 30 MeV and without matching tracks.As training variables we use the number of crystals in the cluster, the ratio of energy deposit in a 3 × 3 matrix of crystals to that in a 5 × 5 matrix, cluster energy, and polar angle.
For the π 0 → γγ candidates we use mass, momentum, and decay angle, which is defined as the angle between the γ momentum measured in the π 0 rest frame and the π 0 boost direction from the laboratory frame.For the K 0 S → π + π − candidates we use mass and a set of parameters describing the displaced vertex of K 0 S .These are the distance of closest approach between the two daughter pions, the impact parameters of the daughter pions, the distance between the IP and the K 0 S vertex, and the angle between the K 0 S momentum and the direction from the IP to the K 0 S vertex; the latter three variables are measured in the plane perpendicular to the beam direction.
The training variables for J/ψ, D, and D * candidates are the signal probability of each daughter (thus the number of variables varies with the channel) and the mass.In case of D, we use also the χ 2 of the mass-vertex fit; if the D decay is a three-body decay, we use invariant masses of pairs of its daughters to take into account signals of ρ, K * , and φ.
For the B meson candidates we use the signal probability of each daughter and χ 2 of the B vertex fit.If there is a D meson in the decay, we include the distance between the B and D vertices, the significance of this distance, and the cosine of the angle between the D momentum and the direction from B to D vertices.In the case when there are several pions or pions and kaons in the decay, we include invariant masses of the combinations in which production of ρ (→ ππ), a 1 (→ 3π) or K * (→ Kπ) is expected.
At the last stage when the training for the B candidates is performed, we include also variables that help to suppress continuum production of light and charm quarks, e + e − → q q.These are the event-shape variable R 2 (the ratio of the second to zeroth Fox-Wolfram moments [13]) and the angle between the thrust axis of the B candidate and the rest of the event.We also include two flags indicating the presence of a muon and an electron, respectively, in the rest of the event.We consider lepton candidates in the c.m. momentum windows 1.0 < p µ < 2.6 GeV/c and 0.8 < p e < 2.6 GeV/c where the contribution of leptons from the semileptonic B decays is enhanced.We require that the leptons are well identified with a likelihood ratio above 0.9 [7].The efficiencies of this requirement are 71% and 76% for muons and electrons, respectively; the probabilities to misidentify hadrons as leptons are 1% and 0.1%, respectively.
Although the training is performed individually for each decay channel of every unstable particle in the decay chain, the signal probability is defined in a universal way so that various channels can be compared.Thus, the signal probability from the classifier is used to rank multiple candidates.At the intermediate stage of the reconstruction, we retain up to 10 best D 0 , D + , and D + s candidates.At the final stage of the B meson reconstruction, we retain only one best candidate selected from all B + and B 0 candidates in the event.
The entire MC statistics corresponds to six times the statistics of real data.We use half of the MC statistics to train the classifier and the other half to determine the efficiency.The efficiency in the part that was used for training is higher by a factor 1.025 ± 0.006, which is an indication of a small overtraining.
The ∆E versus M bc distribution for the Υ(5S) on-resonance data is presented in Fig. 2.
Here we apply a requirement on the B meson signal probability from the FEI classifier of P B > 0.1.In the ∆E projection, the signal events are concentrated near zero.
We optimize the requirements on the P B and ∆E variables individually for each B decay channel.We evaluate the contribution of each channel to an overall figure of merit (FoM) defined as N S / √ N S + N B , where N S and N B are the numbers of signal and background events, respectively.This optimization is performed based on the Υ(5S) on-resonance  The MC simulation at Υ(5S) shows that the background is dominated by the cc production.The B s mesons do not produce a significant contribution in the region of the B ( * ) B( * ) signals.
In the optimization of the channel-dependent selection requirements, N S and N B are determined using the ∆E signal and sideband regions with an additional requirement of 5.27 < M bc < 5.34 GeV/c 2 .The optimization for each channel is performed iteratively by scanning FoM in turn in P B and |∆E |.Since we use the Υ(5S) on-resonance data for the optimization, the measured B meson yield could be biased.We note, however, that the statistics in each B decay channel is high and we use a relatively large step in the P B and |∆E | scanning; therefore, statistical fluctuations in the FoM dependence on P B and |∆E | are small.To further study this issue, we divide the Υ(5S) data sample into two approximately equal parts.We then use part 1 for optimization and reconstruct part 2 using the resulting selection requirements.Similarly, we reconstruct part 1 using selection requirements optimized with part 2. In this way we completely avoid any bias in the yields due to statistical fluctuations in the value of FoM during the optimization.We find a 1.3% smaller ratio of the yields at Υ(5S) and Υ(4S) compared to the default procedure.We consider this value as a symmetric systematic uncertainty due to the optimization procedure.

Analysis of the Υ(4S) data sample
Here our goal is to describe the M bc distribution in the Υ(4S) data in terms of the E cm spread and the B B cross section shape.This experience is essential for the analysis of the Υ(5S) and energy scan data.
The information about the B B cross section shape in the Υ(4S) region is rather limited.Current values for the Υ(4S) mass and width are dominated by the BaBar measurement in 2004 [14].There is also a more precise scan performed by BaBar in 2008 [1], however it was not fitted in the original paper.We attempt to fit the 2008 scan results using the Υ(4S) parameterization from Ref. [14], which is based on the Quark Pair Creation model [15].However, the fit function overestimates the measured values at high-mass side of Υ(4S).Since we need the cross section shape to calculate the signal shape in the M bc distribution, we adopt the following strategy.We perform a simultaneous fit to the M bc distributions and the scan results of BaBar [1].As a suitable model is missing, the cross section is described by a high-order Chebyshev polynomial.

M bc fit function
The signal component of the M bc fit is calculated numerically as a sequence of convolutions.It takes into account the beam-energy spread, the energy dependence of the production cross section, the ISR, and the momentum resolution.
The energy spread of the colliding beams is described by a single Gaussian with the mean E cm0 , which is a nominal c.m. energy.The distribution in E cm is multiplied by the energy dependence of the cross section.We convolve the obtained function with the Kuraev-Fadin radiation kernel [16] to account for the ISR.We then change the argument of the function from E cm to the B meson momentum p 0 .At this step we account for the ISR recoil momentum of the B mesons.
We convolve the distribution in p 0 with the momentum resolution functions.We use three resolution functions to describe the candidates of three types: (1) MC truth-matched candidates in the ∆E signal window, ( 2) not MC truth-matched candidates in the ∆E signal window, and (3) candidates in the ∆E sidebands.The candidates of type 1 correspond to the signal, while candidates of type 2 and 3 correspond to the peaking background.Each of the three resolution functions, f i , is a sum of two Gaussians with special factors that account for the fact that p can not be negative: The special factors are found by considering the momentum resolution function in three dimensions and analytically integrating out all variables other than p.The parameters of the resolution function are determined from the MC simulation (Table 1).The functions The smooth background is described by a threshold function E cm /2 − x multiplied by a third-order Chebyshev polynomial.The shape of the smooth background is the same in the ∆E signal region and sidebands while the normalizations are allowed to float independently (the smooth component in the sidebands is multiplied by a floated parameter r s.b. ).
When fitting the data, we introduce a shift and a width-correction factor for each component of the momentum resolution function, s i and φ i .Using the ∆E distributions we find that for the signal component the shift is negligibly small, while the width-correction factor is φ 1 = 1.187±0.012.We float the shift and the width-correction factor of the peaking background in the sidebands, s 3 and φ 3 , and find that the values are consistent with zero and one, respectively (Table 2 below).Therefore, the shift and the width-correction factor of the peaking background in the signal region, s 2 and φ 2 , that are poorly constrained by the M bc fit, are fixed at s 2 = 0 and φ 2 = 1.
The peaking background components in the ∆E signal region and sidebands are multiplied by a common normalization factor n which is floated in the fit.Thus, the floated parameters related to the signal in the M bc distribution are the signal yield N (integral of the signal component; it does not include the integral of the peaking background), the E cm spread σ Ecm , the peaking background normalization factor n, the shift s 3 , and the width-correction factor φ 3 .

Cross section fit function
We describe the energy dependence of the dressed cross section2 by an 11th order Chebyshev polynomial (in case of the 10th order, χ 2 of the fit is higher by about 20, while in case of the 12th order it is almost unchanged).We require that the cross section is zero at the B B threshold and is never negative.We then apply the ISR correction and convolve with the BaBar energy spread of 4.83 MeV [14].We use 9 BaBar points located between the B B and B B * thresholds.The accuracy of the c.m. energy calibration at BaBar is 1.5 MeV.
We introduce a common shift for all nine points, ∆E BaBar , and float it in the fit with the uncertainty constraint adding a term ∆E 2 BaBar /(1.5 MeV) 2 to χ 2 .

Results of the simultaneous fit
The results of the simultaneous fit to the M bc distributions in the ∆E signal and sideband regions and to the cross section energy dependence are presented in Figs. 4, 5 and Table 2.
We perform the fit at several values of E cm0 , the average c.m. energy of the Belle Υ(4S)   the visible cross section, ∆E cm .We find that ∆E cm is equal to zero at E cm0 = 10.5787GeV and use this E cm0 value in the default fit.If E cm0 is floated, we find E cm0 = 10.5791 ± 0.0003 GeV, which is 1.6 σ away from the constraint.The p-value of the default fit is 1.8%.
To determine the systematic uncertainty, we consider a variation in E cm0 of ±0.5 MeV that corresponds to a variation in ∆E cm of ±1.5 MeV.The decrease of the visible cross section at ∆E cm = ±1.5 MeV is about 1%.This variation produces a negligible change in the yield N , but is a dominant uncertainty for σ Ecm and ∆E BaBar .
We increase the order of the Chebyshev polynomial that describes the cross section shape (11th order to 12th) and which is used in the smooth background component (3rd order to 4th).In both cases we find a negligible change in all fit results.
The normalization factor n of the peaking background is found to be 1.16 (Table 2).This value is determined primarily by the ∆E sidebands.To estimate the systematic uncertainty related to the peaking background in the ∆E signal region, we introduce a separate normalization factor for this component, n 2 , and repeat the fit fixing n 2 at 1.08 and 1.24.This source is dominant for the yield while for σ Ecm and ∆E BaBar the changes are small.The total systematic uncertainty for the yield N , σ Ecm , and ∆E BaBar is presented in Table 2.The value of the yield is used to determine the efficiency at the Υ(4S) energy (Section 7).
The value of the E cm spread at Υ(4S), (5.36 ± 0.19) MeV, and the measurements at other energies are shown in Fig. 6.The value at Υ(5S), (5.36 ± 0.13) MeV, is measured in Ref. [9] using the e + e − → Υ(nS)π + π − (n = 1, 2, 3) processes.To find the values at Υ(1S, 2S, 3S), we use the visible cross sections that are determined based on the event yields and luminosities in Ref. [8].The procedure is the following.We determine the energy dependence of the dressed cross section near the resonance using corresponding parameters: mass, width, and electron width [4].We then apply the ISR correction by performing a convolution with the Kuraev-Fadin radiation kernel [16] and account for the E cm energy spread by performing another convolution with a Gaussian.The integral of the cross section depends on the total and electron widths, while the shape of the cross section is determined by the energy spread, since the Υ(1S, 2S, 3S) resonances are very narrow.Thus, the maximum of the visible cross section is sensitive to σ Ecm .For each Υ(nS) (n = 1, 2, 3) we fit a single point: the measured visible cross section at the resonance maximum.The value of the fit function is the maximum of the calculated cross section.The parameters of the resonance are floated in the fit within the uncertainties of their worldaverage values [4], thus their contribution to the uncertainty in the spread is accounted for.We find that the dominant contribution is the uncertainty in the electron width.The spread values at Υ(1S), Υ(2S), and Υ(3S) are (4.439± 0.157) MeV, (5.19 ± 0.59) MeV, and (5.95 ± 0.80) MeV, respectively.If the bunch current of the positron beam, I e + , is above 0.5 mA, there is a microwave instability that increases the E cm spread [17].The increase factor, f σ , depends linearly on I e + reaching f σ = 1.20 at I e + = 1.0 mA.The average over data taking period value of f σ is determined based on I e + and integrated luminosity of each run.For the Υ(nS) (n = 1, 2, 3, 4, 5) data samples, we find f σ = 1.070, 1.157, 1.223, 1.203, and 1.179, respectively.The uncertainty in f σ is negligibly small compared to the uncertainty in σ Ecm .The corrected values σ cor Ecm = σ Ecm /f σ at various E cm are shown in Fig. 6.The dependence of σ cor Ecm on E cm is consistent with the proportionality hypothesis; from the fit we find: The f σ values in the scan data samples are in the range 1.11 − 1.20.To determine σ Ecm of each scan data sample we use its f σ factor and the σ cor Ecm value from Eq. (4.2).The shift in E cm of BaBar, (−1.75 ± 0.68) MeV, is 1.2 σ away from zero in terms of the BaBar accuracy of 1.5 MeV.This shift could be used in future phenomenological analyses of the BaBar scan results.
As a consistency check we estimate the ISR correction factor to be (1 + δ ISR ) = 0.626 ± 0.012.The uncertainty here is a systematic one due to variation of E cm0 .This value agrees with the result of Ref. [18] of 0.622 ± 0.018.

Analysis of the Υ(5S) data sample
To fit the Υ(5S) data we include also the e + e − → B B * and e + e − → B * B * signals.The decay B * → Bγ leads to additional smearing of the B momentum that we take into account in the fit function by performing additional convolution; relativistic kinematics is used in this calculation.The distribution in the helicity angle of the B * → Bγ decay, defined as the angle between the B momentum measured in the B * rest frame and the boost direction of the B * from the c.m. frame, is expected to be 1+ cos 2 θ h for e + e − → B B * and 1+a h cos 2 θ h with −1 ≤ a h ≤ 1 for e + e − → B * B * .In the fit, we float the parameter a h .
The energy dependence of the e + e − → B ( * ) B( * ) cross sections, that is needed for the fit function, is taken from the measurements described below; the analysis is performed using an iterative procedure.The parameters of the momentum resolution function are determined from the Υ(5S) simulation.They are found to be close to those at Υ(4S).The factor n for the normalization of the peaking background is taken to be the same as in the Υ(4S) data (Table 2).The E cm spread is fixed to the fitted value (Eq.(4.2)) multiplied by the microwave instability correction factor; the result is σ Ecm = (5.44 ± 0.09) MeV.The smooth background is described by a threshold function E cm /2 − x multiplied by a 6th order Chebyshev polynomial; the order is higher than at Υ(4S) because we use a broader fit interval.The result of the fit to the Υ(5S) data is shown in Fig. 7 and Table 3.We  0.998 ± 0.007 report the total yield N total and the fractions of various channels N B ( * ) B( * ) / N total .One of the fractions is not floated as its value is determined from the constraint that the sum of the three fractions is equal to 1.To find its statistical uncertainty, we repeat the fit with a different choice of the fraction, which is not floated.The yield is defined as the integral of the fit function in the region 5.27 < M bc < 5.35 GeV/c 2 .The fit describes the data well, and its p-value is 87%.The M bc distribution in the ∆E signal region is fitted only up to M bc = 5.375 MeV/c 2 to avoid the region near the M bc threshold where the contribution of the e + e − → B B * π and e + e − → B * B * π processes is expected [19].We extrapolate the fit function beyond the fit interval and indeed find an excess of events near the threshold.The study of the three-body channels B ( * ) B( * ) π will be a subject of separate analysis.The value of N total is used to calculate the efficiency in the Υ(5S) data (Section 7), therefore we estimate its systematic uncertainty (Table 4).The shape of the signal depends upon the energy dependence of the B B, B B * , and B * B * cross sections.We find the corresponding uncertainty as described in Section 8. To estimate the uncertainty due to the E cm spread, we vary its value within the uncertainty of ±0.09 MeV.The normalization factor of the peaking background n is varied between 1.08 and 1.24.In the fit to the Υ(4S) data we find that the shift s 3 and the width-correction factor φ 3 of the peaking background in the ∆E sidebands are consistent with zero and one, respectively (Table 2), while in the Υ(5S) data they deviate from the above values with a combined significance of 2.6 σ (Table 3).To estimate the uncertainty associated with the shape of the peaking background, we repeat the fit fixing s 3 = 0 and φ 3 = 1.To account for the uncertainty in the shape of the smooth background we change the order of the polynomial which is used in the parameterization from 6th to 7th and 8th.We account also for the systematic uncertainty due to the optimization procedure of 1.3% (Section 3).The deviations in the yield under the variations of the analysis are considered as systematic uncertainties due to a given source.The total systematic uncertainty (given both in Tables 3 and 4) is estimated as a sum in quadrature of the individual contributions.
We study the distributions of the B B, B B * , and B * B * in the polar angle of the B meson and find that they agree with the expectations.The details are provided in Appendix B.

Fits at various energies
The fit function is the same as at Υ(5S), except for the normalization of the signal.Here the integral of each signal component is not normalized to unity in the range 5.27 < M bc < 5.35 GeV/c 2 .Instead, it is equal to the integral of the ISR kernel [16] multiplied by the relative change of the cross section with energy.This normalization value is equal to the (1 + δ ISR ) correction factor and thus the measured yields include the ISR correction and can be used directly to determine the dressed cross sections.This approach was used in previous energy scan papers [9,20].
We fix the shift s 3 , the width-correction factor φ 3 and the relative normalization of the smooth background r s.b. to the fit results at Υ(5S) (Table 3).The angular distribution parameter a h is floated within the allowed range −1 ≤ a h ≤ 1.In case of the smooth background only the normalization is floated.The fits at various energies are shown in Appendix C.
Measured yields are used to calculate the dressed cross sections as described in Section 8 after we present the determination of the efficiency in the next Section.The values of a h in the scan data samples are poorly constrained by the fit; their uncertainties are close to the size of the allowed range.

Determination of the efficiency
To determine the efficiency at the Υ(4S) energy we use the measured B meson yield N (Υ(4S)) (Table 2) and the total number of the B B pairs in the Υ(4S) SVD2 data, N B B (Υ(4S)) = (619.6± 9.4) × 10 6 [8].This number is obtained by counting hadronic events at Υ(4S) and subtracting the continuum contribution, which is determined using data collected 60 MeV below Υ(4S).The transitions from Υ(4S) to lower bottomonia have total branching fraction of 0.26% [4] and are neglected.The efficiency is calculated as To determine precisely the ratio of the B meson yields in the Υ(4S) and Υ(5S) data samples, we use five final states with low multiplicity for which the distribution over phase space is well known and which can be reliably simulated: The signal-to-background ratio in these final states is high, and they are reconstructed without application of FEI.Selection requirements are taken to be the same as in Ref. [19].To minimize the sensitivity to the peaking background, we fit the ∆E spectra simultaneously in the Υ(4S) and Υ(5S) data samples applying a requirement 5.27 < M bc < 5.35 GeV/c 2 .The fit is performed separately for each channel.The fit function for the Υ(4S) data is a sum of two Gaussians to describe the signal and a first or second order polynomial (depending on the channel) to describe the background.The fit function for the Υ(5S) data is the same except that the ratio of yields at Υ(5S) and Υ(4S), shift, and width-correction factor are introduced for the signal.The background components in the Υ(4S) and Υ(5S) data samples are floated independently.To estimate the systematic uncertainty, we consider variations of the polynomial order and the fit interval.The systematic uncertainty is calculated as the root-mean-square (RMS) of the deviations.We repeat the same fits for the MC samples and find that the efficiency ratio is consistent with one for all the channels.The efficiency-corrected yield ratios are 0.0393 ± 0.0017, 0.0376 ± 0.0023, 0.0399 ± 0.0021, 0.0365 ± 0.0033, and 0.0391 ± 0.0022 for the channels from one to five, respectively; the uncertainties include the statistical and systematic contributions.The average yield ratio is R 5 chan = 0.03882 ± 0.00097.
In this calculation we implicitly assume that the ratio of the B + B − and B 0 B0 production rates, f B + /f B 0 , is the same at Υ(4S) and Υ(5S).Since Υ(5S) is far from the B B threshold, one can expect that isospin violation at Υ(5S) is small and f B + /f B 0 = 1.At Υ(4S), f B + /f B 0 was measured to be 1.058 ± 0.024 [4]; thus, it is shifted from 1 by 2.4 standard deviations.We repeat the calculation taking into account the isospin nonconservation at Υ(4S), and find that R 5 chan increases by 0.53%.Since the change is very small, the isospin non-conservation is neglected.
The efficiency at the Υ(5S) energy is determined from the total B meson yield N total (Table 3): 3) The ratio of the efficiencies at Υ(5S) and Υ(4S), 1.049 ± 0.032, agrees with the MC expectation of 1.028 ± 0.004.
From MC simulation we find that the dependence of the efficiency on the B meson momentum is consistent with being linear.Thus, for all energies and various B ( * ) B( * ) final states we determine the efficiency ε based on the corresponding average momentum and the values ε Υ(4S) and ε Υ(5S) , assuming linear dependence on the B meson momentum.To find the uncertainty in ε, we separate uncertainties in ε Υ(4S) and ε Υ(5S) into common and uncorrelated parts; we then assume that the uncorrelated part varies linearly with the B meson momentum.

Results for the cross sections
The dressed cross sections are calculated as where the ratio N/(1 + δ ISR ) is directly obtained from the fit.The cross sections are presented in Table 5 and in Fig. 8.The cross sections show a non-trivial behaviour with  cross sections measured by Belle [21] and BaBar [1].The sum is compatible with the total b b cross section up to E cm = 10.82GeV; this value is close to the B * s B * s threshold.The deviation at higher energy is presumably due to the contributions of B s mesons, multibody final states B ( * ) B( * ) π(π), and production of bottomonia with light hadrons.
To calculate the shapes of the signals at various energies and to determine the ISR corrections, we need to parameterize the energy dependence of the cross sections.Since currently there is no suitable phenomenological model for the cross section energy dependence, we fit the cross sections using high-order Chebyshev polynomials.The fit to the total b b visible cross section in the Υ(4S) region (Fig. 5) and the analysis in Ref. [18] (Fig. 9) show that the dressed B B cross section goes to zero at the B B * threshold.Thus, for the  respectively.These orders provide sufficient flexibility to describe the available data.The polynomials are constrained to be positive by adding a penalty term to the χ 2 in case the polynomial becomes negative at any energy.To account for the E cm spread, we convolve the polynomials with the Gaussian.The results of the simultaneous fit are presented in Figs. 8 and 10.
The uncertainty in the shape of the cross section energy dependence contributes to the systematic uncertainty in the cross section measurements at various energies.The uncertainty in the shapes originates from the parameterization and from the limited statistical accuracy in the cross section measurements.In addition to the default set of the polynomial orders of (10,17,12), we consider also sets (11,18,13), (12,19,14), (9,16,11), and (8, 15, 10) that provide a conservative estimation of the possible cross section behaviours.Corresponding fit results are shown in Fig. 11.
To estimate the influence of the statistical accuracy, we use toy MC.We generate pseudoexperiments using the fitted cross sections as central values and statistical uncertainties in data as standard deviations.We fit the energy dependence of the cross sections in each pseudoexperiment and, based on the fit results, determine the M bc signal shapes for all energies.We then refit the data using the new shapes and repeat the measurement of the cross sections.The RMS of the deviations is taken as a systematic uncertainty; the parameterization and the statistical accuracy are considered as separate sources.We find that for the three high-statistics Υ(5S) on-resonance points the deviations are strongly correlated, therefore they are accounted for in the correlated uncertainty of these points; for other points the uncertainties are considered to be uncorrelated.
To study the uncertainty due to the shape of the smooth background in the M bc fits, we multiply the corresponding contribution by the Chebyshev polynomial of the first or second order with floated parameters.The RMS of the deviations of the yields are used to calculate the uncertainties.
We vary the E cm values within their uncertainties, the deviations in the yields are found to be negligible.The cross section shape contributions and the smooth background shape contribution are added in quadrature to obtain the total uncorrelated systematic uncertainties.They are found to be small compared to the statistical uncertainties as shown in Fig. 8.
Various contributions to the correlated systematic uncertainty, estimated for the Υ(5S) high-statistics data, are presented in Table 6.The contribution of the cross-section shape is estimated as discussed above.The contributions of the E cm spread and the peaking background yield and shape are estimated as described in Section 5 for the total B meson yield.We account also for the uncertainty in the efficiency (Eq.(7.3)) and the uncertainty in the integrated luminosity of 1.4%.Total correlated uncertainty is estimated as a sum in quadrature of the individual contributions.
The sources of the correlated systematic uncertainty at energies other than Υ(5S) onresonance are the same as listed in Table 6, except for the cross section shape source which is accounted for in the uncorrelated uncertainty.The contributions of the E cm spread and the shape of the peaking background are assumed to be the same as listed in Table 6.The uncertainty in the efficiency varies with the B meson momentum as described in Section 7.
To transform the multiplicative correlated uncertainty in the cross section into the additive one, we use formula: where the symbol ⊕ denotes addition in quadrature.The measured dressed cross sections at various energies with statistical, uncorrelated systematic, and correlated systematic uncertainties are presented in Table 5.

Discussion
Figure 12 presents a comparison of the measured exclusive cross sections with the predictions of the Unitarized Quark Model (UQM) [2].The data confirm the prediction that the cross sections show oscillatory behaviour.Also, there is a rather good agreement in the positions of the minima, that in the UQM are due to zeros in the Υ(4S, 5S, 6S) wave functions.The UQM fails to describe the absolute values of the cross sections.Contrary to the expectations, the cross sections in the minima are not zero, which suggests that the UQM misses some general non-resonant offset.
In the UQM there are narrow structures in all the channels that correspond to the signals of Υ(5S).Data do not show such structures.Thus, in the final states B B, B B * , and B * B * we find no clear Υ(5S) signal.As follows from Fig. 9, the narrow peak in the Υ(5S) region is present in other b b final states, B ( * ) s
The sum of exclusive B The polarization of the B * B * channel is described by three amplitudes, as discussed in Appendix B. To measure these amplitudes the reconstruction of a photon from the B * → Bγ decay would be needed.Such a measurement requires higher statistics than currently available 1 fb −1 at the scan points.
The separation between the points is rather large in the low-energy region.In particular, the model [2] predicts an additional zero in the B B and B B * cross sections which is in the gap between two scan points.More scan data with smaller energy step sizes and larger integrated luminosity in this region are needed to understand reliably the shape of the exclusive cross sections.These data could be collected by the Belle II experiment.
10 Measurement of visible cross sections and event fractions at Υ(5S) As a byproduct, we measure the visible cross sections e + e − → B B X, e + e − → B B, e + e − → B B * , and e + e − → B * B * at Υ(5S), as well as corresponding fractions of events and the fraction of B ( * ) s

B( * )
s events f s .To find the inclusive e + e − → B B X cross section, we use the same method as for the measurement of the efficiency at Υ(5S) (Section 7).The only difference is that the requirement M bc < 5.35 GeV/c 2 is not applied; thus, we consider not only two-body e + e − → B ( * ) B( * ) but also multi-body e + e − → B ( * ) B( * ) π(π) processes.The ratio of the B meson yields in the Υ(5S) and Υ(4S) SVD2 data samples, averaged over the five lowmultiplicity B-decay channels that were used in Section 7, is where the uncertainty includes statistical and systematic contributions.For the cross section we find: where L is the total integrated luminosity of the three high-statistics points in Table 5.The fraction of B B X events is where σ b b = (340±16) pb is the total b b cross section at Υ(5S) [23].Here and in the following we take into account that in the ratio of cross sections the uncertainty due to the integrated luminosity cancels.The remaining events, f non−B B X = 1 − f B B X = 0.249 ± 0.040, contain B s mesons or bottomonia with light hadrons.
To estimate the fraction of bottomonium events, we consider all final states with bottomonium listed in PDG 2020 [4].These are Υ(nS)π + π − (n = 1, 2, 3), Υ(1S)K + K − , Υ(1D)η, h b (nP )π + π − (n = 1, 2), and χ bJ π + π − π 0 .The sum of corresponding fractions is (3.50 +0.40  −0.42 )%, where we assume that the uncertainties in various fractions are uncorrelated.Using isospin relations, we account also for the final states with neutrals: the fractions for the π + π − transitions are multiplied by 1.5, while the fraction for the K + K − transition is multiplied by 2.0.The resulting sum is (4.92 +0.53  −0.56 )%.Belle reported preliminary results on the Υ(nS)η (n = 1, 2) and Υ(1D)π + π − transitions [24] that show that corresponding fractions are not large.There are also rather strict upper limits on fractions of the η b (nS)ω (n = 1, 2) [25] and h b (nP )η (n = 1, 2) [26] transitions.However, there are still channels for which no experimental information is available.Among them are 4π transitions to Υ(nS), h b (nP ) and Υ(1D), as well as Υ(5S) → Z b π → η b (1S)ρπ.To estimate the total fraction of bottomonium, we assume that all the channels that are not in PDG 2020 contribute no more than already measured channels and thus assign a large positive uncertainty: Finally, we estimate the fraction of the events with B s mesons: where the uncertainty includes statistical and systematic contributions.This result is consistent with the PDG 2020 value of 0.201 ± 0.031 [4] and the value that follows from the measurement of the e + e − → B ( * ) s B( * ) s cross section in Ref [27]: To measure the visible e + e − → B ( * ) B( * ) cross sections, we use the formula: where N B ( * ) B( * ) is the B ( * ) B( * ) signal yield in the interval 5.27 < M bc < 5.35 MeV (Table 3).here N B ( * ) B( * ) / N total are the fractions of various channels given in Table 3.
We study the systematic uncertainty in the values of the expression (N B ( * ) B( * ) / N total )× k B ( * ) B( * ) in the same way as described in Section 5 for N total .The results are presented in Table 7.
The contribution of the ISR events in the Υ(4S) region is obtained by integrating the ISR kernel multiplied by the dressed cross section shown in Fig. 5 in the region between the B B and B B * thresholds.The resulting value is (10.53 ± 0.31) pb, where the uncertainty includes the contributions from the statistical (1.2%) and systematic (2.6%) errors in R b [1], and the systematic uncertainty due to the variation of E cm0 (0.7%).The values of the measured visible cross sections are presented in Table 8.We show also corresponding event fractions.Previous Belle measurement of these values [19] did not take into account the ISR tails of the signals.The results shown in Table 8 supersede those in Ref. [19].
We also measure the ratio σ vis (e + e − → B B * )/σ vis (e + e − → B B X), which is needed for the measurement of f s using lepton-charge correlations in the dilepton events [28].Based on Eqs.(10.2) and (10.8), we find  where the uncertainty includes statistical and systematic contributions.While calculating the ratio R all 5 chan /R 5 chan = 0.2962 ± 0.0176, we take into account the correlation of the corresponding uncertainties.

Conclusions
To conclude, we report the first measurement of the energy dependence of the exclusive cross sections, e + e − → B B, e + e − → B B * , and e + e − → B * B * in the region from 10.63 to 11.02 GeV.The results are presented in Table 5 and Fig. 8.The cross sections show non-trivial behavior with several maxima and minima.They can be used in future phenomenological studies to shed light on b b-quark and B ( * ) B( * ) -meson interactions in this energy region.
As a byproduct, we measure at Υ(5S) the fraction of the events containing the B s mesons, f s = 0.200 +0.040 −0.064 , where the error contains statistical and systematic contributions.We measure also the visible cross sections and corresponding fractions of the events for the B B X, B B, B B * and B * B * final states (Table 8 and Eq.(10.9)).

A B and D decay channels
The B-and D-meson decay channels that are used in this analysis are listed in Tables 9  and 10.
B Angular analysis at Υ(5S) As a consistency check, we measure signal yields in various intervals of the B meson polar angle (Fig. 13).From the MC simulation we find that the variation of the efficiency with cos θ can be neglected.The expected distributions for B B and B B * are sin 2 θ and 1+cos 2 θ, respectively.The data agree with these expectations well, the p-values of the corresponding fits are 69% and 24%, respectively.The B * B * pairs can be produced in three states: L = 1, S = 2; L = 3, S = 2; L = 1, S = 0, where L is the orbital angular momentum and S the total spin of the B * B * pair.The expected total polar angle distribution is 1 + b cos 2 θ, with −1 ≤ b ≤ 1.From the fit we find b = −0.20 ± 0.03, the p-value of the fit is 88%.
We measure the parameter a h in various intervals of cos θ and do not find a significant variation.5 (from left to right and from top to bottom).Legend is the same as in Fig. 7.
-32 -   5 (from left to right and from top to bottom).Legend is the same as in Fig. 7.
-33 -    5 (from left to right and from top to bottom).Legend is the same as in Fig. 7.

Figure 1 .
Figure 1.The R b scan results from BaBar [1] (left), and the expected in Unitarized Quark Model contributions of the B B, B B * , and B * B * channels [2] (right).

Figure 2 .
Figure 2. The ∆E versus M bc distribution for the Υ(5S) data.

Figure 3 .
Figure 3.The M bc distributions in the Υ(5S) (left) and Υ(4S) (right) data.Black solid histogram shows the ∆E signal region, while red points with error bars show the normalized ∆E sidebands.

Figure 4 .Figure 5 .
Figure 4.The M bc distributions in the ∆E signal region (left) and sidebands (right).Points with error bars are data, solid red histogram is the result of the simultaneous fit to these distributions and the cross section energy dependence (Fig. 5), dashed red histogram is the smooth background, dotted black histogram is the peaking background in the ∆E signal region.The lower panels show the residuals.

Figure 6 .
Figure 6.Energy dependence of the E cm spread.Black dots with error bars are the measurements, red open dots are σ Ecm corrected for the microwave instability effect, red line is the fit result.

Figure 7 .
Figure 7.The M bc distributions in the ∆E signal (top) and sideband (bottom) regions.Points with error bars are data, solid histogram is the result of the simultaneous fit, dashed histogram is the smooth background.Black dotted histogram indicates the contribution of the B B channel that includes a peak near threshold due to the ISR production of Υ(4S).Vertical red line at 5.375 GeV/c 2 in the top panel indicates the upper boundary of the fit interval.

Figure 8 .
Figure 8. Measured dressed cross sections at various energies for e + e − → B B (left), e + e − → B B * (right), and e + e − → B * B * (bottom).The outer error bars indicate statistical uncertainties and inner red error bars indicate uncorrelated systematic uncertainties.Solid curves show the result of the simultaneous fit to these distributions and the total b b cross section energy dependence (Fig. 9).Dashed curves show the fit function before the convolution to account for the E cm spread.

Figure 9 .
Figure 9. Energy dependence of the total b b dressed cross section obtained in Ref.[18] from the visible cross sections measured by Belle[21] and BaBar[1] (black dots).Open red circles represent the sum of the exclusive B B, B B * , and B * B * cross sections measured in this work.Right panel is a zoom of the low cross section region.

Figure 10 .
Figure 10.Energy dependence of the total b b dressed cross section from Ref. [18] (blue dots).Solid black curve is the result of the simultaneous fit to this distribution and the exclusive B B, B B * , and B * B * cross section energy dependence (Fig. 8).Vertical line at 10.75 GeV indicates the upper boundary of the fit interval; dashed black curve is an extrapolation of the fit function.Also shown are the individual contributions of B B (blue dashed curve), B B * (red dotted curve), and B * B * (green dash-dotted curve).

Figure 11 .
Figure 11.Measured dressed cross sections at various energies for e + e − → B B (left), e + e − → B B * (right), and e + e − → B * B * (bottom).The outer error bars indicate the statistical uncertainties and the inner red error bars indicate the systematic uncertainties due to the cross section parameterization.Solid curves show the fit results for the default set of polynomial orders.Dotted curves show the fit results for the polynomial orders varied by ±1 and ±2.

Figure 12 .
Figure 12.Measured dressed cross sections at various energies for e + e − → B B (left), e + e − → B B * (right), and e + e − → B * B * (bottom).Points and dashed curves are the same as in Fig. 8. Solid curve shows the predictions of the Unitarized Quark Model [2].
The factor k B ( * ) B( * ) = 1.042, 1.120, and 1.089 for B B, B B * , and B * B * , respectively, accounts for the ISR tail in the M bc > 5.35 MeV region.For the B B channel, k B B includes the ISR contribution down to the B B * threshold, while the region between the B B and B B * thresholds (the Υ(4S) region) is accounted for separately.For the B B * and B * B * channels, k B ( * ) B( * ) include the ISR contributions down to the corresponding thresholds.The efficiency ε B ( * ) B( * ) is equal to ε Υ(5S) × r B ( * ) B( * ) , where the factor r B ( * ) B( * ) = 1.010, 1.004, and 0.998 for B B, B B * , and B * B * , respectively, accounts for the small variation of the efficiency with the momentum of the B meson (Section 7).Using Eq. (7.3), we obtain σ vis (e + e − → B ( * ) B( * ) ) = N B ( * ) B( * ) N total k B ( * ) B( * ) N B B (Υ(4S)) × R 5 chan r B ( * ) B( * ) × L , (10.8)

Figure 13 .
Figure 13.The yield of B * B * (top set of points), B B * (middle) and B B (bottom set of points) as a function of the polar angle of the B meson.The histograms show the fit results.

Figure 14 .
Figure 14.The M bc distributions for the points 1 to 6 in Table5(from left to right and from top to bottom).Legend is the same as in Fig.7.

Figure 15 .
Figure 15.The M bc distributions for the points 7 to 12 in Table5(from left to right and from top to bottom).Legend is the same as in Fig.7.

Figure 16 .
Figure 16.The M bc distributions for the points 13 to 18 in Table5(from left to right and from top to bottom).Legend is the same as in Fig.7.

Figure 17 .
Figure 17.The M bc distributions for the points 19 to 23 in Table5(from left to right and from top to bottom).Legend is the same as in Fig.7.

Table 1 .
Parameters of the momentum resolution functions of Eq. (4.1).

Table 2 .
Results of the simultaneous fit to the Belle M [1]distribution and the BaBar cross section scan results[1].The first error is statistical, the second one (if present) is systematic.N(581.2± 1.1 ± 3.2) × 10 3

Table 3 .
Results of the fit to the M bc distribution at Υ(5S).The errors are statistical.

Table 4 .
Systematic uncertainties in the total signal yield N total at Υ(5S) (in %).

Table 5 .
Energies (in MeV), luminosities (in fb −1 ) for various data samples and the results for the dressed cross sections (in pb).The first error in the cross section is statistical, the second is uncorrelated systematic, and the third is correlated systematic.
B B B B * B * B *

Table 9 .
Decay channels of B + and B 0 used in FEI.

Table 10 .
Decay channels of D 0 , D + and D + s used in FEI.