SecDec: A general program for sector decomposition

SecDec: A general program for sector decomposition

Jonathon Carter Gudrun Heinrich IPPP, Department of Physics, University of Durham, Durham DH1 3LE, UK Max-Planck-Institute for Physics, Föhringer Ring 6, 80805 Munich, Germany
Abstract

We present a program for the numerical evaluation of multi-dimensional 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 multi-loop 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, multi-dimensional parameter integrals, numerical integration

IPPP/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, multi-dimensional 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. multi-loop 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 unit-hypercube. Those integrals are evaluated numerically by Monte Carlo integration.
Restrictions: Depending on the complexity of the problem, limited by memory and CPU time. Multi-scale 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 WRITE-UP

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 multi-scale 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 two-loop 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 three-loop and four-loop 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 multi-loop 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 four-particle 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 Frixione-Kunszt-Signer 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 non-linear 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. Half-integer 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 double-indices 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 semi-definite 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 semi-definite 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 light-like, 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 multi-loop 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 so-called 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 tree-like 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 order111To account for half-integer 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

Figure 1: Directory structure of the SecDec 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.

Figure 2: Flowchart showing the main steps the program performs to produce the result files. In each of the subdirectories loop or general, the file Template.m can be used to define the integrand. The produced files are written to a subdirectory created according to the settings given in param.input. By default, a subdirectory with the name of the graph or integrand is created to store the produced functions. This directory will contain subdirectories according to the pole structure of the integrand. The perl scripts (extension .pl) are steering the various steps to be performed by the program.

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 Cuba-2.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 user-defined 3-loop example, and a 4-loop example to be written to the scratch disk is given in Fig. 3.

Figure 3: Example for a directory structure created by running the loop demo programs NPbox, QED, ggtt1, A61. A four-loop example defined by the user to be written to the scratch disk is also shown.

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 book-keeping 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 non-logarithmic poles will disappear upon -expansion. In particular, for half-integer powers, the labelling does not correspond to “true” poles, but rather to terms which can be cast into functions like , which are well-defined 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.

Figure 4: Example for a directory tree corresponding to the pole structure of the graph QED contained in the demo programs.

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, non-standard 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 non-planar massless on-shell two-loop 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

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

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

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

  4. 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 epsilon-expansions and writes the Fortran functions. Depending on the “exe-flag” specified in the parameter file (see below for a detailed explanation of the flag), this script also launches the compilation and the numerical integrations.

  5. 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.

  6. 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: Non-planar massless two-loop box

The non-planar massless two-loop box is a non-trivial example, as the sector decomposition applied to the standard representation, produced by combining all propagators simultaneously with Feynman parameters, exhibits “non-logarithmic 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 re-parametrised such that poles of this type do not occur, because these poles require complicated subtraction terms which slow down the calculation. Non-linear 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 non-planar massless two-loop box in appendix 8.2, and the template to calculate the graph in this way can be found in SecDec/general/demos.

Figure 5: The non-planar two-loop box, called NPbox in example 1.

To obtain results for the non-planar massless two-loop 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 non-planar 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
Table 1: Numerical results for the points and of the massless non-planar double box.

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 non-planar box graphs.

5.1.2 Example 2: Planar two-loop ladder diagram with massive on-shell 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 on-shell particles, depicted by solid (blue) lines in Fig. 6.

Figure 6: Blue (solid) lines denote massive particles.

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 third-last 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
Table 2: Numerical results up to order for the points and of the two-loop ladder diagram shown in Fig. 6. An overall factor of is not included in the numerical result.

5.1.3 Example 3: Non-planar two-loop diagrams with two massive on-shell legs

This example gives results for two non-planar 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.

ggtt2
Figure 7: Non-planar graphs occurring in the calculation of at NNLO. Blue (solid) lines denote massive particles.
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
Table 3: Numerical results for the diagrams shown in Fig. 7. The finite diagram ggtt1 has been calculated up to order . An overall factor of is extracted.

5.1.4 Example 4: A rank one tensor two-loop box

In order to demonstrate how to run the program for integrals with non-trivial numerators, we give the example of a rank one planar massless on-shell two-loop 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
Table 4: Numerical results for the point of the rank one two-loop ladder diagram given by eq. (LABEL:eq:dbr1). An overall factor of has been extracted.

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 three-loop 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.

Figure 8: The three-loop vertex diagram with the dotted propagator raised to the power .
0.16666 1.8334 18.123 125.32 889.96 5325.3
Table 5: Numerical results for the diagram shown in Fig. 8 with the dotted propagator raised to the power . The errors are below one percent.

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
Table 6: Results for the hypergeometric function at . The timings in the last column are the ones for the numerical integration. The time taken for decomposition, subtraction and -expansion was 11 seconds.

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
Table 7: Results for the hypergeometric function at .

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 pseudo-state 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.

Figure 9: Interference of diagrams leading to factors of in the denominator.
(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
Table 8: Results for the integral given by eq. (29). The factor is not included in the numerical result.

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