# Implementing new processes in McMule

In this section we will discuss how new processes can be added to McMule. Not all of the points below might be applicable to any particular process. Further, all points are merely guidelines that could be deviated from if necessary as long as proper precautions are taken.

As an example, we will discuss how Møller scattering \(e^-e^-\to e^-e^-\) could be implemented.

A new process group may need to be created if the process does not fit any of the presently implemented groups. This requires a new folder with a makefile as well as modifications to the main makefile as discussed in Section Creating a new process group.

In our case, \(ee\to ee\) does not fit any of the groups, so we create a new group that we shall call

`ee`

.Calculate the tree-level matrix elements needed at LO and NLO: \(\mathcal{M}_{n}^{(0)}\) and \(\mathcal{M}_{n+1}^{(0)}\). This is relatively straightforward and – crucially – unambiguous as both are finite in \(d=4\). We will come back to an example calculation in Section Example calculations in Mathematica.

A generic matrix element file is needed to store ‘simple’ matrix elements as well as importing more complicated matrix elements. Usually, this file should not contain matrix elements that are longer than a few dozen or so lines. In most cases, this applies to \(\mathcal{M}_{n}^{(0)}\).

After each matrix element, the PID needs to be denoted in a comment. Further, all required masses as well as the centre-of-mass energy, called

`scms`

to avoid collisions with the function \({\tt s(pi,pj)}=2{\tt pi}\cdot{\tt pj}\), need to be calculated in the matrix element to be as localised as possible.In the case of Møller scattering, a file

`ee/ee_mat_el.f95`

will contain \(\mathcal{M}_{n}^{(0)}\). For example, \(\mathcal{M}_{n}^{(0)}\) is implemented there as shown in Listing 21.

```
FUNCTION EE2EE(p1, p2, p3, p4)
!! e-(p1) e-(p2) -> e-(p3) e-(p4)
!! for massive (and massless) electrons
implicit none
real(kind=prec), intent(in) :: p1(4), p2(4), p3(p4), p4(4), ee2ee
real(kind=prec) :: den1, den2, t, scms, m2
t = sq(p1-p3) ; scms = sq(p1+p2) ; m2 = sq(p1)
den1 = sq(p1-p3) ; den2 = sq(p1-p4)
ee2ee=(8**m2**2 - 8*m2*scms + 2*s**2 + 2*scms*t + t**2)/den1**2
ee2ee=ee2ee+2*(12*m2**2 - 8*m2*scms + scms**2) / den1 / den2
ee2ee=ee2ee+(24*m2**2 + scms**2 + t**2 - 8*m2*(s + t))/den2**2
ee2ee = ee2ee * 128*pi**2*alpha**2
END FUNCTION
```

Further, we need an interface file that implements the abstract interface

`partInterface`

(cf. Section Particle string framework). In our case this is called`ee/ee.f95`

, as shown in Listing 22.In most cases, the

`particle`

string framework can automatically define all soft limits and NTS matrix elements. In rare cases, this needs to be done manually in`ee/ee.f95`

. [1]

```
FUNCTION EE2EE_part(p1, p2, p3, p4)
!! e-(p1) e-(p2) -> e-(p3) e-(p4)
!! for massive (and massless) electrons
implicit none
real(kind=prec) :: p1(4), p2(4), p3(p4), p4(4)
type(particles) :: ee2ee_part
ee2ee_part = parts((/part(p1, 1, 1), part(p2, 1, 1), part(p3, 1, -1), part(p4, 1, -1)/))
END FUNCTION
```

Because \(\mathcal{M}_{n+1}^{(0)}\) is border-line large, we will assume that it will be stored in an extra file,

`ee/ee2eeg.f95`

. The required functions are to be imported in`ee/ee_mat_el.f95`

.Calculate the one-loop virtual matrix element \(\mathcal{M}_{n}^{(1)}\), renormalised in the OS scheme. Of course, this could be done in any regularisation scheme. However, results in McMule shall be in the FDH (or equivalently the FDF) scheme. Divergent matrix elements in McMule are implemented as \(c_{-1}\), \(c_0\), and \(c_1\)

\[\mathcal{M}_{n}^{(1)} = \frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)}\Bigg( \frac{c_{-1}}\epsilon + c_0 + c_1\epsilon+\mathcal{O}(\epsilon^2) \Bigg)\,.\]For \(c_{-1}\) and \(c_0\) this is equivalent to the conventions employed by Package-X [20] up to a factor \(1/16\pi^2\). While not strictly necessary, it is generally advisable to also include \(c_{-1}\) in the Fortran code.

For NLO calculations, \(c_1\) does not enter. However, we wish to include Møller scattering up to NNLO and hence will need it sooner rather than later anyway.

In our case, we will create a file

`ee/ee_ee2eel.f95`

, which defines a functionFUNCTION EE2EEl(p1, p2, p3, p4, sing, lin) !! e-(p1) e-(p2) -> e-(p3) e-(p4) !! for massive electrons implicit none real(kind=prec), intent(in) :: p1(4), p2(4), p3(p4), p4(4) real(kind=prec) :: ee2eel real(kind=prec), intent(out), optional :: sing, lin ... END FUNCTION

The function shall return \(c_0\) in

`ee2eel`

and, if`present`

\(c_{-1}\) and \(c_1\) in`sing`

and`lin`

.At this stage, a new subroutine in the program

`test`

with reference values for all three matrix elements should be written to test the Fortran implementation. This is done by generating a few points using an appropriate phase-space routine and comparing to as many digits as possible using the routine`check`

.In our case, we would construct a subroutine

`TESTEEMATEL`

as shown in Listing 23.

```
SUBROUTINE TESTEEMATEL
implicit none
real (kind=prec) :: x(2),y(5)
real (kind=prec) :: single, finite, lin
real (kind=prec) :: weight
integer ido
call blockstart("ee matrix elements")
call initflavour("mu-e", mys=40000._prec)
scms = 40000.
musq = mm**2
x = (/0.75,0.5/)
call psx2(x,q1,me,q2,me,q3,me,q4,me,weight)
call check("ee2ee:",ee2ee (q1,q2,q3,q4),5684.9581122306845)
mat = ee2eel(q1,q2,q3,q4)
call check("ee2eel fin:",mat,221816.84997000592)
call initflavour("mu-e", mys=40000.)
y = (/0.3,0.6,0.8,0.4,0.9/)
call psx3_fks(y,q1,me,q2,me,q3,me,q4,me,q5,weight)
call check("ee2eeg",ee2eeg(q1,q2,q3,q4,q5),434.98382653487744)
call blockend(3)
END SUBROUTINE
SUBROUTINE TESTMEEVEGAS
xinormcut1 = 0.2
xinormcut2 = 0.3
call blockstart("Moller VEGAS test")
call test_INT('ee2ee0', sigma_0, 2,10,10, NaN)
call test_INT('ee2eeF', sigma_0, 2,10,10, NaN)
call test_INT('ee2eeR', sigma_1, 5,10,10, NaN)
call blockend(3)
END SUBROUTINE
```

In addition, McMule provides built-in routines for testing the convergence of real-emission matrix elements to the corresponding soft limits, for ever smaller photon energies.

In our case, we would construct a subroutine

SUBROUTINE TESTEESOFTN1 implicit none real(kind=prec) y0(5) call blockstart("e-e \xi->0") call initflavour("muone") xinormcut1 = 0.3 y0 = (/0.01,0.6,0.8,0.999,0.01/) call test_softlimit(y0, ["ee2eeR"]) END SUBROUTINE

where

`test_soft_limit`

compares the real matrix element (`ee2eeR`

) with its soft limit implemented in`ee/ee.f95`

. The comparison starts at the energy defined by the phase-space point`y0`

and proceeds with ever smaller photon energies. A flavour (`muone`

) as well as`xinormcut1`

are required in order to complete the phase-space generation.Define a default observable in

`user`

for this process. This observable must be defined for any`which_piece`

that might have been defined and test all relevant features of the implementation such as polarisation if applicable.Add the matrix elements to the integrands defined in

`integrands.f95`

. This is done by adding a new`case`

corresponding to the new`which_piece`

in the`initpiece()`

.for a IR-finite, non-radiative piece (i.e. LO but also VP), add

case('eb2eb0') call set_func('00000000', eb2eb) ps => psx2 ; fxn => sigma_0 nparticle = 4 ; ndim = 2 masses(1:4) = (/ Me, Me, Me, Me /)

which adds a

`which_piece`

`ee2ee0`

that is calculated using the matrix element`ee2ee`

. The phase space is generated with`psx4()`

and integrated using`sigma_0()`

(no subtraction). The process involves 4 particles and, since it is a \(2\to2\) process, the integration is two-dimensional. The masses of the involved particles are all`Me`

.for pieces with an IR cancellation between real and virtual corrections, i.e. calculations involving photon loops, we need to specify

`xieik1`

(at one-loop) and/or`xieik2`

(at two-loop)case('eb2ebF') call set_func('00000000', eb2ebf) ps => psx2 ; fxn => sigma_0 nparticle = 4 ; ndim = 2 masses(1:4) = (/ Me, Me, Me, Me /) xieik1 = xinormcut*(1.-(2*me)**2/scms)

One needs to take care that \(\xi_c\) is properly normalised. The user will enter a value from 0 to 1 which needs to be matched to \(\xi_c\) as defined in (12)

for real corrections we need to use a subtracting integrand, i.e.

`sigma_1()`

for single-real and`sigma_2()`

for double-real corrections.case('eb2ebR') call set_func('00000000', eb2ebg) call set_func('00000001', eb2eb) call set_func('11111111', eb2eb_part) ps => psx3_fks ; fxn => sigma_1 nparticle = 5 ; ndim = 5 masses(1:5) = (/ Me, Me, Me, Me, 0._prec /) xicut1 = xinormcut*(1.-(2*me)**2/scms)

Additionally to the real matrix element

`eb2ebg`

, we also specified the reduced matrix element`eb2eb`

and the`particle`

string function`eb2eb_part`

. Note further changes to the number of particles and phase space dimension to accommodate the extra photon.

Additionally to these required parameters, there are number of optional parameters such as

`symmfac`

(which is set to 2 for \(e^-e^-\to e^-e^-\) because of the two indistinguishable final state particles),`polarised`

(whether to consider the process polarised), and`softCut`

and`collCut`

. For a full list of parameters, see Section Optional parameters for integrands.Once

`integrands`

are implemented, a second test routine should be written that runs short integrations against a reference value. Because`test_INT`

uses a fixed random seed, this is expected to be possible very precisely. Unfortunately, COLLIER [6] might produce slightly different results on different machines. Hence, integrands involving complicated loop functions are only required to agree up to \(\mathcal{O}(10^{-8})\).After some short test runs, it should be clear whether new phase-space routines are required. Add those, if need be, to

`phase_space`

as described in Section Phase-space generation.Per default the stringent soft cut, that may be required to stabilise the numerical integration (cf. Section Implementation of FKS schemes), is set to zero. Study what the smallest value is that still permits integration.

Perform very precise \(\xi_{c}\) independence studies. Tips on how to do this can be found in Section Study of \xi_{c} dependence.

At this stage, the NLO calculation is complete and may, after proper integration into McMule and adherence to coding style has been confirmed, be added to the list of McMule processes in a new release. Should NNLO precision be required, the following steps should be taken

Calculate the real-virtual and double-real matrix elements \(\mathcal{M}_{n+1}^{(1)}\) and \(\mathcal{M}_{n+2}^{(0)}\) and add them to the test routines as well as integrands.

Prepare the \(n\)-particle contribution \(\sigma_n^{(2)}\). In a pinch, massified results can be used also for \(\hat{\mathcal{E}}(\xi_{c})\mathcal{M}_{n}^{(1)}\) though of course one should default to the fully massive results.

Study whether the pre-defined phase-space routines are sufficient. Even if it was possible to use an old phase-space at NLO, this might no longer work at NNLO due to the added complexity. Adapt and partition further if necessary, adding more test integrations in the process.

Perform yet more detailed \(\xi_{c}\) and soft cut analyses.

In the following we comment on a few aspects of this procedure such as the \(\xi_{c}\) study (Section Study of \xi_{c} dependence), the calculation of matrix elements (Section Example calculations in Mathematica), and a brief style guide for McMule code (Section Coding style and best practice).

## Creating a new process group

Adding Møller scattering to McMule, the example discussed above, requires the addition of a new process group `ee`

.
For this we create a new folder in McMule called `ee`

containing a meson.build (Listing 24), a `mat_el`

file (`ee/ee_mat_el.f95`

, Listing 25) and a module file (`ee/ee.f95`

, Listing 26).
Finally, the name of the process group needs to be added to the `GROUPS`

and `WGROUPS`

variables of the makefile.

```
# add files that require compiler optimisation here
src_opt = []
# add all other files here
src = [
'ee_ee2eel.f95',
'ee_ee2eeg.f95',
'ee_mat_el.f95',
'ee.f95'
]
lib_opt = static_library(
'ee_opt', src_opt, dependencies : mcmule_pre_dep,
fortran_args : [opt_level1]
)
dep_opt = declare_dependency(
link_with : lib_opt,
include_directories : lib_opt.private_dir_include()
)
groups_lib += static_library(
'ee', [src, lib_opt.extract_all_objects(recursive : true)],
dependencies : [mcmule_pre_dep, dep_opt]
)
groups_incl += lib_opt.private_dir_include()
src_test += '../ee/ee_test.f95'
foreach i : src + src_opt
src_groups += 'ee/' + i
endforeach
```

```
!!!!!!!!!!!!!!!!!!!!!!!!!
MODULE EE_MAT_EL
!!!!!!!!!!!!!!!!!!!!!!!!!
use functions
use ee_ee2eel, only: ee2eel
use ee_ee2eeg, only: ee2eeg
implicit none
contains
FUNCTION EE2EE(p1,p2,q1,q2)
!! e-(p1) e-(p2) -> e-(q1) e-(q2)
!! for massive electrons
...
END FUNCTION EE2EE
FUNCTION EE2EEF(p1,p2,q1,q2)
!! e-(p1) e-(p2) -> e-(q1) e-(q2)
!! massive electrons
real(kind=prec) :: p1(4),p2(4),q1(4),q2(4)
real(kind=prec) :: ee2eef, mat0, Epart
Epart = sqrt(sq(p1+p2))
mat0 = ee2ee(p1,p2,q1,q2)
ee2eef = ee2eel(p1,p2,q1,q2) + alpha / (2 * pi) * mat0 * &
Ieik(xieik1, Epart, parts((/ part(p1, 1, 1), part(p2, 1, 1), &
part(q1, 1, -1), part(q2, 1, -1) /)))
END FUNCTION EE2EEF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE EE_MAT_EL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
```

```
!!!!!!!!!!!!!!!!!!!!!!!!!
MODULE EE
!!!!!!!!!!!!!!!!!!!!!!!!!
use functions
use phase_space, only: ksoft, ksoftA, ksoftB
use ee_mat_el
implicit none
contains
FUNCTION EE2EE_part(p1, p2, p3, p4)
!! e-(p1) e-(p2) --> e-(p3) e-(p4)
!! both massive and massless electrons
real (kind=prec) :: p1(4),p2(4),p3(4),p4(4)
type(particles) :: ee2ee_part
ee2ee_part = parts((/part(p1, 1, 1), part(p2, 1, 1), part(p3, 1, -1), part(p4, 1, -1)/))
END FUNCTION EE2EE_part
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE EE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
```

## Study of \(\xi_{c}\) dependence

When performing calculations with McMule, we need to check that the dependence of the unphysical \(\xi_{c}\) parameter introduced in the FKS scheme (cf. Appendix The FKS^2 scheme) actually drops out at NLO and NNLO. In principle it is sufficient to do this once during the development phase. However, we consider it good practice to also do this (albeit with a reduced range of \(\xi_{c}\)) for production runs.

Because the \(\xi_{c}\) dependence is induced through terms as \(\xi_{c}^{-2\epsilon}/\epsilon\), we know the functional dependence of \(\sigma^{(\ell)}_{n+j}\). For example, at NLO we have

where \(\xi_{c}\) independence of \(\sigma^{(1)}\) of course requires

At NNLO we have

We require

\[a_{0,i} + a_{1,i} + a_{2,i} = 0\]

for \(i=1,2\). However, the IR structure allows for an even stronger statement for the \(a_{j,2}\) terms

\[a_{0,2} = a_{2,2} = -\frac{a_{1,2}}2\,.\]

Of course we cannot directly calculate any of the \(a_{1,i}\) or \(a_{2,i}\) because we use numerical integration to obtain the \(\sigma^{(\ell)}_{n+j}\).
Still, knowing the coefficients can be extremely helpful when debugging the code or to just quantify how well the \(\xi_{c}\) dependence vanishes.
Hence, we use a fitting routine to fit the Monte Carlo results *after* any phase-space partitioning has been undone.
Sometimes non of this is sufficient to pin-point the source of a problem to any one integrand.
However, if the goodness of, for example, \(\sigma^{(2)}_{n+2}(\xi_{c})\) is much worse than the one for \(\sigma^{(2)}_{n+1}(\xi_{c})\), a problem in the double-real corrections can be expected.

## Example calculations in Mathematica

A thorough understanding of one-loop matrix elements is crucial for any higher-order calculation. In McMule, one-loop matrix elements either enter as the virtual contribution to NLO corrections or the real-virtual contribution in NNLO calculations. In any case, a fast numerical routine is required that computes the matrix element.

We perform all one-loop calculations in FDF as this is arguably the simplest scheme available. For theoretical background, we refer to [25] and references therein.

We use Qgraf for the diagram generation.
Using the in-house Mathematica package `qgraf`

we convert Qgraf’s output for manipulation with Package-X [20].
This package is available on request through the McMule collaboration

An example calculation for the one-loop calculation of \(\mu\to\nu\bar\nu e\gamma\) can be found in Listing 27. Of course this example can be made more efficient by, for example, feeding the minimal amount of algebra to the loop integration routine.

When using `qgraf`

for fdf some attention needs to be paid when considering diagrams with closed fermion loops.
By default, `qgraf.wl`

evaluates these traces in \(d\) dimensions.
`RunQGraf`

has an option to keep this from happening.

```
<<qgraf.wl
onshell = {
p.p -> M^2, q.q -> m^2, p.q -> s/2
};
A0 = (4GF/Sqrt[2]) "diag1"/.RunQGraf[{"mum"},{"nu","elm"},0] //. {
line[_, x_] -> x, p1->p, q1->p-q, q2->q, _δZ | δm -> 0
};
A1 = pref /. RunQGraf[{"mum"},{"nu","elm"},1] //. {
line[_, x_] -> x, p1->p, q1->p-q, q2->q, _δZ | δm -> 0
};
M0=Block[{Dim=4},Simplify[Contract[
1/2 Z2[m] Z2[M] FermionSpinSum[
A0 /. DiracPL -> (Dirac1 - Z5 γ5)/2,
A0 /. DiracPL -> (Dirac1 + Z5 γ5)/2
]
]] /. onshell]/.{
Z2[M_] -> 1 + (α/(4π)) (-3/(2ε)-5/2 + 3/2 Log[M^2/Mu^2]),
Z5 -> 1 - (α/(4π))
};
M1=Block[{Dim=4},Simplify[Contract[
1/2 FermionSpinSum[
A1/.γ.k1 -> γ. 4[k1]+I γ5 μ,
A0
]
] /. onshell /. {
μ^n_ /; EvenQ[n] -> μ2^(n/2), μ -> 0
}/.{
4[k1]. 4[k1] -> k1.k1 + μ2, 4[k1] -> k1
}]];
M1bare = Simplify[KallenExpand[LoopRefine[LoopRelease[
Pro2LoopIntegrate[
Coefficient[M1, μ2, 0]/(16 π^2)
]
+ μIntegrate[
Coefficient[M1, μ2, 1]/(64 π^3),
1
],
onshell
]]] /. e -> Sqrt[4 πα]];
```

There is a subtlety here that only arise for complicated matrix elements.
Because the function Package-X uses for box integrals, `ScalarD0IR6()`

, is so complicated, no native Fortran implementation exists in McMule.
Instead, we are defaulting to COLLIER [6] and should directly evaluate the finite part of the `PVD`

function above.
The same holds true for the more complicated triangle functions.
In fact, only the simple `DiscB()`

and `ScalarC0IR6()`

are natively implemented without need for external libraries.
For any other functions, a judgement call is necessary of whether one should `LoopRefine`

the finite part in the first place.
In general, if an integral can be written through logarithms and dilogs of simple arguments (resulting in real answers) or `DiscB()`

and `ScalarC0IR6()`

, it makes sense to do so.
Otherwise, it is often easier to directly link to COLLIER.

## Coding style and best practice

A large-scale code base like McMule cannot live without some basic agreements regarding coding style and operational best practice. These range from a (recommended but not enforced) style guide over the management of the git repository to how to best run McMule in development scenarios. All aspects have been discussed within the McMule collaboration.

Fortran code in McMule is (mostly) written in accordance with the following style guide. If new code is added, compliance would be appreciated but deviation is allowed if necessary. If in doubt, contact any member of the McMule collaboration.

Indentation width is two spaces. In Vim this could be implemented by adding the following to

`.vimrc`

autocmd FileType fortran set tabstop=8 softtabstop=0 expandtab shiftwidth=2 smarttab

Function and subroutine names are in all-upper case.

A function body is not indented beyond its definition.

When specifying floating point literals specify the precision when possible, i.e.

`1._prec`

.Integrands should have

`ndim`

specified.Internal functions should be used where available.

Masses and other kinematic parameters must be calculated in the matrix elements as local variables; using the global parameters

`Mm`

and`Me`

is strictly forbidden.These rules also hold for matrix elements.

For python code, i.e. pymule as well as the analysis code, PEP8 compliance is strongly encouraged with the exception of `E231`

(Missing whitespace after `,`

, `;`

, and `:`

), `E731`

(Do not assign a lambda expression, use a `def`

) as well, in justified cases, i.e. if required by the visual layout, `E272`

(Multiple spaces before keyword), and `E131`

(Continuation line unaligned for hanging indent).

McMule uses a git repositories for version management.
Development usually happens on feature branches that are merged into the `devel`

branch semi-frequently by the McMule collaboration after sufficient vetting was performed.
Finally, once a project has been finished, the `devel`

branch gets merged into the `release`

branch that is to be used by McMule’s users.

In general, developers are encouraged to not commit wrong or unvetted code though this can obviously not be completely avoided in practice. To avoid uncontrollable growth of the git repository, large files movements are strongly discouraged. This also means that matrix elements should not be completely overhauled barring unanimous agreement. Instead, developers are encouraged to add a new matrix element file and link to that instead.

Even when running McMule for development purposes the usage of menu files is strongly encouraged because the code will do its utmost to automatically document the run by storing the git version as well as any modification thereof. This allows for easy and unique reconstruction of what was running. For production runs this is not optional; these must be conducted with menu files after which the run folder must be stored with an analysis script and all data on the AFS as well as the user file library to ensure data retention.