SecDec: A tool for numerical multiloop calculations
Abstract
The version 2.0 of the program SecDec is described, which can be used for the extraction of poles within dimensional regularisation from multiloop integrals as well as phase space integrals. The numerical evaluation of the resulting finite functions is also done by the program in an automated way, with no restriction on the kinematics in the case of loop integrals.
1 Introduction
Nowadays we are in the long awaited situation of being confronted with a wealth of high energy collider physics data, enabling us to explore physics at the TeV scale. However, for the analysis and interpretation of these data, precise theory predictions are mandatory. In most cases, this means that calculations beyond the leading order in perturbation theory are necessary. Such calculations either involve integrations over loop momenta for the virtual corrections, or integrations over phase spaces for the real radiation corrections. In both cases, multidimensional parameter integrals need to be evaluated, which can contain ultraviolet singularities and, in the presence of massless particles, contain soft and/or collinear singularities. These singularities, if regulated by dimensional regularisation, appear as poles in , but factorising the poles from complicated multiloop or multiparameter integrals is a highly nontrivial task. The program SecDec, presented in [1, 2], performs this task in an automated way, based on the algorithm of sector decomposition [3, 4, 5]. Other implementations of sector decomposition in public programs can be found in [6, 7, 8]. However, the latter programs, including SecDec 1.0 [1], are not designed to cope with kinematics where physical thresholds can be present in the integration region. The program SecDec 2.0 [2] is able to achieve this task, by an automated deformation of the integration contour into the complex plane [9, 10, 11, 12]. In this talk the emphasis is on the presentation of the new features of the program SecDec.
2 The algorithm
2.1 General framework
The decomposition algorithm is described in detail in [13], here we only describe the basic concepts. Consider an dimensional parameter integral, for example
(1) 
where are polynomial functions which do not vanish for . However, we would like to factorize the integrand such that the term in square brackets is nonvanishing in the limit . This can be achieved by a decomposition into sectors where and are ordered: we multiply equation (1) with and substitute in the first sector and in the second sector, to arrive at
(2) 
where are functions of the transformed variables. As are nonvanishing at the origin, the terms in square brackets cannot lead to singularities anymore; the singularities are factored out into the monomials of type instead, and subtraction of the poles and expansion in are straightforward in this form. In more complicated cases, the decomposition procedure may have to be iterated to achieve a full factorisation. A detailed description can be found e.g. in [3, 13].
2.2 Loop integrals
A scalar Feynman integral in dimensions at loops with propagators, where the propagators can have arbitrary, not necessarily integer powers , has the following representation in momentum space:
(3) 
where the are linear combinations of external momenta and loop momenta . Introducing Feynman parameters leads to
(4) 
where is a polynomial of degree and is of degree in the Feynman parameters. The Lorentz invariants which can be formed from the external momenta of the diagram, as well as propagator masses, are contained in . As a simple example, consider the function for a massless oneloop box diagram:
(5) 
The functions and can also be constructed directly from the topology of the corresponding Feynman graph [14, 15], and the implementation of this construction in SecDec version 2.0 is one of the new features of the program.
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 . is a positive semidefinite function. A vanishing function is related to the UV subdivergences of the graph. 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. 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 a divergence are given by the Landau equations [16]:
(6)  
(7) 
If all kinematic invariants formed by external momenta are of the same sign, the necessary condition for an IR divergence can only be fulfilled if some of the parameters go to zero. These singularities can be regulated by dimensional regularisation and factored out of the function using sector decomposition. The same holds for dimensionally regulated UV singularities contained in . However, after these singularities, appearing as poles in 1/, have been extracted using sector decomposition, for nonEuclidean kinematics we are still faced with integrable singularities related to kinematic thresholds. How we deal with these singularities will be described in Section 4.4.
2.3 General parameter integrals
The program can also deal with parameter integrals which are more general than the ones related to multiloop integrals, for example phase space integrals involving massless particles, where the regions in phase space corresponding to soft and/or collinear configurations lead to singularities which can be extracted as poles in . 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. We assume that the singularities are regulated by noninteger powers of the integration parameters, where the noninteger part is the of dimensional regularisation or some other regulator. The general form of the integrals is
(8) 
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.
3 The SecDec program
The program consists of two parts, an algebraic part and a numerical part. The algebraic part uses code written in Mathematica [17] 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 or C++ functions forming the coefficient of each term in the Laurent series in are integrated using the Monte Carlo integration programs contained in the Cuba library [18, 19], or Bases [20]. The different subtasks are handled by perl scripts. The flowchart of the program is shown in Fig. 1 for 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 ) and src/util (miscellaneous useful functions). The documentation, created by robodoc [21] 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.
In order to use the program, 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 and propagators, resp. of the topology; optionally numerator, nonstandard propagator powers, spacetime dimensions

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

The program comes with example input and template files in the subdiretories loop/demos respectively general/demos, described in detail in [1].
4 Installation and usage
4.1 Installation
The program can be downloaded from
http://projects.hepforge.org/secdec.
Unpacking the tar archive via tar xzvf SecDec.tar.gz will create a directory called SecDec with the subdirectories as described above. Then change to the SecDec directory and run ./install.
Prerequisites are Mathematica, version 6 or above, perl (installed by default on most Unix/Linux systems), a Fortran compiler (e.g. gfortran, ifort), or a C++ compiler if the C++ option is used.
4.2 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 specify the integrand 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. 
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.

If the jobs have been submitted to a cluster, when all jobs have finished, execute the command ./results.pl [d myworkingdir p myparamfile]. This will create the files containing the final results in the graph subdirectory specified in the input file.


After the calculation and the collection of the results is completed, you can use the shell command ./launchclean[graph] to remove obsolete files.
It should be mentioned 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 submit the jobs which are expected to take longest.
4.3 New features
Version 2.0 of SecDec contains the following new features, some of which will be illustrated by examples in Section 5. More details are given in [2].

Multiscale loop integrals can be evaluated without restricting the kinematics to the Euclidean region.

The possibility to loop over ranges of parameter values is automated, with the option of outputting results in a format suitable for plotting.

The user can define functions at a symbolic level and specify them only later after the integrand has been transformed into a set of finite functions for each order in .

The regulator of the parameter integrals can be different from the dimensional regulator . This is particularly useful to define e.g. measurement functions at a later stage of the calculation.

For scalar multiloop integrals, the integrand can be constructed directly from the topology of the diagram. The user only has to provide the labels of the vertices connected by the propagators and the propagator masses, but does not have to provide the momentum flow.

The files for the numerical integration of multiscale loop integrals are written in C++ rather than Fortran. For integrations in Euclidean space, both the Fortran and the C++ versions are supported.

Both the algebraic and the numerical part allow full parallelisation.
4.4 Implementation of the contour deformation
Unless the function in equation (4) is of definite sign for all possible values of invariants and Feynman parameters, the integrand of a multiloop integral will vanish within the integration region on a hypersurface given by the solutions of the Landau equations (6),(7). However, we can avoid the poles on the real axis by a deformation of the integration contour into the complex plane. As long as the deformation is in accordance with the causal prescription of the Feynman propagators, and no poles are crossed while changing the integration path, we can make use of Cauchy’s theorem to choose an integration contour such that the integration is convergent. The prescription for the Feynman propagators tells us that the contour deformation into the complex plane should be such that the imaginary part of should always be negative. For real masses and Mandelstam invariants , the following Ansatz [9, 10, 11, 12] is convenient:
(9) 
In terms of the new variables, we thus obtain
(10) 
such that acquires a negative imaginary part of order .
The convergence of the numerical integration can be improved significantly by choosing an “optimal” value for . Values of which are too small lead to contours which are too close to the poles on the real axis and therefore lead to bad convergence. Too large values of can modify the real part of the function to an unacceptable extent and could even change the sign of the imaginary part if the terms of order get larger than the terms linear in . This would lead to a wrong result. Therefore we implemented a fourstep procedure to optimize the value of , consisting of

ratio check: To make sure that the terms of order in equation (10) do not spoil the sign of the imaginary part, we evaluate the ratio of the terms linear and cubic in for a quasirandomly chosen set of sample points to determine a maximal allowed .

modulus check: The imaginary part is vital at the points where the real part of is vanishing. In these regions, the deformation should be large enough to avoid large numerical fluctuations due to a highly peaked integrand. Therefore we check the modulus of each subsector function at a number of sample points, and pick the fraction of the value of maximising the function with the minimal modulus, i.e. the value of lambda which keeps furthest from zero.

individual adjustments: If the values of are very different in magnitude, it can be convenient to have an individual parameter for each subsector function and each Feynman parameter .

sign check: After the above adjustments to have been made, the sign of Im is again checked for a number of sample points. If the sign is ever positive, this value of is disallowed.
5 Examples and results
In this example, we will demonstrate three of the new features of the SecDec 2.0 program: the construction of directly from the topology of the graph, the evaluation of the graph in the physical region, and how results for a whole set of different numerical values for the invariants can be produced and plotted in an automated way. We will use the twoloop diagram shown in Fig. 2 as an example. Numerical results for this diagram have been produced in [22, 23], analytical ones in [24, 25].
5.1 Topologybased construction of the integrand
The template file templateP126.m in the demos subdirectory
contains the following lines:
where each entry in corresponds to a propagator of the diagram; the first entry is the mass of the
propagator, and the second entry contains the labels of the two vertices which the propagator connects.
Note that if an external momentum is flowing into the vertex, the vertex also must have the label .
For vertices containing only internal propagators the labeling is arbitrary.
The onshell conditions in the above example state that
.
5.2 Results in the physical region
To run this example, from the loop directory, issue the command ./launch d demos p paramP126.input t templateP126.m. The timings for the finite part and a relative accuracy better than 1%, using Cuba3.0 [19], are about 100 seconds per point on an Intel(R) Core i7 CPU at 2.67GHz, where the timings are very similar far from the threshold and close to threshold.
5.3 Producing data files for sets of numerical values
To loop over a set of numerical values for the invariants and
once the C++ files are created, issue the command
./multinumerics.pl d demos p multiparamP126.input.
This will run the numerical integrations for the values of and specified
in the file demos/multiparamP126.input.
The files containing the results will be found in demos/2loop/P126, and the
files p2.gpdat, p1.gpdat and p0.gpdat will contain the
data files for each point, corresponding to the coefficients of
and respectively.
These files can be used to plot the results against the analytic results using gnuplot.
This will produce the files P126R0.ps, P126I0.ps which will look like Fig. 3.
6 Conclusions
We have presented the program SecDec 2.0, which can be used to factorise dimensionally regulated singularities and numerically calculate multiloop integrals in an automated way. As a new feature of the program, it now can deal with physical kinematics, i.e. is not restricted to the Euclidean region anymore. A new construction of the integrand, based entirely on topological rules, is also included, and the new features are demonstrated for the examples of a massive twoloop threepoint function. In addition, the program can produce numerical results for more general parameter integrals, as they occur for example in phase space integrals for multiparticle production with several unresolved massless particles, and offers the possibility to include symbolic functions which can be used to define measurement functions like e.g. jet algorithms at a later stage. We are looking forward to applications of the program for the calculation of higher order corrections to various observables.
Acknowledgements
J.C. was supported by the British Science and Technology Facilities Council (STFC). We also acknowledge the support of the Research Executive Agency (REA) of the European Union under the Grant Agreement number PITNGA2010264564 (LHCPhenoNet).
References
Footnotes
 Based on a talk given at the 14th International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT), Uxbridge, London, UK, September 2011.
References
 Carter J and Heinrich G 2011 Comput.Phys.Commun. 182 1566–1581 (Preprint 1011.5493)
 Borowka S, Carter J and Heinrich G 2012 (Preprint 1204.4152)
 Binoth T and Heinrich G 2000 Nucl. Phys. B585 741–759 (Preprint hepph/0004013)
 Roth M and Denner A 1996 Nucl. Phys. B479 495–514 (Preprint hepph/9605420)
 Hepp K 1966 Commun. Math. Phys. 2 301–326
 Smirnov A, Smirnov V and Tentyukov M 2011 Comput.Phys.Commun. 182 790–803 (Preprint 0912.0158)
 Bogner C and Weinzierl S 2008 Comput. Phys. Commun. 178 596–610 (Preprint 0709.4092)
 Gluza J, Kajda K, Riemann T and Yundin V 2011 Eur.Phys.J. C71 1516 (Preprint 1010.1667)
 Soper D E 2000 Phys. Rev. D62 014009 (Preprint hepph/9910292)
 Nagy Z and Soper D E 2006 Phys. Rev. D74 093006 (Preprint hepph/0610028)
 Binoth T, Guillet J P, Heinrich G, Pilon E and Schubert C 2005 JHEP 10 015 (Preprint hepph/0504267)
 Anastasiou C, Beerli S and Daleo A 2007 JHEP 05 071 (Preprint hepph/0703282)
 Heinrich G 2008 Int. J. Mod. Phys. A23 1457–1486 (Preprint 0803.4177)
 Nakanishi N 1971 Graph Theory and Feynman Integrals (Gordon and Breach, New York)
 Tarasov O V 1996 Phys. Rev. D54 6479–6490 (Preprint hepth/9606018)
 Landau L D 1959 Nucl. Phys. 13 181–192
 Wolfram S , Mathematica, Copyright by Wolfram Research
 Hahn T 2005 Comput. Phys. Commun. 168 78–95 (Preprint hepph/0404043)
 Agrawal S, Hahn T and Mirabella E 2011 (Preprint 1112.0124)
 Kawabata S 1995 Comp. Phys. Commun. 88 309–326
 Slothouber F and et al http://www.xs4all.nl/ rfsber/Robo/robodoc.html
 Bonciani R, Mastrolia P and Remiddi E 2004 Nucl.Phys. B690 138–176 (Preprint hepph/0311145)
 Ferroglia A, Passera M, Passarino G and Uccirati S 2004 Nucl.Phys. B680 199–270 (Preprint hepph/0311186)
 Fleischer J, Kotikov A and Veretin O 1998 Phys.Lett. B417 163–172 (Preprint hepph/9707492)
 Davydychev A I and Kalmykov M 2004 Nucl.Phys. B699 3–64 (Preprint hepth/0303162)