SecDec: A general program for sector decomposition
Abstract
We present a program for the numerical evaluation of multidimensional polynomial parameter integrals. Singularities regulated by dimensional regularisation are extracted using iterated sector decomposition. The program evaluates the coefficients of a Laurent series in the regularisation parameter. It can be applied to multiloop integrals in Euclidean space as well as other parametric integrals, e.g. phase space integrals.
PACS: 12.38.Bx, 02.60.Jh, 02.70.Wz
keywords:
Perturbation theory, Feynman diagrams, infrared singularities, multidimensional parameter integrals, numerical integrationIPPP/10/94, DCPT/10/188
,
PROGRAM SUMMARY
Manuscript Title: SecDec: A general program for Sector Decomposition
Authors: J. Carter, G. Heinrich
Program Title: SecDec
Journal Reference:
Catalogue identifier: AEIR_v1_0
Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html
Programming language: Wolfram Mathematica, perl, Fortran
Computer: from a single PC to a cluster, depending on the problem
Operating system: Unix, Linux
RAM: depending on the complexity of the problem
Keywords: Perturbation theory, Feynman diagrams, infrared singularities, multidimensional parameter integrals, numerical integration
PACS:
12.38.Bx,
02.60.Jh,
02.70.Wz
Classification:
4.4 Feynman diagrams,
5 Computer Algebra,
11.1 General, High Energy Physics and Computing.
Nature of problem:
Extraction of ultraviolet and infrared singularities from parametric integrals
appearing in higher order perturbative calculations in gauge theories, e.g.
multiloop Feynman integrals, Wilson loops, phase space integrals.
Solution method:
Algebraic extraction of singularities in dimensional regularisation using iterated sector decomposition.
This leads to a Laurent series in the dimensional regularisation parameter ,
where the coefficients are finite integrals over the unithypercube.
Those integrals are evaluated numerically by Monte Carlo integration.
Restrictions: Depending on the complexity of the problem, limited by
memory and CPU time. Multiscale integrals can only be evaluated at Euclidean points.
Running time:
Between a few minutes and several days, depending on the complexity of the problem.
LONG WRITEUP
1 Introduction
Sector decomposition is an algorithmic method to isolate divergences from parameter integrals as they occur for instance in perturbative quantum field theory. Originally it was devised by Hepp [1] in the context of the the proof of the BPHZ theorem in order to disentangle overlapping ultraviolet singularities. Similar ideas, applied to the subtraction of infrared divergences, can be found e.g. in [2]. It was employed later to extract logarithmic mass singularities from massive multiscale integrals in the high energy limit at two loops [3, 4].
In [5], the concept of sector decomposition was elaborated to a general algorithm in the context of dimensional regularisation, allowing the isolation of ultraviolet as well as infrared singularities from Feynman parameter integrals in an automated way. First applications of this algorithm were the numerical evaluation of twoloop box diagrams at certain Euclidean points, see e.g. [5, 6, 7]. More recently, the method has been used to numerically check a number of analytic threeloop and fourloop results [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21], most of them produced by either the public program FIESTA [22, 23] or the code which is described in the present article. Further references about recent applications of sector decomposition to multiloop calculations can be found in [23, 24].
Sector decomposition also has been combined with other methods for a numerical calculation of loop amplitudes, first on a diagrammatic level in Refs. [25, 26], later for whole amplitudes in Refs. [27, 28, 29, 30]. The latter approaches contain a combination of sector decomposition and contour deformation [31, 32, 33, 34, 35], which allows one to integrate the Feynman parameter representation of an amplitude numerically in the physical region.
As phase space integrals in dimensions can be written as dimensionally regularised parameter integrals, sector decomposition can also serve to factorise entangled singularity structures in the case of soft and collinear real radiation. This idea was first presented in [36] and was subsequently applied to calculate all master fourparticle phase space integrals where up to two particles in the final state can become soft and/or collinear [37]. Shortly after, this approach has been extended to be applicable to exclusive final states as well by expressing the functions produced by sector decomposition in terms of distributions [38]. Further elaboration on this approach [39, 40] has lead to differential NNLO results for a number of processes [41, 42, 43, 44, 45, 46, 47, 48, 49]. The combination of the FrixioneKunsztSigner subtraction scheme [50] for soft and collinear real radiation with the decomposition into sectors to treat real radiation at NNLO, as proposed recently in [51], is also promising with regards to the reduction of the number of functions produced by the decomposition. A combination of sector decomposition with nonlinear variable transformations as proposed in [52] can also serve to reduce considerably the number of functions to integrate, but is less straightforward to automate completely.
To date, the method of sector decomposition has been applied successfully to a considerable number of higher order calculations, for a review we refer to [53, 24]. Here we will concentrate on the method of sector decomposition from a programming point of view.
Despite its success in practical applications, for quite some time there was no formal proof for the existence of a strategy for the iterated sector decomposition such that the iteration is always guaranteed to terminate. This gap has been filled in Ref. [54], by mapping the problem to Hironaka’s Polyhedra game [55] and offering three strategies which are proven to terminate. Bogner and Weinzierl also implemented the algorithm in a public computer program for iterated sector decomposition written in C++ [54]. A Mathematica interface to this program, which also allows the calculation of contracted tensor integrals, has recently been published in [56].
A different strategy guaranteed to terminate, leading to less subsectors than the strategies of Ref. [54], was given by A.V. Smirnov and M. Tentyukov, who implemented the algorithm in the public program FIESTA [22]. Based on a detailed analysis of Hepp and Speer sectors in Ref. [57], an alternative strategy, which is based on Speer sectors, has been implemented in FIESTA2 [23]. As the latter strategy also uses information on the topology of the graph, it can perform the decomposition more efficiently in certain cases.
Another group has implemented [58] the sector decomposition algorithm in FORM [59]. Mapping sector decomposition to convex geometry and using algorithms in computational geometry lead to a guaranteed terminating strategy which seems to be optimal with regards to the number of produced subsectors [60, 61].
Compared to the existing packages, the present program for sector decomposition has several new features, the main ones being a sophisticated treatment of (potential) numerical instabilities and the possibility to apply the program to more general functions, like e.g. parameter integrals occurring in phase space integrals over real radiation matrix elements, or functions including symbolic parameters. Halfinteger powers of polynomial functions are also allowed.
The structure of the article is as follows. In section 2, we briefly review the theoretical framework. Then we give an overview of the program structure and explain the individual components. The installation and usage of the program are described in section 4, followed by a number of examples illustrating different aspects of the program in section 5. The appendix contains more details about timings and phase space parametrisations as well as a section with further information for the user.
2 Theoretical background
2.1 Feynman integrals
A general Feynman graph in dimensions at loops with propagators and loop momenta in the numerator, where the propagators can have arbitrary, not necessarily integer powers , has the following representation in momentum space:
(1) 
where the are linear combinations of external momenta and loop momenta . Introducing Feynman parameters leads to
(2)  
where , is a matrix containing Feynman parameters, is an dimensional vector composed of external momenta and Feynman parameters, and contains kinematic invariants and Feynman parameters.
To perform the integration over the loop momenta , we perform the following shift in order to obtain a quadratic form for the term in square brackets in eq. (2):
(3) 
After momentum integration one obtains
where  
(5)  
and denotes the nearest integer less or equal to .
The expression
stands for the sum over all different combinations
of doubleindices distributed to metric tensors and
vectors . The double indices
denote the Lorentz index, belonging to the loop momentum.
As can be seen from Eq. (2.1), the difference between scalar () and tensor () integrals, once the Lorentz structure is extracted, is given by the fact that there are additional polynomials of Feynman parameters in the numerator. These polynomials can simply be included into the sector decomposition procedure, thus treating contracted tensor integrals directly without reduction to scalar integrals.
The functions and also can be constructed from the topology of the corresponding Feynman graph. For more details we refer to [62, 63, 53, 24].
is a positive semidefinite function. Its vanishing is related to the UV subdivergences of the graph. Overall UV divergences, if present, will always be contained in the prefactor . In the region where all invariants formed from external momenta are negative, which we will call the Euclidean region in the following, is also a positive semidefinite function of the Feynman parameters . Its vanishing does not necessarily lead to an IR singularity. Only if some of the invariants are zero, for example if some of the external momenta are lightlike, the vanishing of may induce an IR divergence. Thus it depends on the kinematics and not only on the topology (like in the UV case) whether a zero of leads to a divergence or not. The necessary (but not sufficient) conditions for an IR divergence are given by the Landau equations [64, 65, 66], which, in parameter space, simply mean that the necessary condition for an IR divergence can only be fulfilled if some of the parameters go to zero, provided that all kinematic invariants formed by external momenta are negative.
For a diagram with massless propagators, none of the Feynman parameters occurs quadratically in the function . If massive internal lines are present, gets an additional term . If the power of the Feynman parameters in the polynomial forming is larger than one for at least two different parameters, initially or at a later stage in the iterated decomposition, an infinite recursion can occur. This happens in the example given in section 5.1.2 if the “naive” decomposition strategy is employed. We have implemented a heuristic procedure to change to a different decomposition strategy only in cases where at least two Feynman parameters occur quadratically, which lead to a terminating algorithm without producing a large number of subsectors. We do not claim that this procedure is guaranteed to terminate, but it proved useful for practical purposes.
2.2 General parameter integrals
The program can also deal with parameter integrals which are more general than the ones related to multiloop integrals. The only restrictions are that the integration domain should be the unit hypercube, and the singularities should be only endpoint singularities, i.e. should be located at zero or one. The general form of the integrals is
(6) 
where are polynomial functions of the parameters , which can also contain some symbolic constants . The user can leave the parameters symbolic during the decomposition, specifying numerical values only for the numerical integration step. This way the decomposition and subtraction steps do not have to be redone if the values for the constants are changed. The are powers of the form (with such that the integral is convergent; note that half integer powers are also possible). We would like to point out that most phase space integrals in dimensions over real radiation matrix elements can also be remapped to functions of the type (6). Examples are given in section 5.
2.3 Iterated sector decomposition
Here we will review the basic algorithm of iterated sector decomposition only briefly. For more details we refer to [5, 24].
Our starting point is a function of the form of Eq. (6). Loop integrals in the form of Eq. (2.1), with open Lorentz indices contracted by external momenta or metric tensors, are a special case thereof, distinguished by the presence of a distribution which constrains the sum of the integration parameters. Therefore loop integrals are treated somewhat differently than the more general functions, meaning that the constraint is integrated out in a special way leading to the socalled primary sectors.
I. Generation of primary sectors (loop integrals only)
We split the integration domain into parts and eliminate the distribution in such a way that the remaining integrations are from 0 to 1. To this aim we decompose the integration range into sectors, where in each sector , is largest:
(7) 
The integral is now split into domains corresponding to integrals from which we extract a common factor: . The generated sectors will be called primary sectors in the following. In the integrals we substitute
(8) 
and then integrate out using the constraint. As are homogeneous of degree , , respectively, and factorises completely, we have and and thus, using , we obtain
(9) 
Note that the singular behaviour leading to –poles still comes from regions where a set of parameters goes to zero.
II. Extraction of the singular factors
For functions of the form Eq. (6), we first have to determine which of the Feynman parameters generate singularities at , and which ones can lead to singularities at zero and one. The parameters for which a denominator vanishes at but not at should be remapped by the transformation . If the integrand can become singular at both endpoints of the integration range for a parameter , we split the integration range at 1/2: After the split
(10) 
and the substitution in and in , all endpoint singularities occur at only. This splitting is done automatically by the program; the user only has to define which variables should be split.
At this stage our starting point is a parametric integral where the integrand vanishes if some of the integration parameters go to zero. Our aim is to factorise the singularities, i.e. extract them in terms of overall factors of type . We proceed as follows.
 1.

Determine a minimal set of parameters, say , such that at least one of the functions vanishes if the parameters of are set to zero.
 2.

The corresponding integration range is an cube which is decomposed into subsectors by decomposing unity according to
(11)  3.

Remap the variables to the unit hypercube in each new subsector by the substitution
(12) This gives a Jacobian factor of . By construction factorises from at least one of the functions .
For each subsector the above steps have to be repeated as long as a set can be found such that one of the rescaled functions vanishes if the elements of are set to zero. This way new subsectors are created in each subsector of the previous iteration, resulting in a treelike structure after a certain number of iterations. The iteration stops if the functions contain a constant term, i.e. if they are of the form
(13) 
where are polynomials in the variables , and is a constant, i.e. is nonzero.
The resulting subsector integrals have the general form
(14) 
The singular behaviour of the integrand now can be read off directly from the exponents , for a given subsector integral.
III. Subtraction of the poles
For a particular , the integrand after the factorisation described above, is of the form
(15) 
If , no subtraction is needed and one can go to the next variable . If , one expands into a Taylor series around . Subtracting the Taylor series (to order^{1}^{1}1To account for halfinteger exponents, e.g. , we use , denoting the nearest integer less or equal to . for ) and adding it back in integrated form, we obtain a part where the poles are subtracted and a part exhibiting poles times a function depending only on the remaining integration parameters.
(16) 
For , expanding the above expression in is equivalent to an expansion in “plus distributions” [67, 38]
where  
(17) 
with the integrations over the terms containing already carried out.
After having done the subtractions for each , all poles are extracted, such that the resulting expression can be expanded in . This defines a Laurent series in
(18) 
where the coefficients are finite parameter integrals of dimension for and of dimension for . denotes the leading pole, which can be at most for an loop integral. The finite coefficient functions can be integrated by Monte Carlo integration if the Mandelstam invariants in respectively the numerical constants in a general integrand have been chosen such that the integrand does not vanish in the integration domain.
Improving the numerical stability
For in eq. (2.3), the singularity is of logarithmic nature, i.e. if a lower cutoff for the parameter integral was used. In renormalisable gauge theories, linear () or even higher () poles should not occur. However, they can occur at intermediate stages of a calculation, and as they are formally regulated by dimensional regularisation, a method has been worked out for the program to be able to deal with higher than logarithmic singularities efficiently. This method relies on integration by parts (IBP) in a way which aims at maximal numerical stability.
Let us consider an integrand after subtraction, as the one in eq. (2.3), and focus on only one variable, , and define
(19) 
where, in the case , is as by construction. Integration by parts gives
(20)  
(21) 
As vanishes as , the term in the square bracket in
eq. (20) is zero.
Also notice that as , so can be treated in the same way.
For the case we thus obtain, expanding in
(22) 
For , we can use eq. (21) to iterate this procedure until we reach .
This method is certainly beneficial in the case of numerical instabilities coming from terms of the form or in , leading to differences of large numbers for . The IBP method raises the powers of , trading them for additional logarithms in the integrand, which can be integrated more easily by the Monte Carlo program. Note also that the whole procedure is linear in , so it allows one to split into smaller functions which can be dealt with more easily by the numerical integrator.
Error treatment
We usually integrate the sum of a small set of functions stemming from the expansion of a certain pole structure individually, and afterwards sum all the individual results contributing to a certain pole coefficient. The set of functions to sum before integration is defined by the size of the individual functions: the functions are summed until their sum reaches about two Megabytes.
The errors are calculated by adding the Monte Carlo errors of the individual integrations in quadrature. It is possible that there are large cancellations between different functions contributing to the same pole coefficient. In such cases it is better to first sum all coefficient functions and then integrate. The program offers the possibility to do so as an option which can be specified in the input parameter file.
The user should be aware that for complicated functions containing many subtractions, the Monte Carlo error estimate is not quite appropriate: it is calculated on a purely statistical basis, scaling like if is the number of sampling points. However, this is only a reliable error estimate under the assumption that the sampling has mapped all the important features of the function (i.e. all peaks) sufficiently precisely, and strictly is only valid for square integrable functions. If the function is not square integrable (but integrable), the Monte Carlo estimate for the integral will still converge to the true value, but the error estimate will become unreliable. For more involved integrals, we are faced with functions which have gone through numerous decompositions and subtractions, such that their shape in the unit hypercube is quite complicated, and therefore the naive Monte Carlo error estimate tends to underestimate the “true” error.
Often the main source of underestimated errors in the final result is the fact that there are a large number of integrations to sum, and so adding the errors in quadrature would only give a truly appropriate error estimate if there were no systematic errors in the numerical integration.
We should also note that converting the result to extract a different dependent prefactor may lead to cancellations between different contributions to a certain pole coefficient such that the error estimate may be too optimistic in these cases.
3 Structure of the program
The program consists of two parts, an algebraic part and a numerical part. The algebraic part uses code written in Mathematica [68] and does the decomposition into sectors, the subtraction of the singularities, the expansion in and the generation of the files necessary for the numerical integration. In the numerical part, Fortran functions forming the coefficient of each term in the Laurent series in are integrated using the Monte Carlo integration program BASES, version 5.1 [69], or one of the routines from the CUBA library, version 2.1 [70]. The different subtasks are handled by perl scripts. The directory structure of the program is shown in Fig. 1, while the flowchart in Fig. 2 shows the basic flow of input/output streams.
The directories loop and general have the same global structure, only some of the individual files are specific to loops or to more general parametric functions. The directories contain a number of perl scripts steering the decomposition and the numerical integration. The scripts use perl modules contained in the subdirectory perlsrc.
The Mathematica source files are located in the subdirectories src/deco: files used for the decomposition, src/subexp: files used for the pole subtraction and expansion in , src/util: miscellaneous useful functions. For the translation of the Mathematica expressions to Fortran77 functions we use the package Format.m [71]. The subdirectories basesv5.1 and Cuba2.1 contain the libraries for the numerical integration, taken from [69] and [70], respectively. The documentation, created by robodoc [72] is contained in the subdirectory doc. It contains an index to look up documentation of the source code in html format by loading masterindex.html into a browser.
The intermediate files and the results will be stored in a subdirectory of the working directory whose name mysubdir can be specified by the user (first entry in param.input, leaving this blank is a valid option). A subdirectory of mysubdir with the name of the graph, respectively integral to calculate will be created by default. If the user would like to store the files in a directory which is not the subdirectory of the working directory, for example in /scratch, he can do this by specifying the full path in the second entry in param.input. An example of a directory structure created by running the examples NPbox, QED, ggtt1, A61, a userdefined 3loop example, and a 4loop example to be written to the scratch disk is given in Fig. 3.
The directory created for each graph will contain subdiretories according to the pole structure
of the graph. The labelling for the pole structure is of the form e.g. 2l0h0,
denoting 2 logarithmic poles, no linear and no higher poles.
We should point out that this labelling does not necessarily correspond to the
final pole structure of the integral. It is merely for bookkeeping purposes,
and is based on the counting of the powers of the factorised integration variables.
In more detail, if variables have power , variables have a power
and variables have a power , the labelling will be
, even though the nonlogarithmic poles will disappear
upon expansion. In particular, for halfinteger powers, the labelling does not
correspond to “true” poles, but rather to terms which can be cast into functions like
, which are welldefined in the context of dimensional regularisation,
where can be regarded as an arbitrary
(complex) parameter.
Note also that in the case of a prefactor containing poles multiplying the parameter
integral, the poles which are flagged up at this stage of the program will only correspond to
the poles read off from the integration parameters.
In any case, the final result will be given to the order specified by the user
in param.input.
Each of these “polestructure” directories contains
further subdirectories where the files for a particular power in epsilon are stored.
An example is given in Fig. 4.
The user only has to edit the following two files:

param.input: (text file)
specification of paths, type of integrand, order in , output format, parameters for numerical integration, further options 
Template.m: (Mathematica syntax)

for loop integrals: specification of loop momenta, propagators; optionally numerator, nonstandard propagator powers

for general functions: specification of integration variables, integrand, variables to be split

To give a specific example rather than empty templates, the files param.input and Template.m in the loop subdirectory contain the setup for example 1, the nonplanar massless onshell twoloop box diagram, while those in the general directory contain the setup for example 6, a hypergeometric function of type .
Apart from these default parameter/template files, the program comes with example input and template files in the subdiretories loop/demos respectively general/demos, described in detail in section 5.
The user can choose the numerical integration routine and the settings for the different integrators contained in the Cuba library in the file param.input. The compilation of the chosen integration routine with the corresponding settings will be done automatically by the program.
4 Installation and usage
Installation
The program can be downloaded from http://projects.hepforge.org/secdec/.
Installation is done by unpacking the tar archive, using the command tar xzvf SecDec.tar.gz. This will create a directory called SecDec with the subdirectories as described above. Change to the SecDec directory and run ./install.
Prerequisites are Mathematica, version 6 or above, perl (installed by default on most Unix/Linux systems) and a Fortran compiler (e.g. gfortran, ifort). The install script only checks if Mathematica and perl are installed on the system and inserts the corresponding path into the perl scripts. The install script does not test the existence of a Fortran compiler because the compiler should be specified by the user in param.input. If no compiler is specified, it defaults to gfortran.
Usage

Change to the subdirectory loop or general, depending on whether you would like to calculate a loop integral or a more general parameter integral.

Copy the files param.input and Template.m to create your own parameter and template files myparamfile, mytemplatefile.

Set the desired parameters in myparamfile and define the propagators etc. in mytemplatefile.

Execute the command ./launch p myparamfile t mytemplatefile in the shell.
If you omit the option p myparamfile, the file param.input will be taken as default. Likewise, if you omit the option t mytemplatefile, the file Template.m will be taken as default. If your files myparamfile, mytemplatefile are in a different directory, say, myworkingdir, use the option
d myworkingdir, i.e. the full command then looks like ./launch d myworkingdir p myparamfile t mytemplatefile, executed from the directory SecDec/loop or SecDec/general.
Alternatively, you can call the launch script from any directory if you prepend the path to the launch script, i.e. the command
path_to_launch/launch p myparamfile t mytemplatefile executed from myworkingdir would run the program in the same way. path_to_launch can be either the full or relative path for SecDec/loop or SecDec/general.The ./launch command will launch the following perl scripts:

makeFU.pl: (only for loop integrals) constructs the integrand functions and the numerator function from the propagators and indices given in Template.m.

decompose.pl: launches the iterated sector decomposition

subexp.pl: launches the subtractions and epsilonexpansions and writes the Fortran functions. Depending on the “exeflag” specified in the parameter file (see below for a detailed explanation of the flag), this script also launches the compilation and the numerical integrations.


Collect the results. Depending on whether you have used a single machine or submitted the jobs to a cluster, the following actions will be performed:

If the calculations are done sequentially on a single machine, the results will be collected automatically (via results.pl called by launch). The output file will be displayed with your specified text editor. The results are also saved to the files [graph]_[point]epstothe*.res and [graph]_[point]full.res in the subdirectory subdir/graph (loops) respectively subdir/integrand (general integrands) (name specified in param.input, where you can also specify different names for different numerical points).

If the jobs have been submitted to a cluster: when all jobs have finished, execute the command ./results.pl [d myworkingdir p myparamfile] in a shell from the directory SecDec/loop or SecDec/general to create the file containing the final results.
If the the user needs to change the batch system settings: manually edit perlsrc/makejob.pm and perlsrc/launchjob.pm. This writes the desired syntax to the scripts job[polestructure] in the corresponding subdir/graph or subdir/integrand subdirectory.


After the calculation and the collection of the results is completed, you can use the shell command ./launchclean[graph] to remove obsolete files.
If called with no arguments, the script only removes object files, launch scripts, makefiles and executables, but leaves the Fortran files created by Mathematica, so that different numerical points can be calculated without rerunning the Mathematica code. If called with the argument ’all’ (i.e. ./launchclean[graph] all), it removes everything except the result files displaying the final result and the timings.
The ’exe’ flag contained in param.input offers the possibility to run the program only up to certain intermediate stages. The flag can take values from 0 to 4. The different levels are:
 exe=0:

does the iterated sector decomposition and writes files containing lists of subsector functions (graphsec*.out) for each pole structure to the output subdirectory. Also writes the Mathematica files subandexpand*.m for each pole structure, which serve to do the symbolic subtraction, epsilon expansion and creation of the Fortran files. Also writes the scripts
batch[polestructure] which serve to launch these jobs at a later stage.  exe=1:

launches the scripts batch[polestructure]. This will produce the Fortran functions and write them to individual subdirectories for each pole structure.
 exe=2:

creates all the additional files needed for the numerical integration.
 exe=3:

compilation is launched to make the executables.
 exe=4:

the executables are run.
If the first steps of the calculation, e.g. the decomposition or the creation of the Fortran functions, are already done, the following commands are available to continue the calculation without having to restart from scratch:

finishnumerics.pl [d myworkingdir p myparameterfile]:
if the ’exe’ flag in param.input resp. myworkingdir/myparameterfile is set smaller than four, this will complete the calculation without redoing previous steps. 
justnumerics.pl [d myworkingdir p myparameterfile]:
if you would like to redo just the numerical integration, for example to produce results for a different numerical point or to try out a different number of sampling points, iterations etc. for the Monte Carlo integration: change the values for the numerical point resp. the settings for the Monte Carlo integration and the “name of the numerical point” in the parameter file, and then use the command ./justnumerics.pl [d myworkingdir p myparameterfile] to redo only the numerical integrations (if the Fortran files f*.f have been produced already). Using this option skips the Mathematica subtraction and epsilon expansion step which can be done once and for all, as the variables at this stage are still symbolic. After completion of the numerical integrations, use the command ./results.pl [d myworkingdir p myparameterfile] to collect and display the results as above.
The program tries to detect the path to Mathematica automatically. In case you get the message “path for Mathematica not automatically found”, please insert the path to Mathematica on your system manually for the variable $mathpath in the file perlsrc/mathlaunch.pl.
We also should mention that the code starts working first on the most complicated pole structure, which takes longest. This is because in case the jobs are sent to a cluster, it is advantageous to first send the jobs which are expected to take the most time.
5 Description of Examples
5.1 Loop integrals
The examples described below can be found in the subdirectory loop/demos.
5.1.1 Example 1: Nonplanar massless twoloop box
The nonplanar massless twoloop box is a nontrivial example, as the sector decomposition applied to the standard representation, produced by combining all propagators simultaneously with Feynman parameters, exhibits “nonlogarithmic poles” (i.e. exponents of Feynman parameters ) in the course of the decomposition. We should point out that, even though the program can deal with linear or higher poles in a completely automated way, it is often a good idea to investigate if the integrand can be reparametrised such that poles of this type do not occur, because these poles require complicated subtraction terms which slow down the calculation. Nonlinear transformations as e.g. described in [52] can be useful in this context. Further, integrating out first one loop momentum, and then combine the remaining propagators with the obtained intermediate result using another set of Feynman parameters often leads to a representation where at least one of the parameters can be factorised without sector decomposition, thus speeding up the calculation considerably. This is demonstrated for the nonplanar massless twoloop box in appendix 8.2, and the template to calculate the graph in this way can be found in SecDec/general/demos.
To obtain results for the nonplanar massless twoloop box shown in Fig. 5 without doing any analytical steps, copy the file loop/param.input to a new parameter file, say paramNPbox.input, and specify the desired order in , the numerical point and possible further options. Likewise, copy the file loop/Template.m to a new template file, say templateNPbox.m (this already has been done for the examples described here, see the subdirectory loop/demos). Then, in the loop/demos directory, use the command ../launch p paramNPbox.input t templateNPbox.m. The file templateNPbox.m already has the propagators of the nonplanar double box predefined. In paramNPbox.input, we defined the prefactor such that a factor of is not included in the numerical result: if we define
(23) 
the program should yield the results for , given in Table 1. Note that according to eq. (2.1), we always divide loop integrals by , so this factor is never included in the numerical result. The decomposition produces 384 subsectors.
(s,t,u)  (1,1,1)  (1,2,3) 

1.75006  0.41670  
2.99969 0.00055  0.9313 0.00067  
22.821 0.003  5.8599 0.0035  
113.629 0.013  42.79 0.02  
395.27 0.05  162.730.09 
The result for the graph called NPbox at the numerical point called point in the input file will be written to the file NPbox_[point]full.res in the subdirectory 2loop/NPbox, where 2loop is a subdirectory which has been created by the program, using the directory name the user has specified in the first entry of paramNPbox.input. By default, a subdirectory with the name of the graph is created, but the user can also specify a completely different directory (e.g. scratch) where the results will be written to (second entry in paramNPbox.input).
More information about the decomposition is given in the file NPboxOUT.info.
Information about the numerical integration is contained in the files
[point]intfile.log in
the subdirectories graph/polestructure/epstothe[i],
where “polestructure” is of the form e.g. 2l0h0,
denoting 2 logarithmic poles and 0 linear, 0 higher poles.
It should be emphasized that in param.input, the numbers for the Mandelstam invariants should be defined as the Euclidean values, so the values for should always be negative in param.input. Note also that the condition cannot be fulfilled numerically in the Euclidean region, so it should not be used in onshell={…} in the template file to eliminate from the function in the case of nonplanar box graphs.
5.1.2 Example 2: Planar twoloop ladder diagram with massive onshell legs
The purpose of this example is to show how to deal with diagrams where the decomposition could run into an infinite recursion if the default strategy is applied. The rungs of the ladder are massless particles (e.g. photons), while the remaining lines are massive onshell particles, depicted by solid (blue) lines in Fig. 6.
To run this example, execute the command ../launch p paramQED.input t templateQED.m from the loop/demos directory. Only the primary sectors number one and seven are at risk of running into infinite recursion, therefore they are listed in the thirdlast item of paramQED.input as the ones to be decomposed by a different strategy. The results for the numerical point called point will be written to the file QED_[point]full.res in the subdirectory 2loop/QED. Numerical results for some sample points are given in Table 2. The kinematic points are defined by the mass and the Mandelstam variables . We extracted a prefactor of .
(s,t, m)  (0.2,0.3,1)  (3/2,4/3,1/5) 

1.56161 1.33  2.1817 0.0003  
5.3373 0.0018  1.4701 0.0026  
1.419 0.025  30.191 0.014  
62.46 0.18  140.730.057  
284.76 0.87  450.670.19 
5.1.3 Example 3: Nonplanar twoloop diagrams with two massive onshell legs
This example gives results for two nonplanar graphs occurring in the calculation of at NNLO, shown in Fig. 7. The analytic results for these graphs are not yet available. Numerical results at Euclidean points can be produced by choosing numerical values for the invariants in paramggtt1.input respectively paramggtt2.input and then executing the command ../launch p paramggtt1.input t templateggtt1.m in the loop/demos directory, analogously for ggtt2. Results for two sample points are shown in Table 3.
ggtt1  

(0.5,0.4,0.1,0.17,0.17)  (1.5,0.3,0.2,3,1)  
38.07970.0027  0.19904  
263.22 0.015  0.71466  
936.86 0.06  1.45505 0.0002  
ggtt2  
(0.5,0.4,0.1,0.17,0)  (1.5,0.3,0.2,3,0)  
10.9159 0.0006  0.13678  
43.5213 0.0075  0.2087 0.00024  
165.384 0.048  3.3417 0.0014  
20.8420.268  6.5930.007  
2117.5 1. 57  20.420.04 
5.1.4 Example 4: A rank one tensor twoloop box
In order to demonstrate how to run the program for integrals with nontrivial numerators, we give the example of a rank one planar massless onshell twoloop box, where we contract one loop momentum in the numerator by .
where we omitted the terms in the propagators. The result for the kinematic sample point is shown in Table 4.
(s,t,u)  (3,2,5) 

0.319449  
0.46536  
0.5848 0.0004  
3.3437 0.0013  
1.6991 0.0035 
Note that in this example, we used a positive value for the Mandelstam invariant , which seems to contradict the requirement to have only Euclidean values for the invariants. However, in this case we can do this because the function does not depend on at all. The numerator does depend on , but as a numerator which is not positive definite does not spoil the numerical convergence, we can as well choose a numerical value for such that the relation is fulfilled. This has the advantage that it allows us to use the latter relation to simplify the numerator.
5.1.5 Example 5: A threeloop vertex diagram with dependent propagator powers
This example shows how to calculate diagrams with propagator powers different from one. The results for the graph (notation of Ref. [8]), given in Table 5, can be produced by running ../launch p paramA61.input t templateA61.m from the loop/demos directory.
0.16666  1.8334  18.123  125.32  889.96  5325.3 
The analytical result for this diagram with general propagator powers is given in Ref. [8] and is also given in the file 3loop/A61/A61analytic.m to allow comparisons between analytical and numerical results for arbitrary propagator powers.
5.2 More general polynomial functions
The examples described below can be found in the subdirectory general/demos.
5.2.1 Example 6: Hypergeometric functions
As an example for “general” polynomial functions, we consider the hypergeometric functions , using the integral representation recursively:
(25)  
Considering with the values , we obtain the results shown in Table 6. The “analytic result” has been obtained using HypExp [73, 74].
order  analytic result  numerical result  time taken (secs) 

1  1.0000002  2  
0.189532  0.1895960.00036  21  
2.299043  2.3060.011  124  
55.46902  55.61 0.39  248  
1014.39  1018.45.9  429 
This result can be produced by typing ./launch d demos p param5F4.input t template5F4.m in the subdirectory general, or by typing ../launch p param5F4.input t template5F4.m in the subdirectory general/demos.
The program can also deal with functions containing half integer exponents. Table 7 shows results for with arguments . These results can be produced by the command ../launch p param4F3.input t template4F3.m in the subdirectory general/demos.
order  analytic result  numerical result  time taken (seconds) 

1  0.999997  1.6  
4.27969  4.2810 0.0055  54  
26.6976  26.625 0.121  90 
5.2.2 Example 7: Phase space integrals
Sector decomposition can be useful for the calculation of phase space integrals where infrared divergences are regulated dimensionally. This is particularly the case for double real radiation occurring in NNLO calculations involving massive particles, where analytic methods show their limitations.
Here we give examples of phase space integrals, which should be considered as part of a phase space written in factorised form. We choose particles 3 and 4 to be massless, while is the momentum of a massive state, either a single particle or a pseudostate formed by additional momenta in the final state, i.e. . After all integrations have been mapped to the unit interval, we have integrals of the form
(26)  
(27) 
The derivation is given in the appendix, section 8.3.
The invariants in this parametrisation are given by
where
(28)  
We would like to point out that for the examples below, more convenient parametrisations, i.e. parametrisations where the variables in the denominator factorise, and/or reflect symmetries of the squared matrix element, certainly do exist. However, the purpose of the examples is to illustrate that the code can deal with denominators which are amongst the most complicated ones which do occur in NNLO real radiation involving two (unresolved) massless particles in the final state, where they cannot always be “rotated away” by suitable transformations. A hybrid approach combining sector decomposition with convenient parametrisations/transformations is certainly the method of choice for real radiation at NNLO. The program can be used to evaluate the integrals occurring in such an approach.
Three massless particles in the final state
We first consider a case where is a massless particle, i.e. the limit in eq. (26). If we combine the phase space with the toy matrix element , we have singularities at , and . Such denominators come e.g. from the interference of diagrams as shown in Fig. 9.
(29)  
The term goes to zero for . Although does not lead to a singularity in the above example, for numerical stability reasons, and having in mind the presence of more complicated matrix elements than our toy example, it is preferable to transform this factor to an expression which is finite in the above limits. Splitting the integration at 1/2 and then doing sector decomposition achieves this goal. The program will do this automatically if the template file contains splitlist={3}, to tell the program that the integration over should be split at 1/2. Of course the singularities at will also be extracted automatically.
Using the command ../launch p params23s35.input t templates23s35.m in the subdirectory general/demos, sector decomposition leads to the result given in Table 8.
1.5705 0.0005  4.3530 0.0025  1.712 0.005  31.040 0.014 
Two massless and one massive particles in the final state
The example in this subsection illustrates the program option to exclude certain parts of the integrand from the decomposition, even though they can become zero at certain values of the integration parameters. This can be useful if a particular term is known not to lead to a singularity. Note that terms with powers are excluded from the decomposition by default.
As an example we pick an integral over , where the line singularity has been remapped already (see appendix, section 8.3).
(30)  
Choosing the option “n” for “no decomposition” in the definition of the integrand for the term in square brackets