Bhabha scattering at NNLO with next-to-soft stabilisation

A critical subject in fully differential QED calculations originates from numerical instabilities due to small fermion masses that act as regulators of collinear singularities. At next-to-next-to-leading order (NNLO) a major challenge is therefore to find a stable implementation of numerically delicate real-virtual matrix elements. In the case of Bhabha scattering this has so far prevented the development of a fixed-order Monte Carlo at NNLO accuracy. In this paper we present a new method for stabilising the real-virtual matrix element. It is based on the expansion for soft photon energies including the non-universal subleading term calculated with the method of regions. We have applied this method to Bhabha scattering to obtain a stable and efficient implementation within the McMule framework. We therefore present for the first time fully differential results for the photonic NNLO corrections to Bhabha scattering.


Introduction
Electron-positron or Bhabha scattering is one of the best studied processes in the Standard Model [1]. It is well suited for luminosity measurements at e + e − colliders because of its large cross section and clean signature. Furthermore, for energies well below the electroweak scale the radiative corrections are dominated by quantum electrodynamics (QED) which allows for a very precise theory prediction. As a consequence, much work has been put into the calculation of higher-order matrix elements as well as the development of Monte-Carlo event generators.
The next-to-leading order (NLO) matrix elements have been known in the full Standard Model for quite some time [2][3][4][5]. At next-to-next-to-leading order (NNLO) the situation is different. In the case of the electroweak corrections only logarithmically enhanced terms have been calculated [6][7][8][9]. On the QED side much more is known. The full two-loop matrix element with vanishing electron mass was calculated some time ago [10]. Subsequently, this result was extended to also include leading-order mass effects [11][12][13][14]. The subset of the two-loop matrix element containing closed electron loops has been computed without any approximations [15]. Although the exact mass dependence of the full two-loop contribution is still not known, leading power-suppressed mass effects were recently taken into account in [16]. The one-loop corrections to the radiative matrix element were calculated in [17].
In addition to the work that has been put into the calculation of the matrix elements various Monte-Carlo event generators were developed, combining the matrix elements to physical observables such that non-trivial detector geometries and acceptances can be taken into ac-count [4,[18][19][20][21][22][23][24][25][26][27]. In particular, the BABAYAGA event generator that is based on the matching of the exact NLO results to a parton shower algorithm has achieved a precision of below 0.1% [28]. A detailed analysis of the impact of fixed-order fermionic NNLO contributions was presented in [29].
Even though all necessary ingredients are available, a Monte Carlo that includes also NNLO photonic corrections was missing. The main bottleneck in this regard has been the real-virtual contribution that suffers from numerical instabilities when integrated over the phase space of the emitted photon. The source of these instabilities can be traced back to the disparate scales in the process introduced by the small electron mass that acts as a regulator of collinear divergences. This problem is exacerbated in the presence of soft radiation.
In this paper we present a method to reliably integrate the real-virtual matrix element over the full phase space. It is based on the expansion for small photon energies E γ ≡ ξ × √ s/2 including the non-universal next-to-soft contribution at O(ξ −1 ). To verify our method we have compared with approximate results from BABAYAGA at the cross section as well as at the differential level and found agreement within the expected 0.1% precision. With this method it is therefore possible to make reliable predictions for Bhabha scattering at the differential level including the full set of NNLO QED corrections. This paper is organised as follows: We begin by briefly introducing our calculational framework in Section 2. The main result is presented in Section 3 where we describe how the stabilisation of the real-virtual matrix element was achieved via the next-to-soft approximation. We verify our method in Section 4 and conclude in Section 5.

Overview of the calculation
We consider the scattering process e − (p 1 )e + (p 2 ) → e − (p 3 )e + (p 4 ){γ(p 5 )γ(p 6 )} (1) up to NNLO in QED. As we are mainly interested in establishing the stabilisation method we restrict ourselves to purely photonic corrections, i.e. we do not take into account contributions from closed fermion loops. Ultraviolet (UV) and infrared (IR) divergences are regularised in d = 4 − 2ǫ dimensions and the renormalisation is performed in the on-shell scheme. All tree-level and one-loop matrix elements were calculated with the full electron mass dependence. The corresponding diagrams were generated using QGraf [30] and evaluated with the Mathematica code Package-X [31]. In the case of the numerically delicate real-virtual matrix element this Mathematica calculation serves mostly as a reference calculation. In the bulk of the phase space we instead rely on OpenLoops [32,33]. As we will discuss in Section 3, for small photon energies we switch to a next-to-soft approximation.
As mentioned in the introduction the full mass dependence of the photonic two-loop matrix element is not known. However, for most practical applications we can assume the scale hierarchy m 2 ≪ Q 2 ∈ {s, t, u} with the electron mass m and the Mandelstam invari- For our purposes it is therefore justified to neglect the power-suppressed terms of O(m 2 /Q 2 ). Due to the universal structure of collinear divergences the leading mass effects can be straightforwardly included based on the massless result. This massification procedure was developed in the context of Bhabha scattering [11][12][13] and was recently extended to processes with a heavy mass [34].
The matrix elements are implemented in the integrator McMule, a Monte Carlo for MUons and other LEptons [35]. This framework is based on the FKS ℓ subtraction scheme [36] which is an extension of the original FKS scheme [37,38] beyond NLO for QED. This subtraction scheme allows to consistently remove the singularities arising from soft photon emission in order to calculate observables in a fully differential way. The simplicity of the FKS ℓ subtraction scheme is due to the absence of collinear and the simple structure of soft singularities. All collinear divergences are regulated by finite fermion masses. Following the notation from [36] the soft singularities exponentiate according to the YFS formula [ All soft poles of the ℓ-loop matrix element (squared amplitude) with n final-state particles M with the scaled photon energy ξ = 2E γ / √ s and the eikonal factor E.
The most challenging part of the calculation presented here is to ensure a reliable integration of the real emission contributions in the phase-space region where the photon becomes collinear to the emitting fermion or where it becomes soft. In the former case the smallness of the electron mass acting as an infrared regulator results in large pseudo-collinear singularities. To address this issue we use a dedicated tuning of the phase-space parametrisation to help the vegas integration [40] find and deal with these problematic regions. In the latter case the integrand develops an unregularised soft singularity that is subtracted with the IR counterterm, resulting in a large cancellation. For this cancellation to work the matrix element has to be evaluated with very high precision. Due to the analytical and algebraic complexity of the real-virtual matrix element this is a highly non-trivial task and has presented the main obstacle to a complete, fully differential NNLO calculation of Bhabha scattering in the past. Our solution to this problem is the main result of this paper and is discussed in detail in the next section.

Real-virtual stabilisation via next-tosoft approximation
This section discusses how an implementation of the realvirtual matrix element can be obtained that ensures a stable and efficient integration in the soft phase-space region. As alluded to above, any general-purpose calculation of a one-loop matrix element will run into numerical instabilities at some point. In particular, for processes with an external photon with ever smaller energies, the IRsubtracted matrix element is a typical numerical pitfall whereby two expressions diverging as 1/ξ are combined to obtain an integrable integrand diverging as 1/ √ ξ. The crucial question is whether these instabilities appear only for small enough ξ such that the integration can be done reliably. In this context, QED calculations are particularly delicate since the final states tend to be much less inclusive than jet cross sections computed for hadronic collisions.
We use OpenLoops [33] for the bulk of the phase space since it shows a remarkable numerical stability. In order to test for which values of ξ the instabilities start to appear, we compare OpenLoops to a dedicated computation of the real-virtual matrix element in Mathematica using arbitrary precision arithmetic. In Figure 1 we show the deviation of OpenLoops from the 'exact' Mathematica result. For illustration we use an arbitrary phase-space point as well as one where the photon is emitted nearly collinear to the initial-state electron. For the former, at ξ = 10 −5 the relative error is 10 −8 . In the collinear case the numerical instabilities are strongly enhanced with a relative difference of 10 −1 for ξ = 10 −5 . All numbers that enter Figure 1 including the particle momenta p i of the two phase-space points are publicly available under [41].
An obvious idea is to expand the real-virtual matrix element for small photon energies and to switch to this approximation for sufficiently small ξ. The leading O(ξ −2 ) n . Using this approximation amounts to using the same algebraic expression for both terms in the subtracted integrand, albeit with different kinematics. As can be seen from Figure 1 this approach is insufficient in the collinear region. If an accuracy below 10 −3 is to be aimed at, in this case one has to switch to the expansion at ξ ∼ 10 −3 . However, the exact matrix element is not sufficiently well approximated by the leading soft contribution in this region. To ensure a decent approximation we have therefore to include the non-universal O(ξ −1 ) term in the soft expansion.
The next-to-soft terms of tree-level matrix elements have been considered a long time ago (Low-Burnett-Kroll theorem) [42,43]. Going beyond tree level, it is tempting to try to apply effective-field-theory methods. However, for QED with massive fermions the appropriate effective theory is the QED version of heavy quark effective theory and the genuine one-loop contribution to the next-tosoft effects is expected to be given by the process dependent soft function. From a practicable point of view we have thus decided to directly calculate the non-universal O(ξ −1 ) term in the soft expansion.
To be precise, we have computed the real-virtual matrix element in terms of scalar Passarino-Veltman functions using the Mathematica calculation of the real-virtual matrix element described in the previous section. Next, we have employed the power counting and expanded in the book-keeping parameter λ. The expansion of the rational coefficients and simple Passarino-Veltman functions was performed with Mathematica. More complicated triangle-and box-functions were expanded at the (loop-)integrand level using the method of regions [44]. For the purpose of calculational efficiency we have used its formulation in the parametric representation [45]. In this case, the contributing regions can be easily found using the public code asy.m [46]. Most resulting integrals could be straightforwardly computed. The remaining ones were calculated using Mellin-Barnes techniques [47,48]. For the most involved integrals a twofold Mellin-Barnes representation was necessary. However, the reduction to single contour integrals was possible in this case by resolving the singularity structure with the Mathematica package MBresolve.m [49]. In summary, the main technical difficulties in performing the next-tosoft expansion are the calculation of the integrals and the treatment of large intermediate expressions.
We have checked that the first term of the expansion indeed reproduces the result from (3). The non-universal subleading contribution was verified numerically. This is also shown in Figure 1 where the inclusion of the nextto-soft O(ξ −1 ) term significantly improves the approximation. This allows us to switch to a reliable expansion as early as ξ ∼ 10 −3 . We can therefore conclude that the next-to-soft approach ensures the numerical stability of the real-virtual matrix element for small photon energies which is a prerequisite for the IR subtraction to work. This is further emphasised by comparing integrated results with and without stabilisation. While the nextto-soft stabilisation ensures that results after successive Monte Carlo iterations are in agreement with each other, a drifting mean value is observed otherwise resulting in a significant discrepancy between the two results. Furthermore, the evaluation of the obtained expansion is a few 100 to over a 1000 times faster than OpenLoops, depending on the details of the kinematics. Since vegas tends to sample predominantly in the soft and collinear region this speed-up is noticeable even in the integration over the full phase space. While OpenLoops provides settings to work at higher accuracy this comes at a cost of speed.

Results and verification
To test the next-to-soft approach of stabilising the realvirtual contribution we have compared to BABAYAGA. For all results presented in this section we have switched from OpenLoops to the next-to-soft approximation at ξ = 10 −3 . As mentioned in the introduction, the event generator BABAYAGA is based on a parton shower algorithm matched to the exact NLO result. Contrary to our complete fixed-order calculation it therefore only gives the logarithmically enhanced contributions at NNLO. For the comparison we use set-up (a) of [28] that is tailored to φ factories with a centre-of-mass energy of √ s = 1020 MeV. The detector configuration is approximated with the kinematical cuts E min = 408 MeV, where E min is the minimum energy of the final-state electron/positron, θ − (θ + ) is the scattering angle in the centre-of-mass frame between the incoming and outgoing electron (positron), and ζ max is the maximally allowed The order-by-order contributions, σ (i) , to the integrated cross section, σ 2 = σ (0) + σ (1) + σ (2) , are presented in Table 1. Additionally, we show the corresponding K factors defined as To avoid comparing to contributions from the parton shower beyond NNLO we have not directly compared to [28] but instead were provided truncated results [50]. We find complete agreement for the LO as well as the NLO result. Our fixed-order NNLO correction σ (2) agrees at the level of 17% with the NNLO contribution from the matched parton shower. This translates to an agreement for the total cross section σ 2 of 0.07%, consistent with the 0.1% precision aimed at in [28]. We note that the sizable K factors are a consequence of the cut on the acollinearity that suppresses hard radiation. Figure 2 shows differential results with respect to the electron scattering angle θ − . The differential cross section at LO as well as at NNLO are displayed in the upper panel. In addition, the lower panel shows the differential K factors The comparison with the truncated parton shower at the differential level yields a similar result as for the total cross section, i.e. complete agreement up to NLO and deviations of below 0.1% at NNLO.

Conclusion and outlook
One of the main complications in the calculation of fully differential higher-order corrections in QED is the occurrence of numerical instabilities in real-emission contributions due to finite but small fermion masses acting as collinear regulators. In the case of Bhabha scattering instabilities arising in the real-virtual matrix element have represented the main bottleneck for a fixed-order Monte Carlo at NNLO accuracy.  In this paper we have presented a new method that ensures a stable and efficient integration of these problematic contributions. Since the main instabilities occur for soft photon emission we have expanded the real-virtual matrix element for small photon energies including the non-universal subleading contribution using the method of regions. This then allows for the reliable use of Open-Loops in the bulk of the phase space and to switch to the next-to-soft approximation otherwise. While an analogous approach could also be used for the real and the double-real matrix element, a stable and fast implementation can be obtained without the next-to-soft stabilisation.
We were therefore able to calculate for the first time the photonic NNLO corrections for Bhabha scattering in a fully differential way. This was implemented in the McMule framework. We have cross-checked our exact NNLO results at the level of the total cross section as well as for differential distributions with the logarithmic approximation implemented in the parton shower generator BABAYAGA.
Numerical instabilities of the kind described above are a critical point of higher-order QED calculations. We therefore expect that the next-to-soft method will prove useful in other processes as well. This is obvious in the case of Møller scattering which is related to the Bhabha process via crossing. Corresponding results relevant for the experiment PRad II [51] will be presented in a forthcoming paper [52]. Furthermore, in the context of muon-electron scattering our approach could turn out to be invaluable where a fully differential NNLO calculation is highly desirable [53] and therefore aimed at [54].
Because of the wide range of applicability of the nextto-soft expansion, an investigation of a potential universal structure would be desirable to allow for a more efficient calculation of the expansion. This could be done in the framework of heavy quark effective theory. Furthermore, a similar approach could be pursued in the collinear region. Switching to a leading collinear expansion could result in a significant speed-up of the phase-space integration. This would entail the calculation of the currently unknown one-loop splitting functions for massive fermions.

Acknowledgement
It is a great pleasure to thank C. Carloni Calame for engaging in a detailed comparison with BABAYAGA. We are also grateful to M. Zoller for providing the real-virtual matrix element in OpenLoops as well as for assisting us with its proper usage. Without the impressive numerical stability of OpenLoops this project would not have been possible. Our thanks are further extended to L. Dixon for sharing the massless two-loop matrix element in electronic form. Finally, we would also like to thank T. Becher for helpful discussions. PB and TE acknowledge support by the Swiss National Science Foundation (SNF) under contract 200021 178967. YU acknowledges partial support by a Forschungskredit of the University of Zurich under contract number FK-19-087 as well as by the UK Science and Technology Facilities Council (STFC) under grant ST/T001011/1. PB acknowledges support by the European Union's Horizon 2020 research and innovation program under the Marie Sk lodowska-Curie grant agreement No 701647.