NNLO leptonic and hadronic corrections to Bhabha scattering and luminosity monitoring at meson factories

Virtual fermionic N_f = 1 and N_f = 2 contributions to Bhabha scattering are combined with realistic real corrections at next-to-next-to-leading order in QED. The virtual corrections are determined by the package bha_nnlo_hf, and real corrections with the Monte Carlo generators Bhagen-1Ph, Helac-Phegas and Ekhara. Numerical results are discussed at the energies of and with realistic cuts used at the Phi factory DAFNE, at the B factories PEP-II and KEK, and at the charm/tau factory BEPC II. We compare these complete calculations with the approximate ones realized in the generator BabaYaga@NLO used at meson factories to evaluate their luminosities. For realistic reference event selections we find agreement for the NNLO leptonic and hadronic corrections within 0.07% or better and conclude that they are well accounted for in the generator by comparison with the present experimental accuracy.


Introduction
The Bhabha scattering process e + e − → e + e − (1.1) is an invaluable tool for the luminosity determination in various experiments. Both low energy devices, operating from about 1 GeV to several GeV and high energy devices, planned to operate at hundreds or thousands of GeV, require theoretical predictions for the Bhabha cross section with quite accurate determinations of QED radiative corrections. The latter contain, besides exponentiated leading logarithmic terms, also the complete fixed order contributions, and in particular the complete two-loop QED corrections. Aiming at per mille accuracy or slightly better, the radiative corrections may neglect the constant terms in the electron mass m e . At next-to-leading order (NLO) final states with unresolved photons will contribute, e + e − → e + e − (γ). (1.2) Further, new mass scales start to play a role, but only in the one-loop self-energy insertions. At next-to-next-to-leading order (NNLO), further final states of Bhabha scattering contain unresolved photons, fermion pairs, or hadrons: e + e − → e + e − (γ, γγ), e + e − (e + e − ), e + e − ( f + f − ), e + e − (hadrons). (1.3) At this order of perturbation theory, a variety of Feynman diagrams depend on additional mass parameters, and one may formally distinguish between N f = 1 corrections (with only electrons) and N f = 2 corrections, and the latter ones are technically more complicated due to the additional mass scale.
The NLO corrections, with inclusion of certain leading higher order terms, are known since a while and several Monte Carlo (MC) programs are carefully tuned. A recent comprehensive review on precision predictions for scattering experiments at meson factories contains a detailed discussion of the state of the art [1].
As far as virtual corrections are concerned, in the last few years there has been major progress in the evaluation of the corrections at the NNLO accuracy. In fact, the photonic two-loop QED corrections were first evaluated in the massless case in [2]. The photonic corrections to massive Bhabha scattering with enhancing powers of ln(s/m 2 e ) were soon derived from that [3]. The missing constant term in m e [4] plus the corrections with electron loop insertions [5,6,7,8], called the N f = 1 case, followed not much later. The heavy fermion (or N f = 2) corrections were first derived in the limit m 2 e << m 2 f << s, |t|, |u| [8,9], where m f is the mass of the heavy fermion and s,t, u are the usual Mandelstam variables, and soon after also for m 2 e << m 2 f , s, |t|, |u| [10,11]. Finally, using dispersion relations, also hadronic corrections became known [12,13,14].
A complete collection of all the relevant formulae for the massive virtual NNLO corrections used in this paper can be found in the just mentioned papers and in [15].
When fermion loops or virtual hadronic corrections are taken into account, the question of considering also the real emission of the corresponding particles arises. That was studied for the emission of electron pairs in [16] in the soft limit of electron pair energy and in logarithmic accuracy. It is shown that the leading logarithmic corrections ln 3 (s/m 2 e ) cancel with those from the irreducible two-loop vertex corrections with electron loops. A similar cancellation is expected for the combination of heavy fermion pair emission with irreducible two-loop vertex corrections with a heavy fermion loop. In practice, however, the situation is evidently a bit more involved, especially at smaller energies, when s, |t|, |u| ∼ m 2 f . Then, the logarithms are not numerically dominating and more diagrams get important. Nowadays, MC programs can do that job. In this article, we will perform such a study of reaction (1.3) due to the additional emission of real pairs of leptons with the Fortran packages HELAC-PHEGAS [17,18,19,20], similarly to what was done in the 1990s for small-angle Bhabha scattering at LEP [21,22]. For heavy fermions and hadrons, we present here the corresponding results for the first time. The case of real hadron emission in Bhabha scattering deserves special attention. The only existing event generator contains only pion pair emission and the results obtained for the real pion pair emission serve as an indication of the size of other left over corrections. For a consistent treatment, we replace the hadronic dispersion integrals in the virtual and soft real hadronic corrections as described in [12,13] by the pion pair form factor. The prediction is then combined consistently with real pion pair emission as evaluated with the MC event generator EKHARA [23,24,25,26]. In the calculations we use the pion form factor from [27]. For the virtual and hard photon corrections the full hadronic corrections were obtained using the vacuum polarisation insertions.
The results may be compared with those from the Bhabha generators which are usually applied for experimental simulations of Bhabha scattering; here we look at BabaYaga [28,29,30], in particular to the latest and most accurate version BABAYAGA@NLO [31].
The aim of the article is to put together all the above discussed NNLO corrections to Bhabha scattering taking into account real experimental conditions and examine how well they are accounted for in the event generator BABAYAGA@NLO used at meson factories for their luminosity measurements. So far, at NNLO level, virtual corrections have been checked in detail only for situations where the dependence of soft radiation on the maximum soft photon energy ω (or, equivalently, the minimal hard photon energy) is "switched off" by setting ω = √ s/2 [4,12,13,14]. This was a good way to compare results obtained by different theoretical groups, but certainly has nothing to do with reality. We restrict ourselves here to low energies (meson factories) because presently they are of immediate relevance from the experimental point of view.
We just mention for completeness the last remaining NNLO issue: that of radiative loop corrections, i.e. the NNLO contributions from the interference of photonic bremsstrahlung off one-loop diagrams with lowest order real photon contributions, first studied in [32] with a restriction to the factorising diagrams. The technical complications arise from non-factorising diagrams, the socalled pentagon diagrams. Recent papers on this issue are [33,34,35,36], but so far without explicit numerical results. Sample numbers for the virtual one-loop (plus real soft) QED corrections to the hard-bremsstrahlung process e + e − → e + e − γ are given in [36]. For future measurements, it would be worthwhile to answer the question if (and when) these corrections have to be yet included in the MC event generators employed for simulating Bhabha scattering events at low-energy high-luminosity electron-positron colliders.
The paper is organised as follows. In Section 2 we discuss the exact NNLO massive corrections to Bhabha scattering and present benchmark results for event selections close to the experimental ones. In Section 3 we describe the approximate treatment of these corrections in the BABAYAGA@NLO event generator and derive benchmark results for the same event selections as for the exact results. In Section 4 we show detailed numerical studies of the quality of the approximations used in BABAYAGA@NLO. We draw our conclusions in Section 5.

The NNLO massive corrections
The complete NNLO N f = 1, 2 corrections to Bhabha scattering consist of three parts, each of them with contributions from virtual and real electron pair corrections (N f = 1 case) and corrections due to muon pairs, τ pairs and hadrons (N f = 2 cases): We want to concentrate here on the interplay of virtual and real corrections. For the various pure self-energy corrections in σ NNLO virt we refer to [1,13] and the references quoted therein. This is in accordance with the approach chosen in the MC packages used for the interpretation of experimental results, and we will not include these pure two-loop self-energy corrections in the numerical results discussed below.
As a result, the following contributions will be studied:  • virtual plus soft photonic corrections: • hard photon radiation: Here we take into account realistic experimental phase space cuts. The sum σ NNLO v+s + σ NNLO h now is independent of ω.
• Bhabha scattering with the additional production of fermionic pairs or of hadrons, σ LO real .  will take the exact expressions, and for the other contributions the approximation m 2 e << s, |t|, |u| is applied.
The irreducible vertex diagrams are infrared finite, but the reducible vertices and the box contributions are not. To make the latter two infrared finite, one has to add the corresponding soft real photon emission: Here, the term σ NLO virt × F γ,soft (ω) arises from the interference of (soft) single-photon bremsstrahlung diagrams, where one of the diagrams has a self-energy insertion. The σ NLO virt comes also from a Bhabha cross section due to diagrams with first-order fermionic (hadronic) self-energy insertions, and F γ,soft (ω) is the usual soft photon eikonal factor; for explicit expressions see [8,13]. The parameter ω is the infrared cut-off and has to be adapted such that soft-photon emission has the Born kinematics and the sum of soft and hard photon radiation is numerically independent of ω; typically, it is ω/E beam = 10 −6 · · · 10 −3 . The evaluation of the NNLO virtual corrections has been detailed elsewhere and we may restrict ourselves here to few remarks on the specifics of this article.
In the simplest case of a one-loop self-energy insertion, one may just evaluate Feynman diagrams and gets after renormalisation for the case of a light fermion with mass m, m 2 << s: For the virtual hadronic and heavy lepton pair corrections we use the dispersion approach. Here, the photon propagator is substituted in the Feynman diagrams as follows: with the propagator For heavy leptons with mass m l and charge Q l = −1, it is at one loop (exact in the mass m l ): For hadronic corrections, a natural choice for R had is a representation by experimental data: where σ had (z) ≡ σ (e + e − → γ → hadrons; z). The real hadronic emission can be studied at present only for pion pair production as only for this hadronic final state a MC code is available. Correspondingly, for the pion case we use instead of (2.14) the undressed pion form factor F π : where is the pion velocity. For the cross section ratio R, which we need here, this transforms to: The pion form factor F π (z) has been determined in [27]. The numerical studies using this parameterisation are presented in Section 4.
All this applies to one-loop insertions in reducible diagrams. For irreducible two-loop vertex and box diagrams, one has to perform an additional loop integration, and the result is more involved. The complete virtual NNLO N f = 1 corrections (due to electron self-energy corrections) σ NNLO v+s,e + e − are known exact in m e [7] in terms of harmonic polylogarithms [37]; a re-calculation and also corresponding formulae in the kinematic limit m 2 e << s,t, in terms of ordinary polylogarithms, are given in [8]. The virtual NNLO N f = 2 corrections (due to µ, τ, hadronic self-energy corrections) are determined with dispersion formulas.
At the end of this Section, we shortly comment on the structure of the various contributions to the virtual plus soft photon cross section: The eikonal factor is in the limit of small m e [13]: The virtual contributions are: The sum over M covers m µ , m τ , m π with the corresponding parameterisations of R M (z); for leptons see (2.13). The lower integration bound is M 2 0 = 4m 2 for leptons. For hadrons, where one sums over all the hadronic contributions, the lower bound is M 2 0 = m 2 π 0 , corresponding to π 0 γ, the lightest hadronic final state.
The kinematical factors C (rational functions of s and t) and the kernel functions K are universal. A special role play the irreducible vertex and box diagrams. In these diagrams, the self-energy correction is part of a loop insertion, and due to the additional loop momentum integration the replacement (2.11)-(2.12) leads to more involved kernel functions compared to (2.12); for the vertex [38]: Here Li 2 (x) is the usual dilogarithm and ζ 2 = Li 2 (1) = π 2 /6. From the irreducible box diagrams we have three different, lengthy box kernel functions K box,A (x, y, z), K box,B (x, y, z), K box,C (x, y, z); see for explicit expressions Eqs. (71)-(73) of [13].
For practical reasons, it makes sense to split the σ NNLO virt into two pieces, namely the infrared finite σ NNLO vert of (2.22) and the so-called "rest" σ NNLO rest , We just remind the reader that a third piece, the pure self-energy corrections, is not included in the study. The σ NNLO rest is the sum of all the infrared divergent contributions. In [13], it is detailed in Section VI and in Eq. (87). This sum is infrared finite, but depends on the photon energy cut ω related to the separation of soft and hard photons. In the energy regions of relevance with s > M 2 0 , the net result may adapted from Eq. (93) of [13]: The explicit expressions for F 1 to F 4 are given in (88)-(91) of [13]. They are infrared finite, but depend on 2ω/ √ s. It is well-known that the irreducible vertex corrections from a fermion pair with mass m contribute to the cross section with terms of order ln 3 (s/m 2 ). For electrons, this is a huge enhancement: 27) and the form factor is for m 2 e /x << 1 [39]: The vertex function may be found in [40] exact in m e . For heavy leptons, the logarithmic terms in , but a deviation appears in the constant term which becomes v 2 f [39,8]: It is these logarithmic dependences which make the real pair emission contributions so important, because cancellations of the leading terms appear.
Technically, we evaluate the electron corrections dσ NNLO v+s,e /dΩ with the Mathematica program CROSSSECTION.M [15]. The contribution dσ NNLO v+s,e dΩ is represented there by function NNLOEL. However, it includes also the iterated one-loop self-energies NNLOFE2 and the genuine two-loop self-energy NNLOFE1. Both are calculated separately, but are not included here in the numerics. They have been subtracted from NNLOEL in order to estimate what we call here σ NNLO virt,e . The heavy fermion corrections have been calculated with a Fortran package applying the dispersion technique described here. In order to cover also pion pair corrections σ NNLO v+s,π , the σ NNLO v+s is determined with an updated version BHA_NNLO_HF of the Fortran package BHBHNNLOHF [13,15].

Hard photonic corrections
The NNLO hard photonic corrections σ NNLO h with a self-energy insertion arise from the classes of diagrams shown in Fig. 2.3, with emission of one hard photon. They were calculated with the Fortran program BHAGEN-1PH-VAC [41] based on the generator BHAGEN-1PH [42]. This cross section depends on the soft photon cut-off E min γ = ω and only after adding them to σ NNLO v+s , the sum of the two σ NNLO v+s+h is independent of the cut-off. We calculate here separately the contributions from diagrams with the electron, muon, tau, pion and complete hadronic vacuum polarisation insertions. The dependence on additional cuts is crucial and varies considerably with the experimental set-up. A careful discussion is given in Section 4.
Even if the insertion of the vacuum polarisation corrections to the square of the tree level amplitude [43] is straightforward, for completeness we give below the formulae, which are used in the unpublished program [41]. We consider the process and follow the notation of [43] The differential cross section for the process (2.31) can be written as where E + , E − , E γ are the energies of the final positron, electron and photon, respectively. The quantities X,Y, Z refer to the s-channel annihilation, the t-channel scattering and the interference part of the squared amplitude, respectively. Keeping only the diagrams with a virtual photon exchange (the original formulae [43] contained also the weak Z boson contributions) and with the vacuum polarisation corrections included they read: (2.36) The generation of the phase space was not changed with respect to the original program BHAGEN-1PH [42].

Real electron pair contributions
The most important real fermion pair corrections to Bhabha scattering are, at all energies, the unresolved electron pair corrections. There are 36 diagrams of this kind contributing to σ LO e + e − (e + e − ) , part of the Bhabha cross section (2.1). Sample diagrams are shown in Fig. 2.4.
The e + e − pair corrections fall into three classes: four s-channel diagrams with two e + e − pairs in the final state, 24 diagrams (8 in s-channel and 16 in t-channel) with one e + e − pair, and 8 peripheral t-channel diagrams (no e + e − pair). What is usually considered as the electron pair corrections, are those with two electron pairs in the final state, Fig. 2.4(a). The contribution from such soft electron pairs is known [16], see also [13]. It is, in the limit of small m e and small energy cut-off parameter D of the unresolved e + e − pair, proportional to the lowest order Bhabha Born cross section σ LO e + e − : with: δ e so f t = 1 3 The parameter D has to fulfill: From the sum of (2.37) and (2.27), the compensation of the leading mass singularities (contained here in the L 3 s , L 3 t , L 3 u terms) in the cross section becomes evident. If there are unresolved contributions, which do not fulfill Born kinematics, or if the logarithms are not really big, or if the experimental accuracy is at the per mille level or better, a complete calculation is needed, and this is part of the present study.
The Feynman diagrams are finite as long as the electron mass is assumed to be finite, so that a straightforward Feynman diagram calculation of the 2 → 4 process can be performed without any true singularities.

Real muon and tau pair contributions σ LO
For both e + e − → e + e − µ + µ − and e + e − → e + e − τ + τ − there are 12 diagrams. Samples of them are shown in Fig. 2

.5.
There are four classes of diagrams, to be discussed here for muon corrections: two s-channel diagrams with production of both an e + e − -and a µ + µ − -pair; two s-channel diagrams with an e + e − -pair; six diagrams (two in s-channel, four in t-channel) with a fermion pair; and finally two peripheral t-channel diagrams. Contrary to real electron pair corrections, for heavy lepton pairs it is at the meson factory energies never appropriate to assume that the lepton mass is much smaller than e.g. √ s. So, one has to perform a complete lowest order 2 → 4 Feynman diagram calculation. Again, the cross section is finite as long as all masses are assumed to be finite. The same discussion holds for the process e + e − → e + e − τ + τ − .
The results for the electron, muon and tau pair corrections to the Bhabha scattering process have been obtained in the framework of the HELAC-PHEGAS leading-order MC program [17,20]. The phase space integration was executed with the help of PHEGAS [18], a general purpose multichannel phase space generator. HELAC-PHEGAS generates, in a fully automatic manner, events for all possible parton level processes at hadron and lepton colliders within the Standard Model. More precisely, integrated cross sections and kinematic distributions with arbitrary cuts on particles in the final state and with full spin correlations can be obtained. It has already been extensively used and tested in phenomenological studies, see e.g. [19,44,45,46].
In the present study, the exact QED 2 → 4 matrix elements for e + e − → e + e − + − processes, where ± = e ± , µ ± , τ ± , have been generated, including all Feynman diagrams (36, 12 and 12 respectively) and mass terms. Let us mention here, that in order to generate pure QED contributions in HELAC-PHEGAS the coupling of + − to the Z boson has to be simply set to zero. We have checked that the Z contributions are negligible for event selections used in this paper and do not affect any conclusions.

Real pion pair contributions
The lightest hadronic final state produced in e + e − scattering via the one photon exchange mechanism is the charged pion pair. This final state (i.e. e + e − π + π − ) was investigated in details in [24,25,26] and implemented into the Monte Carlo generator EKHARA [24,23]. Since other hadronic final states produced this way were never implemented in a generator, their studies are impossible at present and results obtained for the reaction e + e − → e + e − π + π − will be used also as a hint towards understanding the importance of other reactions like e + e − → e + e − + (π + π − π 0 , K + K − , K S K L , · · · ). The corresponding set of diagrams consists of 14 diagrams, with their representatives shown in Fig. 2.6. To model the pion-photon interactions we use the vector dominance model with the pion form factor from [27]. These contributions can be seen as: initial state electron pair emission (a), final state electron/pion pair emission (b,c), pion pair emission from the t-channel Bhabha process (d) and γ * − γ * pion pair production (e). The last set of diagrams is small for large electron and positron angles and its modelling is relatively crude, neglecting scalar meson production and the subsequent decays to pion pairs. Figure 2.6: Sample diagrams with real pion pair emission.

Other hadronic corrections
As discussed in the previous Section the lack of an event generator for processes e + e − → e + e − plus hadrons does not allow the computation of the real hadron emission with the exception of charged pion pairs. However contributions coming from virtual hadronic vacuum polarisation insertions can be calculated completely (see Section 2.7). In this Section we show how the vacuum polarisation insertions from pions compare to the full hadronic corrections. It has to be stressed that even if this difference can give some indication of the size of the missing real emission contributions, its actual size depends heavily on the event selection used by a given experiment. Thus a reliable estimation of the missing terms is not possible without Monte Carlo simulations. Specific problems are caused by narrow resonances like J/ψ, ψ(2S), ... and one should devote them a special treatment. Narrow resonances with mass M res and partial width Γ e + e − res can be described approximately by the ansatz R res (z) = 9π Based on this, their contributions to the NNLO Bhabha process can be derived from the general formulae of [13]. We discuss here as an example the contribution from the "rest" (Eq.  We conclude that for experiments performed on top of a narrow resonance, this resonance cannot be treated as a mere correction and more detailed studies have to be performed. These should include examining of finite width effects, beam spread effects, estimation of NNNLO corrections and the accuracy of the vacuum polarisation insertions in a close vicinity of these resonances. Having this in mind we do not present here hadronic contributions for the BES-III experiment running at J/ψ and ψ(2S) energies and we plan to devote to this issue a separate study.
We now come to the net hadronic vacuum polarization effects, i.e. look now at the sum of the so-called "rest" terms and the irreducible vertex corrections. To obtain all numerical results below we use "rest" as given by Eq. 2.26 together with the corresponding formula for the vertex, Eq. (2.22). The recent update of R had valid in the range m 2 π 0 < s < (100 GeV) 2 [47] is applied, and for higher s the R had is taken from [48]. For further details on the implementation of R had see Appendix E of [13].
The numerical integrations for the hadronic virtual and soft contributions were performed by means of the adaptive integration routine VEGAS [49], which works efficiently even for so narrow  In Fig. 2.7 we compare the full result for R had and the contributions coming from pions only. At low energies (Fig. 2.7 a) the biggest difference comes mostly from contributions of three pions and from kaon pairs pronounced at ω and φ resonances, while at high energies the pion contributions vanish rapidly and do not play any significant role (Fig. 2.7 b).
From Table 2.3 it is clear that the pion contributions are not sufficient for an accurate evaluation of the hadronic contributions and the lack of a generator for generic hadronic final states does not allow to draw a final conclusion about the real hadronic corrections with the exception of the charged pion pair emission and event selections which kinematically exclude the hadron production (KLOE energy).
The R enters the results with weight functions, which give more relevance to the low energy range, however the differences are important also there. Partially one can see the effect of the weight functions comparing the full hadronic corrections to the vacuum polarisation with pion pair contributions ( Fig. 2.8), however the complete weight functions are complicated and different for   the virtual and real contributions, so the careful evaluation of the integrals is needed.

Exact NNLO numerical results
In this Section we collect exact NNLO results which can serve as a benchmark for further investigations. The event selections used here are very close to event selections used at meson factories to measure their luminosity and are described in Appendix A. We give separately the contributions from electron, muon, tau and pion pair production as well as the complete hadronic contribution. For the last one, as mentioned already in the previous Section, there exist no generator to give the contributions from the real hadron emission beyond the pion pair production. An educated guess of the size of the missing contributions, based on the fact that the pions are the lightest hadrons produced and that the highest energy of the meson factories is about 10 GeV, is that they should not be much higher than the contributions from pion pairs. Thus based on information from Table 2.7 we can conclude that they should be completely negligible for the event selections used at meson factories for the luminosity measurements.
In Tables 2.4-2.8 the meaning of the different entries is the following: σ B is the Born cross section, σ v+s the cross section with NNLO virtual plus soft photon corrections (Section 2.1), σ h the NNLO cross section with a self energy insertion corrected by the emission of one hard photon (Section 2.2), σ v+s+h the sum of σ v+s and σ h and σ pairs the leading-order cross section with emission of real pairs (Sections 2.3-2.5). The total NNLO massive correction can be obtained by summing σ v+s+h with σ pairs .   (14) As one can see from Tables 2.4-2.8 the discussed corrections cannot be ignored for high precision luminosity measurements. Actually the relative size of electron pair corrections amount to about 0.3% at KLOE, 0.1 − 0.2% at BES and BaBar, and below the 10 −3 level at Belle only. The contribution of muon pair and hadronic corrections is slightly smaller, in the 0.05-0.1% range at all meson factories, while the contribution of tau pair corrections is, not surprisingly, generally Table 2.6: Results for tau pair corrections at different energies, in GeV, for the reference event selection defined in Appendix A. The σ B is the Born cross section within the acceptance cuts. The separation of soft and hard photons is fixed at ω/E beam = 10 −4 , where ω = E min γ (soft photon cut-off). All cross sections are given in nb.    (1) negligible.Therefore, the pair corrections have to be implemented into a MC event generator at least in an approximated form (see the next Section for their implementation in the generator BABAYAGA@NLO). In particular, for the real emission one can conclude that only the reaction e + e − → e + e − e + e − gives significant contributions to the cross section used in the luminosity measurements. When the accuracy of the experiment reaches the level 10 −3 this process has to be considered and its contributions added to the theoretical cross section or alternatively subtracted as a background from the experimental cross section.

The program
BabaYaga is an event generator for precise simulations of the processes e + e − → e + e − , µ + µ − , γγ in QED. It was developed for precision measurements with per mille accuracy of the luminosity of  (6) (1) GeV scale e + e − colliders. It has been adopted and is still presently used for this purpose by KLOE, BES/BES-III, CLEO, BaBar and Belle experiments. BabaYaga is released in two versions, the more accurate BABAYAGA@NLO and BabaYaga 3.5, http://www2.pv.infn.it/∼hepcomplex/babayaga.html. In BabaYaga 3.5 the most relevant QED radiative corrections are included in a pure Parton Shower (PS) approach, even though improved to include radiation interference effects for a more accurate description of radiative events. BABAYAGA@NLO also includes non-logarithmically enhanced O(α) corrections, which were the main source of theoretical error of BabaYaga 3.5. As summarised in Section 3.2, the necessary NLO ingredients are matched in BABAYAGA@NLO with the PS approach, in order to preserve exponentiation of large contributions and ensure normalisation at NLO accuracy.
Concerning the NNLO corrections that are the concern of the present study, the contribution of NNLO massive corrections included in the code is approximate and comes, as detailed in Section 3.3, from a) insertion of self-energy corrections in NLO virtual + soft photon correction; b) insertion of self-energy corrections in NLO hard photon correction.
From a theoretical point of view, the part b) coincides with the contribution to the NNLO corrections described in Section 2.2, although the e + e − → e + e − γ matrix element, phase space and related MC integration are completely independent of the calculation previously discussed. From a numerical point of view, one should expect a priori that the two calculations provide results in agreement within the respective statistical uncertainties, as it will be studied in the following. On the other hand, the contribution a) is just a subset of the full NNLO correction described in Section 2.1, as further detailed in Section 3.3.
Furthermore, the contributions due to real pair emission of the type 2 → 2 + (2), which are included in the complete NNLO benchmark calculation (Sections 2.3-2.5), are presently neglected in the program.

General formulation of BABAYAGA@NLO
As far as photon corrections are concerned, the two basic ingredients combined in BABAYAGA@NLO to guarantee the target theoretical accuracy are 1. exact NLO QED corrections (soft+virtual and hard contributions);

leading QED logarithms due to multiple collinear and soft radiation beyond O(α).
The matching is performed according to the following formula [31] is the Sudakov form factor. It describes universal (process independent) virtual + soft radiation up to the energy fraction ε (soft-hard separator) In the program, the scale Q 2 is chosen such that L = log Q 2 m 2 e = log st um 2 e − 1 ≡ L − 1 where s, t and u are the Mandelstam variables of the process and m e is the electron mass. This choice is dictated by the perturbative NLO calculation of the Bhabha process and ensures that, in addition to initial and final state leading contributions, also initial-final state interference leading effects are resummed to all orders.
2. |M n,LL | 2 : n-photon emission squared amplitude in collinear approximation, with phase space factor dΦ n 3. F SV and F H : residuals of the exact NLO calculation w.r.t. the leading log approximation ensured by the ingredients 1. and 2. above, i.e. The explicit expressions of all the above contributions can be found in [31].
A further necessary ingredient is the correction due to vacuum polarisation. It is included in the Bhabha Born matrix element setting r s = α(s)/α and r t = α(t)/α and rescaling the s and t channel amplitudes as follows In the code, we use the resummed expression α(q 2 ) = α/(1 − ∆α(q 2 )), where ∆α(q 2 ) is the fermionic contribution to the photon self-energy. It is treated analytically for the leptonic and top-quark one-loop contributions, while the non-perturbative five quark (hadronic) contribution, ∆α (5) had , is included according to the latest Jegerlehner [50] and Teubner et al. [47] parameterisations. In order to include an important class of O(α 2 ) factorizable corrections, we insert the vacuum polarisation correction in the NLO cross section too, both in the soft plus virtual contribution and the hard photon matrix element. The inclusion of the vacuum polarisation both in the soft plus virtual and hard photon part guarantees that their sum is independent of the soft-hard separator ε, as we explicitly checked.

NNLO massive corrections in BABAYAGA@NLO
As detailed in [31], it is possible to extract from the matched formula given in Eq. (3.1) the different pieces contributing to the cross section at NNLO. In particular, to explain how NNLO massive corrections are (approximately) taken into account in BABAYAGA@NLO, let us write the expansion up to O(α 2 ) of the cross section with soft plus virtual corrections. It can be derived from the first (n = 0) term of the infinite sum in Eq. (3.1). In order to highlight the s, t and interference contributions, first we define where V = −(2α/π)I + L is the O(α) truncation of the Sudakov form factor, E i and B i have been defined above and r s,t are the vacuum polarisation corrections for the s and t channels, respectively. If we define 1/(1 − ∆α(q 2 )) ≡ 1/(1 − δ q 2 ), the r 2 S , r 2 t and r s r t read r 2 s = 1 + 2δ s + 3δ 2 s , r 2 t = 1 + 2δ t + 3δ 2 t , r s r t = 1 + δ s + δ t + δ 2 s + δ 2 t + δ s δ t . (3.6) The first line of Eq. (3.6) is the Born contribution, the second line is the one soft photon plus one loop virtual correction (notice that it is equal to E s + E t + E st because B s + B t + B st = 1), the third line is the vacuum polarisation correction at O(α) and the remaining lines represent the cross section with O(α 2 ) soft plus virtual corrections. From the latter it is simple to disentangle the NNLO contribution due to the insertion of the vacuum polarisation correction in the O(α) soft plus virtual coefficients. The relevant terms are those containing a δ i factor. Among them, there is a pure two-loop self-energy correction (sixth line in Eq. (3.6)), which we discard for the comparison with the exact NNLO calculation, in accordance with the discussion of Section 2.1. Therefore the formula of interest reduces to Equation (3.7) is used in the present study to validate the approximate treatment of NNLO massive corrections as in BABAYAGA@NLO in the soft plus virtual regime. Equation (3.7) must be added to the contribution obtained by dressing the hard bremsstrahlung cross section with self-energy corrections. The hard photon matrix element is the sum of eight amplitudes where the real photon is attached to a s or t channel-like diagram. As in Eq. (3.2), those amplitudes are rescaled by r s and r t , respectively, to account for the effect of vacuum polarisation. In BABAYAGA@NLO, the squared amplitude for the emission of a real photon |M 1 | 2 is exact as calculated through FORM [51] and cross checked with the output of the ALPHA algorithm [52].
To summarize, the best knowledge of running α QED and all the virtual factorizable NNLO vacuum polarization corrections are included in BABAYAGA@NLO. In particular, within the full set of (reducible and irreducible) NNLO virtual massive corrections present in the exact calculation described in Section 2.1, only the subset of loop-by-loop corrections (see Figure 2.2) is taken into account in the code. In other words, the contributions of purely irreducible NNLO vertex and box corrections is not included. However, since the real pair emission corrections are neglected as well, this means that there is no imbalance of ln 3 (s/m e ) terms in BABAYAGA@NLO, thus respecting the compensation mechanism of leading mass singularities discussed in Section 2.1 and Section 2.3 (see Eq. (2.27) and Eq. (2.37)).

Numerical results of BABAYAGA@NLO
Below we give the benchmark results from the BABAYAGA@NLO MC event generator. Their detailed comparison to the exact results will be given in the next Section.  We repeat here that BABAYAGA@NLO does not generate the real pair emission contribution. This explains why the corresponding predictions for the cross section σ pairs and σ hadrons are not given in Tables 3.1-3.4. However, from Section 2.7 (see the next Section for a more detailed discussion) it is clear that only the reaction e + e − → e + e − e + e − gives a contribution which is relevant and needs a Table 3.3: BABAYAGA@NLO results for tau pair corrections at different energies, in GeV, for the reference event selection defined in Appendix A. The σ B is the Born cross section within the acceptance cuts. The separation of soft and hard photons is fixed at ω/E beam = 10 −4 , where ω = E min γ (soft photon cut-off). All cross sections are given in nb. Comparing the exact NNLO results of Tables 2.4-2.8 with those of BABAYAGA@NLO one can note that, while differences are present as expected, in the soft + virtual cross section, the independent predictions of the two calculations for the hard photonic correction σ h are in agreement within the MC statistical uncertainty.

BABAYAGA@NLO versus the exact NNLO massive corrections
In this Section we would like to answer how well the NNLO massive corrections are accounted for in the BABAYAGA@NLO event generator. The results are summarised in Table 4 In Table 4.1 we show the impact of the exact corrections and the corrections given by BABAYAGA@NLO for the reference event selections at meson factories. One can observe that the leptonic corrections are dominated by the electron corrections. Moreover they are well accounted for in the BABAYAGA@NLO code with the biggest unaccounted correction of the order of 0.5 ‰. In Figs. 4.1-4.4 we have studied the stability of the obtained results against changes of the event , where σ BY is the full (with all radiative corrections included) prediction of BABAYAGA@NLO according to Eq. (3.1). One can observe that even if for very stringent cuts one can have the differences up to 0.9 ‰ (Belle, Fig.4.3), generally the results are not varying rapidly. Comparisons between results for Belle ( Fig.4.3) and BaBar event selections (Fig.4.4) show also that the actual difference is sensitive to the event selection used and it is recommended to study the effect on the NNLO massive corrections, whenever the event selection is changed by a given experiment. Similar effects are observed for the hadronic corrections (with the exception of J/ψ and ψ(2S) energies, which will be examined in a separate paper). The biggest hadronic correction missing in BABAYAGA@NLO is about 0.4 per mille and in addition for KLOE and BES energies the missing hadronic contribution is of the opposite sign of the missing leptonic contribution, thus partly cancelling each other. For B-factories the missing leptonic and hadronic corrections are of the same sign, but even there the sum of the missing parts does not exceed one per mille (0.7 per mille for the reference event selections).

Conclusions
We presented an exact calculation of NNLO massive corrections to Bhabha scattering, beyond approximations previously addressed in the literature, and performed detailed studies of the missing part of these corrections in the BABAYAGA@NLO MC program. The approximate BABAYAGA@NLO formulae were confronted with the exact NNLO results for event selections close to the experimental ones. The stability of the results with changing of the event selections was also examined.  We found that NNLO massive corrections are relevant for precision luminosity measurements with 10 −3 accuracy, their relative total contribution to the Bhabha scattering cross section being of a few per mille. Hierarchically, the electron pair contribution is the largely dominant part of the correction, the muon pair and hadronic correction, which are the next-to-important effect, are quantitatively on the same grounds, while the tau pair contribution is negligible for the energies of meson factories.
Thanks to the exact NNLO calculation we tested the theoretical accuracy of the generator BABAYAGA@NLO which includes such corrections according to an approximate treatment through insertion of self-energy contributions into NLO correction factors. On the grounds of several numerical results, we concluded that the very bulk of the correction is taken into account in the program. For reference realistic event selections the maximum observed difference is at the level of 0.07%. When cuts are varied the sum of the missing pieces can reach 0.1%, but for very tight hadrons:  acollinearity cuts only.
As a possible perspective, it is worth mentioning that the leading logarithmic part of the missing pieces in BABAYAGA@NLO, coming from the interplay between NNLO virtual and real pair corrections, could be accounted for by means of appropriate singlet/non-singlet QED structure functions, as discussed e.g. in [53,54,55] and done in the past for small-angle Bhabha scattering at LEP [56]. However, this improvement of the theoretical formulation of BABAYAGA@NLO would be necessary only whenever the experimental precision of the luminosity measurements should require it and should be accompanied by the inclusion of other sources of NNLO ingredients presently neglected in the code and contributing at an accuracy below the 10 −3 level [1,57]. As already discussed in the paper, open issues of the present work are a more detailed study of hadronic NNLO corrections and of the uncertainty induced by hadronic vacuum polarisation insertions to Bhabha scattering in a close vicinity of narrow resonances. We plan to address the above perspectives in future works.
KEK. This work is a part of the activity of the "Working Group on Radiative Corrections and Monte Carlo Generators for Low Energies" [http://www.lnf.infn.it/wg/sighad/]. M. Worek acknowledges support by the Initiative and Networking Fund of the German Helmholtz Association "Physics at the Terascale", contract HA-101. This work is also supported in part by Sonderforschungsbereich/Transregio SFB/TRR 9 of DFG "Computergestützte Theoretische Teilchenphysik", European Initial Training Network LHCPHENOnet PITN-GA-2010-264564 and by Polish Ministry of Science and Higher Education from budget for science for 2010-2013 under grant number N N202 102638.   • | cos θ | < 0.8, where θ is the polar angle of the electron or positron in the lab system, this corresponds to the barrel region of BES-III detector. Since in BEPC, e + and e − beams are colliding with equal energy but at a 22mrad crossing angle, the lab system is slightly different from the CoM system.
• E e + > 1.0 GeV and E e − > 1.0 GeV, where E is the energy deposited in the electromagnetic calorimeter (EMC). 1 To see how the NNLO corrections depend on the event selection we obtained results also for | cos θ | < 0.7, 0.75, 0.85 and 0.9. Belle runs at an asymmetric e + e − collider, but all criteria are expressed in the CoM frame. The selection of Bhabha events by Belle are: • √ s =10.58 GeV • 50.5 • < θ ± < (180 − 50.5) • • Two charged tracks have momentum > 2.645 GeV • The track with maximum deposited energy in EMC greater than 2 GeV, • The sum of the deposited energies of all tracks in EMC is greater than 4 GeV (both charged and neutral particle can deposit energy in EMC and it is not checked if the particle is charged or neutral) • Acollinearity angle (2D) ξ 2dmax < 10 • • Transverse momentum of any observed charged particle greater than 0.