Two-loop QED corrections to the scattering of four massive leptons

We study two-loop corrections to the scattering amplitude of four massive leptons in quantum electrodynamics. These amplitudes involve previously unknown elliptic Feynman integrals, which we compute analytically using the differential equation method. In doing so, we uncover the details of the elliptic geometry underlying this scattering amplitude and show how to exploit its properties to obtain compact, easy-to-evaluate series expansions that describe the scattering of four massive leptons in QED in the kinematical regions relevant for Bhabha and M{\o}ller scattering processes.

We study two-loop corrections to the scattering amplitude of four massive leptons in quantum electrodynamics.These amplitudes involve previously unknown elliptic Feynman integrals, which we compute analytically using the differential equation method.In doing so, we uncover the details of the elliptic geometry underlying this scattering amplitude and show how to exploit its properties to obtain compact, easy-to-evaluate series expansions that describe the scattering of four massive leptons in QED in the kinematical regions relevant for Bhabha and Møller scattering processes.
In recent years, particle physics has seen several interesting developments in experiments at the low-energy precision frontier.Among these are the discrepancy between theory predictions and the experimental value for the muon anomalous magnetic moment, most recently measured to 0.20 ppm at Fermilab [1], as well as the so called 'proton radius puzzle'.The latter consists in a discrepancy between the proton charge radius as determined in [2,3] compared to previous results [4].The upcoming PRad II experiment [5] will perform an independent measurement to attempt to resolve this inconsistency.This experimental program requires matching efforts on the theoretical side to provide equally precise and reliable predictions.An important part of these efforts is the recent development of the Monte Carlo event generator McMule [6].With the newly developed next-to-soft stabilization [7,8] for real-virtual corrections, McMule has the potential to describe Bhabha [9] (e + e − → e + e − ) and Møller scattering (e − e − → e − e − ) at the fully differential level up to next-to-next-to-leading order (NNLO) in Quantum Electrodynamics (QED).While Bhabha scattering is important for luminosity measurements at lepton colliders, Møller scattering is the main source of systematic uncertainty for the PRad II experiment [5] quoted above.Møller scattering is also important in searches for parity violation and for precise measurements of the weak mixing angle [10].Finally, precision measurements of Møller scattering at very low energies (2.5 MeV) [11] have also been undertaken recently.
The remaining outstanding ingredient to make theoretical studies in NNLO QED at arbitrary energy scales possible, are the two-loop virtual amplitudes for the scattering of four massive leptons, retaining full dependence on the lepton mass.The calculation of these virtual corrections has received much attention in the last decades.In QED with massless leptons, these amplitudes were computed more than two decades ago [12].Full event simulations at leading and next-to-leading order (NLO) [13][14][15], as well as power-suppressed mass effects up to NNLO have also been studied in detail [16][17][18][19][20][21], while fermionic loop corrections were computed with full mass dependence in [22].Moreover, also logarithmically enhanced electroweak contributions have been considered to NNLO [23][24][25][26].However, even though the computation of the relevant two-loop integrals was initiated already more than two decades ago [27][28][29][30][31][32][33][34][35], complete results for the two-loop virtual amplitudes where the full mass-dependence on the lepton mass is retained, are still not available, mostly due to the complexity of the integrals involved.
In this letter, we move an important step towards the exact inclusion of mass effects to Bhabha and Møller scattering up to NNLO in QED, by performing the first fully massive calculation of the two-loop QED corrections to the polarized and unpolarized scattering amplitude of four massive leptons.While retaining the full mass dependence renders the required two-loop integrals considerably more complicated, it will allow us to study the phenomenological impact of so-far neglected mass-effects in extreme regions of phase-space in upcoming phenomenological studies.In addition to their phenomenological relevance, these amplitudes also provide us with an invaluable playground to test recent developments in the theory of elliptic generalizations of multiple polylogarithms [36][37][38][39][40][41][42] and about the generalization of so-called canonical differential equations [43] to arbitrary geometries [44][45][46][47][48][49][50][51][52].

KINEMATICS AND TENSOR DECOMPOSITION
We work in QED with one single type of massive lepton, which for definiteness we refer to as the electron.We study higher-order corrections to the scattering of four electrons where all momenta are outgoing and satisfy the on-shell condition p 2 i = m 2 , i = 1 . . .4, as well as momentum conservation p µ 1 + p µ 2 + p µ 3 + p µ 4 = 0.The corresponding amplitude, A(1 e + , 2 e − , 3 e + , 4 e − ), can be parameterized as a function of the fermion mass m, and three Mandelstam invariants where, due to momentum conservation, s + t + u = 4m 2 .
Following [53,54], we work in 't Hooft-Veltman dimensional regularization scheme [55] (tHV) and decompose the scattering amplitude into eight independent Lorentzcovariant, physical tensors T i and respective scalar form factors F i , We choose the tensors as where the spinor chains t i are defined as i V e (p 1 ) × U e (p 3 ) Γ i V e (p 4 ) , (5) and i } represent the following sets of Dirac matrices In this computational scheme, external momenta and polarizations are considered four-dimensional, while internal states and loop momenta are treated in D dimensions.
One can then show that the number of tensors is equal to the number of independent chirality configurations to all orders in perturbation theory [53,54].In our case there are 2 4 = 16 configurations, of which only half are independent in a parity-invariant theory such as QED, which matches the eight tensors above.Furthermore, we note that the process under consideration is invariant under the simultaneous exchange p 2 ↔ p 3 and p 1 ↔ p 4 .We find that under this transformation two tensor structures are mapped onto each other, i.e. t 2 ↔ t 3 and t 6 ↔ t 7 , cf. (6), which in turn implies that T 7,8 are odd.Accordingly, by symmetry we conclude that the corresponding form factors must be zero F 7 = F 8 = 0 to all orders.Following the standard approach, we compute the form factors in (3), by defining a set of projection operators Here "•" denotes the scalar product between tensors and their dual, which is realized in practice by summation over spins of the external fermions, such that By applying the projectors on the corresponding relevant QED Feynman diagrams, we can express each form factor as a linear combination of scalar Feynman integrals and organize the one-and two-loop integrals into several integral topologies.On the technical level, our computation proceeds as follows.We begin by generating relevant Feynman diagrams with QGRAF [56].Using the computer algebra system FORM [57-60], we insert Feynman rules and apply projectors.We employ the public tool Reduze2 [61] to find mappings onto topologies and to expose their symmetries.Finally, with the help of Kira [62][63][64] we solve the required integration-by-parts (IBP) relations [65,66] and reduce all integrals to 267 master integrals.This is achieved following Laporta's algorithm [67], improved by finite field techniques [68,69].

CANONICAL BASES FOR THE NON-PLANAR FEYNMAN INTEGRALS
While all planar two-loop topologies have been known in fully analytic form for some time [34,35], their nonplanar counterparts have remained elusive due to the appearance of new mathematical functions of elliptic type.In particular, we are interested in the non-planar family of integrals displayed in the left graph of fig (1).We work in dimensional regularization and define the integrals as where γ E denotes the Euler-Mascheroni constant, D = 4−2ϵ is the dimension of space-time, and µ is an auxiliary scale introduced to render Feynman integrals dimensionless.The propagators are given by The integrals in (7) are functions of homogeneous coordinates [ s : t : m 2 ] on CP 2 and, without loss of generality, we may set µ = m, or equivalently work on the patch [ y : z : 1 ] with y = s/m 2 and z = t/m 2 .For definiteness, we will display our formulas in the region s > 4m 2 , t < 0 , though all results can also be easily continued to any other kinematic region.By solving IBP identities, all integrals can be expressed in terms of 52 independent master integrals.The latter fulfil a system of first-order partial-differential equations [70][71][72][73][74] in the kinematical invariants To solve this system, it is useful to search for a basis transformation to a so-called ϵ-factorized form: Such a system can be formally solved by a path-ordered exponential where the path γ connects the initial boundary point (y 0 , z 0 ) to a generic point (y, z).In the polylogarithmic case, if the matrix A can be expressed only through logarithmic differential forms, this matrix is said to be in canonical form, and the new integral candidates ⃗ J are called a canonical basis [43].While the generalization of a canonical basis beyond polylogarithms in not yet fully understood, advances have been made in extending ϵ-factorized bases to arbitrarily complicated geometries [45,46,[49][50][51]75].
For our problem, we achieved an ϵ-factorization by leveraging many of these developments.In particular, for the planar topologies, and for all polylogarithmic sub-sectors of the non-planar topology, we used unitarity cuts and multivariate residue analysis [76] to select integral candidates with unit leading singularities, see also [34,35].Starting from the six-propagator non-planar integrals generalizations of these methods to genus-one geometries become necessary.In fact, it is easy to show that the maximal cut of the irreducible six-propagator non-planar four-point graph (see right panel in fig. 1) in Baikov representation [77,78] can be expressed as .
By further taking the residue at z 2 = 0 in (12), one is left with an integral over a family of elliptic curves with the four roots given by We choose as first period for E 4 the integral where K(λ) is the complete elliptic integral of the first kind and its argument reads In order to arrive at an ϵ-factorized form, we first notice that all integrals corresponding to the graph of I 110111100 are reduced to six independent master integrals (plus subtopologies).We therefore expect two masters integrals which satisfy a coupled differential equation and map to the generators of the first de Rham cohomology group H 1 dR (E 4 ), plus four additional ones corresponding to independent punctures on the elliptic curve [51].Candidates for the first two masters can be found for example starting from the ansatz [79-81] where is the Wronskian of the second-order Picard-Fuchs equation associated to the elliptic curve.The explicit expression of J 46 is immaterial for this discussion, and is given in the supplemental material.The remaining four candidates can be identified by analysing their integrand representation and the structure of the resulting differential equations.As a last step, in order to obtain a fully ϵ-factorized form, one needs to integrate out some inhomogenous entries in the differential equation matrix, which leads to the appearance of additional transcendental integrals.In this way, the final ϵ-factorized system (10), is expressed in terms of 87 distinct one-forms ω i .It is easy to verify that the integrability condition dA = A ∧ A is satisfied and that all ω i are the closed dω i = 0.
The individual differential forms can be simplified by exploiting the underlying geometry of the family of elliptic curves in (13).As an example, consider the following two functions which are among the objects required to express the matrix A in (10).Again, formulas are given assuming y > 4 and z < 0 for definiteness.While the details of the construction are immaterial for this paper and are discussed elsewhere [82,83], it suffices to say that one can parameterize the kinematical variables by where {[x : Y : 1], t 4 } are the canonical coordinates on the moduli space of elliptic curves given by the variety Y 2 = (x 2 − 1)(x 2 − t 4 ).In these coordinates, the period in (15) becomes We can identifying t 4 with a Hauptmodul for the congruence subgroup Γ 1 ⊂ SL 2 (Z) [84].Strikingly, it turns out that by changing variables to the canonical coordinates, one can easily see that the two transcendental integrals in (18) are just combinations of simpler functions where f (t 4 ) is given by and F(x, t 4 ) is the derivative of the Abel map: Other differential forms in the alphabet can also be substantially simplified and all double integrals over the periods can be rewritten in terms of rational functions of T 1 and T 2 .One can then show that all differential forms are given by combinations of the five algebraic functions , and the three transcendental functions {K(t 4 ) , f (t 4 ) , F(x, t 4 )}.We want to stress that the choice of canonical coordinates in (19) is not merely an academic curiosity, and the final, simplified form is essential to efficiently implement the numerical evaluation of the iterated integrals described below.To explicitly solve the integrals, we first expand (11) in ϵ.At each order, the solution of the differential equation is expressed by Chen iterated integrals [85] and we fix all boundary conditions imposing regularities at different phase-space points.In this way we obtain fully analytic results for the non-planar master integrals in terms of Chen iterated integrals.
Currently, there are no public numerical routines to evaluate the special functions that appear in the nonplanar sector.We therefore obtain generalized series expanions for all master integrals.More precisely, we start from the differential equations in ϵ-factorized form in order to algorithmically obtain a small mass expansion for the individual master integrals.In particular, we obtain a generalized power series (including logarithms of the mass), whose coefficients can be expressed in terms of harmonic polylogarithms [86].We obtain results that are valid both for the kinematics relevant for Bhabha and Møller scattering.As a cross check, we compared individual master integrals against a direct numerical evaluation with AMFlow [87], both for Bhabha and Møller scattering kinematics, and found agreement to high precision.Our series expansions allow for fast numerical evaluation, appropriate for phenomenological studies.A precise description of the numerical implementations can be found in the description of the ancillary files along with the arXiv submission of this manuscript.

UV RENORMALIZATION AND IR FACTORIZATION
Using the master integrals calculated above, as well as the planar integrals from [34,35], we can obtain an analytic result for the bare amplitude for both polarized and unpolarized scattering.The UV divergences can then be renormalized according to with the relation between bare and physical quantities Here Z 2 and Z m are on-shell wave function and mass renormalization constants, and Z e refers to coupling constant renormalization either in the MS or on-shell (OS) scheme.The relevant quantities are collected in the supplemental material.As expected [88], after UV renormal-ization we are left with IR poles which are one-loop-exact, where C is the finite remainder function, α is the on-shell electromagnetic coupling, and Z IR 1 is the anomalous dimension which controls the soft singularities of the amplitude to all-orders through exponentiation [89,90].The exact form of Z IR is immaterial for the present discussion and we report it for completeness in the supplemental material.
We performed several checks on our results.First of all, we verified that our two-loop amplitudes have the correct UV and IR behavior, as illustrated above.In addition, we compared both the bare and the finite remainders of our one-loop amplitudes against OpenLoops [91,92] and found perfect agreement.We stress here that the unpolarized finite remainders in Conventional Dimensional Regularization equal those in the tHV scheme, while the bare and UV-renormalized amplitudes in general differ.The equality of the finite remainders provides another check of our calculation.

DISCUSSION AND CONCLUSIONS
Our results for the two-loop amplitudes for Bhabha and Møller scattering are given as generalized series expansion in x = m/E CM .They are provided as computerreadable files in the ancillary material of the arXiv submission for both the polarized and unpolarized scattering amplitudes.We provide sufficiently high orders to obtain reliable predictions for the low-energy experiments mentioned in the introduction, where we expect the mass effects to be the largest.In the following we discuss some of the phenomenological implications of our results.We focus here on unpolarized Møller scattering, but all conclusions equally apply to Bhabha scattering.
Let us start by assessing the accuracy of the smallmass expansion.We begin by noticing that we expect the expansion to become unreliable in the extreme forward or backward regions, where the coefficients of the series in x develop large logarithms in (−t)/s which can invalidate the convergence of the expansion. 1To quantify the region of convergence, we compare the exact results for the oneloop amplitude A 1l exact with the corresponding expansion A 1l 20 to O(x 20 ) and study the ratio δ 1l exact,20 = (A 1l exact − A 1l 20 )/A 1l exact .Depending on the scattering energy E CM = √ s, we find that δ 1l exact,20 ≤ 1% for different ranges of the scattering angle θ: where the energy values are chosen to match those probed at present and future experiments.This shows that at very low energies the expansions must be interpreted with care outside of the central region.To extend this to the two-loop amplitudes, we repeat the same analysis at one and two loops, comparing this time the series expanded to order 20 with the one expanded to order 18.We find that the same applies: for L = 1, 2 (A Ll 20 − A Ll 18 )/A Ll 20 ≤ 1% for the same values of θ as in (27).In fig. 2 we display the various orders of the series for the two-loop amplitude, for different values of the scattering at the intermediate energy of E CM = 32m.We highlight the lack of convergence for θ not in the range [9  x 0 x 4 x 8 x 12 x 16 x 20 FIG. 2: Convergence of the mass expansion.Plotted are the 2-loop finite remainders C †(2) C (0) as functions of scattering angle in degrees, at various truncation orders.
After having assessed the validity of our small-mass expansions, let us comment on the phenomenological relevance of the mass effects.We only discuss here the mass effects in the purely virtual corrections.So far two-loop mass effects had only been included to leading-power, O(x 0 ).We expect that the finite-mass effects are more pronounced for small values of E CM .In fig. 2 we see that, for E CM = 32m, the two-loop leading-power approximation does not capture the full extend of the mass effects for θ ≳ 150 • (for small angles, we are outside the region of ( 27)).We therefore expect that in that region precise NNLO results can only be obtained by including the subleading terms we have computed.The effect is even more pronounced for E CM = 5m: in fig. 3 we show that, even in the range of intermediate angles in (27), the leading-power approximation does not provide a reliable prediction of the finite-mass effects.At the same time, we observe a very nice convergence of the mass expansion, corroborating that we can provide reliable and precise predictions for the two-loop corrections even at such low energies.A full discussion of the size of the NNLO QED corrections will be presented elsewhere.
To conclude, in this letter we have addressed the calculation of the two-loop QED corrections to the scattering of four identical massive leptons, retaining full dependence on the lepton mass.This constitutes the last outstanding ingredient necessary to perform NNLO QED phenomenological studies for standard processes as Bhabha and Møller scattering.In addition to the phenomenological interest behind these calculations, the scattering amplitudes computed in this paper are an important example of physical processes that receive a nontrivial contribution from Feynman integrals of elliptic type.We presented a strategy to compute these amplitudes analytically through the differential equations method and provided a robust numerical implementation.We demonstrated that for low values of E CM , the mass effect can be sizeable and is not captured by the leading-power approximation.We therefore expect that our results will play an important role in making precise predictions for lepton collider experiments possible.
Introducing for convenience Ψ 1 = r t Ψ 0 , the ϵ-factorized basis for the sector corresponding to the right-hand graph of fig. 1   Finally, we also provide the ϵ-factorized basis for the top sector, given in the left graph of fig. 1.The three integrals read J 50 = r s r u I 1111111−10 , arXiv:2311.06385v1 [hep-ph] 10 Nov 2023

4 FIG. 1 :
FIG.1:The non-planar topology (left) and its next-to-top sector (right) with 6 master integrals.Solid lines correspond to massive propagators of mass m, dashed lines correspond to massless propagators.