Susy Theories and QCD: Numerical Approaches
Abstract
We review onshell and unitarity methods and discuss their application to precision predictions for LHC physics. Being universal and numerically robust, these methods are straightforward to automate for nexttoleadingorder computations within Standard Model and beyond. Several stateoftheart results including studies of /+3jet and +4jet production have explicitly demonstrated the effectiveness of the unitarity method for describing multiparton scattering. Here we review central ideas needed to obtain efficient numerical implementations. This includes onshell looplevel recursions, the unitarity method, color management and further refined tricks.
This article is an invited review for a special issue of Journal of Physics A devoted to “Scattering Amplitudes in Gauge Theories”.
pacs:
12.38.Bx, 13.85.t, 13.87.aUCLA11TEP120, NSFKITP11175
1 Introduction
The physics program at the Large Hadron Collider (LHC) relies heavily on the ever increasing theoretical control over modeling highenergy proton collisions. The detailed theoretical understanding not only increases the reach in new physics and particle searches, but also allows to study the fundamental dynamics and properties of particles. Formulating new observables for addressing specific physics questions is a typical task which relies on quantitative reliable theoretical input. In the long run, high statistics measurements combined with precision prediction from theory will allow systematic probes of fundamental particle theory at ever deeper levels.
The theory of Quantum Chromodynamics (QCD) concisely describes the collisions of protons, however, the dominant dynamics differ depending on observables and the regions of phasespace [1, 2, 3, 4, 5]. Here we have in mind proton collisions with large momentum transfer that are typical for the production of heavy particles. These include the Higgs, top, new particles within theories of supersymmetry as well as more conventional Standard Model processes at large scattering angles. To a good approximation such scattering processes factorize into the longdistance dynamics of quarks and gluons within protons, short distance hard interactions between these partons and, finally, the formation of hadrons and observable jets from the emerging partonic states and remaining proton fragments. Monte Carlo event generators deal with all aspects of the multilayered simulation of the proton collisions. For some purposes, i.e. sufficiently inclusive observables, it is accurate to use a simpler, purely partonic description of events. To this end one combines finalstate partons into observable jets, consistently ignoring corrections from showering and hadronization. Numerical methods are commonly used for the evaluation of differential cross sections, being well suited for comparisons to experimental data.
The hard scattering process is a central stage in the simulation of proton collisions. It is described through scattering amplitudes which are accessible through first principle computations in quantum field theory. The complexity of scattering process allows only a perturbative approach, based on the expansion in the strong coupling . A first step towards a theoretical understanding of QCD is the evaluation of cross sections at leading order (LO) in the strong coupling . Many tools [6, 7, 8, 9, 10] are available to generate predictions at leading order. Some of the methods applied incorporate highermultiplicity leadingorder matrix elements into partonshowering programs [11, 12], using matching (or merging) procedures [13, 14].
The truncation of the perturbative expansion introduces an explicit dependence on the unphysical renormalization scale leading to a theoretical uncertainty. QCD cross sections can have strong sensitivity on higherorder corrections, motivating the challenging quest for perturbative corrections. Nexttoleading (NLO) order predictions significantly reduce renormalization and factorizationscale dependence – a feature that becomes more significant with increasing jet multiplicity (see e.g. [15]). In addition, NLO corrections take into account further physics effects including initial state radiation, parton merging into jets and a more complete account of partonic production channels. Fixedorder results at NLO can also be matched to parton showers [16] with the prospect of complete event generation at nexttoleadingorder in the strong coupling. The value of a first principle understanding of scattering processes in addition to the increased quantitative control motivates the quest for cross sections at NLO and sometimes beyond this [17, 18].
In this contribution we will describe looplevel onshell [19, 20, 21] and unitarity methods [22, 23]. Our main focus will be on generalized unitarity [24, 25, 26] and its extensions for the numerical computation [27, 28, 21, 29] of hard scattering matrix elements. In addition, we will discuss numerical onshell looplevel recursions [21]. The central strategy of these approaches is to make maximal use of universal physical principles and mathematical structures in order to describe complex multiparticle processes including quantum corrections. Numerical algorithms based on unitarity and onshell methods are efficient, numerically stable and can be automated for generic scattering processes. The rapid recent developments in this field have already led to many new studies of complex scattering processes at the LHC [30, 31, 32, 33, 34, 35, 36, 37, 38]. Several further implementations of unitarity methods have been reported [39, 40, 41, 38, 42]. For related analytic approaches see the chapter in this review by Britto [43].
For processes at the LHC with complicated final states, computations can be very challenging. In principle, one can use the traditional Feynman diagram representation of amplitudes for numerical evaluation. However, even at leading order, the number of Feynman diagrams grows more than factorial as the number of finalstates increases. Computation times scale accordingly and refined approaches are needed.
At leading order recursive approaches allow to reduce the growth in complexity. Two central strategies are commonly used for numerical computations: Offshell recursion relations [44, 45, 46], based on DysonSchwinger equations, optimise the reuse of recurring groups of Feynman graphs. In contrast, onshell recursions [47, 48] take advantage of the remarkable simplicity of the physical scattering amplitudes (see e.g. [49]). The simplicity arises in part from symmetry properties of tree amplitudes [50, 51, 52] that are present in QCDlike theories; see the chapter in this review by Brandhuber, Spence and Travaglini [53]. For most practical purposes the efficiency of the two approaches is comparable, depending on the explicit realization of the algorithms.
At nexttoleadingorder additional challenges arise, in particular, for the virtual corrections due to the loopmomentum integration. NLO cross sections are built from several ingredients: virtual corrections, computed from the interference of treelevel and oneloop amplitudes; realemission corrections; and a mechanism for isolating and integrating the infrared singularities in the latter. Automated approaches [54, 55, 56] to deal with these issues are available based on subtraction methods [57, 58, 59]. Recursive methods are effective for computations of realemission corrections. Such methods, however, are not directly applicable to virtual corrections. Traditional methods evaluate the loop integrals of Feynman diagrams (see e.g. [60, 61, 62]), and have to overcome two central challenges: growth of the number of Feynman diagram expressions and the evaluation of tensorial loop integrals, while maintaining gauge invariance. Means to deal with tensor integral reductions [63, 64] as well as strategies to recycle substructures have been shown to reduce complexity inherent in Feynman diagram approaches [65, 61]. For a more complete discussion of important NLO computations along these lines see ref. [66].
The unitarity method [22], in contrast, constructs loop amplitudes from onshell tree amplitudes; gauge invariance is built in and maintained throughout computations. In addition, complexity arising from large numbers of Feynman diagrams is avoided by recursive methods for tree evaluations. Similarly, looplevel recursions [19, 20, 67, 21] construct amplitudes efficiently using purely onshell lowerpoint input. Effective numerically stable implementations of these onshell methods have been demonstrated by various groups [21, 68, 35, 40, 39, 41, 38, 42]. Beyond this, unitarity approaches have already been applied to stateoftheart NLO computations [36, 30, 31, 32, 69, 34, 37].
In more detail, numerical unitarity based approaches use a combination of methods. Scattering amplitudes are naturally split into two parts; a part with logarithmic dependence on kinematic invariants and a rational remainder. Typically, unitarity approaches in strictly four dimensions are used for the computation of the logarithmic parts, although, onshell recursions [20] may as well be applied in certain cases. At present, there are three main choices for computing the rational remainder within a processnonspecific numerical program: onshell recursion [21], dimensional unitarity [23, 29], and a refined Feynmandiagram approach [70, 38]. We will discuss here numerical unitarity approaches in four and dimensions following refs. [21, 29] as well as numerical looplevel onshell methods [21].
Several recent developments allow us to use a purely numerical approach at loop level. A key tool is generalized unitarity [24, 25] which imposes multiple unitarity cuts (more then two propagators) and gives a refined system of consistency relations that is easier to solve. In addition, generalized cuts allow for a hierarchical approach; computing coefficients of fourpoint integrals, in turn threepoint and, finally, twopoint etc. integrals from successively cutting four, three and two propagators. (For the related maximal cut technique for multiloop amplitudes see [71].) Analogous approaches may be applied in dimensions [23] (see also [72]), where higher point integrals have to be considered. Further simplifications arise from working directly with the loop integrand as opposed to the integrated loop amplitude. The unitarity method then turns into a purely algebraic approach. The starting point for this approach is a generic representation of the loop integrand [27], whose free, process dependent parameters are to be determined by unitarity relations. A particularly useful parametrisation of the loop integrand was given by [27] (see also ref. [28]) and extended in [29] to dimensions. Importantly, despite the restrictions imposed by onshell conditions in unitarity cuts, sufficient freedom remains in loopmomentum parametrization in order to uniquely determine all integral coefficients for renormalizable theories and beyond.
Modern onshell and unitarity methods may be set up to take advantage of refined physical properties and formal structures of scattering amplitudes. We will discuss the uses of structures like spinor helicity [73, 74], analyticity properties, color decomposition [75, 76, 77] and supersymmetry properties [78, 24] in order to make computations efficient. The scattering amplitudes are then decomposed into a fine set of gauge invariant pieces (primitive amplitudes), which are computed individually and eventually assembled into the full matrix element. This approach leads to excellent numerical stability and can be further exploited [32] for caching and efficiency gains through importance sampling (as used in colorexpansions). In addition to aiming for efficiency it can be helpful to use methods which are easy to automate within existing frameworks [70, 38] or fulfill further computational constraints [40].
Furthermore, a numerical approach to amplitudes requires attention to numerical instabilities induced by roundoff error. A natural way to deal with roundoff errors is to setup a rescue system which monitors numerical precision and invokes an alternative computational approach when checks fail. A convenient rescue strategy [79, 21] for unitarity based approaches is the use of higher precision arithmetic. The advantage of a fine splitup of loop amplitudes into gauge invariant subparts is that one can setup a very targeted and thus efficient rescue system.
The present chapter of this review is organized as follows. In section 2 we discuss a representative example for nexttoleadingorder multijet computations at hadron colliders pointing out the importance of onshell and unitarity methods. In section 3 we discuss the key properties of in matrix elements that can be exploited for the computation of loop amplitudes. Finally, numerical unitarity approaches will be discussed in section 4 and looplevel recursions in section 5. We end with conclusions and an outlook.
2 NLO Predictions for Hadron Colliders
As an example, to point out the key features of NLO QCD predictions we will focus on processes of massive vectorboson production in association with jets. In particular, we focus on recent progress due to the use of numerical unitarity approaches.
The production of massive vector bosons in association with jets at hadron colliders has been the subject of theoretical studies for over three decades [80, 81, 82, 83, 84, 85, 86, 87]. These theoretical studies have for example played an important role in the discovery of the top quark [88]. The oneloop matrix elements for jet and jet production were determined [24] via the unitarity method [22] (see also ref. [89]), and incorporated into the partonlevel MCFM [90] program. Studies of production in association with heavy quarks have also been performed [91, 92, 93, 94, 95, 96].
Beyond this, the numerical unitarity approach allowed to include additional finalstate objects. Studies of jet production can be found in [30, 31, 32, 69] and jet in [34]. The stateoftheart in perturbative QCD for hadron colliders are currently partonlevel nexttoleading order computations with up to five finalstate objects. The first and only such process to be computed so far is jet production [36]. More generally, several important QCD processes with four finalstate objects have been computed to date [97, 35, 98, 37].
Processes of  and boson production in association with jets have a particularly rich phenomenology at the electroweak symmetry breaking scale, being important backgrounds to many searches for new physics and particles, for Higgs physics, and will continue to be important to precision topquark measurements. Decays of the massive vectorbosons into neutrinos mimic missing energy signals and are of particular importance for supersymmetry searches [99, 100, 101]. The clean leptonic decays of the vectorbosons open a high resolution view of underlying QCD dynamics. Inclusive production cross section provides valuable information about parton distribution functions as well as fundamental parameters of the Standard Model. The signal of vectorboson production in association with jets per se includes physics of jetproduction ratios [102] including comparative studies of , and photon production. Experimental studies of vectorboson + jet production at the Tevatron were published by the CDF and D0 collaborations [103, 104, 105] as well as at the LHC by the ATLAS collaboration [106].
2.1 Validation & Prediction
Before turning to the theoretical and technical issues it is useful to assess how good the results are by comparing to experimental data. The quantitative impact of NLO corrections can be directly validated against data from Tevatron collisions.
Fig. 1 compares the distribution of the thirdmost energetic jet in CDF data [103] to the NLO and LO predictions for jet production [30, 32]. (For a similar analysis using a leading color approximation see also [31].) The upper panels of fig. 1 show the distribution itself, while the lower panels show the ratio of the LO value and of the data to the NLO result.The NLO predictions match the data very well, and uniformly with well matching slope. The central values of the LO predictions, in contrast, have different shapes from the data. The change of shape between LO and NLO is attributed to an imprecise scale in the coupling constant, that governs the production of the softest observed jet, in the leading order computation which gets corrected once loop effects are included as discussed in refs. [107, 31].
Scaledependence bands indicate rough estimates of the theoretical error. Those are obtained by varying the renormalization and factorization scales by factors of two around a central scale. For such scale variations, the dependence is of the order of for jet processes at LO. At NLO, the scale dependence shrinks to , and we obtain a quantitatively reliable answer. A more detailed discussion can be found in refs. [30, 32].
As another example, we consider predictions for the LHC. For the inclusive production of jets, a basic quantity to examine is the distribution for the softest observed jet.
Fig. 2 shows the distributions of the softest observed jet in jet () production at LO and NLO, respectively. For details on our analysis setup we refer to [36]. The predictions are normalized to the central NLO prediction in the lower panels. Comparing the distributions successively starting from +1jet production, we observe the reduction in differential cross section of about a factor of from one panel to the next; each observed jet leads to an additional power in the strong coupling. At the same time the lower panels in fig. 2 show an approximately linear increase in the LO scale variation bands, up to about for jets. The scale variation of the NLO result, displayed in the lower panels of fig. 2, is strongly reduced to about for the present setup.
In summary, in the above examples the advantages of NLO computations over the leading order appear through several quantitative improvements. Firstly, NLO predictions show a greatly reduced dependence on unphysical renormalization and factorizationscales as compared to leading order. The second improvement we pointed out concerns the shapes of distributions. Due to inclusion of radiation from an additional parton as well as a more truthful description of the scale dependence shapes of distributions are modeled better at NLO.
The above results validate our understanding of the jet processes for typical StandardModel cuts. It will be interesting, and necessary, to explore the size of corrections for observables and cuts used in newphysics searches. A related process that contributes an irreducible background to certain missing energy signals of new physics is jet production. We expect that the current setup [36] will allow us to compute NLO corrections to jet production, as well as to other complex processes, thereby providing an unprecedented level of theoretical precision for such backgrounds at the LHC.
Partonlevel NLO simulations of this kind are first principle predictions whose outcome directly reflect properties of the underlying theory. Although NLO computations are more challenging, in general they yield results with better reliability and agreement with measurements.
2.2 Setup of Complete Computation
The computation of differential distributions is the end product of combining many important ingredients pulled together in a Monte Carlo program; these include parton distribution functions and couplings, phasespace integration, matrix elements, analysis framework etc. Various tools are available to deal with complete NLO computations. One such tool is MCFM [108], which contains an extensive library of analytic matrix elements for NLO computations. Another approach (see [35] and references) uses tools including Helac [7] and CutTools [79] for a numerical approach [27]. Here we will describe another setup [30, 32, 36] based on onshell and unitarity methods that was used for the computation of jet production in section 2.1.
In addition to LO components of MonteCarlo programs, at NLO the computations rely on further similarly important ingredients. For the jet production [30, 32, 36] the realemission and dipolesubtraction terms [58], are provided by the SHERPA package [12]. SHERPA was used to perform phasespace integration. BlackHat was used to compute the realemission tree amplitudes for jets using onshell recursion relations [47], along with efficient analytic forms extracted from superYangMills theory [109].
In terms of scattering amplitudes we need the input of up to eightpoint oneloop QCD amplitudes as well as up to ninepoint treelevel QCD amplitudes; example Feynman diagrams are shown in fig. 3. The squared matrix elements are summed over all initial and finalstate partons, parton helicities and colororderings. For the present computation the boson is decayed into an observable electron and a neutrino. Amplitudes of this kind can be obtained from QCD amplitudes; with the lepton pair replaced by a quark pair and the boson exchange related to a gluon exchange. Appropriate dressing with coupling constant and propagator terms are needed. A recent analysis of high multiplicity tree amplitudes of this kind can be found in [109].
3 Structure of Oneloop Matrix Elements
The evaluation time of matrix elements is often dominating cross section computations, thus, emphasising the importance of efficient numerical algorithms. Beyond this, matrix elements are objects of fundamental theoretical interest; new physics effects observable at high energy colliders may originate in properties of matrix elements (see e.g. [110]).
Matrix elements are functions of a large number of variables, which characterize particles, polarization states, color quantum numbers, and kinematics. To nexttoleadingorder in the strong coupling the parton level cross sections for resolved finalstate objects, , depend on squared born matrix elements,
(1) 
as well as the interference terms,
(2) 
where , , and are respectively the momentum, helicity (), and color index of the th external gluon or quark. The shorthand refers to the complex conjugate part that has to be included. The efficient management of parton and helicity sums is important. For simplicity, we will consider scattering amplitudes involving quarks and gluons only. Much of the methods can be carried over to more general particle spectra.
As inspired by analytic approaches (see e.g.[111]) we disentangle degrees of freedom in order to arrive at a fine set of gauge invariant objects. To this end several structures are used: Color decomposition into color ordered sub amplitudes disentangles color information and kinematics. Use of spinor helicity notation aligns notation of kinematics and polarization vectors. Spinor variables, in addition, lead to a natural way to work in complex momentum space. This in turn allows to exploit analyticity properties of the basic color ordered scattering amplitudes. Use of a standard basis of integral functions will allow a further fine splitup of the loop amplitude into integral functions and their integral coefficients.
Several features motivate us to disentangle matrix elements into a fine set of primitive amplitudes. First of all, do these have cleaner physical and analytic properties than the full matrix elements, as will be discussed in detail below. This can be exploited for the construction of computational algorithms, as will be discussed in detail below. For numerical approaches, a detailed understanding of physical properties (e.g. IR/UV pole structure of primitive amplitudes, or, tensor rank of integrals) allows one to monitor precision and stability of the computation. Furthermore, caching systems built on primitive objects (here, amplitudes) lead to important efficiently gains through reduction of redundant computations. Finally, during numerical phasespace integration, one can introduce importance sampling, computing computationally expensive but numerically small parts on a subset of phase space points.
3.1 Color Decomposition
At treelevel, amplitudes for gauge theory with external gluons can be decomposed into colorordered partial amplitudes, multiplied by an associated color trace (see e.g. [74]). Summing over all noncyclic permutations gives the full amplitude,
(3) 
the coupling is set to one, and is the set of noncyclic permutations of . The are the set of hermitian generators of the color group. The coefficients of the color structures define the treelevel colorordered partial amplitudes, .
One of the important features of this set of amplitudes is that it forms a closed set under collinear, soft and multiparticle factorization. They have manifest transformation properties under parity transformation and reversal of the ordering of external legs. Similarly, amplitudes with fermions can be decomposed into colorordered subamplitudes [112].
For oneloop amplitudes, one may perform a similar color decomposition to the one at treelevel in eq. (3). In this case, there are up to two traces over color matrices [76],
(4) 
where is the largest integer less than or equal to . The leading colorstructure factor is just times the tree color factor, and the subleading color structures are given by the double trace expressions, . is the set of all permutations of objects, and is the subset leaving invariant.
The leading partial amplitudes take the form (see e.g. [113]),
(5) 
with and being colorordered subamplitudes, or primitive amplitudes. While is fixed to have only gluons propagating in the loop, is restricted to have a Weyl fermion propagating in the loop. The external gluons are colorordered; the diagrams contributing to the loop amplitudes can be generated from colorordered Feynman rules, see e.g. ref. [113].
The coefficients of the subleading color structures; the subleading partial amplitudes, can be expressed in terms of sums [77, 78, 76] of the primitive amplitudes, where different ordering of the external states appear. Beyond the fact that such a decomposition exists, we will not need details of its form here.
The primitive amplitudes,
(6) 
form a generating set of amplitudes, such that given these amplitudes, the full oneloop matrix elements can be computed. For fundamental fermions a similar splitup of partial amplitudes is typically more complicated [78, 24]. In addition to the ordering of the external leg the routing around the loop (left or rightturner) of each of the fermion lines has to be specified. See figure 4 for examples of primitive amplitudes including also external fermion lines.
Primitive amplitudes, have a more transparent analytic structure than full matrix elements, because their external colored legs follow a fixed cyclic ordering. In particular, properties of factorization and branch cut singularities simplify as can be summarized by the following:

Only factorization poles and branch cut singularities in adjacent legs appear.

Primitive amplitudes and colorordered tree amplitudes form a closed set under factorization and unitarity cuts.
With the notation specification of ’adjacent legs’ we refer to the fact that factorization poles and unitarity cuts appear only in a specified subset of kinematic invariants , with an ordering of momenta identical to the one of external gluons. Closure under factorization means that color ordered amplitudes factorize onto color ordered amplitudes. In particular, for unitarity approaches as well as onshell recursions, primitive loop amplitudes can be constructed from colorordered tree amplitudes alone.
3.2 Structure of Loop Amplitude
Any oneloop amplitude can be written as a sum of terms containing branch cuts in kinematic invariants, , and a rational remainder ,
(7) 
The cutcontaining part can in turn be written as a sum over a basis of scalar master integrals [114, 115],
(8) 
The scalar integrals – bubbles, triangles, and boxes – are known functions [116]. We omitted tadpole functions, which in dimensional regularization vanish for massless particles circulating in the loop. The explicit dimension dependence is contained in the integral functions with their coefficient functions strictly four dimensional. Feynman diagrams of the integral functions are shown in fig. 5.
The computation of a oneloop amplitude amounts to determining the rational
coefficient functions and
For the analysis of oneloop amplitudes it is often useful to have two distinct forms of the integrals in mind. One can think of the integral functions as logarithms and polylogarithms of kinematic invariants. As examples we give explicitly a bubble integral function,
(9) 
exposing discontinuities in kinematic invariants through branch cut singularities of the underlying logarithmic functions. For the kinematic invariants the standard definitions, , are used. The constant is defined by, The scale corresponds to the renormalization and factorization scales which, for convenience, are set to be equal here; . The integrals are also viewed as Feynman amplitudes of a scalar field theory,
(10) 
where sums of two momenta are denoted by the shorthand, .
3.3 The Loop Integrand
For explicit computations it is useful to consider loop amplitudes before integration, that is to find a universal parametrization [27, 29] of the loop integrand. In addition to the scalar Feynman diagrams fig. 5, implied from eq. (8), tensorial numerator terms have to included. These tensorial terms describe as well the angular distribution of the virtual particles, which averages out upon integration. The explicit relation between numerator tensors and angular variables will be topic of section 4.3.
A generic form of the loop integrand is given by [27, 29] (see also [28]),
(11)  
where stands for a given discrete dimension and we restricted external momenta and polarizations to be fourdimensional. If we allow dimensional polarization vectors and external momenta, higher polygons have to be added in a natural way. The pentagon terms should be dropped when working in strictly four dimensions, . In the above expression, propagators are denoted by ; for simplicity, of notation we restrict the discussion on massless internal states. Furthermore, we omitted the single propagator terms which drop out of the final results in the absence of massive states. When used with the explicit argument for the loop momentum, denotes the integrand as opposed to the loop amplitude .
The numerators , , and are sums of tensors contracted with the loop momentum . The tensorrank is bounded by powercounting. We will refer to these tensors contracted with loopmomentum somewhat imprecisely as numerator tensors. These numerator tensors can be expanded in terms of a basis of tensors in momentum space multiplied by scalar loopmomentum independent coefficients. The scalar coefficients then characterize the loop amplitude. See below in section 3.3.3 for an explicit representation in terms of a basis of tensors and scalar coefficients.
Integrand parametrizations (11) are common in unitarity approaches; for a discussion in the context of multiloop computations see e.g. refs. [71]. A particularly useful parametrization of the oneloop integrand has been given in [27, 28, 29], as will be discussed further in the following section.
Loop Integration.
With an appropriate representation of the loop integrand the loop integrations can be performed trivially. This is achieved by writing the integrand numerators as a direct sum [27] of terms that integrate to zero and nonvanishing scalar terms. The form of the integrand in eq. (11) then directly relates to the from eq. (8). When we evaluate the known analytic expressions for the basis of integrals, we thus obtain an exact numerical algorithm to go from an offshell integrand to the integrated loop amplitude. An approach of this kind was used in refs. [27, 21, 28, 29]. We will motivate a canonical form for the loop integrand in this section leaving a more complete discussions to the next. Such a choice can be viewed as an implicit integral reduction procedure.
To have an example in mind, consider the box numerators in the form,
(12) 
Here the vector is understood to be orthogonal to the external momenta of the box function. (For an explicit definition and properties of see later in eq. (17).) The coefficients and are the free parameters of the ansatz that have to be determined. The coefficients we eventually need to compute are the coefficients of the scalar term which correspond to the scalar basis integral coefficients in eq. (8). The tensor coefficient , although necessary at intermediate steps of the computation drops out after integration,
(13) 
That is, after integration the numerator loop momentum gets replaces by a linear combination of the external momenta which are orthogonal to ; giving a vanishing tensor integral.
It is instructive to consider another form of the tensor numerator, including a term with an inverse propagator,
(14) 
Clearly the inverse propagator may be cancelled against one of the box propagators turning into the coefficient of a scalar triangle function. If we were to treat this term at the box level, we obtain tensorial contributions that have to be transcribed into scalar integral functions with some care. Assuming for the moment some prior knowledge about generalized unitarity cuts we can make some further observations. That is, this term actually vanishes on the quadruple cut, leaving the coefficient unspecified at first. It may then be fixed using the triplecut equations, although only the sum of the scalar triangle and may be fixed. The numerical unitarity relations eq. (51) are then not triangular. Then box, triangle and bubble coefficients cannot be solved for consecutively. This is of course obvious in the present example, given that we rewrote a scalar triangle coefficient in boxform. This observation emphasises the need for a good basis of numerator tensors.
A less trivial deformation of the numerator tensors would be to mix in a propagator term with the linear tensor,
(15) 
Here again a scalar triangle contribution is pulled back into the box integral. This time the triangular nature of the cutequations (see below in eq. (51) for the explicit form of the cutequations) stays intact and no redundancy is introduced into the numerator tensors. However, one has to pay attention not to drop the coefficient when integrating the tensor box integrals. We can read off the box coefficient directly, , and in addition, we have to include to the related triangle coefficient.
Terms may be moved around between integral functions in this way, effectively introducing a change of basis of integral functions. As a nontrivial application, a cut completion, that is a subtraction of gramdeterminant poles, may be achievable in this way. The form of the numerator tensors has important implication; it allows to keep the unitarity relations triangular, and, keeps the integration of the loop integrand simple.
Numerator Tensors.
In the numerical unitarity approach one is naturally lead to obtain equations for the integrand expression eq. (11). It is then convenient to use an explicit form of the numerators in terms of a basis of tensors. Computing a loop amplitude then amounts to determining the free tensor coefficients.
There are several natural requirements [27] for a good basis of numerator tensors. The first requirement is, that numerator tensors should of course be general enough to parametrize the loopintegrand we are interested in. Typically one uses all tensors up to a given rank, as determined by powercounting. Furthermore, optimally one would like to use a minimal set of tensors. A last requirement is then, that it should be easy to relate the integrand basis back to the integral representation in eq. (8). It turns out that an optimal tensor basis can be found, which satisfies all the above requirements [27, 28, 29].
For numerator tensors in strictly dimensions the tensor basis looks particularly simple. (We will discuss the dimensional generalizations below in section 4.8.) In fact, the result will be a basis of tensors, called spurious numerators in the literature [27], which integrate to zero,
(16) 
where stands for a representative basis tensor. Upon integration, the loopmomentum dependent numerators in eq. (11) may thus be dropped and the remaining scalar (rankzero) terms are directly identified with the integral coefficients in (8).
In order to obtain this basis of tensors it is convenient to introduce the NeervenVermaseren basis [115] for vectors in momentum space; a distinct basis for each of the integral functions. Each integral defines a distinguished set of momenta; the momenta in eq. (16). Momentum space is decomposed into the direct sum of two subspaces; the physical space parametrized by and its complement, spanned by the vectors [28],
(17) 
where we assumed . For the complementary case momentum space is parametrized solely in terms of a linearly independent set of vectors .
The vectors are dual to the external momenta
(18) 
which is naturally related to the metric of momentum space.
A generic numerator tensors can be expressed as tensor products of the vectors (17). A basis of tensors is thus given by,
(19) 
In fact, the set of all numerator tensors needed is given by the symmetric traceless tensors in the transverse space,
(20) 
where the curlbrackets denote the operation of symmetrization and subtraction of traces within transverse space. By trace we mean the contraction,
(21) 
of Lorentz indices with the metric tensor in transverse space, . For symmetric tensors it is sufficient to single out two indices for contraction.
Again the case of a vanishing transverse space is special. For this we are left only with the rankzero scalar numerators. In fact, for this case a further simplification appears; the scalar point integral can be written in terms of lowerpoint integrals. To show this one uses identities implied by inserting vanishing Gram determinants, , into gon integral. Repeated reasoning along these lines leads to reduce gon integrals to gons or lower. We refer to a recent discussion on this in [117] for further details.
Examples of symmetric traceless tensors are then, , which is traceless due to the orthogonality of the vectors . A further example is the tensor, , for which the trace was explicitly subtracted.
The form (20) of the spurious terms can be understood in the following way. To start with, it turns out that tensors (19) with components pointing along the physical space, i.e. , are redundant. For the simplest case (with ) the contraction of with the loop momentum ,
(22) 
gives rise to inverse propagators and , and a scalar term . Although we started with a rankone tensor integrand, after the inverse propagators is cancelled, we obtain lowerpoint integrals and a scalar integral,
(23)  
The tensor integrand we started with is thus redundant, as it can be expresses solely in terms of lowerrank and lowerpoint terms. A similar reasoning, applied recursively, shows the redundancy of tensors with multiple components and enforces in the notation of eq. (19). Tensors including physical directions are thus redundant; they either lead to linear combination of lowerpoint integrals or are expressed as tensors of lower rank. We thus do not need to consider these tensors further, once we account for this.
What remains to be considered are tensors with components purely in the transverse space. Of these only a subset of tensors is linearly independent. In particular, a tracecontaining term in the transverse space is related to a trace in physical space in addition to a metric tensor. Thus, a tracecontaining term yields inverse propagators and loopmomentum independent terms when contracted with loop momentum,
(24) 
where we used equation (18). The contractions of can be transcribed into inverse propagators and terms independent of loop momentum as in eq. (22). Traces thus lead to lowerpoint or lowerrank tensors and are thus linearly dependent. The only choice left are traceless, symmetric tensors in the transverse space.
One might still worry that additional hidden relations can be found to relate integrals with distinct tensor numerators of the from (20). No further relations can in fact be found. The independence of the tensor of the same propagator structures can be argued using an explicit onshell loopmomentum parametrization in eq. (52) as we will discuss further in section 4.3. The independence of tensor integrals (20) with distinct propagator structures is due to their differing factorization properties; e.g. triangle integrals with numerator tensors (20) cannot mimic the quadruple cut divergence of fourpoint integrals.
Tensor integrals with symmetric traceless numerator tensors integrate to zero. Due to Lorentz and parityinvariance, see e.g. [27, 118], a generic tensor integral is written in terms of productions of the vectors and metric tensors ,
(25) 
The integrals simply depend on no other vectors and tensors are excluded by parity invariance. Upon contraction with a symmetric traceless tensor, from (20), the left hand side of (25) turns into a tensor integral (16) while the right hand is easily seen to vanishes; is traceless and the vectors being orthogonal to . We thus verified that symmetric traceless numerator tensors (20) lead to vanishing integrals (16).
In summary, the symmetric traceless tensors (20) fulfill all criteria of an optimal basis as discussed in the beginning of this section.
Tensor Basis.
An explicit form of the numerator tensors in eq. (11) in terms of the vectors (17) was given in [28]. For box coefficients we have,
(26) 
for triangles,
(27)  
and bubbles,
(28)  
respectively. Here we introduced . The vectors differ between the three equations and are defined for each associated propagator structure individually as defined in eq. (17). The coefficients we wish to compute are the , and terms which correspond to the scalar basis integral coefficients. The tensorial expressions vanish upon integration but have to be kept at intermediate steps of the computation.
The above representation is not unique; not only may one chose a different basis for the transverse space and thus different basis vectors . One may also alter the tensor basis used in this parametrization. For example, the above expression uses terms of the form which are not traceless. An alternative, traceless symmetric representation would be instead . The difference of the two approaches amounts to moving loopmomentum tensors between bubble and trianglesfunctions. Given that both forms integrate to zero, either of the above choices leads directly to the same final result. (For a related discussion see also section 3.3.1.)
3.4 Spinor Helicity
Spinor variables give a unified way to express polarization vectors of gluons, fermion helicity states and kinematics of a scattering process. Furthermore, spinor variables lead to a natural way to work with complex momenta. Complexvalued onshell momenta are important in order to fully exploit analyticity properties of amplitudes. We will see examples of this in computations of integral coefficients, sections 4.4 and 4.7, and, later in section 5, when we consider onshell recursions.
We follow the standard spinor helicity notation and conventions as in refs. [73, 74]. As a shorthand notation for the twocomponent (Weyl) spinors we use,
(29) 
Lorentzcovariant spinor products of left and righthanded Weyl spinors can be defined using the antisymmetric tensors and ,
(30) 
These products are antisymmetric, , .
One can reconstruct the momenta from the spinors, using ,
(31) 
Equation (31) shows that a massless momentum vector, written as a bispinor, is simply the product of a lefthanded spinor with a righthanded one. In order to specify onshell momenta we will often use the abbreviation,
(32) 
where spinorial indices are suppressed and the index contractions are indicated by the parenthesis; . Spinor products of the above momentum are then given by, and .
The usual momentum dot products can be constructed from the spinor products using the relation,
(33) 
We will also use the notation,
(34) 
A further important class of quantities are Gram determinants defined by,
(35) 
Gram determinants appear naturally in unitarity cuts; when solving for onshell momenta negative powers of Gram determinants appear. These then enter the computation of integral coefficients.
Gram determinants can be associated to the linearly independent momenta of integralfunctions. The respective integral coefficients typically have inverse powers of these Gram determinants in addition to the ones inherited from reduction of higherpoint tensor integrals. For looplevel onshell recursions, section 5.3.2, we will see that Gram determinants play an important role.
Basic Tree Amplitudes.
When using the spinorhelicity formalism treelevel scattering amplitudes simplify significantly. Further simplifications arises in part from symmetry properties of tree amplitudes [50, 51, 52] that are present in QCDlike theories (see also [53]). Numerical implementations of onshell recursions may recurse all the way to threepoint vertices. More efficiently, they can easily be combined with a library of compact analytic trees. The recursion is stopped and analytic expressions are used whenever available, leading to an efficient numerical algorithm.
One example of simplifications due to the use of spinor helicity variables are the infinite set of vanishing treelevel gluon amplitudes,
(36) 
with all helicities identical, or all but one identical. Parity may of course be used to simultaneously reverse all helicities.
The infinite set of ParkeTaylor amplitudes [119, 120] is another striking example for which the use of spinor helicity formalism yields a particularly simple form,
(37) 
with two negative helicities and the rest positive. The only gluons with negative helicity are in positions and . Helicities are assigned to particles with the convention that they are outgoing. The parity conjugate amplitudes may be obtained by exchanging the left and righthanded spinor products in the amplitude, .
Furthermore, implicit supersymmetry properties [121] allow to relate fermion, gluon and scalar amplitudes of differing spins. For example, in order to replace the gluons and by scalar states,
(38) 
we simply multiply the pure gluon amplitude with an overall factor. Such relations can be used to speed up the sums (1) over finalstate particles when using trees with manifest supersymmetry properties . (See e.g. [49, 122] for trees with manifest supersymmetry properties.)
Onshell Momenta.
For real momenta, and are complex conjugates of each other up to a sign depending on the sign of the energy component. For the degenerate but important case of threepoint kinematics,
(39) 
only the trivial solution can be found for real momenta. For these real solutions all spinor products vanish.
However, for complex momenta, it is possible to choose all three lefthanded spinors to be proportional, , , while the righthanded spinors are not proportional, but obey the relation, , which follows from momentum conservation, . Then,
(40) 
A second branch of solutions to the onshell conditions can be found as the conjugate set of momenta, .
Such degenerate kinematics are important for unitarity cuts associated to integral functions with massless corners. An explicit computation will be discussed in section 4.4 and section 4.7. For such cases threepoint tree amplitudes,
(41) 
have to be evaluated on solutions to the onshell conditions (39). They are nontrivial on one set of complex solutions and vanish as on the other. The above nontrivial solutions involving complex momenta are then necessary in order to exploit generalized unitarity cuts. The general form of this type of onshell conditions is discussed below eq. (55).
3.5 Supersymmetric Decomposition
The supersymmetric decomposition of the amplitudes is particularly useful when considering rational terms of scattering amplitudes. In particular, the rational parts of amplitudes with gluon and fermion degrees of freedom circulating in the loop can be related to often easier to compute scalar ones.
From powercounting arguments we known a priori that the supersymmetric amplitudes and are cut constructible in four dimensions and are free of rational terms [22]. The multiplet and the chiral matter multiplet are built from a particular combination of gluon, fermion and scalar degrees of freedom. For the case of external gluons, the couplings of the matter particles resembles the one of QCD, leading to the relations,
(42) 
between supersymmetric amplitudes (lhs) and basic field theory amplitudes (rhs). The superscripts, and , indicate the states circulating in the loop, and stand for gluon, Weyl fermion and a complex scalar, respectively. Although the above relations are for adjoint fermions in the loop, they can be directly related to massless fundamental quark loops [123, 22].
Inverting the above relations (42) one obtains the amplitudes for QCD via
(43) 
This then implies that the rational terms within and equal the ones from ,
(44) 
With this decomposition we can then compute the cut containing pieces in strictly four dimensions taking into account the full QCD spectrum in the loop. At the same time one may compute the rational part of the QCD amplitudes purely from amplitudes with a complex scalar in the loop.
When computing the rational terms using the dimensional unitarity approach [23], virtual scalars are much more straightforward to deal with as opposed to gluons or fermions. While the kinematics of internal scalars has to be considered beyond four dimensions, we do not have to worry about dimensional extension of polarization states of gluons and fermions. Computations are then very similar to having a massive scalar [23, 124, 72] in the loop, however, where the mass is related to the dimensional momentum.
4 The Unitarity Method
The modern unitarity method [22] in addition to generalized unitarity [24, 25] are the foundation of powerful approaches for loop computations with phenomenological interest. Many recent generalizations [27, 28, 21, 29], in particular, with a numerical application in mind, have helped to established a standard unitarity algorithm. These numerical unitarity methods were first applied to studies of hadron collider physics in [30, 31, 32, 69, 34], and are by now used by many groups [37]. Beyond this, various other implementations of numerical unitarity approaches have been reported [40, 39, 41, 38, 42].
Below we will focus on key developments of the unitarity method with emphasise on numerical aspects. We will follow aspects of the approaches outlined in [28, 21, 29]. For a discussion of analytic unitarity methods we refer to the chapter in the review by Britto [43] and references therein. A more detailed account of the modern unitarity approach as well as its application for multiloop computations may be found in the chapters of this review by Bern and Huang [125] and by Carrasco and Johansson [71].
4.1 Unitarity Relations
In terms of the nonforward part of the Smatrix, , unitarity conditions , imply the nonlinear equations,
(45) 
When combined with analyticity properties, as present in field theory, the unitarity condition (45) relates branch cut discontinuities of scattering amplitudes, to integrals of products of scattering amplitudes. (See e.g. refs. [126] for an early account of unitarity and analyticity.) At oneloop order the unitarity relations may be written as,
(46) 
where the state sum is over all intermediate physical states in the theory. The phasespace integral , is defined over integration contours with the intermediate momenta and onshell; . The notation, , stands for the branch cut discontinuity in the complexified variable . E.g. for a logarithm we have such that the operator picks out the coefficient of the logarithm. Importantly, the nonlinear unitarity relation links onshell amplitudes of different looporder in perturbation theory. For simplicity, we restrict our discussion to color ordered amplitudes.
For fieldtheory amplitudes Cutkosky [127] generalized eq. (46) further, providing a prescription to directly compute more generic discontinuities [128]. An early version of generalized unitarity for generic field theories was demonstrated in [24, 25], including massless states and an arbitrary number of external particles. Specialized to oneloop, the discontinuities are given by phasespace integrals of multiple onshell scattering amplitudes,