Golem95C: A library for one-loop integrals with complex masses
We present a program for the numerical evaluation of scalar integrals and tensor form factors entering the calculation of one-loop amplitudes which supports the use of complex masses in the loop integrals. The program is built on an earlier version of the golem95 library, which performs the reduction to a certain set of basis integrals using a formalism where inverse Gram determinants can be avoided. It can be used to calculate one-loop amplitudes with arbitrary masses in an algebraic approach as well as in the context of unitarity-inspired numerical reconstruction of the integrand.
keywords:NLO computations, One-loop diagrams, Complex masses, Unstable particles
, , , , , ,
NEW VERSION PROGRAM SUMMARY
Manuscript Title: Golem95C: A library for one-loop integrals with complex masses
Authors: G. Cullen, J.-Ph. Guillet, G. Heinrich, T. Kleinschmidt, E. Pilon, T. Reiter, M. Rodgers
Program Title: golem95-1.2.0
Licensing provisions: none
Programming language: Fortran95
Computer: Any computer with a Fortran95 compiler
Operating system: Linux, Unix
RAM: RAM used per integral/form factor is insignificant
Keywords: NLO computations, One-loop diagrams, Complex masses, Unstable particles
Classification: 4.4, 11.1
External routines/libraries: some finite scalar integrals are called from OneLOop [1, 2], the option to call them from LoopTools [3, 4] is also implemented.
Catalogue identifier of previous version: AEEO_v1_0
Journal reference of previous version: Comput. Phys. Commun. 180 (2009) 2317.
Does the new version supersede the previous version?: yes
Nature of problem: Evaluation of one-loop multi-leg integrals occurring in the calculation of next-to-leading order corrections to scattering amplitudes in elementary particle physics. In the presence of massive particles in the loop, propagators going on-shell can cause singularities which should be regulated to allow for a successful evaluation.
Solution method: Complex masses can be used in the loop integrals to stand for a width of an unstable particle, regulating the singularities by moving the poles away from the real axis.
Reasons for the new version: The previous version was restricted to massless particles in the loop.
Summary of revisions: Real and complex masses are supported, a general parameter for the renormalisation scale is introduced, improvements in the caching system and the user interface.
Running time: Depends on the nature of the problem. A single call to a rank 6 six-point form factor at a randomly chosen kinematic point, using complex masses, takes 0.06 seconds on an Intel Core 2 Q9450 2.66 GHz processor.
Collider experiments at the TeV scale, in particular the LHC experiments, are expected to shed light on the mechanism of electroweak symmetry breaking and possibly guide us to a more complete theory of fundamental interactions than the present Standard Model. In order to achieve these goals, predictions for signal as well as background processes should be well under control, necessitating calculations at next-to-leading order (NLO) accuracy or beyond.
Over the last few years, enormous progress has been made to push the calculation of NLO corrections towards a higher number of particles in the final states, i.e. to “multi-leg” amplitudes, in QCD as well as in the electroweak sector. For reviews see e.g. [5, 6].
Nowadays, the efforts are also increasingly being focused on the goals of automating multi-leg one-loop calculations and making them publicly available. Recent public programs with emphasis on multi-leg one-loop calculations are e.g. CutTools , Samurai [8, 9], other public programs which have been optimized for less than four particles in the final state are e.g. FeynArts/FormCalc [3, 10, 11, 4], MCFM , VBFNLO .
An important ingredient for such programs is an integral library containing the one-loop integrals which are the basic building blocks of any one-loop amplitude unless it is calculated purely numerically. Several libraries are available to date: FF [14, 15], Looptools , QCDLoop , OneLOop , golem95 , Hexagon.F . A code for the calculation of one-loop four-point functions with complex masses (D0C) can be found in . The latter has been integrated into the LoopTools library  where the complex version of infrared finite integrals with less than four legs are already implemented. A complete set of scalar four-point integrals, both in dimensional and in mass regularisation and valid also for complex masses can be found in  in analytic form.
The calculation of scalar one-loop integrals has a long tradition of pioneering work, see e.g. [21, 22, 14, 23, 24, 25, 26]. For processes involving unstable particles, these integrals are also required for complex internal masses, in order to be able to work within the so-called “complex-mass scheme” developed in Refs. [27, 28]. For calculations beyond one loop, complex values for invariants derived from external momenta are also required [29, 30], but we will concentrate on one-loop corrections here.
In this article, we present an extension of the tensor and scalar library of Ref.  to integrals with arbitrary masses, in particular also complex masses. Furthermore, we extend the approach which was previously based solely on form factors as building blocks of the amplitude to an approach which is useful in the context of reconstruction of the integrand using -dimensional unitarity.
This article is organized as follows. In Section 2, we review the theoretical background, with particular emphasis on the treatment of potential numerical instabilities. In subsection 2.4 we give an example to demonstrate how the introduction of complex masses can cure Landau singularities stemming from on-shell massive particles in the loop. Subsection 2.5 is dedicated to the new feature of golem95 to be used in the context of a numerical reconstruction of amplitudes at the integrand level. Section 3 gives a brief overview of the software structure, while a detailed description of the individual software components is provided in Section 4. The installation instructions can be found in Section 5, followed by a listing of the examples which are new in this version in Section 6, before we conclude.
2 Theoretical background
The program is an update of the tensor and scalar integral library described in more detail in Ref. , based on the formalism developed in Refs. [31, 32] to reduce tensor integrals to a convenient set of basis integrals. Similar reduction schemes can be found e.g. in Refs. [25, 33, 34, 35, 36, 37, 38]. Here we will describe the theoretical framework only briefly and focus on the new features of the program.
2.1 Form Factors
Tensor integrals can be divided into a part containing the Lorentz structure and a part consisting of scalar quantities, which we call form factors.
We define an -point tensor integral of rank in dimensions as
where , is the loop momentum, and is a combination of external momenta. Using the shift invariant vectors
we can write
The notation stands for the distribution of the Lorentz indices , and the momentum labels to the vectors and metric tensors in all distinguishable ways. denotes an ordered set of propagator labels, corresponding to the momenta forming the kinematic matrix , defined by
The form factors are linear combinations of so-called reduction coefficients derived from the matrix and -point integrals with . The kinematic matrix is related to the Gram matrix ( for ) by
The golem95 program uses the fact that tensor integrals are related to Feynman parameter integrals with Feynman parameters in the numerator. A scalar integral, after Feynman parametrisation, can be written as
where is an integral with Feynman parameters in the numerator. stands for the nearest integer less or equal to and the symbol indicates that powers of the metric tensor are present. Feynman parameter integrals corresponding to diagrams where propagators are omitted, or pinched, with respect to the “maximal” topology can be defined as
The program golem95 reduces the integrals internally to a set of basis integrals, i.e. the endpoints of the reduction (they do not form a basis in the mathematical sense, as some of them are linearly dependent). The choice of the basis integrals can have important effects on the numerical stability in certain kinematic regions, as will be explained below. Our reduction endpoints are 4-point functions in 6 dimensions , which are IR and UV finite, UV divergent 4-point functions in dimensions, and various 2-point and 3-point functions, some of the latter with Feynman parameters in the numerator. This provides us with a convenient separation of IR and UV divergences, as the IR poles are exclusively contained in the triangle functions. Explicitly, our reduction basis is given by integrals of the type
where , as well as with no Feynman parameters in the numerator, and two-point and one-point functions.
Note that and are UV divergent, while can be IR divergent. In the code, the integrals are represented as arrays containing the coefficients of their Laurent expansion in .
We would like to emphasize that the program also can be used as a library for scalar master integrals.
2.3 Treatment of potential numerical instabilities
2.3.1 Spurious singularities due to inverse Gram determinants
Further reduction of the integrals in eqs. (2.2) to integrals with no Feynman parameters in the numerator introduces factors of , i.e. inverse Gram determinants. A particular feature of golem95 is the fact that the above integrals are not reduced to scalar basis integrals in cases where becomes small, thus avoiding problems with small inverse determinants. In these cases, the above integrals are evaluated numerically. As is a dimensionful quantity, the switch to the numerical evaluation of the basis integrals is implemented such that the value of the dimensionless parameter is tested, where
If , the reduction is performed, else the program switches to the direct numerical evaluation of the integral. The default value is . In particular, we use a certain one-dimensional parameter representation here, obtained after performing two integrations analytically. In this way one can use deterministic integration routines, leading to a fast and precise numerical evaluation. The relative error to be achieved in the numerical integration has been set to the default value . If this precision has not been reached, the program will write a message to the file error.txt.
This switch to a direct numerical evaluation will be done automatically for all triangle integrals, and for box integrals with massless propagators and up to three off-shell legs. For box diagrams with massive internal propagators the one-dimensional parameter representation is not yet implemented, but will be provided in a forthcoming version.
In some cases, calculating in double precision Fortran may not be sufficient. The code is designed such that it can be compiled in quadruple precision as well.
2.3.2 Landau Singularities
After Feynman parametrisation and momentum integration, the denominator of any one-loop integral is given by (see eq. (6)). Necessary conditions for a Landau singularity to occur in a one-loop integral thus can be expressed as
The leading Landau singularity, corresponding to and for all , occurs if all particles in the loop go simultaneously on-shell. Sub-leading Landau singularities occur if a sub-matrix of has a vanishing determinant and at least one of the Feynman parameters is zero, corresponding to pinched diagrams. Leading Landau singularities are also called anomalous thresholds in the literature [40, 41], this term stemming from the fact that they only occur in the physically allowed phase space region if a number of rather special kinematical conditions are fulfilled. An example will be given below.
2.4 Complex masses
In order to demonstrate how the introduction of complex masses, standing for a width in the case of unstable particles in the loop, regulates the Landau singularities, we choose the example of a box diagram contributing to the production of a heavy neutral Higgs boson and a pair in gluon fusion in the context of supersymmetry, where the loop contains two squarks (sbottoms) and two neutralinos, as shown in Fig. 1.
Denoting the momentum of the Higgs boson by and the momenta of the two b-quarks by and , the kinematic matrix associated with the four-point diagram is symmetric under (. To make the singularity structure apparent, we will scan the different thresholds and singularities as a function of the invariant . The -quarks are taken to be massless in this example. Similar investigations have been worked out in detail in [42, 43].
Solving the phase space constraint () for leads to the following boundaries for the physically allowed phase space:
The leading Landau singularity is characterised by all particles in the loop going simultaneously on-shell, leading to . The determinant of the kinematic matrix can be written as 
The determinants and correspond to diagrams where the propagators and respectively in Fig. 1 are pinched. The Kaellen function is (minus) the determinant of the kinematic matrix associated with a two-point function where both and are pinched. is the solution of the equation . Choosing the numerical values GeV, GeV, GeV, GeV, fixing to and combining with the phase space constraints of eq. (12) and the requirement , we encounter the following discontinuities, shown in Fig. 2:
a normal threshold at GeV, where vanishes, corresponding to the production of a squark and a neutralino in the cut two-point diagram associated with a pinch of both and in the box diagram.
an anomalous threshold at GeV, corresponding to
an anomalous threshold at GeV, corresponding to
a leading Landau singularity at GeV, corresponding to
Introducing complex masses moves the poles due to propagators going on-shell away from the real axis and therefore regulates the leading Landau singularity, as can be seen from Fig. 3.
We should remark that Landau singularities due to massive (unstable) particles in the loop going on-shell usually behave like as . This is in contrast to cases with several massless particles, where usually has several zero eigenvalues. An interesting case is the six-photon amplitude, where not only , but also its derivative with respect to the invariants involved vanishes at the singular point, meaning that the singularity is not integrable anymore. Introducing an imaginary part as advocated above would certainly not help in this case, but in the six-photon example we are saved by the gauge structure, which leads to numerators taming the singularity structure when the individual contributions are combined to physical helicity amplitudes. This has been worked out in detail in Refs. [44, 45, 5].
2.5 Tensorial Reconstruction of the Integrand
The library golem95 in its original version, when used in amplitude calculations, relies on some user generated code expressing a diagram111Actually, it is not required to organize the calculation in terms of diagrams. The term diagram here refers to a set of terms of an amplitude sharing some common loop propagators. in terms of form factors () as defined in Eq. (3). This requirement of the previous version of golem95 restricted the applicability of the library to algebraic methods for generating the amplitude. In order to be able to use the library in the context of a numerical reconstruction of the integrand, we included new features described in the following.
The general structure of a one-loop diagram can be written as
where . After decomposing the -dimensional vector into its projection onto the physical 4-dimensional Minkowski space and the radial component of the -dimensional orthogonal space, such that , we can express the numerator function in terms of the following tensor structure,
The coefficients can be determined numerically by the algorithm described in Ref. . The contraction of the coefficients with the tensor integrals can be carried out in a process independent and numerical way. Therefore this tensorial reconstruction of the integrand provides a way of processing both diagrams from algebraic constructions and purely numerical input, using the numerator function and the set of denominators as the only common source of information. We have implemented the reconstruction of the coefficients and their contraction with the tensor integrals as a part of golem95. The new interface is described in Section 4.2.
The contraction with the tensor integrals also requires the implementation of integrals with explicit -dependence, which can be computed from the poles of known form factors,
In practice, the above equation simplifies greatly since the only non-zero contributions come from the form factors and . The only non-zero integral with in the numerator is the box
3 Overview of the software structure
The structure of the golem95 program is the following: There are four main directories:
src: the source files of the program
demos: some programs for demonstration
doc: documentation which has been created with robodoc 
test: supplements the demonstration programs, containing files to produce form factors with user-defined kinematics. The user can specify the rank, numerator, numerical point etc. via a steering file.
The subdirectory structure is the same as described in .
4 Description of the individual software components
We focus here on the new features, for more details on the software components which are the same as in version 1.0, we refer to  and to the documentation contained in the program.
4.1 Form factor evaluation
A typical setup for the calculation of form factors is:
…fill matrix …
…evaluate form factors …
The three subroutines act as interfaces to the library. They are used to initialize and manipulate objects and features of the program. They are implemented in the module matrice_s and will be described in more detail in the following. To pass the desired masses and momenta, the user assigns values to the object s_mat, which are the entries of the kinematic matrix given in eq. (4). The form factors are implemented in the respective modules in src/form_factor. We will sketch the internal structure of the library and the distinction between real and complex masses in the following.
4.1.1 The module matrice_s
This module is located in src/kinematic/matrice_s.f90. It is used to reserve and free memory for the kinematic matrix and related objects, as well as for the computation of the inverse matrices. The three macro functions described above are:
- subroutine initgolem95(dim, opt_set)
allocates memory for two arrays s_mat_c and s_mat_r of dimension dim and rank two. These arrays represent internally the complex valued matrix and its real part. Additional memory is reserved for the inverse of and its submatrices as well as for parameters appearing in the reductions of the form factors. A public pointer s_mat is associated with s_mat_c. This object is used to assign the matrix , both in the real mass and complex mass case. An additional argument opt_set can be given, which is an integer array reflecting the numbering of the propagators. The default value is an integer vector of range 1 to dim. This subroutine also initializes the caching system.
- subroutine preparesmatrix()
creates the derived type object s_mat_p. This consists of two pointers to s_mat_c and s_mat_r as well as two integers encoding information about the entries of complex masses and vanishing masses. We describe the usage of this object in more detail in section 4.1.2. The arrays used later for the inverse matrices and sub-matrices are calculated in this routine. If the user defines values for the matrix which are purely real, only the real arrays are needed and calculated.
- subroutine exitgolem95()
deallocates all arrays, nullifies all pointers and clears the cache.
Note that the subroutines initgolem95 and exitgolem95 also have been designed to make the setup more user-friendly, so the calls to the subroutines allocation_s, allocate_cache, init_invs, deallocation_s, clear_cache and the definition of set_ref in version 1.0.0 are now obsolete.
4.1.2 Internal structure
Internally, the library is subdivided into three layers, implementing the reduction from general form factors down to specific integrals. The global derived type object s_mat_p is passed from one layer to the next. The information encoded in the integer bits of s_mat_p, is used to direct the evaluation of the form factors to the specific integrals needed.
Upon a form factor call, reductions are performed in the respective modules located in src/form_factor. These reductions only need information from the inverse matrices and related objects created with preparesmatrix().
From the form factor modules, generic -point functions are called, passing the derived type object s_mat_p. At this level, further reductions might be performed. Up to this point, the only information needed from the kinematic matrix is the position of the zero-mass entries, which is encoded in the integer s_mat_p%b_zero. In the current implementation, we choose massless six-dimensional boxes and massive four-dimensional boxes as our basis of the reduction. A comparison of b_zero with the pinched propagators gives a fast case distinction.
Only in the last step, when the generic -point functions call specific integrals, does the (sub-) matrix need to be passed over. If the latter contains complex masses, which again can be quickly determined by comparing s_mat_p%b_cmplx with the set of remaining propagators, an implementation of the integral for complex masses is called with s_mat_p%pt_cmplx, otherwise the version for real masses is called with s_mat_p%pt_real. This structure of the library allows for an efficient reduction of the form factors to the basis integrals, without a full knowledge of the kinematic matrix .
4.2 Interface for tensorial reconstruction
The tensorial reconstruction described in Section 2.5 has been implemented in two modules, tens_rec and tens_comb. In the simplest case the user only needs to access the function evaluate from the module tens_comb. For more advanced applications of the interface the user will have to use functions and data types of both modules.
4.2.1 The module tens_rec
The module tens_rec is located in src/interface/tens_rec.f90 and contains all routines required for the reconstruction of the tensor coefficients and for the evaluation of reconstructed numerators.
The module tens_rec defines the data types
with components c () to store the coefficients of a numerator of maximum rank . In particular, the component c is a two dimensional complex array which stores the coefficients of the monomials with distinct components of . The first index of the array labels the ways of choosing of the four components of . The second index labels the monomials of the corresponding polynomial. A certain order of the array entries should not be assumed by the user.
Functions and subroutines
In all functions and subroutines which have a parameter numeval it should point to a function with the following signature:
interface function numeval(Q, mu2) use precision_golem, only: ki real(ki), dimension(0:3), intent(in) :: Q real(ki), intent(in) :: mu2 complex(ki) :: numeval end function end interface
For all momenta the order is assumed.
- subroutine reconstruct(numeval,c0,c1,c2)
for reconstructs the tensor coefficients of the numerator function given by numeval. The parameters c0, c1 and c2 are output parameters; c1 and c2 are optional and only available in the subroutines where . The argument c0 is of type coeff_type_ and contains the constant part (with respect to , i.e. in eq. (15)) of the polynomial. The arguments c1 and c2 are of type coeff_type_ and contain the coefficients of the and part of the polynomial222There is no coeff_type_0: complex(ki) is used instead..
- subroutine print_coeffs(coeffs, unit)
is overloaded and can take any derived coefficient type as its first argument. The second argument is optional and has the default value unit=6 (standard output). This routine pretty-prints the coefficients to the given file or device.
- pure function tenseval(Q, coeffs, max_k)
evaluates the polynomial given by the coefficients coeffs for the given real momentum Q. The optional argument max_k is only used internally and should not be assigned by the user.
- pure function ctenseval(Q, coeffs)
evaluates the polynomial given by the coefficients coeffs for the given complex momentum Q.
The subroutines reconstruct and ctenseval can together serve as an implementation of a presampling for programs of unitarity based reductions at the integrand level [7, 8] or as an implementation of a rescue system when the unitarity based method fails due to vanishing Gram determinants. Both options are described in more detail in .
4.2.2 The module tens_comb
The module tens_comb is located in src/interface/tens_comb.f90 and contains all routines for the contraction of reconstructed tensor coefficients with tensor integrals. It also contains the convenience function evaluate which combines the steps required for reconstruction and contraction.
Functions and subroutines
- function evaluate(numeval, momenta, set, rank)
combines tensorial reconstruction and the contraction of the coefficients with the tensor integrals. The argument momenta is of dimension(:,0:3) and contains the momenta . The argument set can either be an integer array or an integer number (generated by the function packb) and denotes the set of pinched propagators; a value of zero can be used if no propagators are pinched. It should be noted that the argument momenta also includes the momenta belonging to pinched propagators. The last argument, rank, is optional. It specifies the rank of the numerator; if omitted it is assumed that .
- function contract_(coeffs, momenta, set)
contracts the coefficients coeffs of type coeff_type_ with the -point rank- tensor integral. The arguments momenta and set are defined as in evaluate with the restriction that set has to be an integer number.
- function contract_s1(coeffs, momenta, set)
contracts the coefficients coeffs of type coeff_type_ with the corresponding -point tensor integral containing in the numerator. The arguments momenta and set are defined as in evaluate with the restriction that set has to be an integer number.
- function contract_s2(coeffs, momenta, set)
contracts the coefficients coeffs of type coeff_type_ with the corresponding -point tensor integral containing in the numerator. The arguments momenta and set are defined as in evaluate with the restriction that set has to be an integer number.
The function evaluate determines from the size of momenta and the number of pinches in set. It then combines the according calls to reconstruct and contract_, contract_s1 and contract_s2. The functions evaluate and contract… require that the matrix s_mat and the cache are set up properly by the sequence which one would have to call when evaluating form factors.
5 Installation instructions
The program can be downloaded as an archive golem95-1.2.0.tar.gz from the following URL: http://projects.hepforge.org/golem/95/. The installation instructions given below also can be found in the Readme file which comes with the code. Information and updates of the program can also be found at http://projects.hepforge.org/golem/trac/wiki/golem95C.
The installation setup is based on autotools .
To install the golem95 library, type the following commands:
./configure [--prefix=mypath] [--precision=quadruple] [FC=compiler] [F77=fortran77compiler]
The --prefix option denotes the installation prefix, under which the directories lib/ and include/ are generated. If no option is given, on a Linux system the configure script would choose --prefix=/usr/local. The argument --precision selects if double or quadruple precision should be used in the library; it should be noted that quadruple precision is not supported by all Fortran compilers and therefore --precision=double is the default value. If the variable FC is not set the first fortran compiler which is automatically detected will be used. Another variable commonly used is FCFLAGS which allows one to pass compiler flags to the Fortran compiler. The fortran77 compiler is only needed to run the demo file demos/demo_LT.f.
In addition, it is possible to call the finite scalar box and triangle integrals
with internal masses from LoopTools. The setup to do so is automated, all the user has to do
for this option is (a) install LoopTools, and (b) use the option
[--with-looptools=your_path_to_libooptools.a] for the configure script, i.e.
type the following commands:
./configure [--prefix=mypath] [--with-looptools=your_path_to_libooptools.a] [--precision=quadruple] [FC=compiler] [F77= fortran77compiler]
Examples can be found in the subdirectory demos. The program demos/demo_cmplx_masses.f90 shows how the library is used to evaluate form factors. For a given matrix with real or complex masses, six-point form factors from rank 0 to rank 6, as well as one- to five-point form factors are evaluated. The results are written to the files test_ff(6)_c/r.txt.
The program demos/demo_tens_rec.f90 together with the module demos/demo_tens_mod.f90 demonstrates the use of the function evaluate. In this example, a toy amplitude is given by a six-point integral where the numerator consists of three propagators which are also present in the denominator. The scalar three-point function resulting from the direct cancellation of three propagators is compared to the results obtained by expanding the numerator into contracted loop and external momenta, leading to rank 6 hexagons if no propagators are cancelled, rank 4 pentagons if one propagator is cancelled, rank 2 boxes if two propagators are cancelled.
The program used to scan the Landau singularity as in Fig. 3 can be found in demos/SusyLandau.f90.
We have presented a program for the numerical evaluation of scalar integrals and tensor form factors entering the calculation of one-loop amplitudes, which is able to provide results for real as well as complex masses in the loop integrals. The program is an extension of on an earlier version of the golem95 library, but now also can be used in the context of a unitarity-inspired numerical reconstruction of the integrand at the tensorial level. Improvements in the caching system and in the user interface also have been made. The program, available at http://projects.hepforge.org/golem/95/, provides a complete library of scalar and tensor integrals up to rank six 6-point functions, including the option of complex masses.
We would like to thank Thomas Hahn for reading the manuscript and useful comments. This research was supported by the UK Science and Technology Facilities Council (STFC) and the Scottish Universities Physics Alliance (SUPA). T.R. has been supported by the Foundation FOM, project FORM 07PR2556.
-  A. van Hameren, C. G. Papadopoulos, and R. Pittau, Automated one-loop calculations: a proof of concept, JHEP 09 (2009) 106, [arXiv:0903.4665].
-  A. van Hameren, OneLOop: for the evaluation of one-loop scalar functions, arXiv:1007.4716.
-  T. Hahn and M. Perez-Victoria, Automatized one-loop calculations in four and D dimensions, Comput. Phys. Commun. 118 (1999) 153–165, [hep-ph/9807565].
-  T. Hahn, Feynman Diagram Calculations with FeynArts, FormCalc, and LoopTools, arXiv:1006.2231.
-  NLO Multileg Working Group Collaboration, Z. Bern et. al., The NLO multileg working group: summary report, arXiv:0803.0494.
-  SM and NLO Multileg Working Group Collaboration, J. R. Andersen et. al., The SM and NLO multileg working group: Summary report, arXiv:1003.1241.
-  G. Ossola, C. G. Papadopoulos, and R. Pittau, CutTools: a program implementing the OPP reduction method to compute one-loop amplitudes, JHEP 03 (2008) 042, [arXiv:0711.3596].
-  P. Mastrolia, G. Ossola, T. Reiter, and F. Tramontano, Scattering AMplitudes from Unitarity-based Reduction Algorithm at the Integrand-level, JHEP 08 (2010) 080, [arXiv:1006.0710].
-  G. Heinrich, G. Ossola, T. Reiter, and F. Tramontano, Tensorial Reconstruction at the Integrand Level, JHEP 10 (2010) 105, [arXiv:1008.2441].
-  T. Hahn, Generating Feynman diagrams and amplitudes with FeynArts 3, Comput. Phys. Commun. 140 (2001) 418–431, [hep-ph/0012260].
-  T. Hahn and M. Rauch, News from FormCalc and LoopTools, Nucl. Phys. Proc. Suppl. 157 (2006) 236–240, [hep-ph/0601248].
-  J. Campbell and R. K. Ellis, Next-to-leading order corrections to W + 2jet and Z + 2jet production at hadron colliders, Phys. Rev. D65 (2002) 113007, [hep-ph/0202176].
-  K. Arnold et. al., VBFNLO: A parton level Monte Carlo for processes with electroweak bosons, Comput. Phys. Commun. 180 (2009) 1661–1670, [arXiv:0811.4559].
-  G. J. van Oldenborgh and J. A. M. Vermaseren, New Algorithms for One Loop Integrals, Z. Phys. C46 (1990) 425–438.
-  G. J. van Oldenborgh, FF: A Package to evaluate one loop Feynman diagrams, Comput. Phys. Commun. 66 (1991) 1–15.
-  R. K. Ellis and G. Zanderighi, Scalar one-loop integrals for QCD, JHEP 02 (2008) 002, [arXiv:0712.1851].
-  T. Binoth, J. P. Guillet, G. Heinrich, E. Pilon, and T. Reiter, Golem95: a numerical program to calculate one-loop tensor integrals with up to six external legs, Comput. Phys. Commun. 180 (2009) 2317–2330, [arXiv:0810.0992].
-  T. Diakonidis, J. Fleischer, T. Riemann, and B. Tausk, A recursive approach to the reduction of tensor Feynman integrals, PoS RADCOR2009 (2010) 033, [arXiv:1002.0529].
-  D. T. Nhung and L. D. Ninh, D0C : A code to calculate scalar one-loop four-point integrals with complex masses, Comput. Phys. Commun. 180 (2009) 2258–2267, [arXiv:0902.0325].
-  A. Denner and S. Dittmaier, Scalar one-loop 4-point integrals, arXiv:1005.2076.
-  G. ’t Hooft and M. J. G. Veltman, Scalar One Loop Integrals, Nucl. Phys. B153 (1979) 365–401.
-  K. Fabricius and I. Schmitt, Calculation Of Dimensionally Regularized Box Graphs In The Zero Mass Case, Z. Phys. C3 (1979) 51–53.
-  W. Beenakker and A. Denner, Infrared Divergent Scalar Box Integrals With Applications In The Electroweak Standard Model, Nucl. Phys. B338 (1990) 349–370.
-  A. Denner, U. Nierste, and R. Scharf, A Compact expression for the scalar one loop four point function, Nucl. Phys. B367 (1991) 637–656.
-  Z. Bern, L. J. Dixon, and D. A. Kosower, Dimensionally regulated one loop integrals, Phys. Lett. B302 (1993) 299–308, [hep-ph/9212308].
-  Z. Bern, L. J. Dixon, and D. A. Kosower, Dimensionally regulated pentagon integrals, Nucl. Phys. B412 (1994) 751–816, [hep-ph/9306240].
-  A. Denner, S. Dittmaier, M. Roth, and D. Wackeroth, Predictions for all processes 4 fermions + gamma, Nucl. Phys. B560 (1999) 33–65, [hep-ph/9904472].
-  A. Denner, S. Dittmaier, M. Roth, and L. H. Wieders, Electroweak corrections to charged-current 4 fermion processes: Technical details and further results, Nucl. Phys. B724 (2005) 247–294, [hep-ph/0505042].
-  S. Actis, G. Passarino, C. Sturm, and S. Uccirati, Two-Loop Threshold Singularities, Unstable Particles and Complex Masses, Phys. Lett. B669 (2008) 62–68, [arXiv:0809.1302].
-  G. Passarino, C. Sturm, and S. Uccirati, Higgs Pseudo-Observables, Second Riemann Sheet and All That, Nucl. Phys. B834 (2010) 77–115, [arXiv:1001.3360].
-  T. Binoth, J. P. Guillet, G. Heinrich, E. Pilon, and C. Schubert, An algebraic/numerical formalism for one-loop multi-leg amplitudes, JHEP 10 (2005) 015, [hep-ph/0504267].
-  T. Binoth, J. P. Guillet, and G. Heinrich, Reduction formalism for dimensionally regulated one-loop N-point integrals, Nucl. Phys. B572 (2000) 361–386, [hep-ph/9911342].
-  G. Duplancic and B. Nizic, Reduction method for dimensionally regulated one-loop N- point Feynman integrals, Eur. Phys. J. C35 (2004) 105–118, [hep-ph/0303184].
-  W. T. Giele and E. W. N. Glover, A calculational formalism for one-loop integrals, JHEP 04 (2004) 029, [hep-ph/0402152].
-  F. del Aguila and R. Pittau, Recursive numerical calculus of one-loop tensor integrals, JHEP 07 (2004) 017, [hep-ph/0404120].
-  A. van Hameren, J. Vollinga, and S. Weinzierl, Automated computation of one-loop integrals in massless theories, Eur. Phys. J. C41 (2005) 361–375, [hep-ph/0502165].
-  A. Denner and S. Dittmaier, Reduction schemes for one-loop tensor integrals, Nucl. Phys. B734 (2006) 62–115, [hep-ph/0509141].
-  J. Fleischer and T. Riemann, A complete algebraic reduction of one-loop tensor Feynman integrals, arXiv:1009.4436.
-  A. I. Davydychev, A Simple formula for reducing Feynman diagrams to scalar integrals, Phys. Lett. B263 (1991) 107–111.
-  J. Bjorken and S. Drell, Relativistic Quantum Field Theory. McGraw-Hill, New York, 1965.
-  S. Goria and G. Passarino, Anomalous Threshold as the Pivot of Feynman Amplitudes, Nucl. Phys. Proc. Suppl. 183 (2008) 320–325, [arXiv:0807.0698].
-  A. Denner, S. Dittmaier, and T. Hahn, Radiative corrections to Z Z Z Z in the electroweak standard model, Phys. Rev. D56 (1997) 117–134, [hep-ph/9612390].
-  F. Boudjema and L. D. Ninh, b anti-b Higgs production at the LHC: Yukawa corrections and the leading Landau singularity, Phys. Rev. D78 (2008) 093005, [arXiv:0806.1498].
-  Z. Nagy and D. E. Soper, Numerical integration of one-loop Feynman diagrams for N- photon amplitudes, Phys. Rev. D74 (2006) 093006, [hep-ph/0610028].
-  C. Bernicot and J. P. Guillet, Six-Photon Amplitudes in Scalar QED, JHEP 01 (2008) 059, [arXiv:0711.4713].
-  http://www.xs4all.nl/rfsber/robo/robodoc.html.
-  G. V. Vaughan, B. Elliston, T. Tromey, and I. L. Taylor, GNU Autoconf, Automake, and Libtool: Expert insight into porting software and building large projects using GNU Autotools. New Riders, Indianapolis, 2000.