Two-Loop Four-Fermion Scattering Amplitude in QED

We present the analytic evaluation of the two-loop corrections to the amplitude for the scattering of four fermions in Quantum Electrodynamics, $f^- + f^+ + F^- + F^+ \to 0$, with $f$ and $F$ representing a massless and a massive lepton, respectively. Dimensional regularization is employed to evaluate the loop integrals. Ultraviolet divergences are removed by renormalizing the coupling constant in the ${\overline{\text{MS}}}$-scheme, and the lepton mass as well as the external fields in the on-shell scheme. The analytic result for the renormalized amplitude is expressed as Laurent series around $d=4$ space-time dimensions, and contains Generalized Polylogarithms with up to weight four. The structure of the residual infrared divergences of the virtual amplitude is in agreement with the prediction of the Soft Collinear Effective Theory. Our analytic results are an essential ingredient for the computation of the scattering cross section for massive fermion-pair production in massless fermion-pair annihilation, i.e. $f^- f^+ \to F^- F^+$, and crossing related processes such as the elastic scattering $f F \to f F$, with up to Next-to-Next to Leading Order accuracy.

Introduction -The Muon g-2 collaboration at Fermilab has recently confirmed [1] that the observed magnetic activity of the muon is compatible with the earlier findings obtained at Brookhaven National Lab [2][3][4]. The anomalous magnetic moment of the muon, (g − 2) µ , shows a 4.2σ deviation from the prediction of the Standard Model of elementary particles (SM) [5]. However, the theoretical determination of this quantity, obtained via dispersive techniques, might be affected by the improper estimation of the hadronic corrections to the muon-photon interaction, which could be responsible of such a discrepancy. Alternative results obtained through lattice QCD calculations point towards a possible mitigation of the tension between theory and experiments [6].
Recently, a novel experiment, MUonE, has been proposed at CERN, with the goal of measuring the running of the effective electromagnetic coupling at low momentum transfer in the space-like region [7]. As proposed in [8], this measurement would provide an independent determination of the leading hadronic contribution to the (g − 2) µ . Such a measurement relies on the precise determination of the angles of the outgoing particles emerging from the elastic muon-electron scattering [7,[9][10][11]. To extract the running of the effective electromagnetic coupling from the experimental data, the pure perturbative electromagnetic contribution to the electron-muon cross section must be controlled at least up to the second order in the fine-structure constant [12].
The scattering of a muon µ off an electron e in Quantum Electrodynamics (QED) is the simplest reaction among fundamental leptons of different flavors, and represents a paradigmatic case of charged particles interaction mediated by a neutral gauge boson. The Leading Order (LO) process is known since the mid 1950's [13], while the Nextto-Leading Order (NLO) radiative corrections were computed in [14][15][16][17][18][19][20], and more recently studied in [21]. The two-loop diagrams contributing to the Next-to-Next-toleading order (NNLO) virtual corrections were evaluated in [22] assuming purely massless fermions. At the energies of the MUonE experiment, the muon mass plays an important role for the description of the radiative pattern and cannot be neglected [12]. Nevertheless, the evaluation of Feynman integrals usually becomes more demanding as the number of massive particles present either in the loops or in the external states increases.
NNLO QED corrections involve the two-loop amplitude along with the real-virtual and the double-real emission terms. While the matrix elements for the last two contributions can be calculated without difficulties using standard techniques, their integration over the corresponding phase spaces is complicated by the presence of infrared (soft and collinear) singularities, as well as the presence of masses in both the initial and final state of the scattering process. In order to obtain predictions for fully differential observables, it is necessary to adopt a subtraction procedure. The Abelian nature of the interaction leads us to believe that the computational techniques already used for other processes at the LHC can be successfully adapted to this purpose [23,24]. Preliminary Monte Carlo simulations for µe scattering have already been performed by arXiv:2106.13179v2 [hep-ph] 7 Apr 2022 including parts of the NNLO corrections [25,26]. These simulations account for a subset of the two-loop graphs, not yet including the four-point diagrams with complete dependence on the lepton masses. The complete two-loop amplitude is then a missing crucial ingredient for the computation of the full NNLO QED corrections.
In this work, we present the first fully analytic evaluation of the renormalized two-loop amplitude for four fermion scattering in QED, f − + f + + F − + F + → 0, with f and F representing a massless and a massive lepton respectively. In the past years, we have developed efficient mathematical techniques for the evaluation of multi-loop integrals in dimensional regularization, such as the adaptive integrand decomposition [27][28][29] and the Magnus exponential method for differential equations [30,31]. The combination of these techniques with the more traditional decomposition through integration-by-parts identities (IBPs) [32,33], allowed us to obtain for the first time a complete analytic formula for the renormalized two-loop amplitude of a 2 → 2 process with a non-vanishing mass in internal and external lines.
The one-and two-loop amplitudes presented in this Letter can be applied, for instance, to the case where the light fermion is an electron, f = e, and the heavy fermion is a muon, F = µ, and can be used in the elastic scattering eµ → eµ, as well as in crossing related processes, such e + e − → µ + µ − . If the elastic scattering is the key process of the MUonE experiment, the muon pair production in e + e − annihilation is a key process for the center-of-mass energy calibration at present and future e + e − colliders, such as BESIII [34], BELLE-II [35], CEPC [36], and FCCee [37]. Therefore, a precise knowledge of the radiative effects would improve the precision of the results obtainable at these machines.
The structure of the infrared (IR) singularities of the massless and massive gauge theory scattering amplitudes has been studied in [38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53]. In this work, the determination of the virtual NNLO corrections is complemented by the investigation of the IR singularities of scattering amplitudes in QED, which involve massive particles, and whose universal structure can be determined within Soft Collinear Effective Theory (SCET), along the lines of the method presented in [46,53]. The agreement of the residual IR poles coming from the direct diagrammatic calculus of the renormalized amplitude with the IR poles predicted within SCET is an important validation of the diagrammatic calculation. We explicitly verify this agreement in the case of f − f + → F − F + process.
Additionally, let us observe that the two-loop diagrams considered here, also appear in the (color stripped) Abelian subset of graphs contributing to heavy-quark pair production in Quantum Chromodynamics (QCD) [54][55][56][57][58]. The similarities of the analytic structure of the two-loop amplitude between qq → tt in QCD and f − f + → F − F + in QED, where q and f are treated as massless, is exploited to test the structure of the singularities coming from QED diagrams through a tuned comparison to the Abelian part of known results in QCD.
Recently, the evaluation of integrals coming from planar diagrams [59][60][61] indicates that the computation of four-fermion scattering amplitudes at two loops in QED, by keeping full dependence on the masses of all the involved leptons, might become the subject of near-future investigation.
Scattering Amplitude -We consider the four-fermion scattering process involving a mass-less and a massive lepton pair, with m f = 0 and m F = M = 0. The Mandelstam invariants, defined as s = ( The four-point bare amplitude A b admits a perturbative expansion in the bare coupling constant α b ≡ e 2 b /4π, which, up to the inclusion of the second-order corrections, reads where indicates the n-loop bare amplitude, S ≡ (4πe −γ E ) and µ is the 't Hooft mass scale. The Leading Order (LO) term A (0) b , referred to as Born term, receives contribution from a single tree-level Feynman diagram, shown in the upper row of Fig. 1. The squared LO amplitude, summed over the final spins and averaged over the initial states, reads, for external states treated in d = 4 − 2 space-time dimensions according to the conventional dimensional regularization (CDR) scheme [62], that we use throughout the whole computation. The interferences of one-and two-loop bare amplitudes with the Born amplitude read Analytic Evaluation -The analytic evaluation of M and M (2) b is completely automated, within an in-house software, which can be applied to generic one-and twoloop amplitudes. The Mathematica package Fey-nArts [63] is used to generate Feynman diagrams contributing to the one-and two-loop corrections to the scattering amplitudes as well as the counter-term diagrams required for the renormalization: 6 diagrams and 3 counter-term diagrams at one loop; 69 diagrams (12 of which vanish because of Furry's theorem) and 55 counterterm diagrams at two loops. Representative one-and two-loop diagrams are shown in the second and third row of Fig. 1, respectively. The spin sums and the algebraic manipulation to simplify the Dirac-γ algebra are carried out by means of the FeynCalc [64][65][66] package. Each n-loop graph G (interfered with the Born amplitude) corresponds to an integrand written in terms of scalar products between external, p ν i , and internal, k ν i , momenta. Therefore, Eq.(4) can be generically written as, where: N G = N G (p i , k i ) indicates the numerator, and D σ = D σ (p i , k i , M ) are the denominators corresponding to the internal lines of G.
Integrands are simplified by employing the adaptive integrand decomposition method, implemented in the Aida framework [29]. The intermediate results emerging from the integrand decomposition can be further simplified by means of the IBP identities [32,33]. Our software is interfaced with the publicly available codes Reduze [67] and Kira [68], and, for each diagram, it produces the files for the automated generation of the IBP relations. After the decomposition phase, the interference terms M are written as linear combination of a set of independent integrals, say I (n) , called master integrals (MIs), where C (n) is a vector of coefficients, depending on and the kinematic variables, s, t, M 2 . In particular, M b and M (2) b are conveniently expressed, in terms of 12 and 264 MIs, respectively, analytically computed: two-and threepoint functions have been known since long [69][70][71], while planar and non-planar four-point integrals were computed in [72,73], using the differential equation method via Magnus exponential, and independently in [55,56,74]. The analytic expressions of M can be written as a Laurent series around d = 4 space-time dimensions ( = 0), with coefficients that contain Generalized Polylogarithms (GPLs) [75], defined as iterated integrals, through the recursive formula with G(w 1 ; t) ≡ log(1 − t/w 1 ). The arguments w i are known as letters, and their number, corresponding to the number of nested integrations, is called weight. The two-loop interference term contains 4063 GPLs with up to weight four, whose arguments are written in terms of 18 letters, w i = w i (x, y, z), which depend on the Mandelstam variables through the relations, [72,73] for more details).
Renormalization -The one-and two-loop diagrams contributing to M b contain infrared (IR) and ultraviolet (UV) divergences. To remove the UV divergences, the bare lepton fields (ψ , with = f, F , for massless and massive leptons, respectively) and photon field (A µ ), as well as the bare mass of the massive lepton are renormalized as follows, where, to simplify the notation, the label in the lepton fields is understood and restored when required. The renormalization of the QED interaction vertex, can then be entirely fixed using the QED Ward identity, that implies Z 1 = Z 2 . In particular, this leads to a simple relation between the renormalized charge and the bare charge (obtained by applying Eq. (8) to the bare interaction term and comparing the two renormalized expressions) e Z 1 = e b Z 2 √ Z 3 , therefore, one has e = e b √ Z 3 . The lepton wave functions and the mass of the massive lepton are renormalized in the on-shell scheme, namely, The coupling constant is renormalized in the MS scheme at the scale µ 2 , with Z MS α = 1/Z MS 3 . The renormalized amplitude is obtained by multiplying the bare amplitude with a factor Z 2, for any external lepton , hence, following, these are simply indicated as Z j , with j = {α, f, F, M }, respectively. The renormalization constants admit a perturbative expansions in α, and their expressions can be obtained (either directly or after abelianization) from [57,[76][77][78]. After substituting in Eq. (11) the expansions of the bare amplitude, given in Eq. (2), and the ones of the renormalization constants, given in Eq. (12), the UV renormalized two-loop amplitude reads up to second order corrections in α. The n-loop coefficients A (n) are given in terms of the ones appearing in the bare amplitude as The last term in Eq. (14c) contains the extra contribution of one-loop diagrams having an insertion of the mass counter-term in the massive propagators in all possible ways, as depicted in Fig. 2. The bare coupling α b and the bare amplitudes A (n = 0, 1, 2), appearing in Eqs. (42) and (4), can be replaced by the corresponding renormalized quantities α and A (n) , to build the Born term, M (0) , and the renormalized interference terms, at one loop, M (1) , and at two loops, M (2) . The latter two quantities constitute the main results of this Letter.
Infrared Structure -The IR poles appearing in the twoloop corrections after UV renormalization can independently be obtained starting from the tree-level and the one-loop amplitudes, by following the same procedure employed to study the infrared structure of QCD amplitudes [46,53]. The structure of the IR poles is governed by an anomalous dimension Γ that has the following structure, where the γ i (i ∈ {cusp; cusp,M; h; ψ}) coefficients up to O(α 2 ) are extracted in analogy to the QCD case [46,53,79]. We compute the analytic expression of the two-loop amplitude M (2) for the process f − f + → F − F + both in the non-physical region s < 0, t < 0 as well as directly in the production region. In this physical region, the imaginary part of the anomalous dimension in Eq. (15) is computed by adding an infinitesimal positive imaginary part to s. One can then introduce the IR renormalization factor Z IR , where Γ i , Γ i and β i are the coefficients of the expansion of Γ, its derivative w.r.t. ln µ, and the QED beta function, respectively. The IR poles of the n th -order term M (n) can be calculated using Z IR and the lower order contributions, M (0) , . . . , M (n−1) . In particular, the IR pole structures at one and two loops are found to be, The IR poles structure in Eqs. (17), reconstructed starting from the tree-level and one-loop amplitudes, is in perfect agreement with the one obtained starting from Eq. (14c) and directly calculating the two-loop diagrams. This provides a non trivial test of the complete two-loop calculation.
Results -The analytic results of the interference contributions M (1) and M (2) are given as Laurent series in The analytical expression of M (1) is computed both in the non-physical region, and in the pair production region, s > 4M 2 , t < 0. The latter is required to predict the Here IR poles of M (2) directly in the production region; the analytical expression of M (2) is computed in the nonphysical region, s < 0, t < 0, and its analytic continuation is performed numerically. The renormalized one-and twoloop interference terms are conveniently decomposed in gauge-invariant components, labeled by the number of massless (n l ) and massive (n h ) closed fermion loops In Fig. 3, we plot the finite part of one-and two-loop renormalized amplitudes M (i) 0 , i = 1, 2 in the physical region. The threshold singularity is clearly visible and well reproduced up to very small c.m.e., showing full control of the numerical stability. The complete formula for the analytic expression of the renormalized two-loop amplitude is rather large (∼ 60 MB) and cannot be reported here. The figures are obtained by evaluating this formula with high precision on 10,500 evenly spaced grid points, by employing HandyG [80] and Ginac [81] (via the package PolyLogtools [82]) for the numerical evaluation of GPLs. Each evaluation required from seconds CPU time in the almost flat region to up about 1,500 s CPU time for the configurations approaching the threshold singularity. These grids are available from the authors upon request.
Other tests -The master integrals for the Abelian diagrams in QED can be employed to construct the analytic expressions of some gauge-invariant contributions to the two-loop amplitude of the process qq → tt in QCD [54][55][56][57]: in particular, our results (evaluated in the region of heavy-lepton pair production, and properly accounting for the color factors) agree with the numerical coefficients E q l , E q h , F q l , F q lh , F q h provided in the Table 1 of Ref. [54,55,57], which receive contributions from Abelian diagrams only; the agreement on the poles of the above mentioned color coefficients, at other phasespace points, has been verified using the formula for the IR poles of two-loop amplitudes in QCD, given in Ref. [83].
Conclusion -We presented the first fully analytic evaluation of the amplitude for the scattering of four fermions in Quantum Electrodynamics, involving two different types of leptons, one of which is treated as massless, up to the second order corrections in the electromagnetic coupling constant. The calculations were carried out within the dimensional regularization scheme, and the infrared pole structure of the renormalized amplitude is found to obey the universal behaviour predicted by the Soft Collinear Effective Theory. Our result constitutes the first example of a complete scattering amplitude for 2 → 2 processes, with massless and massive particles in the loops as well as in the external states, involving planar and non-planar diagrams at two loops, analytically evaluated.
Our analytic results can be directly applied to the study, at NNLO accuracy, of massive lepton pair production in massless lepton annihilation, and, upon analytic continuation, to the study of the elastic scattering of massive and massless fermions in QED and QCD. Notably, the virtual corrections presented here are relevant for the recently proposed MUonE experiment at CERN. This experiment is devoted to extraction of the hadronic contribution to the (g − 2) µ from the µe scattering.
The MUonE experiment analysis relies on the knowledge of the pure NNLO QED correction to the µe scattering process, which will be the subject of a dedicated study in the near future.

SUPPLEMENTAL MATERIAL
In this supplemental material, we provide further details on the renormalization constants to perform the UV renormalization, and the IR renormalization factor for the predictions of the IR poles of the oneand two-loop four-fermion scattering amplitude in QED, f − + f + + F − + F + → 0, with f and F , representing massless and a massive leptons, respectively.
Renormalization Constants -The renormalization constants for the wave functions of the massive and massless leptons as well as the mass renormalization constant of the massive lepton admit perturbative expansions in α, which can be taken (either directly or after abelianization) from Refs. [57,[76][77][78], and read as, where the individual perturbative coefficients in the onshell scheme read as,

δZ
(1) δZ (2) δZ (1) with L µ ≡ ln µ 2 /M 2 . Additionally, the renormalization constant for the electromagnetic coupling up to second order in the MS scheme reads, Anomalous dimensions -The structure of the IR poles is governed by an anomalous dimension Γ whose structure reads as, and by the IR renormalization factor Z IR , defined by exponentiation of the following expression, where Γ i , Γ i and β i are the coefficient of the expansion of Γ, its derivative w.r.t. ln µ, and the QED beta function, respectively. We hereby present the coefficients appearing in the above formulas, up to the needed order in α: The cusp anomalous dimensions for massless leptons have the coefficients [79], whereas, for massive leptons, withx defined through the relation, The coefficients of the factors related to the massless and massive leptons are [79], The QED beta function has the expansion in which the only needed coefficient for the present calculation is β 0 , The quantity Γ appearing in Eq. (29) is defined as, with the relevant coefficients, Lastly, the expansion of the renormalization factor Z IR in Eq. (29) is decomposed as, with, Furthermore, in order to implement the inverse decoupling transformation for the massive leptons, in such a way that one works with n l + n h active leptons, one needs to include an additional term proportional to n h in Z IR 2 : h .
In Figs 1 and 2, we plot the finite parts of the individual form factors appearing in the decomposition of the one-and two-loop amplitudes, given in Eqs. (44). Finally, in Table I, we showcase the numerical values of the coefficients A, B, C, D, E, F in the massive-fermion pair production region at a particular phase-space point. Note that, upon accounting for a different definition of the Mandelstam variables, this benchmark point corresponds to the one used in Table 1 of [57].