H1jet, a fast program to compute transverse momentum distributions

We present H1jet, a fast code that computes the total cross section and differential distribution in the transverse momentum of a colour singlet. In its current version, the program implements only leading-order $2\to 1$ and $2\to 2$ processes, but could be extended to higher orders. We discuss the processes implemented in H1jet, give detailed instructions on how to implement new processes, and perform comparisons to existing codes. This tool, mainly designed for theorists, can be fruitfully used to assess deviations of selected new physics models from the Standard Model behaviour, as well as to quickly obtain distributions of relevance for Standard Model phenomenology.


Introduction
After the discovery of the Higgs boson in 2012 [2,3], one of the most urgent tasks of the Large Hadron Collider (LHC) is the characterisation of the Higgs sector, in order to shed light on the exact mechanism for electroweak symmetry breaking. In particular, Higgs production data, either in the form of signal strengths [4,5] or cross section measurements [6,7], offer powerful constraints on Higgs anomalous couplings.
However, it has been pointed out that inclusive Higgs production through gluon fusion, the one with the largest rate, is not able to discriminate effectively between the Standard Model (SM) and another theory giving the same effective coupling between the gluons and the Higgs. In fact, top quarks running in loops give a dimension-6 effective interaction between the incoming gluons and the Higgs, with exact top-mass effects giving tiny corrections. Therefore, in many theories, the strength of dimension-6 contact gluon-gluon-Higgs interactions can conspire with an anomalous top Higgs Yukawa coupling to give exactly the same cross section for Higgs production as the SM [8,9,10].
There are essentially two ways of solving this problem. One is to put a direct constraint on the top Yukawa coupling by the observation of top quarks in association with the Higgs [11,12,13]. The other is to break the top loop by looking e.g. at Higgs production at large transverse momentum, where the Higgs recoils against a hard jet [8,9,10]. Both are indirect probes of new physics effects. The latter is more difficult experimentally, in that it relies on appreciating small deviations from the SM in the shape of the Higgs transverse momentum distribution, in a region where the phase space closes. It nevertheless can give a direct access to new physics coupling the gluons to the Higgs through loops.
Higgs sector aside, production of colour singlets at high transverse momentum is commonly used as a probe of new physics. A relevant example is the production of monojets, which can recoil either against dark matter, or against a SM particle decaying into invisible particles (see e.g. refs. [14,15]).
Theoretical predictions for the transverse momentum distribution of a colour singlet, both in the SM and beyond, can be currently obtained with Monte Carlo programs, such as MadGraph5_aMC@NLO [16] or SusHi [17,18]. These codes, although general, have the drawback of being quite slow. Also, interference terms between new physics and the SM, which carry information on the strength of new interactions, are difficult to extract from Monte Carlo event generators because they are not positive definite. The aim of this paper is to describe a method to obtain the transverse momentum spectrum of a colour singlet in a second or less, and its concrete implementation in the program H1jet. This program makes it possible to predict the effects of several models in a short amount of time. This in turn opens the way to devising more refined cut-based search strategies only for the models showing the largest deviations with respect to the SM.
More precisely, H1jet predicts the transverse momentum distribution of a colour singlet fully integrated over rapidity, and completely inclusive with respect to all coloured particles, i.e. the recoiling jets. Such an approximation is not too unrealistic, because the higher the transverse momentum, the more the colour singlet is central, and the more its decay products will be likely to pass the detector acceptance cuts.
The program is based on elementary analytic manipulations on the expression for the transverse momentum distribution. These make it possible to write the spectrum as a one dimensional integral, whose integrand is the product of an amplitude squared, which can be provided by ourselves or by the user, and a parton luminosity, which we extract from an external program. The relevant amplitudes can be either hard coded, or computed automatically and embedded in the program via a simple user interface.
H1jet already comes with a number of hard-coded processes and models. The main process is where the initial state consists of gluons and light quarks, and Higgs production procedes via quark loops. This process can be calculated in H1jet for different physics models including the SM, a CP-odd Higgs, a simplified SUSY model, and composite Higgs models with a single or multiple top-partners. In addition, the bb → H + jet and pp → Z + jet processes for the SM are implemented. Moreover, H1jet is very flexible and can be easily interfaced to use a custom user-specified process. The paper is organised as follows. In section 2 we briefly describe the method underlying H1jet. In section 3 we describe in detail how H1jet works in practice. In particular, we show how it can be installed and run, and present the features currently implemented. In section 4 we present a detailed comparison with the existing program SusHi for Higgs production both in the SM and beyond. In section 5 we explain how a user can implement a model of new physics inside H1jet. We choose axion-like-particle (ALP) production, giving rise to a monojet. We then describe how to obtain ALP transverse momentum distributions from the generation of the Feynman rules with FeynRules, to the calculation of the amplitude with FeynCalc and its subsequent interface with H1jet to obtain the ALP transverse momentum spectrum. Last, section 6 presents our conclusions.

The Method
Before explaining how the H1jet method works, it is instructive to consider first how to compute the Born cross section for producing a colour singlet X, e.g. a Higgs, of mass m X . This will also allow us to set the notation for the rest of the paper. We consider the 2 → 1 process p 1 p 2 → X, where p 1 and p 2 are the two incoming partons and p X is the momentum of the considered colour singlet. From momentum conservation we haveŝ There can be various partonic subprocesses that contribute to the production of the particle X. Let us denote with M ij the amplitude for the subprocess ij → X (e.g. gg → H), with i, j = g, q f ,qf , where f,f denotes quark or antiquark flavours. The Born partonic cross section for each subprocess iŝ The corresponding hadronic cross section is given by where L ij (τ, µ F ) is the partonic luminosity If we are able to obtain the luminosity L ij m 2 X /s, µ F , we are then able to obtain a numerical prediction for the cross section through a simple multiplication. There are indeed numerical tools that are able to compute, tabulate and interpolate luminosities with incredible efficiency, for instance the program HOPPET [19]. Through an interface with HOPPET, we are able to compute the Born cross section given the amplitudes M ij . This procedure is the same adopted in the program JetVHeto [20], that computes cross sections for colour singlets with a veto on additional jets.
A similar strategy can be devised to obtain a fast calculation of distributions in the transverse momentum of particle X. A non-zero transverse momentum for X is obtained via a generic 2 → 2 partonic process p 1 p 2 → p 3 X, where p 1 , p 2 , and p 3 are massless partons, and p X is the momentum of the colour singlet X. We wish to compute dσ/dp T , where p T is the transverse momentum of p X with respect to the beam axis. At Born level only, p T is also the transverse momentum of the recoiling jet originated by p 3 . The partonic subprocesses contributing to dσ/dp T are gg → gX, q fqf → gX, Without loss of generality, in the centre-of-mass frame of the partonic collision, we can parameterise momenta as follows where η is the rapidity of parton p 3 in the centre-of-mass frame. The partonic p T spectrum for the process initiated by partons ij is given by where E X = m 2 X + p 2 T cosh 2 η is the energy of the colour-singlet particle p X . The above equation selects two values of η, as follows The corresponding hadronic cross section reads Since eq. (9) gives two monotonic functions ofŝ for s > p T + m 2 X + p 2 T , varyingŝ in the allowed range spans all possible values of η in the range −η M < η < η M with This allows us to perform the η integration last, and obtain, after some manipulations, where again L ij is the partonic luminosity for the ij incoming channel as defined in eq. (5). If we are able to obtain the partonic luminosity L ij , say, from HOPPET, we can obtain the transverse momentum spectrum with a one-dimensional integration, which can be performed extremely quickly with a Gaussian numerical integrator. Summarising, by interfacing HOPPET with a code that provides amplitudes for 2 → 1 and 2 → 2 partonic subprocesses producing a colour singlet X, we are able to perform fast computations of total cross sections and transverse momentum spectra for X. In the following sections we describe our implementation of the method for Born processes. Note that, if one were able to perform the analytic integration over the phase space of final-state partons, the method can also be applied to higher-order cross sections and differential spectra.

User's Manual
This section describes the most important technical details of H1jet, including its installation and usage.

Installation
The source code of H1jet can be obtained from ref. [1]. The source code consists of a main directory H1jet with the following subdirectories: bin : contains the executable program h1jet after compilation, as well as the Python 3 helper scripts PlotH1jet.py and DressUserAmpCode.py.
src : source files.
The README.md file contains information on installation and usage.
In the main directory, the user needs to run the configure script: ./configure [options] It will attempt to find a Fortran compiler (gfortran or ifort), as well as the dependencies on the user's machine. A specific compiler and/or compiler flags can be selected with the options ./configure FC=<compiler> and ./configure FFLAGS=<flags>. H1jet has a number of external dependencies which it must be linked to: • LHAPDF [21]: Provides the PDF sets for H1jet.
• CHAPLIN [22]: For complex harmonic polylogarithms used to represent scalar integrals in loopinduced processes.
For the CHAPLIN library, it may be necessary to explicitly state the path to the library files with: ./configure LDFLAGS=-L/path/to/chaplin/lib To compile with a custom user interface: ./configure USERFILE=/path/to/custom/user_interface.f90 See Section 5 below for the implementation of custom user-specified amplitudes.
To install in a specific location: ./configure --prefix=/path/to/installation The default installation path is the main H1jet-directory.
The configure script will generate the Makefile.
To compile H1jet with the generated Makefile, run:

make [options]
This command takes the following options: make clean will delete all module and object files; make distclean will delete all module and object files as well as the executable h1jet.
After compilation, the bin-directory can then be added to the user's PATH environment variable. Alternatively, if the user has specified an installation directory with the --prefix option, the executable can be installed with: The executable h1jet can then be found in the bin-directory at the path specified by --prefix.

Usage
After compilation, H1jet can be run from the bin-directory with: H1jet will print out a brief summary of the settings and parameters used, as well as the Born cross section σ 0 , followed by the dσ/dp T and the integrated cross section σ(p T ) with a lower bound in p T for each p T bin. The following standard UNIX options are available: Display the help message along with all possible options.
-v, --version Display the version of the installed H1jet.
H1jet will display the requested information and then terminate. The output can be directed to a file with the option: -o, --out <file> Direct the output to <file>.
Default: standard output.
The physics process can be selected with: --proc <arg> Specify the process. Arguments: user User specified process. See Section 5 below for details on the implementation.
Depending on the process selected, there exists different relevant options.

General Options
The options listed here apply to all processes.

--pdf_name <arg>
Specify the PDF set name from LHAPDF. The specified PDF set must be available in the local installation of LHAPDF.
--scale_strategy <arg> Set the scale strategy, i.e. the dynamic µ = µ R = µ F value. Arguments: The mass M is given by option --mH for processes H and bbH, option --mZ for process Z, and option --mass for process user.
where µ is determined by the choice from --scale_strategy. Default: 0.5.

--nbins <value>
Number of histogram bins in the output of the transverse momentum distribution. Default: 400.

--log
Enables logarithmic x-axis of the histogram, i.e. logarithmic bins in p T . The option --ptmin must be set to a non-zero value in order to use this option, otherwise the program will quit with an error.
Note that the Yukawa couplings are given by where κ q are the dimensionless factors specified by the options --yt and --yb above, and v/ √ 2 is the vacuum expectation value of the Higgs field.
Note also that H1jet uses the G µ scheme for the all electroweak parameters [23]. Hence, the Higgs vacuum expectation value is given by To consider a CP-odd Higgs instead, it is necessary to select the following option: --cpodd Toggle for calculation of CP-odd Higgs.
The interaction between the CP-odd Higgs H and a SM quark q is: where the implementation in H1jet uses by defaultκ t = 1 andκ b = 0. Both parameters can be changed with the options --yt and --yb.
Here, both CP-even and CP-odd Higgs production are loop-induced processes. The amplitudes for 2 → 1 processes are taken from ref. [24]. For CP-even Higgs production in 2 → 2, the amplitudes are taken from ref. [25], and their interface is adapted from HERWIG 6 [26]. We have taken the CP-odd 2 → 2 amplitudes from ref. [9].
Top-partner. H1jet allows the calculation of Higgs production via loops of top partners in addition to top loops. To include a top-partner T in the quark loops, it is necessary to set the top-partner mass m T to a non-zero value by using the --mtp option.
The SM top Yukawa factor can be modified by the mixing angle, where κ t is the Yukawa factor set by option --yt.
The top-partner Yukawa factor will likewise be modified with κ T set by --ytp.
The top-partner specific options are: Default: 0 GeV.
The above is for a simplified composite Higgs model, where the compositeness scale f is set to infinity. The top-partner can also be considered in the explicit composite Higgs models of ref. [27], all with finite f . Four different models are implemented, M1 5 , M1 14 , M4 5 , and M4 14 , which modify the Yukawa coupling factors in the following way: where theκ's are the CP-odd couplings, and, For M1 5 and M1 14 , the option --sth2 sets the mixing angle θ L , while for M4 5 and M4 14 , the same option sets the angle θ R . The reason for this is that we want to reproduce the f → ∞ limit, where θ T = θ L , θ R depending on the chosen model. When needed, the angles θ L and θ R are derived one from the other by using the relation The composite Higgs top-partner model specific options are: --model <arg> Specify the top-partner model. Arguments: --imc1 <value> Imaginary part of the c 1 coefficient, Im(c 1 ). Default: 0.
If the option is not set, all couplings will be automatically computed in the limit f → ∞. If the option is set but no value is provided, the program will return a floating point exception.
Multiple top-partners models. H1jet makes it possible to include multiple top-partners in the particle loops. To do that, it will be necessary to specify an input file with the masses and Yukawa coupling factors for each particle running in the loop, including SM quarks. This can be done with the following option: The dimensionless Yukawa coupling factors κ q andκ q are respectively the CP-even and CP-odd couplings for a quark q, with the following Lagrangian: where m q is the mass of the quark.
The integer value specifying the loop approximations can take the following values: 0 Small mass limit for fermions.
1 Full mass effects for fermions.

2
Large mass limit for fermions.

3
Full mass effects for scalars. 4 Large mass limit for scalars.
See Section 3.2.6 for more information on the loop approximations. Note that, for implemented processes, using an input file is the only way to change the approximation in which loops are computed.

SUSY.
H1jet includes the simplified SUSY model with two stopst 1 andt 2 considered in refs. [28] and [29]. To include the SUSY stopst 1 andt 2 in the quark loops, it will be necessary to set the first stop mass mt 1 to a non-zero value by using the --mst option. The second stop mass is then given by where ∆m is set with the --delta option. The stop Yukawa coupling factors will be given by: where Note that m t , m Z , and sin 2 θ W can be set with the --mt, --mZ, and --sinwsq options respectively, while sin 2 θt and tan β are SUSY specific options and can be set with the --sth2 and --tbeta options. All of the SUSY specific options are: --mst <value> SUSY stop mass, mt 1 [GeV]. Default: 0 GeV.
Note that the top partner mass m T and SUSY stop mass mt 1 can not both be set non-zero at the same time via command-line options. However, if one uses an input file, one can explicitly specify masses, couplings and loop approximations for an arbitrary number of fermions and scalars. This would also allow a user to implement a specific SUSY model with more supersymmetric partners, each with the appropriate coupling.

Relevant Options for Process: bbH
If process bbH is selected, i.e. bb → H + jet, then the following options are relevant:

Relevant Options for Process: Z
If process Z is selected, i.e. pp/pp → Z + jet, then the following options are relevant:

Relevant Options for Process: user
If process user is selected, i.e. a custom user-specified process, any of the above physics options may be relevant if they are used in the custom amplitude code. The code will have to be inspected to determine this. The only built-in process-relevant option is: Additional options may be added depending on the custom process/amplitude. See Section 5 below for more details on the implementation of a custom process.

Loop Approximations
Small and large mass limits can be used as approximations for the quarks in the loop calculations. This requires some knowledge of the meaning of the approximations, hence this needs to be set by the user at compile time.
In the file input.f90 located in the src-directory, the subroutine reset_iloop_array can be found. This subroutine can be used by the user to set the iloop_array, which is an array specifying the approximation used for each loop particle. The approximations that can be used are: iloop_sm_fermion Small mass limit for fermions.
iloop_fm_fermion Full mass effects for fermions.
iloop_lm_fermion Large mass limit for fermions.
iloop_fm_scalar Full mass effects for scalars.
iloop_lm_scalar Large mass limit for scalars.

Output
The helper script PlotH1jet.py facilitates easy and quick plotting of the output from H1jet. The script requires Python 3 installed in order to run. The user needs to simply pipe the output of H1jet to the script: ./bin/h1jet [options] | python PlotH1jet.py Alternatively, the plotting script can run on an output file from H1jet: ./bin/h1jet [options] -o result.out python PlotH1jet.py result.out A resulting example plot with default settings in H1jet is shown in Figure 1.
A comparison between the various built-in models is shown in Figure 2. Default SM parameters has been used with mt 1 = 600 GeV, ∆m = 200 GeV, tan β = 5, m T = 1.7 TeV, sin 2 θt = sin 2 θ T = 0.1, and f = 900 GeV, and considering the M4 5 model as the explicit top-partner model.

Benchmarking
The various processes implemented in H1jet have been compared to those of SusHi [17,18], and have all been found to be in agreement. The relative ratio between the H1jet result and the SusHi result for the p T distribution for the CP-odd Higgs are shown in Figure 3, and is found to be in agreement within the Monte Carlo error of SusHi for a large range of p T values. Overall the agreement with SusHi is within 3 × 10 −4 . Note that the largest discrepancies were observed in the low p T region. To validate the H1jet results we have compared them to the approximate expression valid at low p T dσ dp T where σ 0 is the total Born cross section for gg → H. In Figure 4, we show p T σ 0 dσ dp T with the first term of eq. (24) subtracted, as a function of ln p T m H . For p T → 0 this goes nicely towards a constant as expected.   The relative ratio between the H1jet and SusHi results for the SUSY are shown in Figure 5 and is within 2 × 10 −4 . Again the low p T behaviour can be checked by comparing to the resummed expression in Figure 6.
Note that our numerical accuracy crucially depends not only on the accuracy of the numerical integration, but of that of the auxiliary programs used to compute the PDF evolution (HOPPET) and the scalar integrals (CHAPLIN). We have modified various internal parameters of the two libraries, and we obtained differences that are less than permille level. So, a conservative estimate of the numerical uncertainty of H1jet is 1 × 10 −3 .

Adding New Processes to H1jet
H1jet can be interfaced to use the squared matrix element evaluated from a custom Fortran code. The implementation may be most easily explained with a specific example. This section should be read very carefully before attempting to use the interface.

Example: Axion-Like-Particle (ALP) Effective Theory
We will present here a specific example of adding to H1jet the production of a light axion-like-particle (ALP), a, along with a jet. For simplicity, we only consider the gluon-fusion channel, This is a tree-level process due to an effective ALP-gluon coupling, where G a µν is the gluon field strength tensor andG a µν = (1/2) µνρσ G aµν its dual. The model and the FeynRules [30] model files are described and provided in ref. [31]. We use FeynCalc [32,33,34] to evaluate the amplitude from the model, so we have to convert the FeynRules model to a The resulting FeynArts model files are written to a new directory ALP_linear_FA, which needs to be moved to the FeynArts/Models directory. Note that in the FeynArts model, the ALP field is called S [4] and the gluon fields are called V [4].
In a new Mathematica session, we load FeynCalc with FeynArts: This ensures that the model files works with FeynCalc.
Then we create the tree-level 2 → 2 topologies and insert the relevant fields for our process: While not strictly necessary, it is recommended to enable the SMP option. Any additional substitutions in the amplitude can be specified with the FinalSubstitutions option. We then set up the kinematics:

In[13]:= FCClearScalarProducts[];
In [14]:= SetMandelstam[s, t, u, k1, k2, -k3, -k4, 0, 0, 0, mA]; We introduce here a parameter mA for the ALP mass m a . We then square the amplitude: Setting the SUNNToCACF option in SUNSimplify[] to False is not necessary, nor is it necessary to fix SUNN to 3. This can be handled by the dressing script and H1jet.
Finally, we write the amplitude as Fortran code to a file: Note here that we specify the gluon-gluon channel with the gg = ampsquared input to the function. This is required for the subsequent dressing script to work properly. It is important to specify the 2-particle initial state by using combinations of g, u, d, c, s, b, ubar, dbar, cbar, sbar, and bbar. One can also use q and qbar for all the light quarks and antiquarks respectively, i.e. u, d, c, and s. For example, bbbar will be the bb channel.
The helper script provides a help message which can be called with -h or --help. The name of the output file can be specified with the -o option. Multiple input Fortran files can be given as arguments to the helper script. The full usage is: The provided input Fortran code files does not necessarily have to be generated with FeynCalc. They can be generated by any other program or even be written by hand by providing the appropriate expression for gg or the desired channel.
To use the new dressed custom Fortran code with H1jet, it is necessary to recompile H1jet with the custom Fortran code: ./configure USERFILE=/path/to/custom/user_interface.f90 make clean make Running ./h1jet --help we see that three new additional options have been added:
The leading c_ in the name stands for "custom" and is automatically added in order to avoid naming issues in the code.
The result from the ALP implementation in H1jet is shown in Figure 7 and can be compared to the H1jet result for the CP-odd Higgs by using a single top quark in the loop with an infinite mass limit, resulting in an effective coupling between the CP-odd Higgs and the gluons. In fact, the respective ALP and CP-odd couplings are then related as such, The comparison is shown in Figure 8, where we see agreement within 4 × 10 −6 .
The result can also be compared to the same FeynRules model used with MadGraph5_aMC@NLO [16], where our code takes one second to run, while MadGraph can take up to several hours depending on the number of events, due to MadGraph running a full Monte Carlo integration. We have found that H1jet agrees with MadGraph5_aMC@NLO within Monte Carlo errors. We have also seen that MadGraph5_aMC@NLO runs into numerical instabilities at low p T , while H1jet has by construction the correct behaviour.

The Total Cross Section
While not strictly necessary for the user interface to run, it is still recommended to add the code for the evaluation of the total cross section to the custom user interface. This is easy to do as well. We will here show it for the ALP model. We start with considering the gluon-fusion ALP production, gg → a. In Mathematica, create a tree-level 2 → 1 topology, and insert the fields: Then we set up the amplitude: We again specify the gluon-gluon channel with the xsgg = xsec, but this time indicate with the leading xs that the code is for the Born cross section. Otherwise, the same rules apply. It is crucial to make sure not to save the Born cross section in the same file as the squared amplitude code for the transverse momentum distribution.
The new generated code ALP_xsec.f90 is provided to DressUserAmpCode.py along with the squared amplitude code: python DressUserAmpCode.py ALP_amp.f90 ALP_xsec.f90 And H1jet can be recompiled to include the new source code: ./configure USERFILE=/path/to/custom/user_interface.f90 make clean make After which H1jet will calculate the Born cross section for the custom process.

Conclusions
We have presented a method that allows a fast computation of the transverse momentum distribution of a colour singlet. The method is implemented in the program H1jet, which returns a transverse momentum spectrum for the specified colour singlet in about a second. H1jet is similar in spirit to SusHi, but is incomparably faster.
The program implements various processes, including Higgs production both in gluon fusion and bottom-antibottom annihilation, as well as Z production. Loop-induced Higgs production is implemented not only in the SM, but also in attractive BSM scenarios, such as SUSY or composite Higgs. For SUSY, we implement a simplified model with two stops, as done in ref. [29]. For composite Higgs, we implement both the simplified model of ref. [10], as well as some explicit models with one or more top partners [27]. Loop integrals can be computed either exactly or in the infinite-mass limit. The latter limit implements in practice a dimension-6 contact interaction between the Higgs and the gluon field. The program is very flexible, and the only process-dependent input is the corresponding amplitude in terms of Mandelstam invariants. This can be computed by the user either manually, or with the use of automated programs such as FeynCalc [34], and connected to the program via a simple interface. As an example, we have included in the package the calculation of the transverse momentum distribution of an ALP starting from the general Feynman rules of ref. [31]. Note that the possibility of automatically implementing a new model inside the program is a feature that is not available in SusHi. We also stress that it is also possible to take advantage of input files to obtain results for an arbitrary number of fermions and scalars in loops, with appropriate couplings. This could be used, for instance, to implement the MSSM instead of the provided simplified SUSY model.
We stress that H1jet is not a replacement for a proper Monte Carlo analysis implementing realistic experimental cuts. However, we believe it will be invaluable for BSM experts to assess whether a given model gives sizeable deviations from the SM. In fact, due to its fast implementation, H1jet makes it possible to perform parameter scans in seconds, and to take into account mass effects in specific models. Also, due to the fact that H1jet is not based on a Monte Carlo integration, one can separate interference between different contributions very precisely, something which is very difficult to achieve with Monte Carlo event generators.
H1jet can be also useful to precision phenomenology. In fact, it makes it possible to easily perform theoretical studies of the transverse momentum distribution of a colour singlet, especially those involving the matching of resummed calculations with exact fixed order. Also, although the implemented cross sections are computed at the lowest order in QCD, nothing prevents the inclusion of higher orders, provided one integrates over all coloured particles.

Bubbles. The bubble integral is defined as
where z = 1 2 The argument of the logarithm in eq. (28) has a different form according to the value of s: Note that the only case in which one needs a small imaginary part is the case s > 4m 2 . This imaginary part has the opposite convention as in CHAPLIN. As a solution, we invert the argument of the logarithm and use the identity ln z = − ln(1/z). In practice, after an appropriate analytic continuation of the square root, we define and implement the bubble as follows: Note that a logarithm of a negative number can also be correctly analytically continued by using the default Fortran implementation of the complex logarithm. As for CHAPLIN, Fortran assumes that a negative number has a small positive imaginary part. Therefore, in case we have a small negative imaginary part, we can still use the relation ln z = − ln(1/z), which gives the correct analytic continuation.
Triangles. The triangle integral C 0 (s) is defined as C 0 (s) = 1 2s where z is given in eq. (29). Again, for s > 4m 2 , the argument of the logarithm has the opposite sign with respect to what is implicitly assumed by CHAPLIN. Therefore, we invert again the argument of the logarithm, and using the definition of z in eq. (31), we implement the triangle as follows: s C 0 (s) = H 1, 1; 1 z . (34) Boxes. The scalar four-point function with three massless (the gluons) and one massive (the Higgs boson) external lines is given by [25], which can be expressed in terms of complex dilogarithms by using the exact result 1 st where are real numbers, with x + > 1 and x − < 0, and acquires an imaginary part according to the value of v. In particular, keeping track of the imaginary part of y yields From the above, we see that w for 0 < v < 4m 2 we can use the dilogarithms as given by CHAPLIN. For v < 0, x − /(x − − y) and x + /(y − x − ) acquire a small positive imaginary part, whereas x + /(x + − y) and x − /(y − x + ) a small negative imaginary part. The reverse happens for v > 4m 2 . Therefore, we need to perform some formal manipulations to use the harmonic polylogarithms provided by CHAPLIN. In practice, whenever the argument z of the dilogarithm is complex, we just use the definitory relation Li 2 (z) = H(0, 1; z). When z = x + iε, with x real, we use Li 2 (x + iε) = H(0, 1; x), with H(0, 1; x) the complex number provided by CHAPLIN. If z = x − iε, we use the identities and we select the one that gives the smallest imaginary part. This of course give numerically indistinguishable results when the imaginary part is large, but is of crucial importance when the imaginary part should be zero but it is not because of the specific numerical methods employed by CHAPLIN.