Hadronic cross section of $e^+e^-$ annihilation at bottomonium energy region

The Born cross section and dressed cross section of $e^+e^-$ to $b\bar{b}$ and the total hadronic cross section in $e^+e^-$ annihilation in the bottomonium energy region are calculated based on the Rb values measured by the BaBar and Belle experiments. The data are used to calculate the vacuum polarization factors in the bottomonium energy region, and to determine the resonant parameters of the vector bottomonium(-like) states, the Y(10750), Upsilon(5S), and Upsilon(6S).


I. INTRODUCTION
The cross section of e + e − annihilation into hadrons is essential information for QED as it is related to the vacuum polarization (VP) of the photon propagator. The measurement of these cross sections is one of the important topics in various e + e − colliders from low energy to high, and the precisions of the measurements have been improved from time to time since the running of the first generation e + e − colliders [1]. The data have been used in many calculations involving the photon propagator, especially in the high precision calculations of the anomalous magnetic moment of the µ, a µ , and the running of the fine structure function, α(s), where s is the center-of-mass (CM) energy squared [2,3].
The cross section of e + e − annihilation into hadrons is often reported in terms of R value, defined as where σ B (e + e − → µ + µ − ) = 4πα 2 (0) 3s , is the Born cross section of e + e − → µ + µ − . The experimental measurements of the R values are compiled in Ref. [1]. There are many data at low energies ( √ s < 2 GeV) with precision at 1% level; while the measurements are sparse and less precise at higher energies, for example, the charmonium (3.7 < √ s < 5.0 GeV) and the bottomonium (10.5 < √ s < 11.2 GeV) energy regions. One of the reasons of less measurements at high energy is the smaller contribution to the VP, and another reason is the fact that fewer experiments were designed in these energy regions. The cross sections of e + e − → bb were measured in much higher precision by BaBar [4] and Belle [5] experiments in the bottomonium energy region, i.e., √ s = 10.5 to 11.2 GeV, than by CUSB [6] and CLEO [7] experiments more than 30 years ago. However, neither BaBar nor Belle (let alone CUSB and CLEO) did radiative correction to the measured cross sections, so the data cannot be used directly for many calculations where the Born cross sections are needed as input.
In this paper, we describe how to get the Born cross section based on the published data from the BaBar and Belle experiments with some reasonable assumptions. We report the Born cross sections from these experiments and discuss the usage of the data samples in the calculation of the VP factors especially in the bottomonium energy region, and the fit to the dressed cross sections to extract the resonant parameters of the vector bottomonium states. We also discuss a possible determination of the VP directly by measuring e + e − → µ + µ − cross sections with high luminosity data at the Belle or Belle II experiment, and a strategy to search for the production of invisible particles in e + e − annihilation.

II. RADIATIVE CORRECTION
The experimentally observed cross section (σ obs ) is related to the Born cross section via where σ B is the Born cross section, F (x, s) has been calculated in Refs. [8,9] and 1 |1−Π(s)| 2 is the VP factor; the upper limit of the integration x m = 1 − s m /s, where √ s m is the experimentally required minimum invariant mass of the final state f after losing energy to multi-photon emission.
In this paper, √ s m corresponds to the BB mass threshold which is 10.5585 GeV.
The radiator F (x, s) is usually expressed as [8] F (x, s) = x β−1 β · (1 + δ ′ ) − β(1 − 1 2 x) with and Here the conversion of soft photons into real e + e − pairs is included. The Born cross section is thus calculated from where (1 + δ(s)) is the initial state radiation (ISR) correction factor. It is obvious that both (1 + δ(s)) and 1 |1−Π(s)| 2 depend on the Born cross section from threshold up to the CM energy under study, while the Born cross section is what we want to measure. These two factors can only be obtained using the measured quantities with an iteration procedure.
The pure ISR correction factor (1 + δ(s)) depends only on the line shape of e + e − → bb cross section, while 1 |1−Π(s)| 2 depends also on the R values in the full energy range, we use a two-step procedure to get the Born cross sections.

A. ISR correction factor
The ISR correction factor is obtained with an iterative method via with σ B 0 (s) |1−Π(s)| 2 = σ obs (s). The iteration is continued until the difference between the two consecutive results is smaller than a given upper limit. The result from the last iteration, denoted by (1 + δ f (s)), is regarded as the final ISR correction factor.

B. Vacuum polarization factor
A similar procedure is used to calculate the VP factor in the bottomonium energy region. In this calculation, however, the total hadronic cross section is used rather than that of e + e − → bb only. Instead of depending on the hadronic cross sections in the bottomonium energy region, the VP factor depends on the R values in the full energy region. In addition, there is also contribution from leptons.
The VP factor includes two terms [10] The first term is the contribution from the leptonic loops with for lepton with mass m. For 0 ≤ s < 4m 2 , define a = (4m 2 /s − 1) 1/2 , while for s ≥ 4m 2 , define a = (1 − 4m 2 /s) 1/2 and b = (1 − a)/(1 + a), The second term in Eq. (10) is the contribution from the hadronic loops. This quantity Π h (s) is related to the total cross section σ(s) of e + e − → hadrons in the one-photon exchange approximation through a dispersion relation Using the identity 1 We follow the procedure in Ref. [11] to calculate the first term in the above equation. First, the integration is performed analytically for narrow resonances J/ψ, ψ(3686), Υ(1S), Υ(2S), and Υ(3S). Second, for the high energy part, it is assumed that R(s) = R(s 1 ) is a constant above a certain value s 1 . And third, the integral between threshold and s 1 is carried out numerically after separation of the principle value part. Thus we have where Γ j , Γ j e + e − , and M j denote total width, partial width to e + e − pair, and mass of the resonance j, respectively. Here, σ nr (s) is the σ(s) in Eq. (15) with the contributions from narrow resonances subtracted.

For
where R EW (s) = 3Σ q e 2 q is the purely electroweak contribution neglecting finite-quark-mass corrections with e q the electric charges of the quarks; the QCD correction factor with parameters defined in Refs. [13,16].
In the bottomonium energy region, the dressed cross section of e + e − → bb is denoted by where σ obs (e + e − → bb) is the observed cross section provided by the Belle and BaBar collaborations [4,5], and the Born cross section of e + e − → u, d, s, c-quarks from the pQCD calculation is denoted by σ B udsc (s). Then σ B 0 (s) = σ B udsc (s) + σ dre b (s) is taken as zeroth order approximation of the Born cross section of e + e − → hadrons. Together with the Born cross sections in other energy regions we obtain the first order approximation of the VP factor,

C. Born cross section
The final Born cross section of e + e − → bb can then be calculated with Eq. (6) with the ISR correction factor and VP factor calculated above, i.e., III. THE DATA Both BaBar [4] and Belle [5] experiments measured R b in the bottomonium energy region: where the denominator is the Born cross section of e + e − → µ + µ − . In both experiments, neither ISR correction, nor the VP correction was considered, so the reported R b corresponds to the observed cross section. In both experiments, the contribution of the narrow Υ states from the initial state radiation, i.e., Υ(1S), Υ(2S), and Υ(3S) states can be removed from the data supplied in the papers. The BaBar measurement [4] was based on data collected between March 28 and April 7, 2008 at CM energies from 10.54 to 11.20 GeV. First, an energy scan over the whole range in 5 MeV steps, collecting approximately 25 pb −1 per step for a total of about 3.3 fb −1 , was performed. It was then followed by a 600 pb −1 scan in the range of CM energy from 10.96 to 11.10 GeV, in 8 steps with non-regular energy spacing, performed in order to investigate the Υ(6S) region. Altogether, there are 136 energy points [4]. In the BaBar paper, the ISR produced narrow Υ states, i.e., Υ(1S), Υ(2S), and Υ(3S) states were included in R b , but in the data file supplied, their contribution is listed and can be removed from the data.
The Belle measurement [5] was done with the scan data samples above 10.63 GeV at total 78 data points. The data consist of one data point of 1.747 fb −1 at the peak √ s = 10.869 GeV; approximately 1 fb −1 at each of the 16 energy points between 10.63 and 11.02 GeV; and 50 pb −1 at each of 61 points taken in 5 MeV steps between 10.75 and 11.05 GeV. The non-resonant qq continuum (q ∈ {u, d, s, c}) background is obtained using a 1.03 fb −1 data sample taken at √ s = 10.52 GeV. Belle experiment supplied a data file of R b with the ISR produced narrow Υ states removed (defined as R ′ b in Belle paper [5]). The BaBar and Belle measurements [4,5] are shown in Fig. 1. Notice that the definitions of R b are different in these papers. After removing the ISR contribution of the narrow Υ states from BaBar results, the R b values and the comparison between the two experiments are shown in Fig. 2. In the following analysis, R b refers to the results after removing the ISR contribution of the narrow Υ states, and R dre b refers to the dressed cross section after the ISR correction is applied, and R B b refers to the Born cross section after the ISR and VP corrections are applied.
We can see from Fig. 2 that Belle results are systematically larger than BaBar measurements. To get the size of the systematic difference, we calculate the ratio between the Belle and BaBar measurements in the energy region covered by both experiments. Figure 2 shows the ratio of the R b between Belle and BaBar measurements, the ratios are fitted with a constant with a good fit  quality, χ 2 /ndf = 56/77, where ndf is the number of degrees of freedom. It indicates that the Belle and BaBar measurements differ by a factor of which is more than 7σ from one if they are the same.

A. Combination of Belle and BaBar data
The BaBar experiment measured the R b above 11.1 GeV which is very flat. This indicates that the bottomonium resonance region has passed and the flat continuum region has reached. At CM energy well above the open-bottom threshold, R values (the total cross section of e + e − annihilation) and R b can be calculated with pQCD with five different flavors of quarks. And this can be compared with the BaBar measurement. If we assume the difference in R b between Belle and BaBar can be extrapolated to energy region above 11.1 GeV, by comparing the expected Belle measurements in this energy region and the pQCD expectation, we can check the normalization of the Belle data.
To compare R b with pQCD calculation, ISR correction and VP correction should be applied to the Belle and BaBar measurements since pQCD calculates the Born cross sections. Using the ISR correction factors (point by point correction, average correction factor 1 + δ ≈ 0.901 above 11.1 GeV) and VP factors ( 1 |1−Π| 2 ≈ 1.076) calculated below, a fit to the R B b from BaBar experiment for CM energies between 11.10 and 11.21 GeV yields R B b = 0.316 ± 0.011 with the error dominated by the common systematic error.
Assuming Eq. (21) applies to R B b at CM energy above 11.1 GeV for Belle measurement, we extrapolate the Belle measurement to this energy region so that we would expect Calculating R B b and the total continuum R values from udsc-quarks in pQCD according to Eq. (18), we find that R B b is almost a constant for CM energy between 11.10 GeV and 11.21 GeV, which is 0.351 with an uncertainty negligible compared with the experimental measurement; and R B udsc can be well parameterized as a linear function of CM energy ( √ s in GeV) between 10 and 12 GeV, i.e., Figure 3 shows the comparison among BaBar measurements, Belle expected, and the pQCD calculated R B b , we can find that Belle data agree with pQCD reasonably well (within about 1σ, the common error of the Belle measurements at high energy is about ±0.011, similar to the BaBar measurements) while the BaBar measurements are about 3σ lower than pQCD calculation. As a consequence, we assume BaBar measurement suffers from a normalization bias, and the Belle measurement is normalized properly. In the analysis below, we increase the BaBar measurements by the factor f in Eq. (21) and combine it with the Belle measurements and treat them as a single data set. The normalized and combined data are shown in Fig. 4.
In the rest of this work, the threshold of open-bottom production is set to be 10.5585 GeV, larger than the first two energy points in BaBar experiment. Therefore, these two data were omitted in our analysis.

B. Parametrization of R b
To calculate the ISR correction factors, the measured R b will be used as input. In order to avoid the point to point statistical fluctuation, one may parameterize the line shape with a smooth curve. There is no known function describing the line shape a priori, so one may parameterize the line shape with any possible combination of smooth curves.
We use the "robust locally weighted regression" or "LOWESS" method to smooth the experimental measurements. The principal routine LOWESS computes the smoothed values using the method described in Ref. [17]. The method works very well only for slowly varying data, which makes the procedure at the Υ(4S) region work improperly. As a consequence, we use the data points directly for √ s < 10.66 GeV and use the smoothed data for the other data points. Figure 5 shows the smoothed R b which looks very reasonable. In the following analysis, we use a straight line to connect two neighboring points. As these points are after smoothing and the step is not big, there is no big jump between neighboring points, we do not expect significant difference between a straight line and a smooth curve.

A. Calculation of ISR correction factors
We follow the procedure defined in Eqs. (7), (8), and (9) to calculate the ISR correction factors. In doing this for experimental data, we assumed the detection efficiencies for bb events without ISR and those with different energy of ISR photons are the same. This assumption is reasonable since there is no selection criteria related to the momentum of the final state B mesons in either BaBar [4] or Belle experiment [5]. The iteration is continued until the difference between two consecutive results is less than 1% of the statistical error of the observed R b .
In the energy region where the cross section varies smoothly, the ISR correction factors become stable after a few iterations while in the Υ(4S) energy region, due to the rapid change of the cross section in narrow energy region, the ISR correction factors only converge to within 1% after more than ten iterations. We iterate 20 times, and the maximum difference is less than 0.5% within the full energy region. Figure 6 shows the final ISR factors as well as the corrected R dre b values.

B. Calculation of vacuum polarization factors
Taking the ISR corrected R b , R dre b , (assuming it is R B b ), adding the pQCD calculation of the udsc-quark contribution to R values (refer to Eq. (23)), we calculate the VP factors in the bottomonium energy region. After obtained the VP factor 1 |1−Π| 2 , we use R dre b / 1 |1−Π| 2 as input to calculate 1 |1−Π| 2 again and we iterate this process. It turns out that after three iterations, the VP factor 1 |1−Π| 2 becomes stable so we take the values from this round as the final results, and we obtain R B b with Eq. (20). Figure 7 shows the VP factors from this calculation.

C. Estimation of errors
In the previous two subsections we obtain the ISR correction factor and VP factor and in turn the Born R B b via Eq. (20). During the calculation, however, only the central values of the observed R b are used. Instead of just scaling the errors of the original measurements by the obtained ISR correction factors and the VP factors, we perform a toy Monte Carlo sampling to investigate how the errors (both statistical and systematic) of the original measurements impact the obtained R B b . At each energy point where the Belle or BaBar measurement [4,5] was performed, we perform 10,000 samplings of the observed R b according to a Gaussian distribution whose mean value and standard deviation are the central value and statistical error of the observed R b , respectively. Besides, the uncommon systematic errors are also added to the samples in the same way. For the common systematic error, the same error is added to each sample at all energy points in Belle or BaBar experiment. While at the energy region 0.36 < √ s < 2.0 GeV and 3.7 < √ s < 5.0 GeV, the data from PDG compilation [1,13] and the BES collaboration [14,15] are assumed completely correlated when we perform the sampling. With each sample as input, we do the calculation described in the previous two subsections again to get the ISR correction factor, VP factor, and R B b for this sample. Finally a distribution of R B b at each energy point is observed, so are the ISR correction factor and the VP factor. We find that these distributions also satisfy Gaussian distribution well so the fitted mean and standard deviation are taken as the central value and error of the corresponding quantities, respectively. The covariances of the distributions of Born and dressed cross sections at different energy points are also available in the supplemental material [18]. After all the above operations, we obtain the R B b as well as its total uncertainty from the combined BaBar and Belle measurements as shown in Fig. 8 and the Table in the Appendix. We can find that the R B b values are very different from the R b reported from the original publications [4,5] and the differences are energy dependent. A common feature is that the peaks are even higher and the valleys become deeper, the two dips at the BB * + c.c. and B * B * thresholds are more significant, the peaks corresponding to the Υ(5S) and Υ(6S) increased significantly, and there is a prominent dip at 10.75 GeV.
The total R value corresponding to the production of udscb quarks can be obtained directly by adding the R B b to the udsc-quark contribution calculated from pQCD, as indicated in Eq. (23).

VI. SUMMARY AND DISCUSSIONS
From the BaBar and Belle measurements of the observed e + e − → bb cross sections, we do ISR correction to obtain the dressed e + e − → bb cross sections from 10.56 to 11.21 GeV. These dressed cross sections are the right ones to be used to determine the resonant parameters of the vector bottomonium states. Together with the R values measured at other energy points and the R values calculated with pQCD, we calculate the VP factors. By applying VP correction, we obtain the Born cross section of e + e − → bb from threshold to 11.21 GeV. These cross sections can be used for all the calculations related to the photon propagator, such as a µ , the µ anomalous magnetic moment, and α(s), the running coupling constant of QED [2,3].
In the following parts of this section, we discuss the usage of the data obtained in this study.

A. Vacuum polarization
The VP factors have been calculated by many groups [11,[19][20][21][22], with the experimental data and various theoretical inputs when the data are not available or less precise. Different techniques on how to handle the discrete data points and how to correct possible bias in data were developed. All these different treatments yield very similar results on hadronic contribution to a µ and on the running of the α at M 2 Z 0 , which indicates that the methods are all essentially applicable with current precision of data.
Previous calculations of the VP factors in the bottomonium energy region used either the resonant parameters of the Υ(4S), Υ(5S), and Υ(6S) reported by previous experiments [6,7] which are very crude [20,21] or the experimental data from previous experiments [6,7] which gave the observed cross sections [22]. We recalculate the VP factors by using R B b obtained in this analysis based on high precision data from BaBar and Belle experiments [4,5], with the ISR correction and VP factors properly considered. Although these new data have small effect on the VP factors far from the bottomonium energy region, they do change the VP factors in the bottomonium energy region as is shown in Fig. 9. The difference between this and the previous calculation [20] is significant comparing with the small errors of the results.  [20]. The solid lines are the central values and the error bars or bands show the uncertainties.

B. The bottomonium spectroscopy
There are very clear structures in R B b distribution shown in Fig. 8. From low energy to high, we identify the Υ(4S) at 10.58 GeV, dips due to BB * + c.c. and B * B * thresholds at 10.61 and 10.65 GeV, respectively, a dip at 10.75 GeV which may correspond to the Y (10750) [23], and the Υ(5S) and Υ(6S) at 10.89 and 11.02 GeV, respectively.
The observed R b values were used to extract the resonant parameters of the Υ(5S) and Υ(6S) in the BaBar [4] and Belle [5] publications. As the ISR correction effect is significant and is energy dependent, this suggests that the fit results are not reliable. To avoid the dip at around 10.75 GeV, both BaBar and Belle fitted data above 10.80 GeV only. A recent study of e + e − → π + π − Υ revealed a new state, the Y (10750), with a mass of (10752.7 ± 5.9 +0.7 where m, Γ, Γ e + e − , and φ are the mass, total width, electronic partial width of the resonance, and the relative phase between the resonance and the real continuum amplitude, and they are all free parameters in the fits. Eight sets of solutions are found from the fit [24], with identical total fit curve, identical fit quality (χ 2 = 274 with 188 data points and 13 free parameters), and identical masses and widths for the same resonance, but with significantly different Γ e + e − and φ. Figure 10 shows one of the solutions of the fit, and Table I lists the resonant parameters from 8 solutions of the fit. The masses and widths of the resonances agree with those from Ref. [23] but with improved precision because of the much better measurements used in this work compared with those in exclusive e + e − → π + π − Υ analyses. The Γ e + e − 's determined from this study allow us to extract the branching fractions of π + π − Υ of these resonances by combining the information reported in Ref. [23], and to understand the nature of these vector states [25][26][27][28][29].
In this analysis, we assumed that all the resonances are Breit-Wigner functions with constant widths and the continuum term is a smooth curve in the full energy region and they interfere with each other completely. In fact, the R B b or the total cross section has contributions from different modes, including open bottom and hidden bottom final states, the parametrization of the line shape should be very complicated due to the coupled-channel effect [30] and the presence of many open bottom thresholds, BB, BB * + c.c., B * B * , B sBs , B sB * s + c.c., B * sB * s , BB 1 + c.c., B * B 0 + c.c., B * B 1 + c.c., BB 2 + c.c., and so on, and even πZ b (10610), πZ b (10650). The situation becomes somewhat simpler in a single final state like π + π − Υ(nS) (n = 1, 2, 3) [5] and π + π − h b (mP ) (m = 1, 2) [31] although the intermediate structure in the three-body final state is also complicated. The results from these fits may change dramatically with more information on each exclusive mode included.
We also try to add one more Breit-Wigner function to fit the cross sections, the fit quality improves slightly with a state at m = (10848 ± 9) MeV/c 2 with a width of (28 ± 14) MeV, a state  Table I. The magnitudes of these components are different in different solutions.  1 Γ e + e − (eV) 10.7 ± 0.9 21.3 ± 1.0 9.8 ± 0.5 φ (degree) 260 ± 3 144 ± 2 34 ± 3 2 Γ e + e − (eV) 11.1 ± 0.9 24.8 ± 1.3 at m = (10831 ± 1) MeV/c 2 with a width of (14 ± 5) MeV, or a state at m = (11065 ± 23) MeV/c 2 with a width of (73 ± 37) MeV. In all these cases, the significance of the additional state is less than 4σ. The data obtained in this analysis can be used to extract resonant parameters of these states if a better parametrization of the cross sections is developed.

C. Search for the production of invisible particles
The experimentally observed e + e − → µ + µ − cross section, with radiative correction is expressed as where F (x, s) is expressed in Eq. (3), x m = 1 − s m /s, with √ s m the required minimum invariant mass of the µ pair in the event selection. Here the mass of µ is neglected for charm and beauty factories. In Eq. (24), the cross section is calculated by QED without any ambiguity except the vacuum polarization Π(s) which is expressed by Eqs. (10)(11)(12)(13)(14)(15), with the hadronic contribution depending on the experimental measured data as input.
If the e + e − → µ + µ − cross section is measured to a high precision, then Π(s) can be obtained. This way Π(s) is measured from experiment directly. It can then be compared with what is calculated by Eqs. (10)(11)(12)(13)(14)(15). This provides a test of QED at high luminosity flavor factories. In Eq. (10), the leptonic term Π l (s, m 2 ) is expressed in terms of QED fine structure constant α(0) and lepton masses which have all been known to very high accuracy; while the hadronic term Π h (s) must be evaluated with Eq. (15) with the input of experimental measured hadronic cross sections. It is seen in Eq. (15) that Π h (s) is most sensitive to the hadronic cross section at energies close to s, which can be measured with the same experiment. So such test of QED can be performed with two sets of data, e + e − → µ + µ − and e + e − → hadrons, collected within the same experiment. Any discrepancy would mean that there are missing hadronic final states, or even more, there could be final state which escape the detection, or is invisible by current detection technology. This provides a test of QED and search for new physics. We propose it to be included into the physics goals in future high luminosity frontier physics.