LHC phenomenology at next-to-leading order QCD: theoretical progress and new results
In this talk I report on recent developments and results relevant for LHC phenomenology at next-to-leading order QCD. Feynman diagrammatic and unitarity based methods have both seen considerable improvements and new ideas recently. Current approaches point towards automated evaluation of one-loop multi-particle amplitudes. Many results for notoriously difficult processes are under construction by several groups and, given the enormous recent progress, it can be concluded that LHC phenomenology at full next-to-leading order level will become the standard approximation soon.
With the start of the Large Hadron Collider (LHC) at CERN we enter exiting times in particle physics. Once the machine is fully operating we will be able to explore an energy range which has never been reached in a laboratory. Many people in this audience are involved in the enormous challenge to master the huge amount of data which will be produced during collisions. I will focus in this talk on the theoretical challenge of predicting LHC phenomena as precisely as possible using our good old Standard Model. We all very much hope for a discrepancy between well-understood data sets and sufficiently precise theoretical predictions. Such a discrepancy will be a signal of new and maybe completely unexpected phenomena which will help us to push our understanding of nature.
What should be expected from the LHC? The first and most important issue is the test of the Higgs mechanism within the Standard Model which predicts the prominent Higgs boson. The LEP experiments set a lower bound of 114.4 GeV on its mass from direct searches. As a Higgs boson with a mass higher than about 200 GeV would induce radiative corrections which are in conflict with electroweak precision measurements, the Standard Model Higgs boson is pretty much nailed down to a low mass window. Ironically, this window is not as easy to close as the mass range above the two Z boson threshold of about 180 GeV, where two muon pairs would provide a gold-plated discovery mode for the Higgs boson. The identification of any signal with the Standard Model Higgs boson will also entail detailed measurements of quantum numbers like spin and CP properties [1, 2]. This necessary homework might take even longer than the time we need to find new signals (if we are lucky, of course!). The Standard Model has some deficits, e.g. it does not contain a natural dark matter candidate, it has no direct relation to gravity and it is sensitive to quadratic renormalisation effects. Although the latter is in principle not a fundamental problem for a renormalisable theory like the Standard Model, a popular extension is supersymmetry (susy) which would not only cure the problem above a sensible susy breaking scale of about 1 TeV, but would even allow to view the SM as a low energy limit of some string theory111If it is ever possible to verify phenomenologically whether string theory is realised in nature is of course a different question. It is allowed to be sceptical here.. Supersymmetric extensions of the Standard Model have been extensively studied in the last years, for an overview see [3, 4, 5]. Other extensions embed the SM in higher dimensions [6, 7, 8] or view the Higgs boson as a Goldstone boson to explain its relatively low mass, i.e. little Higgs models . A further possibility is a strongly interacting theory which might explain the mass of the W and Z-bosons even without a fundamental Higgs boson. There is of course another disturbing logical option: we see neither a hint of the Higgs boson nor any sign of BSM physics. Although this would be a nightmare for experimentalists it is actually an exciting option for theoreticians. Such a scenario can only be realised if something ”invisible” disguises a low mass Higgs boson or if there are some strange interactions active in or beyond the electroweak sector which fake a light Higgs boson in the precision measurements by a yet unknown effect. Either we have to add an invisible sector to the Standard Model or we need to question our quantum field theoretical description thereof, both would be truly exciting!
The discovery potential of the Higgs boson has been investigated thoroughly in the experimental studies of ATLAS and CMS [10, 11] and signal significance plots like the one in Fig. 1 show the importance and relevance of different search channels depending on the Higgs mass. Note that in most experimental studies K-factors, i.e. the effect of higher order corrections, have not been taken into account yet because higher order corrections to many background processes were either not known or the computer codes of known calculations are not flexible enough. Higher order corrections have been considered where available in recent experimental studies of the LHC physics potential .
Preferably discoveries at the LHC should not depend on the theoretical status of theoretical simulations. For example the diphoton decay of the Higgs boson should lead to a clear peak over the falling continuous background in the invariant mass distribution. This background should be measurable from the sidebands. Of course, a precise prediction is still needed for realistic studies before measurement [15, 16] and apart from that, who knows, whether the diphoton channel is not further contaminated by new physics? Other discovery modes, e.g. the discovery mode , has no clear signal and background separation and background subtraction has to rely on extrapolations of background to signal regions which do depend on the precision of theoretical input.
In what concerns the evaluation of leading order estimates for LHC cross sections based on tree-level matrix elements the situation is quite satisfying as many tools have been developed in the pre-LHC era, e.g. Alpgen , COMPHEP , Amegic++ , COMIX , FeynArts/FormCalc , GRACE , HELAC [24, 25, 26], Madgraph/Madevent [27, 28], Whizard , and most of them are publicly available. However, any partonic fixed order calculation is of limited use for experimental studies which need realistic hadronic final states. This is the reason why merging of partonic amplitudes with parton showers has been an active field recently. The standard Monte-Carlo tools which contain partonic showering and include a hadronisation model are currently Pythia , Herwig  and Sherpa . Evidently, the more hard matrix elements are included (or better can be recursively evaluated) in these tools, the more accurate is the description of shapes of distributions and jet structure. In this respect Sherpa is especially promising, as it is based on the matrix element generator Amegic++ . On the other hand, the overall normalisation is predicted only very unreliably as long as virtual higher order corrections are not taken into account. This brings us to the main topic of the talk which is the computation of these higher order corrections.
2 Framework for QCD computations
The evaluation of production cross sections in hadronic collisions is based on the factorisation property of QCD. Whenever the hard scale of the process is significantly higher than the hadronisation scale of around one GeV a differential cross section can be obtained by summing incoherently over the initial partons inside the hadron and can be decomposed into long and short distant contributions.
Here stands for some partonic -particle final state. The parton distribution functions, , which stand for the probability to find a certain parton in an incoming hadron, have to be determined experimentally at a certain scale. Their scale dependence on the other hand is determined perturbatively through the DGLAP-evolution equations.
The factorisation scale () dependence is stemming from the infrared (IR) structure of initial state singularities. In contrast, the renormalisation scale () dependence of the strong coupling constant is an ultra-violet (UV) effect. A hadronic cross section in leading order (LO) perturbation theory is thus plagued by a logarithmic dependence on these scales. This dependence is aggravated when the number of coloured partons increases due to higher powers of . A LO computation in QCD has thus only limited predictive power. By including next-to-leading order (NLO) corrections, the leading logarithmic dependence cancels and one obtains a much milder dependence on the unphysical scales.
The evaluation of NLO corrections contains a loop and real emission part. Schematically
The loop contribution has the same particle content and kinematics as the LO contribution. The real-emission contribution is tree-like but contains an additional parton. Apart from initial state collinear divergences, which are absorbed by the parton distribution functions, the IR divergences cancel between the virtual corrections and the integration of the real emission contributions over the soft/collinear phase space regions. In actual calculations the IR divergences are typically treated by subtracting terms which have the same IR limit as the real emission corrections and do allow for the analytic soft/collinear phase space integration. Adding back the integrated terms allows to represent the (NLO virtual) and (NLO real) part as independently finite contributions, which is mandatory in numerical implementations.
The most wildly used variant of these subtraction methods is the so-called dipole subtraction formalism developed by Catani and Seymour  which has been generalised to the case of massive partons in . We see that the ingredients of a NLO QCD calculation consist of (i) tree level amplitudes, (ii) one-loop amplitudes and (iii) dipole subtractions terms, . As was pointed out above, many automated tools are available for the evaluation of the tree-level contribution. Moreover, different groups have recently implemented the dipole formalism in a process independent way [35, 36, 37, 38], and public codes exist which help to perform this labour intensive part of an NLO computation. The real bottleneck of an automated approach for NLO QCD computations is thus the evaluation of the one-loop matrix elements. This will be the main subject of the remainder of this talk.
3 New methods to compute one-loop amplitudes
3.1 Overview of achieved results
The complexity of loop calculations grows rapidly with the number of external legs. What concerns computations with two initial and two final particles, basically all relevant LHC processes are evaluated meanwhile at NLO in QCD and public codes, like for example MCFM [39, 40], the PHOX programs for observables including photons, and [email protected] for weak boson fusion processes , are available for experimental studies on the partonic level. Methods to combine fixed order calculations with parton showers have been worked out and implemented for many important processes, see for example [email protected] , POWHEG , GRACE  and [46, 47]. Also the SANC collaboration has promoted their framework to include QCD processes and progress has also been reported during this conference [48, 49, 50]. A relevant gap of processes has been closed recently. The off-shell vector boson pair production cross sections for are available now as public codes GG2WW [51, 52] and GG2ZZ . These processes describe a relevant background component to also below threshold. These processes are induced by loop diagrams but are conceptually of LO type. Technically one has to square one-loop amplitudes which is numerically more delicate than computing interference terms between tree and loop contributions.
Adding an additional final state particle increases the complexity of one-loop computations already considerably, and before 2005 only a limited number of processes had been evaluated and implemented in computer programs. NLO predictions have been existing for , [55, 56], in gluon fusion  and weak boson fusion (WBF) [58, 59], and [60, 61]. Between 2005 and this date most of the relevant processes for LHC have been completed and the list of accomplished tasks contains also the processes [62, 63], in WBF222The six external particles are connected via at most 5-point one-loop topologies. [64, 65, 66], , [68, 69], [70, 71], [72, 73] and [74, 75, 76]. For the important process also electroweak [77, 78] and sub-leading interference terms have been evaluated [79, 80, 81, 82]. A discussion of the phenomenological importance of the different processes can be found in . Note that the evaluation of the relevant amplitudes is only about half of the cake. The implementation of a full NLO cross section into a reliable computer code is also a considerable effort due to numerical issues which will be discussed below. Further note that most of the corresponding codes are not public.
The time-scale of all these calculations were of the order of person-years and it is no wonder that only a very limited number of results which deal with NLO LHC processes can be found in the literature[84, 85, 86, 87, 88]. In the remainder of this section I will focus on the recent developments in this direction and report on progress which has been achieved by different groups and different methods recently. The two main approaches are based on either Feynman diagrams or on unitarity cuts. Before discussing these two approaches separately let us shortly comment on two important commonly used techniques.
3.2 Colour and helicity management
When evaluating gauge boson amplitudes or dealing with massless fermions it is very useful to work with spinor helicity methods, for an introduction see . It is well known that tree-level amplitudes are most efficiently represented in this way. The same is true for loop amplitudes. A massless Dirac spinor is already defined by two helicity states defined by
Helicity amplitudes can be written in terms of spinor products which are complex numbers:
Massless gauge bosons like gluons and photons can be expressed using the same building blocks
By construction one works with the two physical degrees of freedom. The auxiliary vector defines an axial gauge. In full amplitudes its dependence drops out due to gauge invariance. The generalisation to massive particles is also well-understood. A public implementation of spinor helicity methods can be found in .
Once the helicity amplitudes are evaluated, the squared amplitude can be obtained numerically. If there are helicity amplitudes labelled by one only has to evaluate products
otherwise terms have to be evaluated. To give an example, the -gluon amplitude can be decomposed into helicity amplitudes. Due to parity invariance only have to be evaluated. Nonetheless we see an exponential growth of terms in N. Another source for exponential growth in N-point amplitudes is the colour structure. A widely used colour decomposition is based on the following two elementary rules valid for an algebra with generators in the fundamental representation:
Any colour structure can be mapped to simple Kronecker deltas which indicate the colour flow. This colour flow decomposition, described for example in [91, 92], is going back to tHooft who applied a double line notation for gluons in the context of the expansion. Any amplitude with fermion lines and gluons can be decomposed in terms of products of colour flow lines. The different elements are labelled by the permutation group acting on colour line indices.
For gluon amplitudes this leads to an over-counting of independent colour states, as some singlet contributions are actually forbidden. The N-gluon amplitude has independent colour states. Still, asymptotically we notice a factorial growth of terms for large . A way to fight the exponential growth both in helicity and colour components may be to resort to Monte Carlo sampling over these degrees of freedom, see for example .
3.3 The unitarity based approach
The idea to use the analytic structure of scattering amplitudes to determine their explicit form is very old [94, 95]. The well known non-linear relation between the transfer matrix, , implied by the unitarity of the S-matrix, ,
leads to a relation between the imaginary part of one-loop amplitudes and sewed tree-level diagrams, schematically:
Here the right hand side stands for the sum over all two-particle cuts. The latter correspond to squared tree-level amplitudes integrated over the respective two-particle phase space. Using the standard Feynman diagrammatic approach it can be shown that any one-loop amplitude with massless internal particles can be decomposed in terms of known scalar integrals with two to four external legs333If masses are present one also has to consider one-point integrals. , and a rational part, , which does not contribute to the imaginary part of the amplitude. The coefficients of the integrals, , and the rational term are rational polynomials in terms of spinor products, see Eq. (3.2), and/or Mandelstam variables and masses. The imaginary parts of the different scalar integrals can be uniquely attributed to a given integral. Subsequently unitarity implies that the knowledge of the imaginary part of the amplitude defines the coefficients of the scalar integrals.
Note that there are many different scalar two-, three- and four-point functions present in a given process. This is not indicated in this schematic equation. We see that tree amplitudes and phase space integrals in D=4 dimensions already fix a large part of the amplitude without ever being forced to evaluate any one-loop diagram. By turning the argument around we learn that one-loop diagrams are to a great extent determined by the pole part of the internal propagators which are nothing but cut propagators. The rational part is actually stemming from the UV behaviour of the one-loop amplitude, it comes from terms . In amplitudes with an improved UV behaviour, such as susy amplitudes, these terms are absent, and this explains why these amplitudes are fully determined by the D=4 cut structure. This fact has seen many applications in the last years, see e.g. . As the rational part is not determined by the D=4 cut structure, one has to resort to more sophisticated methods, like for example on-shell recursion relations  or D-dimensional unitarity .
This method has been developed and applied to a great extent by Bern, Dixon, Dunbar and Kosower in the 90s [99, 100], for a recent review see . It has seen a renaissance after the Santa Barbara workshop in 2003 where Witten pointed out a relation between certain string theories and QCD amplitudes using twistor methods . The revived interest lead to new efficient recursive evaluation procedures for multi-parton amplitudes [102, 103, 104]. A main feature of recent variants of the method is that multiple cuts are used to determine integral coefficients, which goes under the name of generalised unitarity [102, 103].
Cutting four lines in an N-point topology amounts to putting the corresponding four propagators on-shell. This procedure fixes the associated loop momentum completely and the coefficient of the related box diagram, see Fig. 2(a), is given as a product of tree diagrams.
One has to sum over the two complex solutions of quadratic on-shell conditions. This fixes the loop momentum. An important point is that the tree level amplitudes already incorporate gauge invariance manifestly. In Feynman diagram computations many graphs have to be combined to result in a gauge invariant expression. The power of this method has been demonstrated by the analytic evaluation of the cut-constructable part of the six-gluon amplitude . Combined efforts of many groups were needed to compute the different pieces of this amplitude, see [105, 106] for a collection of all formulae. The rational part has been first provided by a Feynman diagrammatic computation [107, 108, 109]. Note that the six-gluon one-loop amplitude is part of a full NLO computation of four jet production at the LHC and thus of phenomenological interest. There exists also a very compact result of the six-photon amplitude which was obtained using cutting methods [110, 111]. As has been shown in , the rational part is zero in this case.
The unitarity method, using multiple cuts, has recently been implemented in a numerical code by the Blackhat collaboration [113, 114] and by Ellis et. al.  (discussed below). To this moment the Blackhat collaboration has provided two remarkable applications, the colour ordered part of the (N=6,7,8)-gluon amplitudes in a special helicity configuration, and the leading colour contribution of (this includes crossing related amplitudes) [116, 87].
Important information on the numerical behaviour of the method could be gained in . The large number of algebraic operations in these complex amplitudes typically lead to round up errors in certain exceptional phase space regions. Any automated method has to come up with a reliable fail-safe procedure if the numerical cancellations are spoilt by the finite number of computed digits. In Fig. 3 a comparison of the MHV 6,7,8-gluon amplitudes between the numerical implementation and the known analytical result is shown.
The logarithmic relative error, , is plotted for the double, single-pole and the finite part of these amplitudes for 100.000 phase space points in Fig. 3. While the bulk of all points is evaluated reliably using standard double precision, it can be observed that the distribution has a tail which for a few percent of phase space points indicates a precision loss of almost all digits. Such a situation typically leads to an unstable behaviour of adaptive Monte Carlo event generators, which tend to sample points exactly in regions which induce a large variance of the result. In the given plot one observes that the tail is cut at a precision of around to . Beyond that point multi-precision libraries  are used to avoid a dangerous loss of precision.
The evaluation of the loop amplitude to is the necessary ingredient for a prediction of the important Standard Model background . Such multi-jet plus lepton plus missing energy signals occur for example in supersymmetric extensions of the Standard Model. The Blackhat collaboration has confronted their leading colour prediction for  with Tevatron data, see Fig. 4.
As expected the stability of the prediction under scale variations is largely improved. The leading colour approximation is expected to be phenomenologically justified within a few percent for Tevatron kinematics. Technically this approximation simplifies the calculation considerably, as many topologies are suppressed and can be neglected in this way.
In 2006 Ossola, Papadopoulos and Pittau (OPP) have proposed another unitarity approach which was based on the question how to reconstruct one-loop amplitudes from the integrand of Feynman integrals . Let us consider a sub-amplitude where all external particles are ordered. The full amplitude can be obtained by adding up different permutations. The leg-ordered amplitude can be obtained after summing the corresponding Feynman diagrams written such that the sum is under the integral and has a common denominator of propagators with .
The fact that general one-loop amplitudes can be written in terms of at most 4-point scalar integrals implies that the numerator can be written schematically in the following way
In this expression the coefficients , , , correspond to the actual amplitude coefficients of the various box, triangle and bubble integrals defined above, see Eq. (5). While these coefficients are independent of the loop momentum, the numerator function also contains spurious k-dependent pieces, , , , which have to vanish upon integration over the loop momentum . OPP asked the question how these spurious coefficients can be defined without doing the actual integration. They concluded that one simply has to evaluate the on-shell conditions of the various propagators which is in one-to-one relation to the unitarity cut methods using multiple cuts. This approach defines the coefficients of the terms ,, , as polynomials in . The coefficients of the box-, triangle-, bubble terms are defined by respectively 2,7,9 -polynomial coefficients. OPP proposed to extract these coefficients by numerical interpolation of the various polynomials in . They also define a procedure to access the rational part of the amplitudes [119, 120]. This method works also for individual one-loop Feynman diagrams and it has been implemented in a public computer code, CutTools . The method has been applied successfully to the evaluation of six-photon amplitudes with massless and massive internal fermions . The method has also been applied to the QCD corrections of triple vector boson production at the LHC. Here the amplitude has been integrated to obtain the total and differential cross sections for relevant observables .
A variant of these numerical unitarity based methods has also been investigated in  and has been extended to a D-dimensional approach [122, 123]. In the latter case also 5-point functions have to be included in the function basis. The key observation is that by adding higher dimensional information the rational part of the amplitude is also determined. A D-dimensional numerator function of a leg ordered amplitude is of the form
This approach guarantees that the evaluation of the rational part is now on the same footing as the cut constructable part. The method relies on the evaluation of amplitudes in dimensions different from four. It can be viewed as a D-dimensional generalisation of the OPP approach. The authors observed that the knowledge of tree amplitudes evaluated for complex momenta which are defined through the on-shell conditions of propagators, or equivalently through multiple cuts, defines the complexity of the algorithm modulo the number of different possibilities to cut a given N-point topology. The latter scales like , as any set of five propagators have to be put on-shell in a given diagram. The recursive evaluation of tree-amplitudes scales approximately like . Thus it can be expected that a colour ordered multi-parton amplitude can be evaluated by such a polynomial complexity algorithm, or short -algorithm, [122, 123]. As was pointed out above, for full amplitudes one needs to sum over different helicities and colour structures which unavoidably turns on an exponential growth and thus one ends up with an exponential or -algorithm for full amplitudes. Again Monte-Carlo sampling over colour and helicities might help to single out the relevant contributions numerically in a more efficient way. If a polynomial algorithm can be realised in this way remains to be shown.
The potential of this variant of the unitarity method has been demonstrated by evaluating single phase space points of different helicity amplitudes of the colour ordered -gluon amplitude for up to legs . Sampling over a large number of phase space points shows a good numerical performance concerning speed.
In Fig. 5 the polynomial behaviour of tree and one-loop amplitudes is shown as a function of .
The polynomial behaviour of the tree, , and the one-loop amplitudes, , is confirmed for the multi-gluon amplitudes. At the moment it seems not feasible to produce such a plot based on Feynman diagrams ever, as gauge cancellations and colour decomposition are not manifest and one needs the evaluation of the order of form factors. An independent implementation of the same method and amplitude, confirming the result, has been presented recently [125, 126]. Other applications of the method are the evaluation of the amplitudes  and the leading colour contributions to , [86, 127]. The latter are also implemented in the Rocket code  which handles numerically problematic phase space points by using higher precision libraries. Again, ultimately multi-precision libraries are in use to guarantee numerical reliability .
From the fast progress made in the last two years it can be expected that further cross section calculations for LHC processes with four or more partons/particles in the final state will become available using the discussed unitarity based methods. Note that apart from the multiple cut methods described so far, also single cut techniques are under development .
3.4 The Feynman diagrammatic approach
The traditional method of performing amplitude calculations is by starting from Feynman diagrams. It is based on representing a one-loop amplitude as a sum of diagrams sorted by colour structures
The kinematical functions can be represented by tensor one-loop integrals. The latter can be evaluated recursively in momentum space using the Passarino-Veltman method [129, 130]. Alternatively one can map the momentum integrals to Feynman parameter integrals which also can be reduced recursively [131, 133, 134, 136, 137].
Over the years quite a few groups have gathered a lot of experience in efficiently implementing Feynman diagram calculations. Cross section calculations, with up to five external particles, can generally be mastered now, as is demonstrated by the long list of accomplished tasks discussed above. For processes only a view NLO computations have been accomplished using Feynman diagrams, , , and . Processes of this complexity relevant for LHC phenomenology are under construction right now. Progress on the computation of the important LHC process was reported in 2008 , this group presented the full cross section computation of the quark induced subprocess at NLO in . The amplitude can be written in terms of rank three 6-point tensor integrals
In this calculation the algorithm described in  was used for the tensor reduction of 5- and 6-point integrals, otherwise Passarino-Veltman reduction was used.
Note that the gluon induced subprocess has a more complicated tensor structure, nonetheless a full result for this important Standard Model background can be expected in the near future.
Another collaboration is currently pressing for an automated, Feynman diagrammatic evaluation of multi-leg one-loop amplitudes. The GOLEM collaboration, where the acronym stands for General One-Loop Evaluator of Matrix elements, set up an automated reduction framework which takes special care of numerical issues by providing alternatives for the evaluation of one-loop integrals in critical phase space regions [135, 147, 141]. In this approach rank R N-point functions are reduced first to 6-point tensor integrals and the latter are then expressed by Feynman parameter integrals like
Spurious numerical problems due to so called inverse Gram determinants are avoided by having the option to evaluate the form factors in two different ways, either by numerical means or by reducing the form factors numerically to a scalar integral basis. The method is designed to allow for an efficient isolation of IR and UV divergences from the result and as such very well adapted for QCD calculations. The method can be also used to evaluate the rational part of an amplitude only . All 6-point form factors are coded in a fortran 95 code golem95, which is the first public library for such form factors . The code golem95 relies on the evaluation of scalar integrals. In the present version only integrals with zero internal masses are implemented. Another implementation of 6-point form factors has been presented recently [148, 149].
A full amplitude evaluation relies on diagrammatic input. The GOLEM collaboration uses public tools like QGRAF  and FeynArts  for this step. The colour algebra and the helicity management is done as outlined above. Subsequently two independent strategies are followed. Firstly, the diagrammatic input is converted to a form factor representation
exported to a fortran code and then linked to the form factor library golem95. Secondly a fully symbolic representation of the scalar integral coefficients is generated using FORM  and MAPLE
Using these implementations several computations of relevance for the LHC have been performed [152, 52, 63, 81]. The viability of the approach for processes was shown by the evaluation of the 6-photon amplitude mediated by a massless electron loop . This amplitude has been used as testing ground for various methods and perfect agreement has been obtained in all cases [153, 111, 154, 155].
Currently the NLO corrections for the process are under construction using the developed computational tools. Because the 4b-quark background needs to be known as precisely as possible in the context of two Higgs doublet models, this process was included in the Les Houches “experimentalist wish-list” . The amplitude consists of two different initial states and . The first can be represented by about 250 Feynman diagrams among which 25 pentagon and 8 hexagon diagrams can be found, see Fig. 7 for a selection of LO and NLO diagrams. Note that only the pentagon and hexagon diagrams are computationally challenging. The other topologies are relatively simple and do not present a problem apart from book keeping which is not an issue in automated approaches. The total number of diagrams is thus not a good measure of the complexity of a calculation.
All one-loop Feynman diagrams have been evaluated in two completely independent ways, by symbolic reduction to scalar integrals and by using the form factor decomposition in combination with the golem95 library as discussed above. The virtual correction is a interference term between the LO tree and NLO loop amplitude. After UV renormalisation, IR divergences in form of single and double poles remain. The IR structure of one-loop amplitudes is well-known and a finite function can be obtained after adding an adequate IR compensation term, as defined for example by the Catani-Seymour insertion operator . In this way one can define a finite local K-factor function
To investigate the numerical performance of the approach, 200.000 random phase space points have been evaluated using double and quadruple precision, see Fig. 8. In the figure the numerical precision of the evaluation of the local K-factor and the cancellation of the IR single/double pole are compared.
The result indicates that the size of the finite K-factor is a good indication of the numerical accuracy of the evaluation. Large K-factors are mostly induced through numerical round-up errors. The cancellation of single and double poles, which is used by other collaborations, can also be applied to decide about the quality of the evaluation but seems to be less indicative. In Fig. 9 the improvement of the scale dependence is shown, if the subtracted virtual correction term, as defined in Eq. (11), is included. The cross sections have been evaluated using the experimental cuts , , GeV and the NLO pdf set CTEQ6.5. A further improvement of the behaviour is to be expected, because the logarithmic factorisation scale dependence is not yet compensated.
To obtain the full NLO result one has to add the dipole subtracted NLO real emission corrections. For the full process, , the gluon induced sub-process also has to be included. These parts of the calculation are in progress.
Note that the evaluation strategy of the virtual NLO corrections differs from the conventional approach of integrating a NLO cross section directly. In the given example the LO cross section has been evaluated first to produce an unweighted event sample. The subtracted virtual correction term defines a finite K-factor function which is simply used to reweight these LO events. In this way interference effects between adaptive Monte-Carlo integration and numerical round-up errors are completely avoided. For a more detailed discussion see . A similar strategy has also been applied in [67, 72].
The goal of the GOLEM collaboration is to provide an automated framework for one-loop amplitude computations. The GOLEM implementation is sketched in Fig. 10.
The input files which contain process and model dependent information are processed with a python script which controls QGRAF and FORM codes. The diagrammatic output, including UV and IR subtraction, is defined in terms of form factors and written to fortran 95 code. The latter is linked to the form factor library golem95. At the moment the code can be compiled in double or quadruple precision but in the future a more dynamical solution will be adopted. The resulting computer program provides a function to reweight LO event samples stemming from tree level event generators. For the shown results Whizard  has been used.
3.5 Is there a preferred method?
We have seen that both the unitarity based and also the Feynman diagrammatic method are sufficiently well developed meanwhile to allow for cross section calculations with more than five external legs. While the unitarity based methods were initially used to obtain compact analytical expressions for helicity amplitudes, the quest for automation lead recently to various numerical implementations. Higher precision libraries play a prominent role in the related computer programs which shows that any given method seems to be doomed to deal with exceptional phase space regions in one way or the other. Higher precision libraries seem to be a natural cure for numerical problems. As the number of problematic phase space points is typically only a few percent of all points, increased evaluation times are not a severe issue444As discussed during this conference a hardware implementation of quadruple precision in processors is not yet in sight although technical standards have been already defined for such a step. The particle physics community would surely be grateful for such a development. Unfortunately fundamental science issues do not play the role they deserve in business plans of computer chip producing companies, although most technological developments are essentially based on progress in fundamental science!.
Different methods will soon provide cross section predictions for partonic processes relevant for LHC phenomenology. As was explained above the scaling properties of one-loop N-point amplitude computations seem to indicate that the unitarity based methods are clearly preferable. The factorial growth in diagrammatic computations has to be compared to polynomial algorithms now. For example, the evaluation of a colour ordered multi-gluon helicity amplitudes scales like whereas a Feynman diagrammatic calculation would involve about leg ordered diagrams with roughly form factor evaluations each. However, let us ask the question, whether the asymptotic behaviour gives a guideline which method is preferable for phenomenological applications. The logarithmic ratio of the asymptotic behaviour, is plotted in Fig. 11.
The figure indicates that for phenomenological relevant multiplicities, say , the asymptotic scaling behaviour is not a good measure for the question which method is preferable. Furthermore the true complexity of a method should be measured with full processes, taking into account all helicity amplitudes, relevant colour structures and sub-processes. It is important to note in this respect that especially for background cross sections severe experimental cuts will be applied. Experience shows that any approximation or assumption might be invalidated in that way. For the full task no -algorithm exists.
For , implementation issues, like efficient caching and the re-use of recursive information can still be improved in Feynman diagram calculations. Together with the increasing computer power555We learnt at this conference that the “multi-core” computer era will change to the “many-core” era were the number of idle and available processors will increase enormously. it can be expected that both methods will deliver what they promise, namely cross section predictions for LHC phenomenology at the next-to-leading order level.
4 Conclusion and Outlook
In this talk I have argued that the LHC needs and deserves an effort to predict prominent signal and background processes at the next-to-leading order level in QCD. Especially multi-particle amplitudes have a high power of which induce large scale variations. Absolute rates thus can not be predicted reliably with leading order Monte-Carlo tools and eventually this will hamper the understanding of LHC data and the discovery of new physics. Many relevant Standard Model processes with three final state particles are meanwhile available in the literature beyond the leading order. Unfortunately the results are not always available as a public code which would be most beneficial for the experimental community. The theoretical community is at the moment enormously active to go beyond the given level of complexity and NLO predictions for processes like , , , , , , which are highly relevant for various Higgs and BSM search channels, are under construction. Many different methods have been developed and are applied for the virtual corrections. Unitarity based methods look very promising in this respect but for the required level of complexity both, the Feynman diagram and the unitarity based approaches, will provide valuable results for phenomenology. The avenue to fully numerical approaches to one-loop amplitudes is still not fully explored, but given the ever growing computer power this might be a promising direction avoiding any kind of algebraic manipulations of initial expressions.
Most groups move at the moment towards automated approaches which will provide platforms to do many computations with the same set-up. To cure the problem of numerical instabilities multi-precision libraries seem to be the accepted way out now. It would be highly appreciable, if the community could agree on a standardisation of in- and output formats for NLO codes and would make computer programs publicly available. The structure of one-loop computations is indeed very modular and a minimal agreement on colour and helicity management and on passing IR subtraction terms, which are basically universal anyhow, could result in transportable modules for virtual corrections. This would allow to use computations of different groups interchangeably by combining them with public tree level matrix element generators which should of course also contain the functionality of IR subtractions. In a next step the combination with parton showers could be obtained. Here the inclusion of NLO precision is also well understood, but again a focus on modularity and transportability would help to use synergies between different groups.
The main conclusion of this talk is that the developments of the last few years are spectacular and point towards Monte Carlo tools at full next-to-leading order QCD level, at least as long the number of final state particles is not excessively high. The dominant and most relevant processes for Higgs and new physisc searches will certainly be available in form of flexible and reliable public computer programs and the analysis of LHC data beyond the leading order in will eventually become the standard.
I would like to thank the organisers of the ACAT 2008 workshop for the opportunity to present my view on the progress of the field during the fantastic meeting in Erice, Italy. I also would like to thank my colleagues of the GOLEM collaboration, A. Guffanti, J. Ph. Guillet, G. Heinrich, S. Karg, N. Kauer, T. Reiter, J. Reuter and G. Sanguinetti for many stimulating discussions and their cooperation, which helped to shape this contribution in a substantial way.
-  N. E. Adam et al., arXiv:0803.1154 [hep-ph].
-  S. P. Martin, arXiv:hep-ph/9709356.
-  A. Djouadi, Phys. Rept. 459 (2008) 1 [arXiv:hep-ph/0503173].
-  B. C. Allanach et al., arXiv:hep-ph/0602198.
-  R. Sundrum, (TASI 2004), arXiv:hep-th/0508134.
-  C. Csaki, J. Hubisz and P. Meade, (TASI 2005) arXiv:hep-ph/0510275.
-  G. D. Kribs, (TASI 2006) arXiv:hep-ph/0605325.
-  M. Schmaltz and D. Tucker-Smith, Ann. Rev. Nucl. Part. Sci. 55 (2005) 229 [arXiv:hep-ph/0502182].
-  K. Cranmer, B. Mellado, W. Quayle, S.L. Wu, “Statistical Methods to Assess the Combined Sensitivity of the ATLAS Detector to the Higgs Boson in the Standard Model”, ATL-PHYS-2005-034.
-  CMS Collaboration, “CMS Physics Technical Design Report Volume II: Physics Performance”, CERN/LHCC/2006-021, CMS TDR 8.2.
-  K. Jakobs [ATLAS Collaboration], Int. J. Mod. Phys. A 23, 5093 (2008).
-  V. Buescher and K. Jakobs, Int. J. Mod. Phys. A 20, 2523 (2005) [arXiv:hep-ph/0504099].
-  G. Aad et al. [The ATLAS Collaboration], arXiv:0901.0512.
-  T. Binoth, J. P. Guillet, E. Pilon and M. Werlen, Eur. Phys. J. C 16 (2000) 311 [arXiv:hep-ph/9911340].
-  Z. Bern, L. J. Dixon and C. Schmidt, Phys. Rev. D 66 (2002) 074018 [arXiv:hep-ph/0206194].
-  N. Kauer, T. Plehn, D. L. Rainwater and D. Zeppenfeld, Phys. Lett. B 503, 113 (2001) [arXiv:hep-ph/0012351].
-  M. L. Mangano, M. Moretti, F. Piccinini, R. Pittau and A. D. Polosa, JHEP 0307 (2003) 001 [arXiv:hep-ph/0206293].
-  A. Pukhov et al., arXiv:hep-ph/9908288.
-  F. Krauss, R. Kuhn and G. Soff, JHEP 0202 (2002) 044 [arXiv:hep-ph/0109036].
-  T. Gleisberg and S. Hoche, JHEP 0812 (2008) 039 [arXiv:0808.3674 [hep-ph]].
-  T. Hahn and M. Perez-Victoria, Comput. Phys. Commun. 118, 153 (1999) [arXiv:hep-ph/9807565].
-  A. Kanaki and C. G. Papadopoulos, Comput. Phys. Commun. 132, 306 (2000) [arXiv:hep-ph/0002082].
-  C. G. Papadopoulos and M. Worek, arXiv:hep-ph/0606320.
-  A. Cafarella, C. G. Papadopoulos and M. Worek, arXiv:0710.2427 [hep-ph].
-  F. Maltoni and T. Stelzer, JHEP 0302, 027 (2003) [arXiv:hep-ph/0208156].
-  T. Stelzer and W. F. Long, Comput. Phys. Commun. 81, 357 (1994) [arXiv:hep-ph/9401258].
-  W. Kilian, T. Ohl and J. Reuter, arXiv:0708.4233 [hep-ph].
-  T. Sjostrand, S. Mrenna and P. Skands, JHEP 0605, 026 (2006) [arXiv:hep-ph/0603175].
-  M. Bahr et al., Eur. Phys. J. C 58, 639 (2008) [arXiv:0803.0883 [hep-ph]].
-  T. Gleisberg, S. Hoche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert and J. Winter, JHEP 0902, 007 (2009) [arXiv:0811.4622 [hep-ph]].
-  S. Catani and M. H. Seymour, Nucl. Phys. B 485, 291 (1997) [Erratum-ibid. B 510, 503 (1998)] [arXiv:hep-ph/9605323].
-  S. Catani, S. Dittmaier, M. H. Seymour and Z. Trocsanyi, Nucl. Phys. B 627, 189 (2002) [arXiv:hep-ph/0201036].
-  M. H. Seymour and C. Tevlin, arXiv:0803.2231 [hep-ph].
-  T. Gleisberg and F. Krauss, Eur. Phys. J. C 53, 501 (2008) [arXiv:0709.2881 [hep-ph]].
-  K. Hasegawa, S. Moch and P. Uwer, Nucl. Phys. Proc. Suppl. 183, 268 (2008) [arXiv:0807.3701 [hep-ph]].
-  R. Frederix, T. Gehrmann and N. Greiner, JHEP 0809, 122 (2008) [arXiv:0808.2128 [hep-ph]].
-  J. M. Campbell and R. K. Ellis, Phys. Rev. D 62 (2000) 114012 [arXiv:hep-ph/0006304].
-  T. Binoth, M. Ciccolini, N. Kauer and M. Kramer, JHEP 0503, 065 (2005) [arXiv:hep-ph/0503094].
-  T. Binoth, M. Ciccolini, N. Kauer and M. Kramer, JHEP 0612, 046 (2006) [arXiv:hep-ph/0611170].
-  T. Binoth, N. Kauer and P. Mertsch, arXiv:0807.0024 [hep-ph].
-  Z. Nagy, Phys. Rev. D 68, 094002 (2003) [arXiv:hep-ph/0307268].
-  D. de Florian and Z. Kunszt, Phys. Lett. B 460, 184 (1999) [arXiv:hep-ph/9905283].
-  V. Del Duca, F. Maltoni, Z. Nagy and Z. Trocsanyi, JHEP 0304, 059 (2003) [arXiv:hep-ph/0303012].
-  J. M. Campbell, R. K. Ellis and G. Zanderighi, JHEP 0610, 028 (2006) [arXiv:hep-ph/0608194].
-  T. Han and S. Willenbrock, Phys. Lett. B 273 (1991) 167.
-  T. Figy, C. Oleari and D. Zeppenfeld, Phys. Rev. D 68 (2003) 073005 [arXiv:hep-ph/0306109].
-  W. Beenakker, S. Dittmaier, M. Kramer, B. Plumper, M. Spira and P. M. Zerwas, Nucl. Phys. B 653 (2003) 151 [arXiv:hep-ph/0211352].
-  S. Dawson, C. B. Jackson, L. H. Orr, L. Reina and D. Wackeroth, Eur. Phys. J. C 33, S451 (2004) [arXiv:hep-ph/0311216].
-  T. Plehn and M. Rauch, Phys. Rev. D 72 (2005) 053008 [arXiv:hep-ph/0507321].
-  T. Binoth, S. Karg, N. Kauer and R. Ruckl, Phys. Rev. D 74 (2006) 113008 [arXiv:hep-ph/0608057].
-  B. Jager, C. Oleari and D. Zeppenfeld, JHEP 0607 (2006) 015 [arXiv:hep-ph/0603177].
-  B. Jager, C. Oleari and D. Zeppenfeld, Phys. Rev. D 73 (2006) 113006 [arXiv:hep-ph/0604200].
-  G. Bozzi, B. Jager, C. Oleari and D. Zeppenfeld, Phys. Rev. D 75 (2007) 073004 [arXiv:hep-ph/0701105].
-  A. Lazopoulos, K. Melnikov and F. Petriello, Phys. Rev. D 76 (2007) 014001 [arXiv:hep-ph/0703273].
-  S. Dittmaier, P. Uwer and S. Weinzierl, Phys. Rev. Lett. 98 (2007) 262002 [arXiv:hep-ph/0703120].
-  S. Dittmaier, P. Uwer and S. Weinzierl, arXiv:0810.0452 [hep-ph].
-  S. Dittmaier, S. Kallweit and P. Uwer, arXiv:0710.1577 [hep-ph].
-  J. M. Campbell, R. Keith Ellis and G. Zanderighi, JHEP 0712 (2007) 056 [arXiv:0710.1832 [hep-ph]].
-  T. Binoth, G. Ossola, C. G. Papadopoulos and R. Pittau, JHEP 0806, 082 (2008) [arXiv:0804.0350 [hep-ph]].
-  V. Hankele and D. Zeppenfeld, Phys. Lett. B 661 (2008) 103 [arXiv:0712.3544 [hep-ph]].
-  A. Lazopoulos, K. Melnikov and F. J. Petriello, Phys. Rev. D 77, 034021 (2008) [arXiv:0709.4044 [hep-ph]].
-  A. Lazopoulos, T. McElmurry, K. Melnikov and F. Petriello, Phys. Lett. B 666 (2008) 62 [arXiv:0804.2220 [hep-ph]].
-  F. Febres Cordero, L. Reina and D. Wackeroth, Phys. Rev. D 78, 074014 (2008) [arXiv:0806.0808 [hep-ph]].
-  M. Ciccolini, A. Denner and S. Dittmaier, arXiv:0710.4749 [hep-ph].
-  M. Ciccolini, A. Denner and S. Dittmaier, Phys. Rev. Lett. 99 (2007) 161803 [arXiv:0707.0381 [hep-ph]].
-  M. M. Weber, Nucl. Phys. Proc. Suppl. 160 (2006) 200.
-  J. R. Andersen and J. M. Smillie, Phys. Rev. D 75, 037301 (2007) [arXiv:hep-ph/0611281].
-  J. R. Andersen, T. Binoth, G. Heinrich and J. M. Smillie, arXiv:0709.3513 [hep-ph].
-  A. Bredenstein, K. Hagiwara and B. Jager, Phys. Rev. D 77 (2008) 073004 [arXiv:0801.4231 [hep-ph]].
-  Z. Bern et al. [NLO Multileg Working Group], arXiv:0803.0494 [hep-ph].
-  A. Bredenstein, A. Denner, S. Dittmaier and S. Pozzorini, JHEP 0808 (2008) 108 [arXiv:0807.1248 [hep-ph]].
-  T. Binoth et al., arXiv:0807.0605 [hep-ph].
-  R. K. Ellis, W. T. Giele, Z. Kunszt, K. Melnikov and G. Zanderighi, JHEP 0901 (2009) 012 [arXiv:0810.2762 [hep-ph]].
-  C. F. Berger et al., arXiv:0902.2760 [hep-ph].
-  T. Reiter, arXiv:0903.0947 [hep-ph].
-  L. J. Dixon, (TASI 1995), arXiv:hep-ph/9601359.
-  D. Maitre and P. Mastrolia, arXiv:0710.5559 [hep-ph].
-  A. Kanaki and C. G. Papadopoulos, arXiv:hep-ph/0012004.
-  F. Maltoni, K. Paul, T. Stelzer and S. Willenbrock, Phys. Rev. D 67 (2003) 014026 [arXiv:hep-ph/0209271].
-  C. G. Papadopoulos and M. Worek, Eur. Phys. J. C 50, 843 (2007) [arXiv:hep-ph/0512150].
-  R. E. Cutkosky, J. Math. Phys. 1 (1960) 429.
-  R. J. Eden, P. V. Landshoff, D. I. Olive, J. C. Polkinghorne, (Cambridge University Press, 1966).
-  Z. Bern, L. J. Dixon and D. A. Kosower, Annals Phys. 322 (2007) 1587 [arXiv:0704.2798 [hep-ph]].
-  C. F. Berger, Z. Bern, L. J. Dixon, D. Forde and D. A. Kosower, Phys. Rev. D 74 (2006) 036009 [arXiv:hep-ph/0604195].
-  C. Anastasiou, R. Britto, B. Feng, Z. Kunszt and P. Mastrolia, Phys. Lett. B 645 (2007) 213 [arXiv:hep-ph/0609191].
-  Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, Nucl. Phys. B 435, 59 (1995) [arXiv:hep-ph/9409265].
-  Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, Nucl. Phys. B 425, 217 (1994) [arXiv:hep-ph/9403226].
-  E. Witten, Commun. Math. Phys. 252 (2004) 189 [arXiv:hep-th/0312171].
-  R. Britto, F. Cachazo and B. Feng, Nucl. Phys. B 725, 275 (2005) [arXiv:hep-th/0412103].
-  A. Brandhuber, S. McNamara, B. J. Spence and G. Travaglini, JHEP 0510, 011 (2005) [arXiv:hep-th/0506068].
-  R. Britto, B. Feng and P. Mastrolia, Phys. Rev. D 73, 105004 (2006) [arXiv:hep-ph/0602178].
-  D. C. Dunbar, Nucl. Phys. Proc. Suppl. 183, 122 (2008).
-  D. C. Dunbar, arXiv:0901.1202 [hep-ph].
-  Z. Xiao, G. Yang and C. J. Zhu, Nucl. Phys. B 758, 53 (2006) [arXiv:hep-ph/0607017].
-  X. Su, Z. Xiao, G. Yang and C. J. Zhu, Nucl. Phys. B 758, 35 (2006) [arXiv:hep-ph/0607016].
-  Z. Xiao, G. Yang and C. J. Zhu, Nucl. Phys. B 758, 1 (2006) [arXiv:hep-ph/0607015].
-  T. Binoth, G. Heinrich, T. Gehrmann and P. Mastrolia, Phys. Lett. B 649 (2007) 422 [arXiv:hep-ph/0703311].
-  G. Ossola, C. G. Papadopoulos and R. Pittau, JHEP 0707 (2007) 085 [arXiv:0704.1271 [hep-ph]].
-  T. Binoth, J. P. Guillet and G. Heinrich, JHEP 0702 (2007) 013 [arXiv:hep-ph/0609054].
-  C. F. Berger et al., Phys. Rev. D 78 (2008) 036003 [arXiv:0803.4180 [hep-ph]].
-  C. F. Berger et al., Nucl. Phys. Proc. Suppl. 183 (2008) 313 [arXiv:0807.3705 [hep-ph]].
-  R. K. Ellis, W. T. Giele and Z. Kunszt, JHEP 0803, 003 (2008) [arXiv:0708.2398 [hep-ph]].
-  C. F. Berger et al., arXiv:0808.0941 [hep-ph].
D. H. Bailey, these proceedings.
-  G. Ossola, C. G. Papadopoulos and R. Pittau, Nucl. Phys. B 763 (2007) 147 [arXiv:hep-ph/0609007].
-  G. Ossola, C. G. Papadopoulos and R. Pittau, JHEP 0805, 004 (2008) [arXiv:0802.1876 [hep-ph]].
-  P. Draggiotis, M. V. Garzelli, C. G. Papadopoulos and R. Pittau, arXiv:0903.0356 [hep-ph].
-  G. Ossola, C. G. Papadopoulos and R. Pittau, JHEP 0803, 042 (2008) [arXiv:0711.3596 [hep-ph]].
-  R. K. Ellis, W. T. Giele and Z. Kunszt, PoS RADCOR2007 (2007) 020 [arXiv:0802.4227 [hep-ph]].
-  R. K. Ellis, W. T. Giele, Z. Kunszt and K. Melnikov, arXiv:0806.3467 [hep-ph].
-  W. T. Giele and G. Zanderighi, arXiv:0805.2152 [hep-ph].
-  A. Lazopoulos, arXiv:0812.2998 [hep-ph].
-  J. C. Winter and W. T. Giele, arXiv:0902.0094 [hep-ph].
-  R. K. Ellis, K. Melnikov and G. Zanderighi, arXiv:0901.4101 [hep-ph].
-  S. Catani, T. Gleisberg, F. Krauss, G. Rodrigo and J. C. Winter, JHEP 0809 (2008) 065 [arXiv:0804.3170 [hep-ph]].
-  G. Passarino and M. J. G. Veltman, Nucl. Phys. B 160, 151 (1979).
-  A. Denner, Fortsch. Phys. 41 (1993) 307 [arXiv:0709.1075 [hep-ph]].
-  A. I. Davydychev, Phys. Lett. B 263 (1991) 107.
-  Z. Bern, L. J. Dixon and D. A. Kosower, Phys. Lett. B 302 (1993) 299 [Erratum-ibid. B 318 (1993) 649] [arXiv:hep-ph/9212308].
-  Z. Bern, L. J. Dixon and D. A. Kosower, Nucl. Phys. B 412 (1994) 751 [arXiv:hep-ph/9306240].
-  O. V. Tarasov, Phys. Rev. D 54 (1996) 6479 [arXiv:hep-th/9606018].
-  T. Binoth, J. P. Guillet and G. Heinrich, Nucl. Phys. B 572, 361 (2000) [arXiv:hep-ph/9911342].
-  G. Duplancic and B. Nizic, Eur. Phys. J. C 35 (2004) 105 [arXiv:hep-ph/0303184].
-  W. T. Giele and E. W. N. Glover, JHEP 0404 (2004) 029 [arXiv:hep-ph/0402152].
-  R. K. Ellis and G. Zanderighi, JHEP 0802 (2008) 002 [arXiv:0712.1851 [hep-ph]].
-  G. J. van Oldenborgh, Comput. Phys. Commun. 66, 1 (1991).
-  A. van Hameren, J. Vollinga and S. Weinzierl, Eur. Phys. J. C 41, 361 (2005) [arXiv:hep-ph/0502165].
-  T. Binoth, J. P. Guillet, G. Heinrich, E. Pilon and T. Reiter, arXiv:0810.0992 [hep-ph].
-  T. Hahn, Acta Phys. Polon. B 30, 3469 (1999) [arXiv:hep-ph/9910227].
-  A. Denner, S. Dittmaier, M. Roth and L. H. Wieders, Nucl. Phys. B 724, 247 (2005) [arXiv:hep-ph/0505042].
-  F. Boudjema et al., In the Proceedings of 2005 International Linear Collider Workshop (LCWS 2005), Stanford, California, 18-22 Mar 2005, pp 0601 [arXiv:hep-ph/0510184].
-  G. Lei, M. Wen-Gan, H. Liang, Z. Ren-You and J. Yi, arXiv:0708.2951 [hep-ph].
-  A. Denner and S. Dittmaier, Nucl. Phys. B 734 (2006) 62 [arXiv:hep-ph/0509141].
-  T. Binoth, J. P. Guillet, G. Heinrich, E. Pilon and C. Schubert, JHEP 0510 (2005) 015 [arXiv:hep-ph/0504267].
-  T. Diakonidis, J. Fleischer, J. Gluza, K. Kajda, T. Riemann and J. B. Tausk, arXiv:0812.2134 [hep-ph].
-  T. Diakonidis, arXiv:0901.4455 [hep-ph].
-  P. Nogueira, J. Comput. Phys. 105, 279 (1993).
-  J. A. M. Vermaseren, arXiv:math-ph/0010025.
-  T. Binoth, J. P. Guillet and F. Mahmoudi, JHEP 0402, 057 (2004) [arXiv:hep-ph/0312334].
-  Z. Nagy and D. E. Soper, Phys. Rev. D 74 (2006) 093006 [arXiv:hep-ph/0610028].
-  C. Bernicot and J. P. Guillet, JHEP 0801 (2008) 059 [arXiv:0711.4713 [hep-ph]].
-  W. Gong, Z. Nagy and D. E. Soper, Phys. Rev. D 79, 033005 (2009) [arXiv:0812.3686 [hep-ph]].
-  A. Ferroglia, M. Passera, G. Passarino and S. Uccirati, Nucl. Phys. B 650 (2003) 162 [arXiv:hep-ph/0209219].
-  T. Binoth, G. Heinrich and N. Kauer, Nucl. Phys. B 654 (2003) 277 [arXiv:hep-ph/0210023].
-  Y. Kurihara and T. Kaneko, Comput. Phys. Commun. 174 (2006) 530 [arXiv:hep-ph/0503003].
-  C. Anastasiou, S. Beerli and A. Daleo, JHEP 0705 (2007) 071 [arXiv:hep-ph/0703282].