Analytic forms for the e + e − annihilation cross sections around a resonance including initial state radiation

For the first time, the exact analytic integration of the multiplication of the Kuraev-Fadin radiative function and the Born-level cross section has been achieved for e + e − annihilation around a narrow resonance. The analytic result is in a perfect agreement with the one obtained through direct numerical integration. The extraction of resonance parameters is demonstrated through the unfolding of the initial state radiation effect. The analytic integration significantly accelerates the regression process by more than 160 times compared to numerical integration. Moreover, it facilitates the determination of experimental measurements for physical parameters, such as the two-dimensional confidence region of the branching fraction and the relative phase between the strong and electromagnetic amplitudes, within just one hour instead of one week typically required for numerical integration.


I. INTRODUCTION
The e + e − annihilation experiments provide cleaner experimental environments for quarkonium decays than pp and hadron decays experiments.The measured quarkonium decay widths, or the corresponding branching fractions, can serve as inputs in phenomenological models or be used for comparison with theoretical predictions to test our understanding of quantum chromodynamics (QCD) [1].There is an unavoidable background from the continuum process, which directly produces the final state via e + e − annihilation, i.e., e + e − → γ * → f [2].It has been shown that the cross section from the interference term will lead to imprecise branching fractions for the broad resonance above open heavy flavor threshold, such as ψ(3770) [3].
A novel study shows that the ratio of the cross section from the interference term with respect to the resonance is surprisingly large compared to the precision of the current experiments even for the narrow resonances below the open heavy flavor threshold, such as J/ψ and ψ(3686) [4].Therefore, to achieve precision measurements of the quarkonium decay widths or the corresponding branching fractions, the cross sections around the resonance of at least three energies must be measured [4].
While the Born cross section of the final state is of interest, it is the experimentally observed cross section data that are measured which contain the initial state radiation (ISR) effect [5,6].ISR is an essential quantum electrodynamics (QED) correction for the precise measurement of cross sections in e + e − annihilations.In the ISR process, the electron or positron emits one or more photons before annihilation, thereby reducing the center-of-mass (c.m.) energy.One well-known method to take into account the ISR effect is the structure function method which is an integral transformation, i.e. convolution, of the kernel cross section by the radiative function [7].Therefore, obtaining the Born cross section is mainly an unfolding process of the ISR effect.
One commonly unfolding method involves employing an iterative procedure to obtain the discrete Born cross sections from the experimentally observed cross section data and by using the Monte Carlo method to obtain the uncertainties [5,8].Recently, a new method, named as the naive method, can obtain both the discrete Born cross sections and the uncertainties simultaneously by solving numerical integral equations [6].The advantage for the above methods is that their results are convenient to compare with results obtained in other experiments.With the obtained discrete Born cross sections, the parameters, such as decay width or branching fraction can be determined by fitting with a certain model.However, the above methods cannot apply to the e + e − annihilation processes around a narrow resonance, such as J/ψ and Υ(1S), because the Born cross section varies sharply and the c.m. energy spread is much greater than the resonance width and the distances between the c.m. energy points of data around the resonance.One of the most frequently used methods is to fit the experimentally observed cross section data by the Born cross sections specified in a certain model with the effects of both ISR and c.m. energy spread [9].
The effect of the c.m. energy spread is taken into account as an additional Gaussian convolution.The integral method is straightforward but time-consuming, as it requires numerical integrations in the regression iterations.
The integral method will be significantly accelerated if the analytic form of the ISRcorrected cross sections is obtained [10].It is well known that the radiative function proposed by Kuraev and Fadin (KF) formulism consists of an exponentiated part and a finite-order leading-logarithmic part, accounting for soft multiphoton emission and hard collinear brems-strahlung, respectively [7].In 1987, R. N. Cahn was the first to derive the analytic form with the exponentiated part in the KF radiative function [11].The approximated analytic form, including the leading-logarithmic parts in the KF radiative function and the upper limit correction by the exponential expansions, has been developed in subsequent works [9,10,12,13].
However, it is shown below that for the ISR-corrected cross sections around J/ψ, the accuracy of the approximated analytic form is within a few percent and propagates in the same order to the estimated value by regression for the hadronic branching fraction.In a sense, the approximated analytic form introduces an imperceptible yet non-negligible systematic uncertainty, which is comparable to the statistical or even the total systematic uncertainty for experiments.
In this paper, we provide the formalism of the ISR corrected Born cross sections and validate the precision of the exact analytic integrations in Sections II, III and IV introduce the formalism that incorporate the vacuum polarization effect and the c.m. energy spread effect, respectively, demonstrating the corresponding precision.Section V presents the fit performance tests of the toy Monte Carlo (MC) samples in the vicinity of J/ψ.Section VI gives discussions and conclusions.Appendix A gives the expressions of the Born cross sections around a resonance.Appendix B provides the corresponding exact analytic integrations of the ISR-corrected cross sections incorporating the KF radiative function, referred as the KF analytic forms.

II. ISR CORRECTED BORN CROSS SECTIONS
In the vicinity of a resonance, the amplitude A f tot. of the final state f in e + e − colliders is the coherent sum of both resonance A f R and continuum amplitudes A f C .The Born cross section can be written as [4] where ϕ is the relative phase between the continuum amplitude A f c and the resonance amplitude A f R .Therefore, the Born cross section is a sum of three parts as where } denote the Born cross sections from continuum, resonance and interference contributions, respectively.
When the kernel cross section is the Born cross section, the ISR corrected cross section ISR is an integral of the Born cross section σ f Born times the radiation function [4], that is where W is the c.m. energy of e + e − annihilation and W min is the threshold energy equaled to the invariant mass of the final states or the experimental cutoff energy.The KF radiative function F (x, W ) has the form [7] and β = 2α π 2 log W me − 1 .The ISR corrected cross section integral σ f ISR is a sum of three parts, which is where and are the integrals of the continuum, resonance and interference contributions, respectively.
For the process e + e − → µ + µ − , in the vicinity of J/ψ, the Born cross section is [9] where α is the fine structure constant, M and Γ are the mass and total decay width of J/ψ, Γ ee and Γ µµ are the decay widths of J/ψ → e + e − and µ + µ − , respectively.For the processes e + e − → 2(π + π − )π 0 and ηπ + π − with η → π + π − π 0 in the vicinity of J/ψ, the Born cross section can be written in a general form as [9] where A W 2 is the form factor, C 1 and C 2 are the ratios of amplitudes, and Φ is the phase between the strong and electromagnetic decays from J/ψ.For the process of e + e − → 2(π + π − )π 0 , C 1 and ϕ are assumed to be 1 and 0, respectively [9].
It can be derived for the cross sections from continuum, resonance and interference contributions.For example, the Born cross section from the continuum contribution of The Born cross sections are composed by rational functions, while the KF radiative function contains not only exponential functions but also logarithmic functions, such as log x, which makes the ISR corrected cross section integral an integrable singularity at the lower limit x = 0 in Eq. ( 3).Fortunately, the integral still has an analytic form.For example, the analytic form for the improper integral of log where Li 2 (z) is Spence's function which is defined as Since the analytic forms are tedious, we put them into Appendix B.

III. VACUUM POLARIZATION EFFECT
When the vacuum polarization (VP) effect is considered, the Born cross section in the vicinity of a resonance is corrected as [13,15] where Π 0 (W ) is the nonresonant VP factor [15].
The ISR corrected cross section with VP effect is also a sum of three parts which is where (PDG) world averages [14] and W min = 2.0 GeV.In the top row, black lines represent KF method results, blue dashed lines represent approximated analytic form results and red dotted lines represent NI method results.In the middle row, red dashed lines represent the results of In the bottom row, black lines represent the results of and are integrals of the resonance, continuum, and interference contributions, respectively.Approximations are used in Eqs. ( 18) and ( 20) compared to Eqs. ( 17) and ( 19) because the nonresonance VP factors Π 0 (W ) are numerically obtained [16].

IV. C.M. ENERGY SPREAD EFFECT
For the narrow resonances, such as J/ψ whose decay width is 92.6 keV, the c.m. energy spreads of e + e − are much larger than the resonance widths.For example, the c.m. energy spread around J/ψ in BEPCII is less than 1 MeV [9].Therefore, the effect for the c.m. energy spread must be considered in the experimentally observed cross section σ f exp (W ) by a Gaussian convolution with σ f ISR−VP (W ) which is In Fig. 3, we compare ISR corrected cross sections including both the VP and c.m. energy spread effects σ µ + µ − NI−exp obtained through the NI with σ µ + µ − KF−exp obtained using the KF analytic form and σ µ + µ − Approx−exp obtained using the approximated analytic form.The relative deviation between σ µ + µ − NI−exp and σ µ + µ − KF−exp remains at the order of 10 −4 .Moreover, the relative deviation at ϕ = 0 • , 90 • and 180 • .The Born cross section parameters are set to be PDG world averages [14] and W min = 2.0 GeV.In the top row, black lines represent KF method results, blue dashed lines represent approximated analytic form results and red dotted lines represent NI method results.In the bottom row, red dashed lines represent the results of − 1 and black lines represent the results of between σ µ + µ − KF−exp and σ µ + µ − Approx−exp is still 0.3% at Φ = 90 • and 3% at Φ = 0 • and Φ = 180 • .It is evident that the Gaussian convolution of the c.m. energy spread does not impact the precision.

V. FIT RESULTS USING TOY MC SAMPLES
We use the toy MC samples for e + e − → µ + µ − and e + e − → 2(π + π − )π 0 around J/ψ to test the KF analytic form and compare it to the approximated one.Firstly, the toy MC samples for the ISR corrected cross sections of e + e − → µ + µ − are generated with 1% uncertainty at 20 c.m. energies around J/ψ.The minimized χ 2 µ + µ − is performed with the free parameters M , S E and ϕ.The χ 2 µ + µ − is expressed as parameters are set to be PDG world averages [14] and W min = 2.0 GeV.In the top row, black lines represent KF method results, blue dashed lines represent approximated analytic form results and red dotted lines represent NI method results.In the bottom row, red dashed lines represent the results of − 1 and black lines represent the results of where σ µ + µ − i is the cross section at every energy point W i and ∆σ µ + µ − i is the corresponding uncertainty.The estimators M and S E and their corresponding uncertainties σ M and σ S E are taken as the information of the e + e − collider.
Secondly, the toy MC samples for the ISR corrected cross sections of e + e − → 2(π + π − )π 0 are then generated with 1% uncertainty at the same 20 c.m. energies around J/ψ.The minimized χ 2 5π is performed with the free parameters A, C 2 and Φ, and the nuisance parameters M and S E .Therefore, the estimators A, C 2 and Φ can be obtained.The χ 2 5π is expressed as The branching fraction for J/ψ → 2(π Finally, the above toy MC sample generations and χ 2 minimizations are repeated 10,000 times.The fit results of e + e − → µ + µ − are shown in Fig. 4. The estimators of M and S E are centered around the MC truth values for both KF and approximated analytic forms.
Only the estimator ϕ of the KF analytic form is centered around the MC truth value, while the distribution of the estimator ϕ of the approximated analytic form deviates more than 2σ from the MC truth value.The χ 2 µ + µ − of the KF analytic form is better than that of the approximated analytic form.
The fit results for the process e + e − → 2(π + π − )π 0 are depicted in Fig. 5. Two solutions, labeled as Solution I and II, are obtained.The distributions of χ 2 5π are identical for both Solution I and II.The χ 2 5π of the KF analytic form is better than that of the approximated analytic form.Only in Solution I the estimator Φ is centered around the MC truth value of 90 • , while in Solution II, Φ is centered around −90 • .Furthermore, in Solution I, only the estimator Br of the KF analytic form is centered around the MC truth value, whereas the estimator of the approximated analytic form deviates significantly.The toy MC test shows that the KF analytic form avoids a subtle yet non-negligible uncertainty that arises from employing the approximated analytic form.
Table I shows the computing time for the fitting processes.The fitting processes took less than 10 seconds for e + e − → µ + µ − and approximately 6 seconds for J/ψ → 2(π + π − )π 0 with the KF analytic form, while it required about 20 and 16 minutes with NI.Consequently, the KF analytic form can improve the fitting speed by more than 160 times compared to NI.
As it is shown in Fig. 1, the relative derivation between the analytic result and the numerical integration can be pushed down to the level of 10 −12 .Moreover, the different contributions, such as VP effect, can be switched off and on in the case of numerical integration and then test physical effects.As it is shown in Table I, for one typical toy MC sample , S E , ϕ and χ 2 µ + µ − , respectively.Histograms with black edges represent the KF analytic form, blue edges represent approximated analytic form and red lines are the values of MC truth.set, the time consumption of the fitting for J/ψ → 2(π + π − )π 0 is only about 16 minutes.It seems that there is no real practical reason to avoid numerical integration, as it is reliable and fast enough with modern computers.However, the current data analysis no longer solely focuses on reporting the confidence region of the branching fraction Br or the phase angle Φ separately.Equation (24) shows that the branching fraction for J/ψ → 2(π + π − )π 0 is correlated to the phase angle Φ.It is better to provide the confidence regions of both branching fraction Br and phase angle Φ.The 1σ, 2σ, and 3σ two-dimensional confidence regions are estimated using the ∆χ 2 value of 2.3, 6.18 and 11.83 relative to the best fit, respectively, as shown in Fig. 6.The computation time for obtaining the two-dimensional confidence regions is only 55 minutes using the KF analytic form.Given that the KF analytic form is over 160 times faster than the NI method with similar NFC, the estimated computation time with the NI method exceeds 6 days.Consequently, the impracticality of the NI method is evident when considering the computation time for two-dimensional confidence regions.

VI. DISCUSSIONS AND CONCLUSIONS
The ISR correction is essential for precise cross section measurements in e + e − annihilation.For the first time, we present the exact analytic forms for the ISR corrected cross sections around a narrow resonance with the QED structure function formalism by KF.The KF analytic forms can accelerate the regression procedure used to extract physical parameters, such as the branching fraction, by over 160 times.Utilizing the toy MC samples, we reveal a non-negligible few percent systematic uncertainty caused by the approximated analytic form.
The KF formalism is well known [7].It has been applied for description of QED radiative corrections to many observables in high-energy physics, including cross sections of e + e − annihilation processes at different energies.Moreover, the KF formalism was realized in several MC codes.In particular, BabaYaga [17] and MCGPJ [18] programs are widely used at e + e − colliders.It is also coded within the framework of BesEvtGen for the BESIII experiment [19].Since there is a certain freedom in the exponentiation prescription of the QED radiative corrections [20], there are many other exponentiation recipes of the QED structure function formalism [21].Although the uncertainty caused by VP and c.m. energy spread effects is below 0.1%, it does not mean that the theoretical uncertainty of the KF approach can reach the same order.As one can learn from Refs.[17,18,22], such a high precision cannot be provided without application of the complete one-loop corrections and the leading effects in orders higher than O(α 2 ) [23,24].Moreover, certain observables, especially around resonances, require calculation of next-to-leading logarithmic corrections of the order O(α 2 L 1 ) which appear beyond the original formalism of Kuraev and Fadin.Comparisons with the results of alternative realizations, e.g., with BabaYaga [17] and MCGPJ [18], should be performed to estimate the resulting theoretical precision.The difference should be taken as the systematic uncertainties of event selection efficiencies for the experimentally observed cross section data.
For the currently running e + e − collision experiments, the BESIII experiment has accumulated 10 billion J/ψ events [25] and 3 billion ψ(3686) events [26] which are at least 6 times larger than those used in previous BESIII measurements.Furthermore, the Belle II experiment plans to take about 500 fb −1 data for each vector bottomonium state [27] which are at least ten or hundreds of times larger than the Belle experiment.The KF analytic forms will be helpful to handle the ISR correction properly in currently running experiments at BESIII and Belle II, as well as in planned ones, such as the super-tau-charm factories [28] and the super-J/ψ factory [29].
the resonance part is and the interference part is

FIG. 5 .
FIG. 5.The histograms of the estimators obtained from toy MC samples for e + e − → 2(π + π − )π 0 .(a), (b), (c) and (g) show the distributions of A, C 2 , Φ and Br for Solution I, respectively; (d), (e), (f) and (h) show the distributions of A, C 2 , Φ and Br for Solution II, respectively; (i) shows χ 2 5π distributions for both Solution I and II.Histograms with black edges represent the KF analytic form, blue edges represent approximated analytic form and red lines are the values of MC truth.

FIG. 6 .
FIG. 6. Confidence regions of Φ and branching fraction Br from the fit of one typical toy MC sample set for e + e − → 2(π + π − )π 0 .The 1σ, 2σ and 3σ two-dimensional confidence regions are estimated using ∆χ 2 values of 2.30 (red), 6.18 (green) and 11.83 (blue) relative to the best fit.The left and right upper panels provide the one-dimensional ∆χ 2 for Φ obtained by profiling branching fraction (blue line), the dashed lines mark the corresponding 1σ, 2σ, and 3σ intervals; and the red dotted line represents MC truth value.The left and right panels are the same, but for branching fraction, with Φ profiled.The black points mark the best estimates, and the error bars display their one-dimensional 1σ confidence intervals.The red dot and dash-dotted lines mark the MC truth values.All the fit results took 55 minutes.

TABLE I .
The comparison of the computing times and the number of function calls (NFC) for fitting with the typical MC samples.Fitting codes are written in C++ and run on a single core of an Intel 3 GHz CPU.