Lepton-pair scattering with an off-shell and an on-shell photon at two loops in massless QED

We compute the two-loop QED helicity amplitudes for the scattering of a lepton pair with an off-shell and an on-shell photon, 0 → ℓℓ¯γγ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \ell \overline{\ell}\gamma \gamma $$\end{document}*, using the approximation of massless leptons. We express all master integrals relevant for the scattering of four massless particles with a single external off-shell leg up to two loops in a basis of algebraically independent multiple polylogarithm, which guarantees an efficient numerical evaluation and compact analytic representations of the amplitudes. Analytic forms of the amplitudes are reconstructed from numerical evaluations over finite fields. Our results complete the amplitude-level ingredients contributing to the N3LO predictions of electron-muon scattering eμ → eμ, which are required to meet the precision goal of the future MUonE experiment.

1 Introduction The MUonE experiment [1][2][3][4] will measure the hadronic running of the electromagnetic coupling α using low-energy elastic electron-muon scattering, eµ → eµ.This will enable a new and precise determination of the hadronic vacuum polarisation (HVP) contribution a HVP µ [5,6] to the muon anomalous magnetic moment a µ .This is required in light of the recent tensions between experimental [7], Standard Model (SM) data-driven [8], and lattice quantum chromodynamics (QCD) [9] results for a µ .Increasing the precision of the theoretical predictions for eµ → eµ scattering is a high priority for the planned MUonE experiment [10,11] and has seen good progress in the last few years [12][13][14][15][16].The recent completion of full next-to-next-to-leading order (NNLO) quantum electrodynamics (QED) corrections [17] indicates that next-to-next-to-next-to-leading order (N 3 LO) corrections in differential distributions are required to meet MUonE's precision goal of 10 parts per million.Electron-line corrections, meaning corrections to the subprocess with the muon line stripped off (e → eγ * ), are the dominant corrections [17], and a collaborative project was started to perform their fixed-order calculation at N 3 LO [18].With the triple-virtual corrections now available [19][20][21], the main missing ingredient is the real-double-virtual (RVV) matrix element (e → eγγ * ) at two loops.While these contributions could be extracted from amplitudes in the literature [22][23][24], our direct computation provides the massless RVV contribution in a complete and compact form.
Another application of the 0 → ℓ lγγ * amplitudes is in electron-positron annihilation experiments [25].They are required for initial-state corrections in predictions of the ratio of hadron-to-muon production in e + e − collisions, which is an important input for existing SM predictions of a HVP µ [26].The two-loop amplitudes contribute to RVV corrections to e + e + → γ * in direct scan measurements, while radiative return measurements concern corrections to e + e − → γγ * [8].In the latter configuration, the e + e − beam has a fixed centre-of-mass energy of a few GeV and the on-shell photon originates from initial state radiation (ISR).The energy lost to the ISR photon is used to effectively scan over the energies of the decay of the off-shell photon.A differential cross section of, for example, γ * → hadrons with respect to the centre-of-mass energy of the decay, dσ/ds, can be extracted from measurements of the differential cross section with respect to the energy of the ISR photon, dσ/dE γ .State-of-theart predictions for these measurements are currently at next-to-leading order (NLO) [26].We provide the two-loop e + e − → γγ * amplitudes required for the double-virtual (VV) corrections at NNLO, although the bottleneck remains in the hadronic decay.
Our amplitudes are calculated in the approximation of massless leptons.In the NNLO massive eµ → eµ cross section calculation [17], the authors obtain photonic corrections (those with no closed fermion loops), using a small-mass expansion [27][28][29] applied to the two-loop amplitudes with massless electrons for the VV corrections.This approximation relies on the electron mass being much smaller than any other scale, which is valid in the bulk of phase space.Further splitting the photonic corrections, they take the subset of electronline corrections and find that the relative difference to the true massive NNLO differential cross section is generally around 10 −3 α 2 , where α is the fine-structure constant, which is negligible compared to the 10 −5 precision goal.The approximation breaks down in soft and collinear regions, where they treat the amplitudes using infrared (IR) factorisation [30][31][32], and is not used for contributions including closed fermion loops [29,33].Our amplitudes can be used analogously for the RVV corrections at N 3 LO.
Our computation uses the modern technology developed for QCD amplitudes with many scales.The high-multiplicity amplitude frontier in massless QCD lies with twoloop five-particle processes, with leading-colour [34][35][36][37][38][39] and full-colour [40][41][42][43] results in a form ready for phenomenological application becoming available over the past few years.Recently, the first single-external-mass calculations are also appearing [44][45][46][47].These computations have made extensive use of finite-field arithmetic to sidestep large intermediate expressions.This technology has had a considerable impact for solutions of systems of integration-by-parts (IBP) identities [48][49][50] but also applies more widely to scattering amplitude computations [51,52].Motivated by the improved algorithms, we choose to implement a complete finite-field based reduction for the 2 → 2 processes with an off-shell leg.Since the kinematics are relatively simple in comparison with other high-multiplicity configurations, this technology is not essential.It does, however, provide an opportunity to review the new techniques for readers who are not familiar with them.
A key ingredient for computing the scattering amplitudes are analytic expressions for the required Feynman integrals.Complete analytic results up to two loops are already available in the literature [53][54][55][56][57]. Expansions of these integrals up to higher orders in the dimensional regularisation parameter ϵ have also been reconsidered recently [57], in view of their usage for N 3 LO corrections to 2 → 2 processes in QCD [58,59].The state of the art for integrals with this kinematic configuration has reached three loops [60][61][62][63].We revisit the computation of the one-and two-loop integrals following the approach of refs.[45,[64][65][66][67] based on the construction of a basis of independent special functions, which gives a unique and uniform representation of all the required Feynman integrals up to transcendental weight four.This enables a more efficient computation of the amplitudes using the modern workflow based on finite-field arithmetic, and leads to more compact expressions.We give explicit expressions for the basis functions in terms of multiple polylogarithms (MPLs) which can be evaluated in an efficient and stable way throughout the physical phase space.We compute all crossings of all massless one-and two-loop four-particle Feynman integrals with an external off-shell leg, so that our results for the integrals may be of use for any scattering process with these kinematics.
Our paper is organised as follows.In section 2, we describe our decomposition of the helicity amplitudes and detail how we express the off-shell currents.In section 3, we discuss our computation of analytic amplitudes by numerical evaluations over finite fields.In section 4, we present the computation of the Feynman integrals in terms of a basis of special functions.We draw our conclusions in section 5. We provide useful technical details in appendices.We define the relevant families of Feynman integrals in appendix A. In appendix B, we discuss in detail how we handle permutations of the integral families in the IBP reduction.In appendix C, we describe our rational parametrisation of the kinematics.Appendix D is devoted to the ultraviolet (UV) renormalisation and IR factorisation which determine the pole structure of the amplitudes.In appendix E we discuss the analytic continuation of the special functions to the physical kinematic regions.

Structure of the amplitude
We calculate the one-and two-loop QED corrections to the process which we call 0 → ℓ lγγ * for short.Here, ℓ denotes an on-shell massless lepton and γ (γ * ) an on-shell (off-shell) photon, while h i and p i are the helicity and momentum of the i th particle.We take the external momenta p i to be all outgoing.They satisfy the following momentum-conservation and on-shell conditions: 2) The single-off-shell four-particle phase space is described by three independent scalar invariants, which we choose as where s i...j := (p i + . . .+ p j ) 2 .We use dimensional regularisation in the 't Hooft-Veltman scheme [68], with D = 4 − 2ϵ spacetime dimensions (where ϵ is the dimensional regulator) and four-dimensional external momenta.
Because of the off-shell photon in the process, the helicity amplitudes A µ (1 ℓ , 2l, 3 γ , 4 γ * ) are actually off-shell currents carrying a free Lorentz index.We consider the perturbative QED expansion of the helicity amplitudes, with prefactor n ϵ = i(4π) ϵ e −ϵγ E , electromagnetic coupling g e , and α = g 2 e /(4π).We truncate the expansion at L = 2 loops.We set the renormalisation scale µ R to one throughout the computation and restore the dependence on it in the final analytic result by dimensional analysis.For the bare amplitudes we have that There are two independent helicity configurations (h 1 , h 2 , h 3 ), which we take as We derive the analytic expressions for these helicity amplitudes.We obtain the remaining helicity configurations, {+ − + , + − −}, through parity transformation (see appendix C of ref. [42]).We decompose the loop-level helicity amplitudes A (L) µ into gauge-invariant subamplitudes A (L) i,j µ , where the subscript i counts the number of closed massless fermion loops and j the number of external photons attached to closed fermion loops.The non-zero contributions are where n l denotes the number of charged lepton flavours running in the loops.Representative Feynman diagrams contributing to these subamplitudes are illustrated in figure 1. Amplitudes with a closed fermion loop attached to an odd number of photons vanish by Furry's theorem.We decompose the amplitude and subamplitude currents as ) The off-shell external leg is indicated by a bold line.
using the following basis written with the spinor-helicity formalism: (2.9) Readers not familiar with the spinor-helicity formalism may like to consult one of the many good reviews on the subject [69][70][71].Note that q 4 is orthogonal to the momenta p i by construction; one can in fact show that q µ 4 ∝ ε µνρσ q 1ν q 2ρ q 3σ .The subamplitude coefficients a (L) i,j;k can be related to the amplitude ones a (L) k through eq.(2.7).
The scattering amplitudes M (L) for fully on-shell processes (for instance, for 0 → e − e + γµ − µ + ) are obtained by contracting the amplitude currents A (L) µ (for 0 → e − e + γγ * ) with a suitable decay current V µ (in this example, γ * → µ − µ + ), as (2.10) In this manner, the on-shell amplitudes M (L) are given by the scalar product between the vector of coefficients (a 4 ), and that of decay-vector contractions (q 1 •V, . . ., q 4 •V).The coefficients a (L) k depend on the helicities of the three on-shell particles in eq.(2.1), while the decay vector V µ depends on the helicities of the particles the off-shell photon decays to.The helicity-summed interference between the L 1 -loop and the L 2 -loop matrix elements is then given by where the subscripts ⃗ h indicates the helicities of all on-shell particles -that is, including the decay products of the off-shell photon -and the overall constant factor averages over the helicities of the incoming particles.
The output of the computation described in section 3 is the set of four projections A (L) i,j • q k for each helicity configuration listed in eq.(2.6).From these, we determine the subamplitude coefficients a (L) i,j;k by inverting eq.(2.8), as where G is the Gram matrix of the vectors q i , that is, the matrix of entries G ij := q i • q j for i, j = 1, . . ., 4. At loop level, we write the subamplitude coefficients as where mon r (F ) are monomials of special functions F (see section 4), and the coefficients c r,w are rational functions of the kinematics.We drop the dependence on i, j, k, and L on the right-hand side of eq.(2.13) for compactness.We truncate the Laurent expansion around ϵ = 0 to the orders required for computing NNLO predictions.We express the coefficients c r,w as Q-linear combinations of a smaller set of linearly-independent coefficients (see section 3).The analytic expressions of the latter are given explicitly in terms of momentum twistor variables (see appendix C).We simplify these expressions through a multivariate partial fraction decomposition using MultivariateApart [72], and by collecting the common factors.
In the ancillary files [73], the directory amplitudes/ contains Mathematica files describing the bare helicity subamplitude currents A (L) i,j µ by their coefficients a (L) i,j;k in the form of eq.(2.13).The Mathematica script current.m is a reference implementation of the numerical evaluation of the bare amplitude coefficients a (L) k in eq.(2.13), including summation of subamplitudes in eq.(2.7), treatment of dependent helicities, and renormalisation scale restoration in eq.(2.5).The Mathematica script evaluation.wldemonstrates the construction of the five-particle on-shell amplitudes in eq.(2.10) for the process 0 → e − e + γµ − µ + , and their helicity-summation to obtain the squared matrix elements in eq.(2.11).The results of the script are checked against a reference point included in reference_point.json.
We perform the following checks of our amplitudes.

Ward identity
We verify the gauge invariance of the subamplitudes A (L) i,j µ by checking that they vanish on replacing the on-shell photon's polarisation vector with its momentum.

Finite remainder
We verify that the ϵ-poles of the bare amplitudes have the structure predicted by UV renormalisation and IR factorisation [76][77][78][79][80].We then subtract the expected poles and define finite remainders at one and two loops as where the square brackets separate renormalisation of UV poles from subtraction of IR poles.We present the derivation of these formulae in appendix D.

Setup of the calculation
In this section, we outline the workflow we used to calculate our amplitudes.Firstly, we generate all Feynman diagrams contributing to eq. (2.1) using QGRAF [81].Each diagram is then replaced with the corresponding Feynman rules for vertices, propagators, and external states, leading to a collection of D-dimensional Feynman integrals.Next, we filter the integrals according to eqs.(2.4) and (2.7) using a collection of Mathematica and FORM scripts [82,83].Within each subamplitude A (L) i,j µ , we then collect the integrals according to their topology, by which we mean a unique set of denominators.For example, the diagrams in figures 1d and 1e belong to different topologies, but those in figures 1c and 1f belong to the same topology (under the assumption of massless lepton propagators).At this point, the subamplitudes are sums of Feynman integrals over distinct integral topologies, with the numerators given by linear combinations of monomials that depend on the loop as well as the external momenta.To work with the projected helicity subamplitudes A (L) i,j • q k , we specify the polarisations of external particles according to eq. (2.6), as well as the projector q µ k of the off-shell photon from eq. (2.9).It is natural to express helicity-dependent objects using the spinor-helicity formalism.Then, the monomials of loop momenta contain the following scalar products and spinor strings: Their coefficients, on the other hand, are composed of the same type of objects, but do not contain any dependence on loop momenta k i .We express these coefficients using the rational parametrisation of the kinematics discussed in appendix C.This marks the start of our finite-field sampling procedure [51].The goal of this approach is to sidestep the algebraic complexity which typically plagues the intermediate stages of symbolic computations by evaluating numerically all rational coefficients.Using integers modulo some large prime number -which constitute a finite field -for the numerical evaluation allows us to avoid the loss of accuracy inherent to floating-point numbers, as well as the expensive arbitrary-precision arithmetic required by rational numbers.Manipulations needed to further process the rational coefficients are a completely separate problem from the calculation of the integrals or special functions that these coefficients multiply.In fact, they can be implemented as a series of rational transformations over finite fields.We stress that this is the methodology we follow at each step of the computation described below.In particular, we use the package FiniteFlow [52], which is conveniently interfaced to Mathematica.The analytic form of the coefficients is not known at any intermediate step.It is reconstructed from the finite-field samples only at the very end of the workflow.
Firstly, we note that not all integral topologies are independent: some of them can be written as subtopologies of others.For this reason, we define the set of maximal topologies, i.e. topologies with the maximum number of propagators allowed for L-loop, n-particle diagrams.In figure 2, we present the maximal topologies for the process under consideration in an arbitrary ordering of the external momenta (we give their explicit definitions in appendix A).Several orderings of the external momenta are relevant for the amplitudes, and we treat them as distinct families.Next, we map all topologies present so far onto one of these maximal topologies.The loop momenta dependent objects of eq.(3.1) are then expressed through the nine inverse propagators and irreducible scalar products (ISPs) associated with the chosen maximal topology.In this way, each subamplitude is now a sum of integrals compatible with IBP reduction [84,85], while their coefficients depend on the external kinematics and ϵ.We generate the required IBP relations using LiteRed [86].The resulting IBP system is then solved using the Laporta algorithm [87] with FiniteFlow's linear solver to yield the reduction of all the integrals present within our maximal topologies onto a much smaller subset of master integrals (MIs).It is beneficial to choose the MIs such that they satisfy differential equations (DEs) in canonical form [88] (see section 4.1).This property allows us to make the most of the optimisation measures we discuss below.We stress that the IBP reduction is also done numerically over finite fields, since the coefficients of the IBP relations are rational functions of external kinematics and the dimensional regulator ϵ.This is an important simplification, since analytic IBP reduction often proves to be the bottleneck of amplitude computations.For many amplitude applications, multiple permutations of the ordered topologies can appear.We outline a strategy to optimise the reduction in such situations in appendix B.
At this point, each projected helicity subamplitude A (L) i,j • q k is written as a linear combination of MIs multiplied by rational coefficients of ϵ and the kinematic variables.We now write the MIs in terms of a basis of special functions up to the required order in ϵ (see section 4).Finally, we Laurent expand the amplitude around ϵ = 0, the deepest pole being 1/ϵ 2L at L loops.The only task left is to reconstruct the rational coefficients of the special-function monomials from their samples over finite fields.In general, this might be a daunting challenge and its complexity stems from two separate factors.The workflow described so far is a series of rational operations chained together within a so-called dataflow graph [52].As such, we essentially have a black-box algorithm which takes numerical values of the kinematic variables as input, and returns the corresponding numerical values of the rational coefficients of the special-function monomials.The first factor is that several sample points are necessary to infer the analytic expression of these coefficients from their values in the finite fields.The required number is correlated with the polynomial degrees of the rational functions viewed as ratios of polynomials: the higher the degree, the more sample points are required.The second factor affecting the reconstruction complexity is the time it takes to obtain the values of the coefficients at each sample point.The more complicated the dataflow graph is, i.e. the more operations are chained together and the more difficult each operation is, the longer it will take to run the black-box algorithm.The most expensive operation in this regard is the evaluation of the solution to the IBP system.The total reconstruction time can thus be estimated as: reconstruction time ≈ (number of sample points) × (evaluation time per point) .(3.2) We emphasise that the sample evaluations can be run in parallel.For a detailed discussion of various strategies to improve the reconstruction time, see section 4 of ref. [41] and section 3.5 of ref. [47].Here, we give a brief overview of the tools that proved sufficient for this work.
First, we look for Q-linear relations among the rational coefficients of each helicity subamplitude.This typically requires few sample points with respect to the full reconstruction.We then solve these linear relations to express all coefficients in terms of a minimal subset of independent ones.Only the latter need to be reconstructed.Choosing them so that they have the lowest degrees often leads to a decrease in the complexity of the reconstruction.
The second strategy we employ is to match the rational coefficients with factorised ansätze informed by the singularity structure of Feynman integrals.The singularities of Feynman integrals can in fact be read off from the DEs they satisfy.For each coefficient we then write an ansatz made of the following factors: for all i, j, k = 1, 2, 3 such that i ̸ = j ̸ = k.This list includes denominator factors of the DEs satisfied by the MIs (listed by eq.(4.2)), as well as spinor structures aimed at capturing the phase information of helicity amplitudes.We then determine the exponents of the ansätze by comparing them to the coefficients reconstructed on a univariate slice of the kinematic variables [34], which are very cheap to obtain.We find that with this ansatz it is possible to determine all denominator factors -which indeed are linked to the singularity structure of the amplitude -and sometimes also some numerator factors.As a result, the undetermined functions yet to be reconstructed are of lower degree and require fewer sample points.We reconstruct the analytic form of the remaining rational functions using FiniteFlow's built-in multivariate functional reconstruction algorithm.Finally, we note that, for more computationally demanding processes, further ansatzbased techniques -for instance based on partial fraction decompositions or informed by the singularity structure of the amplitudes -may be used to optimise the functional reconstruction; see, for example, refs.[41-43, 46, 47, 89, 90].

Computation of the master integrals
The MIs for the relevant integral families were first computed analytically in refs.[53,54].The basis of special functions for four-particle kinematics with an off-shell leg has been subsequently studied in great detail with attention to the compactness of representations in terms of polylogarithms and numerical evaluation across a complete phase-space [56,57,91] (see also ref. [55] for a thorough discussion of the analytic continuation).We revisit this computation to obtain expressions for the MIs which are well suited for the amplitudecomputation workflow discussed in section 3. To this end, we compute the MIs for all permutations of the external legs in terms of a basis of special functions, following the approach of refs.[45,[64][65][66][67].In other words, we express all the Feynman integrals contributing to the amplitudes in terms of a set of special functions which are (algebraically) independent.Having such a unified and unique representation for all permutations of the integral families allows for simplifications and cancellations among different permutations of the Feynman integrals.This leads to a simpler expression of the amplitudes and to a more efficient functional reconstruction in the finite-field setup presented in section 3.
We emphasise that our results cover all MIs required for computing any two-loop fourparticle amplitude with a single external off-shell leg, and not just the ones required for the amplitudes presented in this work.
We discuss the construction of the basis in section 4.1, and how we express it in terms of multiple polylogarithms (MPLs) in section 4.2.Finally, we give some details about the numerical evaluation and the checks we performed in section 4.3.

Construction of the special function basis
We follow the strategy presented in ref. [67].The starting point are the DEs satisfied by the MIs for each family [92][93][94][95][96].Let τ label an integral family, e.g. the double-box in figure 2b for an arbitrary permutation of the external massless momenta.We choose a basis of pure MIs ⃗ g τ , that is, a basis which satisfies DEs in the canonical form [88] Here, d is the total differential, df are constant n τ × n τ matrices, with n τ the number of MIs of the family τ , and are called letters.We emphasise that this alphabet covers all permutations of the relevant integrals.Specific integral families may contain only subsets of it.Canonical bases for the relevant integral families are already available in the literature [57,62,63].Given that, by today's standards, finding canonical bases for these integral families is simple, we re-derived them using a mixture of methods: the package DlogBasis [97], the analysis of results in the literature for related integral families (massless two-loop five-point planar integrals [98] and two-loop four-point integrals with two massive external legs [99,100]), and a set of heuristic rules (see e.g.ref. [101]).We normalise the MIs such that their expansion around ϵ = 0 starts from order ϵ 0 , For the purpose of computing two-loop scattering amplitudes up to their finite part (i.e., up to order ϵ 0 ), it suffices to restrict our attention to w ≤ 4. Since the MIs satisfy canonical DEs (4.1), the ϵ-order of the MI coefficients ⃗ g (w) τ (⃗ s) equals their transcendental weight [88].We compute the derivatives of the MIs using FiniteFlow [52] and LiteRed [86].We do so only for the integral families with the ordering of the external momenta shown in figure 2, and obtain those for all other orderings of the external massless legs by permutation.We provide the definition of the pure MIs and the corresponding DEs for all one-and twoloop four-point one-mass families in figure 2 in the folder pure_mi_bases/ of our ancillary files [73].
In order to solve the DEs (4.1) we need boundary values, i.e., values of all MIs up to order ϵ 4 at a phase-space point.Due to the simplicity -by today's standards -of the integrals under consideration, an arbitrary (non-singular) phase-space point would do.Nonetheless, we make a more refined choice following some of the criteria of refs.[65,66].We choose the following point in the s 12 channel (see appendix E), motivated by two principles: that it is symmetric under the permutations which preserve the s 12 channel (i.e., swapping p 1 ↔ p 2 ), and that it contains few distinct prime factors.
The first condition reduces the number of permuted integral families we need to evaluate in order to obtain the boundary values.The second condition reduces the number of independent transcendental constants appearing in the boundary values, which simplifies the construction of the basis of special functions.The order-ϵ 0 boundary values ⃗ g (0) τ are rational constants.We obtain them up to their overall normalisation by solving the 'firstentry conditions' [102], i.e., by requiring the absence of unphysical branch cuts in the solutions.We fix the overall normalisation and the higher-order boundary values ⃗ g (w) τ (⃗ s 0 ) (for 1 ≤ w ≤ 4) by evaluating all MIs with AMFlow [103] (interfaced to FiniteFlow [52] and LiteRed [86]) at ⃗ s 0 with at least 60-digit precision.We anticipate from section 4.2 that, although we use floating-point boundary values, our results in terms of MPLs are fully analytic.
The canonical DEs (4.1) and the boundary values for all integral families are the input for the algorithm of ref. [67] for constructing a basis of special functions.We refer to the original work for a thorough discussion.Out of all MI coefficients up to transcendental weight 4, the algorithm selects a subset, denoted F := {F (w) i (⃗ s)}, which satisfy two constraints.First, they are algebraically independent, that is, there are no polynomial functional relations among them.Second, the MI coefficients of all families (including all permutations of the external massless legs) up to transcendental weight 4 are expressed as polynomials in the {F (w) i (⃗ s)} and the zeta values ζ 2 = π 2 /6 and ζ 3 .For example, an arbitrary weight-2 MI coefficient g (2) (⃗ s) has the general form with c i , d ij , e ∈ Q.This special subset of MI coefficients, {F (⃗ s)}, constitutes our special function basis.We give the number of functions in the basis in table 1.Note that there is freedom in the choice of which MI coefficients make up the basis.We make use of this freedom to choose as many basis-elements as possible from the one-loop family, then complement them with coefficients from the planar two-loop families, and finally complete them with coefficients from the non-planar two-loop families.In this way no two-loop MI coefficients appear in the one-loop amplitudes, and no non-planar two-loop MI coefficients appear in those amplitudes where only planar diagrams contribute (as is often the case in the leading colour approximation of QCD).} in the basis weight by weight.
The folder mi2func/ of our ancillary files [73] contains the expression of all MI coefficients (for all one-and two-loop integral families in all permutations of the external massless legs) up to weight 4 in terms of our special function basis.This result enables the efficient amplitude-computation strategy based on finite-field arithmetic discussed in section 3.However, at this stage the basis functions {F (w) i } are expressed in terms of Chen iterated integrals [104] and numerical boundary values ⃗ g (w) (⃗ s 0 ).This representation is excellent for investigating the analytic properties of Feynman integrals and amplitudes, but it is not readily suitable for an efficient numerical evaluation.In the next section we discuss how we construct a representation of the function basis in terms of MPLs and zeta values, which is well suited for an efficient and stable numerical evaluation.

Expression in terms of multiple polylogarithms
In this section we construct a representation of our function basis in terms of MPLs.Note that, while the basis functions {F (w) i } are by construction algebraically independent, the MPLs appearing in their representation constructed in this section may not be.Nonetheless, we take a number of simple measures to reduce their number and optimise the representation.The weight-n MPL of indices {a 1 , . . ., a n } and argument x is defined recursively as G(a 1 , a 2 , . . ., a n ; x) := We refer to ref. [105] for a thorough discussion.Since the letters in eq. ( 4.2) are rational and linear in all variables, we can solve the canonical DEs in eq. ( 4.1) algorithmically in terms of MPLs.Order by order in ϵ, the solution is given by starting from the constant weight-0 boundary values ⃗ g τ determined in the previous subsection.Here, γ is a path connecting an arbitrary base-point ⃗ s base to the end-point ⃗ s. τ (⃗ s base ).For ⃗ s base we may use the boundary point ⃗ s 0 in eq.(4.4), so that the constants ⃗ b (w) τ coincide with the boundary values determined numerically in the previous section.We follow a different approach, which allows us to trade all numerical constants in the expressions for zeta values.
We find it convenient to change variables from (s 12 , s 23 , s 4 ) to (z 1 , z 2 , s 4 ), with This way, there is only one dimensionful variable, s 4 , the dependence on which is fixed as an overall factor by dimensional analysis.We then integrate the canonical DEs as in eq.(4.8) along the following piece-wise path in the (z 1 , z 2 , s 4 ) space: (0, 0, 0) Since the Feynman integrals are divergent at the chosen base-point, the latter is understood in a regularised sense (we refer to section 4 of ref. [106] for a thorough discussion).Choosing (0, 0, 0) as base-point has the important benefit of removing spurious transcendental numbers that would pollute the solution were we to choose a base-point where the integrals are finite.As we will see below, only zeta values appear.Roughly speaking, we define regularised, finite values ⃗ b τ (⃗ s base ) by introducing a regulator and formally setting to 0 the (divergent) logarithms of the regulator.Since the integrals are finite at a generic end-point ⃗ s, the divergences at the base-point must cancel out with divergences arising in the integration.We can thus drop all these divergences.Provided that we do it consistently between the integration and the base-point values ⃗ b  as symbols and integrate the canonical DEs as in eq.(4.8) along the path in eq.(4.10) up to weight 4. We parameterise each piece of the path in eq.(4.10) linearly.For instance, γ 2 (t) = (z 1 , t, 0), with t ∈ [0, z 2 ].
• The γ 1 integration leads to MPLs with indices in {0, 1} and argument z 1 .
• The γ 3 integration leads to powers of log(−s 4 ), fixed by dimensional analysis.
Once we have obtained expressions for all MIs in terms of MPLs and symbolic constants ⃗ b (w) τ , we equate them to the numerical boundary values at ⃗ s 0 , and solve for the ⃗ b (w) τ .We use GiNaC [105,107] to evaluate the MPLs numerically.Finally, we use the PSLQ algorithm [108] to express the ensuing values of ⃗ b Contrary to the functions in the basis {F (w) i }, the MPLs in their representation satisfy functional relations.We make use of this freedom to optimise our expressions in view of their numerical evaluation by reducing the number of distinct MPLs that need to be evaluated.First, we use the shuffle algebra of MPLs to push all trailing zeros into logarithms through eq.(4.7) [105].Next, we employ the scaling relation which holds for x, a n ̸ = 0.As a result, all MPLs have argument 1 and indices Finally, we decompose the MPLs to Lyndon words [109] using PolyLogTools [110]; we refer to the latter work for a thorough explanation, and give here only a simple example.This procedure requires that we choose a symbolic ordering of the MPL indices.We choose l 0 ≺ l 1 ≺ l 2 ≺ l 3 ≺ l 4 , meaning that l 1 is greater than l 0 , and so on.Consider the MPL G(l 1 , l 0 ; 1), whose indices are not sorted according to the ordering above, since l 1 ≻ l 0 .We can use the shuffle algebra of MPLs to rewrite it in terms of MPLs whose indices are sorted according to the chosen ordering, as Doing this consistently throughout all expressions reduces the number of higher-weight MPLs in favour of products of lower-weight ones, which are cheaper to evaluate numerically.
To maximise the impact in this sense, we tested all possible orderings of the indices and selected the one -given above -which minimises the number of weight-4 MPLs We write the latter in terms of logarithms rather than MPLs as they play an important role in the factorisation of the IR divergences in the scattering amplitudes (see appendix D for the IR structure of the amplitudes we compute here).We stress that log(−s 4 ) is the only function of a dimensionful argument in our representation of the function basis.We provide in the folder mi2func/ of our ancillary files [73] the expression of the basis functions {F It is important to stress that the MPLs are multi-valued functions.For unit argument, there is a pole on the integration contour whenever one of the indices lies between 0 and 1.In this case the contour must be deformed in the complex plane, either above or below the pole, leading to different branches.Our MPLs are thus well-defined only in the kinematic region where all MPL indices in eq.(4.12) are either less than 0 or greater than 1, and we need s 4 < 0 for the argument of all logarithms in eq.(4.14) to be positive.We analytically continue the MPLs to the kinematic regions of interest by adding infinitesimal imaginary parts to the indices in the numerical evaluation, as we discuss in appendix E. The analytic continuation of these functions has been analysed in great detail in ref. [55], where it is achieved by suitable changes of variables so that the imaginary parts are extracted explicitly, and all MPLs are real and single valued.This approach could lead to faster evaluations although it was unnecessary for the applications considered here.

Performance and validation
We validated our results for the MIs of all families by crosschecking them against values obtained with AMFlow [103] at several random points in all the physical kinematic regions discussed in appendix E. Furthermore, we find agreement with the results of ref. [57].We employ GiNaC [105,107] to evaluate the MPLs.
Our results allow for an efficient and stable evaluation of the MIs, and are thus ready for immediate deployment in phenomenology.Indeed, the amplitudes we computed in this work have already been implemented in McMule [74,75] to provide the real-double-virtual electron-line corrections to eµ → eµ scattering.The evaluation is efficient, running at ≈ 130 events per second in the bulk of the phase space [111] using handyG [112] for the evaluation of the MPLs.

Conclusions
In this article, we calculated analytically the two-loop QED helicity amplitudes for the process 0 → ℓ lγγ * in terms of a basis of multiple polylogarithms that are suitable for fast and stable numerical evaluation.We employed modern finite-field evaluation techniques to reconstruct the amplitudes directly in terms of the special function basis, sidestepping the symbolic computation in all intermediate stages.As a by-product we have recomputed all two-loop master integrals for four-point functions with an off-shell leg up to transcendental weight four, and provide all the necessary ingredients needed to use them in amplitude computations with the same kinematics.
We hope these new results will now open the path to N 3 LO predictions that can be used for the future MUonE experiment.
Note added: Another computation of these amplitudes [113] became available on the arXiv shortly after ours.
B Optimised IBP reduction procedure for amplitudes with many permuted integral families An amplitude will in general have contributions from permutations of the ordered integral families shown in figure 2. To reduce the tensor integrals in the amplitude, IBP identities must be generated for all the permutations of these ordered families.This can lead to a very large IBP system.The performance of the reduction setup is extremely sensitive to the number of IBP identities required so, to minimise the memory consumption, we choose to generate IBP identities only for the ordered families.Next, we obtain the reduction for any permutation of these families by permuting the 'ordered' reduction numerically over finite fields.The result is then given in terms of MIs of each family permutation, but it is missing the symmetry relations that can be found between subsectors of different families.To express the final result in terms of a minimal set of MIs, we find such relations from a separate computation.One may account for integral symmetries using automated tools such as LiteRed [86].Since we use a pure basis of MIs, the symmetry relations amongst them will have rational numbers as coefficients.This is because the presence of any kinematic invariant would spoil the purity of the canonical DEs (see section 4), and would mean that such a symmetry relation in fact involves non-pure integrals.Therefore, the computation of the missing symmetry relations can be performed with all kinematic invariants set to numeric values, which significantly lowers the complexity of this task.Finally, we note that even if symmetries amongst the MIs were missed, a representation of the integrals in terms of a basis of special functions -as we construct in section 4 -would automatically incorporate the extra simplifications and so the same final result would be obtained.Nonetheless, in practice we do find it useful to include these symmetry relations, as they reduce the number of independent coefficients that have to be processed further.
The procedure can be summarised as follows: 1. Generate (analytic) IBPs for the six ordered families.
2. Compute the mappings between permutations of the MIs of the system above.
3. Take the tensor integrals in the amplitudes for each permutation of these families and solve the linear system over finite fields.
4. Apply the symmetry mappings between the MIs of each family permutation to find the minimal set for the full system.
Since there are a few additional bits of terminology, we can consider a concrete example to clarify everything.At one-loop, a four-point process with a single off-shell leg can be described by a single independent integral family which is simply the box topology (see appendix A for its explicit definition).Following the Laporta reduction algorithm leads to a basis of four MIs, MI box = {j box (1, 1, 1, 1), j box (1, 0, 1, 0), j box (0, 1, 0, 1), j box (1, 0, 0, 1)} , (B.1) which are the scalar box and scalar bubble integrals in channels s 12 , s 23 and s 4 respectively.An amplitude will, in general, be written in terms of three permutations of this family.Let us denote these permutations as j box,1234 , j box,2314 , and j box,3124 , where j box,1234 = j box as above and the additional superscript indices refer to the order of the external legs.Following our procedure we would load one set of IBP relations generated for j box .These identities can then be permuted numerically, for example as FiniteFlow graphs, to reduce tensor integrals in each of the three permuted families.The result is now in terms of twelve MIs: three boxes and nine bubbles.While the amplitude is already in a minimal basis of box integrals, there is clearly an over-complete set of bubbles.The independent bubbles are in the channels s 12 , s 23 , s 13 , and s 4 , so the five additional symmetry mappings are j box,2314 (1, 0, 1, 0) = j box,1234 (0, 1, 0, 1) , j box,3124 (1, 0, 1, 0) = j box,2314 (0, 1, 0, 1) , j box,3124 (0, 1, 0, 1) = j box,1234 (1, 0, 1, 0) , j box,2314 (1, 0, 0, 1) = j box,1234 (1, 0, 0, 1) , j box,3124 (1, 0, 0, 1) = j box,1234 (1, 0, 0, 1) .

(B.2)
After applying these identities we arrive at the final result with seven MIs which cover all permutations of the integral families.This approach would not lead to any significant performance enhancements in this simple example of course, but it can be particularly important when considering high-multiplicity examples where the number of permutations is high.

C Rational parametrisation of the kinematics
Since we are applying finite-field techniques to helicity amplitudes, we employ a rational parametrisation of the external kinematics using Hodges's momentum twistor formalism [115].While this is not essential to combat the algebraic complexity for the kinematics considered here, it does provide a convenient parametrisation of the spinor products.The single-off-shell four-particle phase space p is obtained from a massless five-particle parametrisation q (defined in appendix A of ref. [41] with {x 2 ↔ x 4 , x 3 ↔ x 5 }) under The momentum twistor variables x i for p are then related to the scalar invariants ⃗ s through Momentum twistors allow us to express any spinor expression as a rational function in the variables x i .In this representation the helicity scaling is however obscured, as we have fixed the spinor phases in order to achieve a parameterisation in terms of the minimal number of variables (see e.g.ref. [116]).Therefore, we need to manually restore the phase information at the end of the computation.This can be achieved by multiplying the momentum twistor expression by an arbitrary factor Φ with the same helicity scaling as the helicity amplitude under consideration, divided by that factor written in terms of momentum twistor variables.For example, for the helicity configurations of eq.(2.6), we can use the phase factors which in our momentum twistor parameterisation are given by (C.4) We refer to appendix C of ref. [42] for a thorough discussion of how to restore the phase information in a momentum twistor parameterisation.

D Renormalisation and infrared structure
We renormalise the coupling constant by trading the bare coupling α bare for the renormalised one α R through with S ϵ = (4π) −ϵ e ϵγ E .The renormalisation factor Z α in the MS scheme is [117,118] Z The β-function is defined from the renormalised coupling as and expanded as The photon wavefunction renormalisation factor is Z A = Z α , which we include due to the external off-shell photon.The complete renormalisation procedure then is where α bare is expressed in terms of α R through eq.(D.1).
The IR poles of the renormalised amplitude factorise as [76-80] so that Z(α R ) captures all IR poles and F µ is a finite remainder.We obtain the explicit two-loop expression of the IR factor Z(α R ) by choosing QED parameters (C A = 0, C F = 1, and T F = 1) in the non-abelian gauge-theory expressions of ref. [80].We expand it as The coefficients Z (L) are expressed in terms of the anomalous dimension Γ = γ cusp ln −s 12 µ 2 + 2γ l + γ A , (D.9) and its derivative Here, γ cusp is the cusp anomalous dimension, while γ l and γ A are the lepton's and the photon's collinear anomalous dimensions, respectively.We expand all anomalous dimensions y ∈ {Γ, γ i } as with coefficients Finally, the coefficients of the IR factor Z up to two loop are given by Putting together the subtraction of UV and IR poles, and expanding the resulting finite remainder F µ (α R ) in α R leads to the definitions in eq.(2.14).

E Analytic continuation
We analytically continue the MPLs by adding an infinitesimal positive (or negative) imaginary part to the MPL indices l i in eq.(4.12) whenever they fall between 0 and 1.The imaginary part of each index prescribes how to deform the integration contour around the pole associated with it.We do similarly for the logarithms in eq.(4.14).To this end, following ref.[55], we change variables from (s 12 , s 23 , s 4 ) to (s 12 , s 23 , s 13 ), with s 4 = s 12 +s 23 +s 13 .We then add a small positive imaginary part to the latter variables, as where c 1 , c 2 and c 3 are arbitrary positive constants, and δ is a positive infinitesimal.Finally, we check whether this substitution gives a positive or negative imaginary part to each MPL index l i .This depends on the domain of the kinematic variables.We focus on three kinematic regions which are of phenomenological interest.The analytic continuation for any other region may be obtained similarly.
Electron-line corrections to e − µ − → e − µ − γ.To define the domain of the kinematic variables relevant for this application, we embed the four-particle off-shell process of eq.(2.1) in the five-particle process e − µ − → e − µ − γ.We then determine the kinematic constraints for the five-particle process (see e.g.appendix A of ref. [66]), and from them derive the constraints on the four-point off-shell kinematics.The result is The imaginary part of l 1 may be either negative or positive in P eµ→eµγ .However, it is strictly negative in the subregion of P eµ→eµγ where 0 < l 1 < 1.We therefore assign a negative imaginary part to l 1 whenever 0 < l 1 < 1 in P eµ→eµγ .The analysis of the other indices follows similarly, and is summarised in table 2. The arguments of the three logarithms in eq.(4.14) are positive in P eµ→eµγ .
Corrections to e − e + → γγ * .The relevant domain of the kinematic variables in this case can be derived directly for the four-point kinematics, and is typically named the s 12 channel.It is given by P eē→γγ * := {⃗ s : s 23 < 0 ∧ s 13 < 0 ∧ s 12 > −s 23 − s 13 } . (E.4) The MPL indices l 2 , l 3 and l 4 can never fall between 0 and 1 in P eē→γγ * , and hence require no analytic continuation.We instead need to add a positive imaginary part to l 1 .In this region also the logarithms in eq.(4.14) need to be analytically continued.The argument of log(s 12 /s 4 ) is positive in P eē→γγ * .By adding imaginary parts to the arguments of the other logarithms and studying them where the arguments are negative in P eē→γγ * , we determine that the analytic continuation is achieved through the following replacements: All MPL indices l i in eq.(4.12) are either l i < 0 or l i > 1, hence no analytic continuation is required.The same holds for the first two logarithms in eq.(4.14), whose arguments are positive.The only function which needs to be analytically continued is log(−s 4 ).We achieve this by replacing log(−s 4 ) −→ log(s 4 ) − iπ . (E.7) The information about the imaginary parts of the MPL indices can be fed into the publicly available libraries for evaluating these functions numerically, such as FastGPL [119], GiNaC [105,107], and handyG [112].This typically leads to longer evaluation times with respect to MPLs which do not need analytic continuation.We find that this is not an issue for the planned applications of our results (see section 4.3).Nonetheless, we note that a more performant evaluation may be achieved by tailoring the representation to the kinematic region of interest in such a way that no MPLs require analytic continuation.We refer to ref. [53][54][55][56][57] for a detailed discussion.

16 A 18 C Rational parametrisation of the kinematics 19 D Renormalisation and infrared structure 20 E
Definition of the Feynman integral families 16 B Optimised IBP reduction procedure for amplitudes with many permuted integral families Analytic continuation 21

, 1 µFigure 1 :
Figure 1: Representative Feynman diagrams for the subamplitudes defined in eq.(2.7).The off-shell external leg is indicated by a bold line.

Figure 2 :
Figure 2: The two-loop amplitudes include six ordered integral families: three pentatriangles, a double-box, and two crossed double-boxes.All are planar except for the crossed double-boxes.The off-shell external leg is indicated by a bold line.External legs have outgoing momenta.

τ
, this leads to a finite and unique result.In practice, we fix the finite base-point values ⃗ b s) evaluated at the boundary point ⃗ s 0 against the boundary values discussed in the previous subsection.We therefore keep the ⃗ b (w) τ ζ 2 and ζ 3 .As a result, we obtain a fully analytic representation of all MIs -and thus of our special function basis {F (w) i } -in terms of MPLs and zeta values, up to weight 4.

Table 1 :
Number of functions {F