McMule: Precision QED Framework

• McMule is a Fortran-based precision QED framework that delivers fully differential, NNLO and beyond predictions for low-energy lepton scattering and decay processes.
• It employs techniques like subtraction-based infrared handling, massification, and next-to-soft stabilization to address finite mass effects and numerical instabilities.
• The framework supports diverse applications including muon-electron scattering, Møller and Bhabha processes, evolving from a fixed-order integrator to an NNLO event generator.

McMule, short for Monte Carlo for MUons and other LEptons, is a Fortran-based framework for fully differential higher-order QED calculations for low-energy scattering and decay processes involving leptons. It was developed as a precision-theory infrastructure for experiments in which finite lepton-mass effects, fiducial cuts, and radiative corrections are numerically decisive rather than perturbative afterthoughts. In its mature form, McMule combines exact-mass kinematics, subtraction-based infrared handling, semi-analytic and automated loop technology, and process-dependent numerical stabilisation; it began as a fixed-order Monte Carlo integrator and later gained NNLO event-generation capability (Banerjee et al., 2020, Ulrich, 7 Jan 2025).

1. Scope and scientific role

McMule is designed for the low-energy precision frontier: muon decay, lepton–lepton scattering, and lepton–proton scattering in regimes where observables are differential, often not collinear safe, and sensitive to logarithms regulated by the physical electron and muon masses. The framework is explicitly positioned for programs such as MUonE, MUSE, P2, PRad, PRad II, MEG II, Mu3e, and Belle II, and later work extends that role to MOLLER and other electroweak-sensitive low-energy observables (Banerjee et al., 2020, Ulrich, 7 Jan 2025, Kollatzsch, 10 Mar 2026).

A concise summary of process coverage stated in the literature is as follows.

Process family	Stated accuracy	Representative context
$e\mu\to e\mu$	NNLO+	MUonE
$\ell p\to \ell p$	NNLO	P2, MUSE, PRad, QWeak
$e^-e^-\to e^-e^-$	NNLO	PRad II, MOLLER
$e^+e^-\to e^+e^-$	NNLO	luminosity monitoring
$ee\to \ell\ell$	NNLO$\pm$	$R$-ratio and tau studies
$\mu\to \nu\bar\nu e$	NNLO+	MEG, DUNE
$e^+e^-\to \gamma\gamma$	NNLO	low-energy collider luminosity

In this notation, $+$ means future improvements are expected, and $\ell p\to \ell p$0 means the quoted order is not known exactly, as stated in the framework overview (Ulrich, 7 Jan 2025). The process list is broader than a single muon-decay code or a single scattering package: McMule was conceived as a reusable fixed-order environment in which amplitudes, subtraction, and observable definitions are modular rather than process-specific (Banerjee et al., 2020).

2. Perturbative organization and infrared treatment

The framework uses the standard fixed-order decomposition

$\ell p\to \ell p$1

with NNLO split into double-virtual, real-virtual, and double-real sectors. McMule’s defining theoretical choice is to keep physical lepton masses finite, so collinear divergences are regulated by $\ell p\to \ell p$2 and $\ell p\to \ell p$3, while soft singularities are handled with dimensional regularisation and the all-order subtraction framework FKS$\ell p\to \ell p$4 (Banerjee et al., 2020, Engel, 2022).

This QED-specific strategy is central. Because the observables of interest are often not collinear safe, the massless approximation is not merely inconvenient but physically inadequate. McMule therefore subtracts only the soft singular structure, exploiting Yennie–Frautschi–Suura-type exponentiation and organising the perturbative expansion into finite $\ell p\to \ell p$5-, $\ell p\to \ell p$6-, and $\ell p\to \ell p$7-body contributions. The framework’s published derivations emphasize that this is simpler than massless QCD subtraction but numerically delicate because the remaining mass logarithms are large and experimentally relevant (Banerjee et al., 2020, Ulrich, 2023).

A second core method is massification, used when exact fully massive two-loop amplitudes are unavailable. In its schematic form,

$\ell p\to \ell p$8

so the leading small-$\ell p\to \ell p$9 behaviour is reconstructed from a massless hard amplitude, universal collinear factors, and a soft factor (Broggio et al., 2022). This method is used, for example, in mixed photonic two-loop sectors of $e^-e^-\to e^-e^-$0 scattering and in other $e^-e^-\to e^-e^-$1 QED reactions where the exact massive result is not known (Broggio et al., 2022, Engel, 2022).

A third pillar is next-to-soft stabilisation. Real-virtual one-loop amplitudes become numerically unstable when the emitted photon is soft or pseudo-collinear, particularly in small-$e^-e^-\to e^-e^-$2 kinematics. McMule addresses this by replacing the full matrix element in dangerous regions with the leading and next-to-leading soft expansion, a strategy first deployed for NNLO Bhabha scattering and later generalised to other processes, including polarised ones (Banerjee et al., 2021, Broggio et al., 2022, Kollatzsch et al., 2022). This is one of the framework’s most characteristic numerical innovations: it is not merely a rescue system, but a process-independent methodology anchored in the one-loop extension of the Low–Burnett–Kroll theorem for massive fermions (Engel, 2022).

3. Software architecture and computational model

McMule is implemented in Fortran 95, with Python support tools and a modular architecture that mirrors its perturbative decomposition. The codebase contains modules for global definitions, utility functions, scalar-integral interfaces, phase-space generation, matrix elements, user-defined observables, VEGAS integration, and finite integrands for LO, NLO, and NNLO sectors (Banerjee et al., 2020). The process library is organised into process groups, generic processes, and perturbative pieces, so a given physical channel is selected by both its flavour assignment and the requested perturbative component (Banerjee et al., 2020).

From the user perspective, the framework is built around measurement functions rather than hard-coded observables. Cuts, histogram definitions, and derived quantities are supplied in the user layer, while the framework handles phase-space sampling, subtraction, and perturbative bookkeeping (Banerjee et al., 2020). This is why McMule papers repeatedly describe it as a tool for “arbitrary infrared-safe observables” rather than a catalog of pre-tabulated distributions (Broggio et al., 2022, Kollatzsch et al., 2022).

For a substantial period, McMule was explicitly a Monte Carlo integrator, not a full event generator. That distinction matters: its output was binned differential cross sections obtained by direct integration over finite subtracted integrands (Rocco, 2023). Later work reports that McMule “gained the ability to generate events at NNLO rather than just pre-defined differential distributions,” using cellular resampling during generation to mitigate subtraction-induced negative weights (Ulrich, 7 Jan 2025). In the demonstrated muon-decay application, this reduces the negative-weight fraction from about $e^-e^-\to e^-e^-$3 to $e^-e^-\to e^-e^-$4, and the worst negative-weight outlier relative to the average event weight from about $e^-e^-\to e^-e^-$5 to $e^-e^-\to e^-e^-$6 (Ulrich, 7 Jan 2025). This marks a genuine transition from pure parton-level integration toward detector-facing simulation workflows.

4. Process implementations and benchmark calculations

The framework’s best-known benchmark is muon–electron elastic scattering for MUonE. McMule provides the first complete NNLO QED calculation for $e^-e^-\to e^-e^-$7, including leptonic, hadronic, and photonic contributions, with full mass dependence in the fermionic and line-specific photonic sectors and massification for the mixed double-virtual four-point topologies (Broggio et al., 2022). The associated implementation was described as the first NNLO result for a $e^-e^-\to e^-e^-$8 QED process with two distinct nonzero external masses (Rocco, 2023). Phenomenologically, the MUonE-like studies show NLO corrections up to about $e^-e^-\to e^-e^-$9 and NNLO corrections around $e^+e^-\to e^+e^-$0, with the latter already comparable to or larger than the hadronic signal in parts of phase space; this established NNLO as necessary but not yet sufficient for a $e^+e^-\to e^+e^-$1 ppm extraction (Rocco, 2023).

McMule also supplies full NNLO QED for unpolarised Møller scattering. In the PRad II application, the implementation includes photonic, leptonic, and non-perturbative hadronic vacuum-polarisation corrections, hard-photon radiation, and electron-mass effects throughout, with only the photonic two-loop amplitude treated in a controlled small-$e^+e^-\to e^+e^-$2 approximation (Banerjee et al., 2021). The main phenomenological conclusion is that integrated NNLO effects are modest, but differential NNLO corrections become essential in bins where the LO contribution vanishes and radiative channels provide the first nontrivial support (Banerjee et al., 2021).

For Bhabha scattering, McMule’s importance is methodological as much as phenomenological. The next-to-soft stabilisation program was developed precisely because NNLO Bhabha real-virtual contributions were too unstable for a fully differential fixed-order Monte Carlo without it (Banerjee et al., 2021). Once implemented, McMule delivered the first fully differential photonic NNLO Bhabha results, with NLO and NNLO corrections benchmarked against BABAYAGA-like predictions at the $e^+e^-\to e^+e^-$3 level for the observables studied (Banerjee et al., 2021).

For annihilation channels, McMule now covers lepton pair production, diphoton production, and related luminosity-standard processes. The $e^+e^-\to e^+e^-$4 implementation includes full NLO QED, dominant NNLO QED ISR terms, NLO electroweak effects, beam polarisation, full lepton masses, and hadronic vacuum polarisation, with Belle II-oriented studies of $e^+e^-\to e^+e^-$5-pair production and spin asymmetries (Kollatzsch et al., 2022, Gogniat et al., 14 May 2025). The later $e^+e^-\to e^+e^-$6 implementation is explicitly described as completing McMule’s NNLO coverage of the most important low-energy $e^+e^-\to e^+e^-$7 processes, with applications to luminosity measurements up to a few GeV (Engel et al., 28 Dec 2025).

5. Experiment-facing applications

A major experiment-facing deployment is the precision simulation of ordinary polarized muon decay for endpoint searches at MEG II and Mu3e. In this context McMule computes both the signal $e^+e^-\to e^+e^-$8 and the Standard-Model background $e^+e^-\to e^+e^-$9 fully differentially, with NLO QED corrections for the signal and, for the background, leading weak and hadronic effects, full NNLO QED, higher-order collinear logarithms, and NNLL soft-photon resummation in the endpoint region (Gurgone, 2023). The results were used to implement a new positron event generator in the MEG II software, interfaced to Geant4 for detector transport and full reconstruction (Gurgone et al., 2022). The associated sensitivity study quotes estimated branching-ratio reach of about $ee\to \ell\ell$0 for MEG II and $ee\to \ell\ell$1 for Mu3e for a massless ALP with $ee\to \ell\ell$2 coupling, and states that using only an NLO background prediction would have made the search theory-limited rather than statistics-limited (Gurgone, 2023).

McMule has also become a benchmark tool in comparisons of low-energy $ee\to \ell\ell$3 Monte Carlo programs. In the review of tools for low-energy hadronic cross sections, it is presented as the high-accuracy fixed-order reference against which AfkQed, BabaYaga@NLO, KKMC, MCGPJ, Phokhara, and Sherpa are compared (Aliberti et al., 2024). The review makes the distinction explicit: McMule is a parton-level integrator with NNLO accuracy for several non-radiative $ee\to \ell\ell$4 channels and NLO accuracy for associated radiative processes, rather than a resummed detector-level event generator (Aliberti et al., 2024). In leptonic channels this role is particularly strong; in pion channels the same review emphasizes that McMule currently includes ISC and VPC with a customisable pion form factor, but not the full structure-dependent final-state and mixed corrections (Aliberti et al., 2024).

For Belle II, McMule underpins two different programs. One concerns $ee\to \ell\ell$5, where it is used to evaluate radiative corrections to spin-dependent asymmetries relevant for $ee\to \ell\ell$6, including a cut-dependent modification of the asymmetry coefficient $ee\to \ell\ell$7 that restores cancellation of the dominant $ee\to \ell\ell$8 term in realistic acceptances (Gogniat et al., 14 May 2025). The other is the broader $ee\to \ell\ell$9 implementation with dominant NNLO QED and NLO electroweak effects, where the framework’s polarised next-to-soft formalism enables stable real-virtual integration in Belle II kinematics (Kollatzsch et al., 2022).

6. Limitations, approximations, and post-QED extensions

McMule’s literature is unusually explicit about its limitations. The most persistent is that fixed-order NNLO is not always enough. In MUonE kinematics, the small-angle region remains dominated by logarithmically enhanced radiation, so beyond-NNLO effects, resummation, or a parton shower are described as necessary next steps (Rocco, 2023, Broggio et al., 2022). In pair-production asymmetries, the available “dominant NNLO” implementations omit some FSR and IFI sectors, so fixed-order predictions are incomplete for the most precision-sensitive asymmetries (Kollatzsch et al., 2022). In pion channels, current McMule implementations do not yet include all hadronic-structure effects beyond one-photon exchange contributions (Aliberti et al., 2024).

There are also controlled approximation layers internal to the framework. In $\pm$0 scattering, the only approximation in the quoted NNLO prediction is the neglect of terms polynomially suppressed in the electron-mass expansion of the mixed double-virtual contribution (Broggio et al., 2022). In Møller and Bhabha scattering, massification similarly drops $\pm$1-suppressed photonic terms while retaining the logarithmically enhanced and constant pieces (Banerjee et al., 2021, Banerjee et al., 2021). These are not ad hoc shortcuts but explicit design decisions, typically supported by dedicated validation and uncertainty studies.

Recent work extends McMule beyond pure QED in two directions. One is the inclusion of low-energy electroweak effects within LEFT, used for parity-violating Møller scattering and the MOLLER asymmetry $\pm$2; in that setup McMule combines NNLO QED with NLO electroweak corrections and LL Wilson-coefficient running (Kollatzsch et al., 23 Jul 2025, Kollatzsch, 10 Mar 2026). The other is disperon QED, which introduces dispersive effective propagators to incorporate hadronic vacuum polarisation and pion form factors directly inside loop amplitudes, with an implementation for $\pm$3 in McMule (Fang et al., 11 Dec 2025, Kollatzsch, 10 Mar 2026). A notable result of the electroweak extension is that the model dependence of non-perturbative $\pm$4-$\pm$5 mixing in MOLLER becomes sub-percent for the preferred low-scale RGE-improved setup, so it is not expected to be the dominant theory limitation there (Kollatzsch, 10 Mar 2026).

Taken together, these developments define McMule less as a single program than as a precision-computation platform: originally a fixed-order QED framework for massive leptons, then an NNLO event generator, and now an infrastructure into which low-energy electroweak EFT and data-driven hadronic loop effects are being incorporated (Ulrich, 7 Jan 2025, Kollatzsch, 10 Mar 2026). The central scientific role has remained stable throughout: to provide experimentally deployable, fully differential predictions in regimes where finite masses, radiative tails, and fiducial cuts are inseparable from the physics.