Lepton pair production at NNLO in QED with EW effects

We present a fully differential calculation of lepton pair production, taking into account the dominant next-to-next-to-leading order QED corrections as well as next-to-leading order electroweak and polarisation effects. We include all lepton masses, hard photon emission, as well as non-perturbative hadronic corrections. The corresponding matrix elements are implemented in the Monte Carlo framework McMule. In order to obtain a numerically stable implementation, we extend next-to-soft stabilisation, a universal technique based on a next-to-leading-power expansion, to calculations with polarised leptons. As an example, we show results tailored to the Belle II detector with the current setup as well as a potential future configuration that includes polarised beams.


Introduction
Thanks to its high luminosity, Belle II is expected to produce about 45 billion τ τ events over its lifetime [1], roughly fifty times more than Belle I [2] and a hundred times more than at BaBar [3].This increase in statistics will allow for precision measurements of very rare Standard Model (SM) decays such as τ → ν νℓℓ ′ ℓ ′ or τ → ν νℓγ as well as put bounds on charged-lepton-flavour-violating decays such as τ → ℓℓ ′ ℓ ′ .For the SM decays, even differential measurements in terms of Michel parameters will be possible.In this case, spinspin correlations between the two taus of ee → τ τ can be exploited [4].Further, it was recently proposed to measure the anomalous magnetic moment of the tau to 10 −6 by exploiting transverse and longitudinal asymmetries [5].
Hence, the production cross section for ee → τ τ needs to be known as precisely as possible.For the centre-of-mass (CMS) energy used at Belle II ( √ s ≈ 10.5 GeV), QED effects still dominate even though electroweak (EW) effects start to become relevant.With EW effects, we mean all contributions due to EW interactions but without contributions due to pure QED.The Monte Carlo code KKMC [6][7][8] combines a parton shower with fixed-order EW corrections at next-to-leading order (NLO).It has been very useful for many experimental studies [1].The EW effects were further studied with the SANC program, accounting for the polarisation of the incoming leptons [9].
However, the improvements expected from Belle II warrant a renewed theoretical effort.Currently no NNLO-QED calculation for ee → ℓℓ (with ℓ = e) exists as the necessary two-loop integrals are not yet known with full mass dependence.However, a recent theoretical interest in e-µ scattering [10] was inspired by the MUonE experiment [11][12][13].As part of this effort, the necessary integrals were computed in the limit of vanishing electron mass m e → 0 [14,15].This was very recently used to assemble to the full two-loop matrix element (squared) for ee → ℓℓ with m e = 0 [16] which is an important part of the full NNLO-QED.Assuming m e is small compared to all other scales of the process, this matrix element can be used to obtain the full matrix element up to terms suppressed by O(m 2 /Q 2 ).This massification procedure was first developed in [17][18][19] and later extended [20] to the case of a second, heavy mass.
However, the smallness of the electron mass means, that for a first estimate of the NNLO-QED correction, it is sufficient to just consider the electronic corrections, i.e. those due to the electron, and ignoring the more complicated mixed contributions.This can be done in a gauge-invariant manner by assigning different charges for each lepton family and only take contributions proportional to Q2 ℓ Q 6 e .This was demonstrated at NLO-QED [21] and then exploited to calculate the dominant NNLO-QED contributions to eµ → eµ [22,23].Note that these Q 2 ℓ Q 6 e corrections can be calculated exactly in the electron mass m e without approximation as the virtual corrections are just the heavy-quark form factor [24].
In this paper, we use the McMule framework to extend our previous calculation [22] of eµ → eµ to cover the electronic, or initial-state radiation (ISR), NNLO-QED corrections to ee → ℓℓ.This means that our calculation includes the full NLO-QED (incl.the Q 3 ℓ Q 3 e box contributions) but only the leading, i.e.Q 2 ℓ Q 6 e , NNLO-QED corrections.In particular, we do not include final-state radiation (FSR, ) since these are suppressed. 1 We further treat the incoming electrons as polarised and include EW effects at NLO since these are becoming relevant at the energy at which Belle II operates.
Since the CMS energy of Belle II ( √ s ≈ 10.5 GeV) is significantly higher than for MUonE ( √ s ≈ 0.4 GeV), new numerical problems arise in the real-virtual matrix element, especially in the case of soft emission.These can be efficiently handled using next-to-soft (NTS) stabilisation [26], i.e. using a next-to-leading power (NLP) expansion if the emitted photon becomes soft.
This paper is organised as follows: in Section 2, we briefly summarise the calculation as implemented in the McMule framework.Next, we explain how NTS stabilisation changes when considering polarised particles in Section 3. Finally, we present some results for tau production cross sections and asymmetries, both in general and tailored to Belle II in Section 4 before concluding in Section 5.

Overview of the calculation
We consider the scattering process taking into account the full NLO-EW corrections and the electronic NNLO-QED corrections (Q 2 τ Q 6 e ) but drop all remaining NNLO-QED terms, i.e.FSR (Q 6 τ Q 2 e ) and IFI (Q ) as discussed before.Since we are well below the Z peak, we can expand the NLO-EW corrections by considering the masses of the EW bosons (M 2 Z , M 2 W , and M 2 H ) much larger than all other scales of the process (m 2 e , m 2 τ , s = (p 1 + p 2 ) 2 , and t = (p 1 − p 3 ) 2 ).We then expand in the ratio of the heavy scale to the light scale, taking the first two terms of the expansion, i.e. keeping all terms O {s, t, m 2 e , m 2 τ } {M 2 Z , M 2 W , M 2 H } .For simplicity, we write this as an expansion in 1/M Z .This procedure corresponds to how one expands in an effective field theory in 1/Λ while taking into account all effects up to and including dimension six.In this view, the base theory is QED and the underlying theory is not New Physics but rather the full SM.The resulting theory is a subset of what is often referred to as low-energy effective field theory (LEFT) [27][28][29].
Considering only the electronic, i.e.Q 2 τ Q 6 e , contributions at NNLO-QED means that we have exactly calculated the main source of logarithms in the electron mass, i.e. those terms containing α 2 log 2 (m 2 e /Q 2 ) (where e is some other scale).Since we perform this calculation with full m e dependence we do include also power-suppressed terms that are ∝ Q 2 τ Q 6 e .However, since such terms are also contained in the mixed IFI corrections, we have not included all logarithms and power-suppressed terms involving m e .Similar logarithms in the tau mass do appear and those ∝ Q 6 τ Q 2 e (FSR) could be trivially calculated.However, since m 2 e ≪ m 2 τ ∼ Q 2 these logarithms are not expected to be overly large.Diagrams were generated with FeynArts [30] and QGraf [31] and calculated using Package-X [32] with full electron and tau mass dependence.Ultraviolet (UV) and infrared (IR) divergences are regularised in d = 4 − 2ǫ dimensions and the renormalisation is performed in the on-shell scheme up to NLO-EW and NNLO-QED.For the EW corrections, this means that we would like to use e, M W , M Z , M H , and m ℓ as input parameters [33,34].However, it is beneficial to use G F instead of M W as it has much higher precision ( For technical reason, the matrix elements in McMule are expressed through s2 W rather than G F .s 2 W = 0.226202 was obtained through where the one-loop expression of ∆r was taken from [34] to all orders in 1/M Z but without including leading higher-order contributions. The calculation is split into fermionic contributions that are due to fermionic vacuum polarisation (VP) effects (Section 2.1) and bosonic ones (Section 2.2).
Our calculation is performed with longitudinally polarised electrons (see for example [40,41]).We introduce a polarisation vector n i along the beam direction for each particle that takes the form in the particle's rest frame.Of course n i • p i = 0 in any frame.|P i | ≤ 1 is the degree of polarisation that can be chosen as required by the beam parameters.To implement this, we modify the completeness relation of the spinors to Alternatively, it may be simpler to calculate the matrix element for fully polarised initial states, i.e. with This is for example done in OpenLoops [42,43] which we will be using in Section 2.2.However, for most phenomenological applications we are interested in partial polarisation.For the (parity conserving) QED part, we can recover the general result as where −1 ≤ P i ≤ +1 can be arbitrary.

Fermionic corrections
Up to NNLO, all fermionic corrections to our process are due to closed fermion bubbles.They can be split into leptonic contributions (electrons, muons, and taus) and non-perturbative hadronic effects (HVP). 2 At one-loop, these can be written in terms of the photonic and Z currents for the lepton flavour ℓ = e, τ j (ℓ,γ) with the momenta properly chosen.The (Zℓ l)-couplings are flavour-universal The (renormalised) one-loop amplitude for the fermionic vacuum polarisation can be divided into three parts where the transversal fermionic self-energies Σ renorm.

ij,f
are renormalised in the on-shell scheme with the conditions [34] 3 Fermionic contributions due to other particles in the EW sector (such as the Higgs) are suppressed by at least O(1/M 4 Z ) and hence already dropped.Corrections due to boson loops are included in Section 2.2.For the ZZ term, we extract the part that is that arises from the renormalisation of M Z .Hence, it is determined through the fermionic part of the (unrenormalised) self-energy Z where it can be calculated perturbatively as it has no kinematic dependence From the renormalisation conditions (11) we also find the explicit expressions for the renormalised selfenergies [45] By defining the fermionic VP function Πij [45,46] the γγ and γZ amplitudes can be written as Using the definitions of alphaQED [46] ΠγZ where the 3 refers to the third component of the isospin current and γ to the QED current, we arrive at the final expression Here, all non-perturbative Πij (including Πij (M 2 Z )) are taken from alphaQED and can be obtained more or less directly from R ratio data.However, Π3γ is sensitive to fermions flavour in a different way than Πγγ , necessitating a flavour recombination.The simplest strategy for this is assuming that SU(3) f is an exact symmetry which results in Π3γ = 1  2 Πγγ .alphaQEDc19 uses a more complicated strategy based on vector meson dominance (VMD) [46].The perturbative, leptonic contributions to Πij are included in the usual way by calculating the corresponding one-loop diagrams.
Beyond NLO, we have to account for QED-VP insertions into loop diagrams.This is done using the hyperspherical method [47,48] that was used for µ-e scattering [49] and implemented in McMule for µ-e, ℓ-p, and Møller scattering [22,50].There are further two types of leptonic self-energy corrections at two loop in QED: the product of two one-loop self-energy bubbles ("bubble chain") and the genuine two-loop self-energy that can be obtained from [51].Since these are fully perturbative, their inclusion is straightforward.

Bosonic corrections
At NLO, the bosonic corrections are given by two separately divergent types of contributions: real (R) and virtual (V).
The virtual contribution includes in total about 250 one-loop diagrams.The majority of them is due to self-energy corrections arising from purely bosonic loops such as W corrections to the photon propagator (all NLO-EW).The remaining diagrams can be divided into vertex correction diagrams for the electronic (Q 2 τ Q 4 e ) as well as for the tauonic (Q 4 τ Q 2 e ) part and box contributions (Q 3 τ Q 3 e ) involving box diagrams with two photons (NLO-QED), one heavy boson (Z, W, H, or Goldstone bosons) and one photon, or two heavy bosons (all NLO-EW).The latter can still contribute at O(1/M 2 Z ) even though two heavy boson propagators enter the calculation.Once the calculation and expansion of the one-loop diagrams is completed, we renormalise as discussed at the beginning of Section 2.
For the electronic corrections, the VV can be obtained from the on-shell renormalised heavy-quark form factor [24] which is written in terms of harmonic polylogarithms (HPLs) [52].This allows for trivial analytic continuation into the time-like region we are interested in.
We use OpenLoops [42,43] in its standard mode for the RV corrections.While OpenLoops is extremely stable, its standard mode may not be sufficient for soft or collinear emission, especially in the case of small fermion masses.To address this issue, we use NTS stabilisation [26].The basic idea of this method is to switch to an expanded matrix element if the (rescaled) photon energy ξ = 2E γ / √ s drops below a certain cut-off.This cut-off is usually varied between 10 −5 and 10 −2 to ensure that the final result does not depend on its value.
Below the cut-off, we use a matrix element that is expanded for small photon energies up to NLP.To do this, we use an extension of the LBK theorem [53,54] to one loop [55,56] that we will discuss in the next section.

LBK theorem for polarised particles
Following [55,56], we will extend the LBK theorem to polarised cross sections at one loop.We use ξ as the expansion parameter and write the photon momentum as p γ = ξk.
To better understand what happens at one loop, let us first review the changes in the case of polarised particles in the proof at tree-level where similar results have been obtained before [57,58].By splitting A (0) n+1 into contributions from internal and external legs and using gauge invariance, the NTS contribution can be written as [53-56] with the LBK operator When squaring A (0) n+1 , we have to consider the interference between the leading-power (LP, O(1/ξ) at amplitude-level) and NLP (O(ξ 0 ) at amplitude-level) terms.To do this in the unpolarised case, we would use the identity since u(p i )ū(p i ) = / p i + m i .In the polarised case, we have to use ( 5) and ( 22) gets modified accordingly Hence, the matrix element is finally obtained after summing over the polarisation of the photon Thus, the calculation of the NTS term at tree-level remains straightforward as we just need to also calculate the derivatives w.r.t. the polarisation vector.
To extend this discussion to the one-loop level, we use the method of regions [59].It was shown in [55], that the amplitude of the soft contribution is The function S(p i , p j , k) ∼ k can be calculated universally and is presented in [55].The hard contribution is closely related to the LBK theorem (24) with one important subtlety related to the following external leg corrections [56] A The vertex correction to the soft photon emission spoils the basic assumption of the LBK theorem that diagrams with internal emission do not contain any 1/p γ poles.Further, the self-energy correction is technically an external correction and could be expanded using the normal LBK theorem.Hence, these contributions do not reduce to the non-radiative amplitude.Instead, one can show that (26) results in an extra contribution of the form When interfering this contribution with the LP term of ( 20) we find In the unpolarised case, this contribution vanishes after some Dirac algebra.However, if n i = 0, it does not and instead results in Here, we need a modified version of the tree-level matrix element We can now write down a version of the LBK theorem that is valid both at one-loop and in the case of polarised external particles The new term M (0),′ n , while easy to calculate, has severe consequences for the structure of the NTS approximation.Every other term in (32) is directly related to the reduced process, either at one-loop (M n , on the other hand, is a new structure that spoils the elegance of the LBK theorem and its extensions. We have numerically verified that ( 32) is correct by taking the limit ξ → 0 of the real-virtual matrix element relevant for this process (as shown in Figure 1) and also for µ → ν νeγ.

Results
To validate our calculation, we have crossed it to eµ → eµ and compared the NLO-EW with [21] and the NNLO-QED with [22,23] at the level of the differential distributions.We found full agreement in both cases.
In the following, we present some results for e + e − → τ + τ − at √ s = 10.5830052GeV both without any cuts and tailored to Belle II.We stress that these are just examples and that McMule can calculate any IR-safe observable.
We write the total cross section as which is divided into the pure QED and the EW part.The former are further split into LO (σ QED ), and NNLO (σ EW as they turn out to be similar in size.As explained in the next section, this is the result of a cancellation within the LO-EW contributions.
In the interest of Open Science, all raw data, analysis pipelines, and plots can be found at [60] https://mule-tools.gitlab.io/user-library/dilepton/belle

Cross section without cuts
We begin by considering the cross section integrated over all of phase space without any cuts.The relevant K factors are defined as δK (1) = σ (1) QED σ (0) QED , δK (2) = σ We further consider the forward-backward asymmetry in the CMS frame which we denote with a * Following Belle's convention [61], the angles θ τ ± are defined w.r.t. the incoming positron. 4At a given order, A FB is defined to also contain all contributions below it, i.e.
The cross sections and asymmetries are shown in Table 1 for the unpolarised case and the cases where both electrons are polarised parallel (+) or anti-parallel (−) w.r.t.their direction of flight with a degree of polarisation of 70% in their rest frames.Note that in the QED case, parity invariance implies that there are only two independent configurations since In the EW case, parity is violated but CP is still conserved.Hence, implying three independent configurations.The angular distributions used for A FB are shown in Figure 2 for the unpolarised case.We note that even though the NNLO-QED and EW corrections are similar in size at the level of dσ/dθ * , the latter are much smaller for the integrated cross section.This is because the NNLO-QED corrections are symmetric while the EW corrections are largely antisymmetric as can be clearly seen in Figure 2. The dominant antisymmetry of the EW corrections is due to the coupling structure.Calculation at tree-level (unpolarised) with m = M = 0 shows that the leading symmetric contribution ∼ cos 2 θ is suppressed by As a result, the integrated cross section σ The values in ( 40) and ( 41) are given for the unpolarised case.If polarisation effects are taken into account σ EW and σ EW are similar in size.Let us also consider the effects of taking further terms in the 1/M Z expansion.If (40) is extended by dimension-six-squared5 , and dimension-eight contributions, we find   ).When the electrons are polarised, the degree of polarisation is 70% in their respective rest frames.Unless otherwise indicated, all digits are significant.
From this it is clear that at the differential level the expansion is perfectly justified.However, due to the aforementioned antisymmetry of the dimension-six term, the total LO-EW cross section receives sizeable corrections from the dimension-six-squared term.Hence, to avoid these effects, we include the LO-EW contribution without expansion in 1/M Z .Note that the zero crossing of the EW corrections in Figure 2 does not happen at exactly 90 • but is slightly offset due to the small symmetric part in (39).
At LO in QED, A FB is exactly zero as expected, while at NLO in QED, the mixed tauonic-electronic contribution induces a finite but small asymmetry.This is similar to the NLO-QCD effects resulting in a non-zero A FB for the hadronic tt production [62].In principle, this would continue at NNLO in QED.However, the purely electronic contributions considered here are perfectly symmetric w.r.t.θ * τ ± and therefore do not contribute to A FB .The EW corrections are almost perfectly antisymmetric but much smaller than the QED corrections.This means that the full NNLO-QED corrections will be required to give a meaningful result for A FB . 6However, unless calculation of the QED two-loop matrix elements with full m e dependence becomes available, it will not be possible to do this for the (±±)-polarisation configuration which requires a helicity flip that cannot be obtained with massification alone.

Predictions for Belle II
Next, we tailor our calculation to Belle II.The detector is asymmetric since the electron beam's energy is higher than the positron beam's We approximate the detector's geometric acceptance by requiring that the tau leptons are produced within the geometric acceptance [61] 17 The angular distribution of the outgoing taus is shown in Figure 3.The LO distribution vanishes below ≈ 53 • because of the cut on the other particle.However, once real emission is allowed, the angle can become much smaller.Once again, the EW corrections are similar in size to the NNLO-QED ones.The SuperKEKB beams are currently unpolarised which is reflected in our calculation.However, recent proposals suggest that this could be changed in the future, aiming for 70% polarisation [63].To study this QED curves for θ τ + and θ τ − are identical.case, we consider the ratio between the polarised and unpolarised angular distributions of the τ − , both in the lab frame with cuts (θ τ − ) and the CMS frame without (θ * τ − ) Note that the first term in (45) is not centred around one but instead around σ(±+)/σ(00).Hence, we subtract this overall shift to centre R(±+) around zero to make the comparison between R(++) and R(−+) easier.R * and R are shown in Figure 4. Let us first consider the simpler R * (Figure 4a).Since the outgoing taus are not polarised, R * = 0 at LO.At and beyond NLO, this is no longer true because hard-collinear ISR causes a helicity flip in the emitter. 7With (45) normalised as it is, adding all polarisations results once again in a flat line, meaning that we have recovered the unpolarised result.Boosting to the lab frame (Figure 4b) stretches the distributions for forward emissions (53 •  θ τ − 120 • ) and squeezes them for backward (120 • θ τ − ≤ 150 • ).Further, since the cuts (44) mean that hard emission is required for θ τ − 53 • , the effect is significantly enhanced.

Conclusion
We have presented a fully differential calculation of the dominant NNLO-QED and NLO-EW corrections for di-lepton productions, including fermionic and bosonic corrections.We find that the EW corrections are of similar size to the NNLO-QED ones for √ s ≈ 10.5 GeV meaning they are vital for Belle II.To perform this calculation, we have extended the strategy of next-to-soft stabilisation to polarised observables, which, combined with OpenLoops, allows for a fast and stable evaluation of the real-virtual matrix element.
All matrix elements were implemented in the parton-level Monte Carlo code McMule which allows the user to calculate arbitrary IR-safe observables.As a first example, we have calculated differential predictions for Belle II, both for polarised and unpolarised initial states.
Since this calculation only includes the dominant contribution, a natural next step would be the inclusion of the full set of NNLO-QED corrections, especially when considering asymmetries such as A FB .The relevant two-loop matrix elements are known in the unpolarised case for m e = 0 [16] but would need to be extended, in a first step, to the polarised case.Finally, this work will need to be combined with [25] to properly include the numerically delicate real-virtual corrections.Work to that end is currently in progress.
Should even higher precision be required, resummation is required.Currently, McMule is calculating strictly at fixed order which means that some important logarithmically-enhanced contributions are not considered beyond NNLO-QED.These can be resummed to all orders using fragmentation functions (for final state) or parton distribution functions (for initial state).For a recent review on this topic, see [64] and references therein.However, this was not done in the present study as any analytic resummation limits what observables can be calculated.Further, there is an effort within the McMule collaboration to include a YFS parton shower that is similar to PHOTONS++ [65] and matched to NNLO-QED.The ratios R ( * ) (++) in orange and R ( * ) (−+) in blue, both in the CMS frame and in the lab frame.Note that this observable is very insensitive to EW effects and that the shape of the plots originates from the QED corrections.This means that we only need to show two out of the four different polarisation configurations.

Figure 1 :
Figure 1: Convergence of the soft limit at LP and NLP of the dominant one-loop corrections to e + e − → ℓ + ℓ − γ.The reference value M exact is calculated with arbitrary precision in Mathematica.

Figure 2 :Figure 3 :
Figure2: The angular distribution of the two taus in the final state for an unpolarised initial state in the CMS frame.Note that the scale changes from linear to logarithmic at +0.06.
R * in the CMS frame without cuts at NLO-QED (upper panel) and for the full result (lower panel).R in the lab frame with cuts for the full result.Note that the lower panel is zoomed in onto the θ τ − region that is allowed at tree-level.

Figure 4 :
Figure4: The ratios R ( * ) (++) in orange and R ( * ) (−+) in blue, both in the CMS frame and in the lab frame.Note that this observable is very insensitive to EW effects and that the shape of the plots originates from the QED corrections.This means that we only need to show two out of the four different polarisation configurations.

Table 1 :
Cross sections and asymmetries for e + e − → τ + τ − up to NNLO-QED and NLO-EW for all polarisation configurations (e + e −