Susy Theories and QCD: Numerical Approaches.

# Susy Theories and QCD: Numerical Approaches

## Abstract

We review on-shell and unitarity methods and discuss their application to precision predictions for LHC physics. Being universal and numerically robust, these methods are straight-forward to automate for next-to-leading-order computations within Standard Model and beyond. Several state-of-the-art results including studies of /+3-jet and +4-jet production have explicitly demonstrated the effectiveness of the unitarity method for describing multi-parton scattering. Here we review central ideas needed to obtain efficient numerical implementations. This includes on-shell loop-level recursions, the unitarity method, color management and further refined tricks.

###### pacs:
12.38.Bx, 13.85.-t, 13.87.-a

UCLA-11-TEP-120, NSF-KITP-11-175

## 1 Introduction

The physics program at the Large Hadron Collider (LHC) relies heavily on the ever increasing theoretical control over modeling high-energy 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 phase-space [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 long-distance 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 multi-layered 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 final-state 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 higher-multiplicity leading-order matrix elements into parton-showering 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 higher-order corrections, motivating the challenging quest for perturbative corrections. Next-to-leading (NLO) order predictions significantly reduce renormalization- and factorization-scale 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. Fixed-order results at NLO can also be matched to parton showers [16] with the prospect of complete event generation at next-to-leading-order 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 loop-level on-shell [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 on-shell loop-level 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 multi-particle processes including quantum corrections. Numerical algorithms based on unitarity and on-shell 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 final-states 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: Off-shell recursion relations [44, 45, 46], based on Dyson-Schwinger equations, optimise the reuse of recurring groups of Feynman graphs. In contrast, on-shell 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 QCD-like 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 next-to-leading-order additional challenges arise, in particular, for the virtual corrections due to the loop-momentum integration. NLO cross sections are built from several ingredients: virtual corrections, computed from the interference of tree-level and one-loop amplitudes; real-emission 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 real-emission 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 on-shell 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, loop-level recursions [19, 20, 67, 21] construct amplitudes efficiently using purely on-shell lower-point input. Effective numerically stable implementations of these on-shell methods have been demonstrated by various groups [21, 68, 35, 40, 39, 41, 38, 42]. Beyond this, unitarity approaches have already been applied to state-of-the-art 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, on-shell recursions [20] may as well be applied in certain cases. At present, there are three main choices for computing the rational remainder within a process-nonspecific numerical program: on-shell recursion [21], -dimensional unitarity [23, 29], and a refined Feynman-diagram approach [70, 38]. We will discuss here numerical unitarity approaches in four and dimensions following refs. [21, 29] as well as numerical loop-level on-shell 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 four-point integrals, in turn three-point and, finally, two-point etc. integrals from successively cutting four, three and two propagators. (For the related maximal cut technique for multi-loop 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 on-shell conditions in unitarity cuts, sufficient freedom remains in loop-momentum parametrization in order to uniquely determine all integral coefficients for renormalizable theories and beyond.

Modern on-shell 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 color-expansions). 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 round-off error. A natural way to deal with round-off 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 split-up 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 next-to-leading-order multijet computations at hadron colliders pointing out the importance of on-shell 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 loop-level 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 vector-boson 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 one-loop matrix elements for -jet and -jet production were determined [24] via the unitarity method [22] (see also ref. [89]), and incorporated into the parton-level 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 final-state objects. Studies of -jet production can be found in [30, 31, 32, 69] and -jet in [34]. The state-of-the-art in perturbative QCD for hadron colliders are currently parton-level next-to-leading order computations with up to five final-state objects. The first and only such process to be computed so far is -jet production [36]. More generally, several important QCD processes with four final-state 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 electro-weak 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 top-quark measurements. Decays of the massive vector-bosons into neutrinos mimic missing energy signals and are of particular importance for supersymmetry searches [99, 100, 101]. The clean leptonic decays of the vector-bosons 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 vector-boson production in association with jets per se includes physics of jet-production ratios [102] including comparative studies of , and photon production. Experimental studies of vector-boson + 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 third-most 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].

Scale-dependence 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 +1-jet 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 factorization-scales 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 Standard-Model cuts. It will be interesting, and necessary, to explore the size of corrections for observables and cuts used in new-physics 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.

Parton-level 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, phase-space 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 on-shell and unitarity methods that was used for the computation of -jet production in section 2.1.

In addition to LO components of Monte-Carlo programs, at NLO the computations rely on further similarly important ingredients. For the -jet production [30, 32, 36] the real-emission and dipole-subtraction terms [58], are provided by the SHERPA package [12]. SHERPA was used to perform phase-space integration. BlackHat was used to compute the real-emission tree amplitudes for jets using on-shell recursion relations [47], along with efficient analytic forms extracted from super-Yang-Mills theory [109].

In terms of scattering amplitudes we need the input of up to eight-point one-loop QCD amplitudes as well as up to nine-point tree-level QCD amplitudes; example Feynman diagrams are shown in fig. 3. The squared matrix elements are summed over all initial and final-state partons, parton helicities and color-orderings. 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 One-loop 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 next-to-leading-order in the strong coupling the parton level cross sections for resolved final-state objects, , depend on squared born matrix elements,

 ∑ai,hi|Atreen({ki,hi,ai}|2,n=N,N+1, (1)

as well as the interference terms,

 ∑ai,hiAtree∗n({ki,hi,ai})A1-loopn({ki,hi,ai})+c.c.,n=N, (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 split-up 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 phase-space 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 tree-level, amplitudes for gauge theory with external gluons can be decomposed into color-ordered partial amplitudes, multiplied by an associated color trace (see e.g. [74]). Summing over all non-cyclic permutations gives the full amplitude,

 Atreen({ki,hi,ai})=∑σ∈Sn/ZnTr(Taσ(1)⋯Taσ(n)) Atreen(khσ(1)σ(1),…,khσ(n)σ(n)) , (3)

the coupling is set to one, and is the set of non-cyclic permutations of . The are the set of hermitian generators of the color group. The coefficients of the color structures define the tree-level color-ordered partial amplitudes, .

One of the important features of this set of amplitudes is that it forms a closed set under collinear, soft and multi-particle 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 color-ordered sub-amplitudes [112].

For one-loop amplitudes, one may perform a similar color decomposition to the one at tree-level in eq. (3). In this case, there are up to two traces over color matrices [76],

 A1-loopn({ki,hi,ai})=⌊n/2⌋+1∑c=1∑σ∈Sn/Sn;cGrn;c(σ)An;c(σ), (4)

where is the largest integer less than or equal to . The leading color-structure factor is just times the tree color factor, and the sub-leading 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]),

 A1-loopn;1(1,2,3,⋯,n)=Agn(1,2,3,⋯,n)+nfNcAfn(1,2,3,⋯,n), (5)

with and being color-ordered sub-amplitudes, 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 color-ordered; the diagrams contributing to the loop amplitudes can be generated from color-ordered Feynman rules, see e.g. ref. [113].

The coefficients of the sub-leading color structures; the sub-leading 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,

 Agn(1,2,3,⋯,n)andAfn(1,2,3,⋯,n), (6)

form a generating set of amplitudes, such that given these amplitudes, the full one-loop matrix elements can be computed. For fundamental fermions a similar split-up 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 right-turner) 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:

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

2. Primitive amplitudes and color-ordered 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 on-shell recursions, primitive loop amplitudes can be constructed from color-ordered tree amplitudes alone.

For a more complete description of color decomposition we refer the reader to previous reviews [78, 24, 74].

### 3.2 Structure of Loop Amplitude

Any one-loop amplitude can be written as a sum of terms containing branch cuts in kinematic invariants, , and a rational remainder ,

 An=Cn+Rn. (7)

The cut-containing part can in turn be written as a sum over a basis of scalar master integrals [114, 115],

 Cn=∑i1

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 one-loop amplitude amounts to determining the rational coefficient functions and 1 in addition to the remainder . Following the spinor-helicity method [73, 81], we can then re-express all external momenta in terms of spinors. The coefficients of these integrals, , and , as well as the remainder , are then all rational functions of appropriate spinor and momentum variables.

For the analysis of one-loop 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,

 I2(s12,^μ2) = Missing or unrecognized delimiter for \biggl (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,

 I2(s12,^μ2) = ^μ2ϵ∫d4−2ϵℓ(2π)4−2ϵ1ℓ2(ℓ−K12)2, (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]),

 A1-loop,dn(ℓ) = ∑i1<…

where stands for a given discrete dimension and we restricted external momenta and polarizations to be four-dimensional. 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 tensor-rank is bounded by power-counting. We will refer to these tensors contracted with loop-momentum somewhat imprecisely as numerator tensors. These numerator tensors can be expanded in terms of a basis of tensors in momentum space multiplied by scalar loop-momentum 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 multi-loop computations see e.g. refs. [71]. A particularly useful parametrization of the one-loop 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 non-vanishing 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 off-shell 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,

 ¯di1i2i3i4(ℓ) = d0i1i2i3i4+d1i1i2i3i4(n1⋅ℓ). (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,

 ∫d4ℓ(2π)4−2ϵ(n1⋅ℓ)Di1Di2Di3Di4=0. (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,

 ¯di1i2i3i4(ℓ) = d0i1i2i3i4+d1i1i2i3i4(n1⋅ℓ)+d2i1i2i3i4Di1(ℓ). (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 triple-cut 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 box-form. 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,

 ¯di1i2i3i4(ℓ) = Missing or unrecognized delimiter for \biggl (15)

Here again a scalar triangle contribution is pulled back into the box integral. This time the triangular nature of the cut-equations (see below in eq. (51) for the explicit form of the cut-equations) 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 non-trivial application, a cut completion, that is a subtraction of gram-determinant 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 loop-integrand we are interested in. Typically one uses all tensors up to a given rank, as determined by power-counting. 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,

 0 = ^μ2ε∫d4−2εℓ(2π)4−2εnμ1⋯μkℓμ1⋯ℓμkℓ2(ℓ−K1)2⋯(ℓ−K1⋯−Kn−1)2,k>0, (16)

where stands for a representative basis tensor. Upon integration, the loop-momentum dependent numerators in eq. (11) may thus be dropped and the remaining scalar (rank-zero) terms are directly identified with the integral coefficients in (8).

In order to obtain this basis of tensors it is convenient to introduce the Neerven-Vermaseren 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],

 Kifori∈{1,⋯,n−1},nifori∈{1,⋯,d−(n−1)}, ni⋅nj=δij,ni⋅Kj=0, vifori∈{1,⋯,n−1},Ki⋅vj=δji,ni⋅vj=0, (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 2 and are part of the physical space, defined by the external momenta . The vectors are an orthonormal basis in transverse space and are orthogonal to the physical space. Depending on the signature that the transverse space inherits from momentum space, the vectors have to be chosen purely real or imaginary. A further useful quantity is the metric of the transverse space,

 (g⊥)μ1μ2≡d−n+1∑i=1niμ1niμ2=gμν−n−1∑i=1viμKνi, (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,

 nμ1⋯μk=Ki1μ1…Kilμlnj1μ1…njmμm,l+m=k. (19)

In fact, the set of all numerator tensors needed is given by the symmetric traceless tensors in the transverse space,

 nμ1⋯μk=ni1{μ1ni2μ2⋯nikμk}, (20)

where the curl-brackets denote the operation of symmetrization and subtraction of traces within transverse space. By trace we mean the contraction,

 nμ1μ2⋯μngμ1μ2⊥=0, (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 rank-zero scalar numerators. In fact, for this case a further simplification appears; the scalar -point integral can be written in terms of lower-point 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 ,

 ℓ⋅K1=12[ℓ2−(ℓ−K1)2+K21], (22)

gives rise to inverse propagators and , and a scalar term . Although we started with a rank-one tensor integrand, after the inverse propagators is cancelled, we obtain lower-point integrals and a scalar integral,

 K1⋅ℓℓ2(ℓ−K1)2D2⋯Dn = 12[1(ℓ−K1)2D2⋯Dn−1ℓ2D2⋯Dn+ (23) +K21ℓ2(ℓ−K1)2D2⋯Dn].

The tensor integrand we started with is thus redundant, as it can be expresses solely in terms of lower-rank and lower-point 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 lower-point 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 trace-containing term in the transverse space is related to a trace in physical space in addition to a metric tensor. Thus, a trace-containing term yields inverse propagators and loop-momentum independent terms when contracted with loop momentum,

 d−n+1∑i=1niμniνℓμℓν = ℓ2−n−1∑i=1(vi⋅ℓ)(Ki⋅ℓ), (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 lower-point or lower-rank 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 on-shell loop-momentum 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 four-point integrals.

Tensor integrals with symmetric traceless numerator tensors integrate to zero. Due to Lorentz- and parity-invariance, see e.g. [27, 118], a generic tensor integral is written in terms of productions of the vectors and metric tensors ,

 ^μ2ε∫d4−2εℓ(2π)4−2εℓμ1⋯ℓμkℓ2(ℓ−K1)2⋯(ℓ−K1⋯−Kn−1)2= =K{μ1i1…Kμlilgμl+1μl+2…gμk−1μk}f(sij)+⋯. (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,

 ¯di1i2i3i4(ℓ) = d0i1i2i3i4+d1i1i2i3i4t1, (26)

for triangles,

 ¯ci1i2i3(ℓ) = c0i1i2i3+c1i1i2i3t1+c2i1i2i3t2+c3i1i2i3(t1t1−t2t2) (27) +t1t2(c4i1i2i3+c5i1i2i3t1+c6i1i2i3t2),

and bubbles,

 ¯bi1i2(ℓ) = b0i1i2+b1i1i2t1+b2i1i2t2+c3i1i2t3 (28) +b4i1i2(t1t1−t3t3)+b5i1i2(t2t2−t3t3) +b6i1i2t1t2+b7i1i2t1t3+b8i1i2t2t3,

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 loop-momentum tensors between bubble- and triangles-functions. 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. Complex-valued on-shell 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 on-shell recursions.

We follow the standard spinor helicity notation and conventions as in refs. [73, 74]. As a shorthand notation for the two-component (Weyl) spinors we use,

 (λi)α≡[u+(ki)]α,(~λi)˙α≡[u−(ki)]˙α,k2i=0. (29)

Lorentz-covariant spinor products of left- and right-handed Weyl spinors can be defined using the antisymmetric tensors and ,

 ⟨jl⟩=ϵαβ(λj)α(λl)β,[jl]=ϵ˙α˙β(~λj)˙α(~λl)˙β. (30)

These products are antisymmetric, , .

One can reconstruct the momenta from the spinors, using ,

 kμi(σμ)α˙α=(to7.499886pt\hss/\hssto0.0pt$k$i)α˙α=(λi)α(~λi)˙α. (31)

Equation (31) shows that a massless momentum vector, written as a bi-spinor, is simply the product of a left-handed spinor with a right-handed one. In order to specify on-shell momenta we will often use the abbreviation,

 ℓ=λi(to7.499886pt\hss/\hssto0.0pt$K$λj), (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,

 ⟨lj⟩[jl]=2kj⋅kl=sjl. (33)

We will also use the notation,

 ⟨a|Ki1…im|b] = m∑k=1⟨aik⟩[ikb],andsi1⋯im=K2i1⋯im. (34)

A further important class of quantities are Gram determinants defined by,

 Δn(K1,⋯,Kn)≡det(2Kj⋅Kl). (35)

Gram determinants appear naturally in unitarity cuts; when solving for on-shell 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 integral-functions. The respective integral coefficients typically have inverse powers of these Gram determinants in addition to the ones inherited from reduction of higher-point tensor integrals. For loop-level on-shell recursions, section 5.3.2, we will see that Gram determinants play an important role.

#### Basic Tree Amplitudes.

When using the spinor-helicity formalism tree-level scattering amplitudes simplify significantly. Further simplifications arises in part from symmetry properties of tree amplitudes [50, 51, 52] that are present in QCD-like theories (see also [53]). Numerical implementations of on-shell recursions may recurse all the way to three-point 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 tree-level gluon amplitudes,

 Atreen(1±,2+,…,n+)=0,n>3, (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 Parke-Taylor amplitudes [119, 120] is another striking example for which the use of spinor helicity formalism yields a particularly simple form,

 Missing or unrecognized delimiter for \right (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 right-handed 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,

 Missing or unrecognized delimiter for \left (38)

we simply multiply the pure gluon amplitude with an overall factor. Such relations can be used to speed up the sums (1) over final-state particles when using trees with manifest supersymmetry properties . (See e.g. [49, 122] for trees with manifest supersymmetry properties.)

#### On-shell 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 three-point kinematics,

 k!+k2+k3=0,withk2i=0, (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 left-handed spinors to be proportional, , , while the right-handed spinors are not proportional, but obey the relation, , which follows from momentum conservation, . Then,

 ⟨ij⟩≠0, but[ij]=0. (40)

A second branch of solutions to the on-shell 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 three-point tree amplitudes,

 Atree3(1−,2−,3+)=i⟨12⟩4⟨12⟩⟨23⟩⟨31⟩, (41)

have to be evaluated on solutions to the on-shell conditions (39). They are non-trivial on one set of complex solutions and vanish as on the other. The above non-trivial solutions involving complex momenta are then necessary in order to exploit generalized unitarity cuts. The general form of this type of on-shell 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 power-counting 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,

 AN=4n ≡Agn+4Afn+3Asn, AN=1chiraln ≡Afn+Asn, (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

 Agn = AN=4n−4AN=1chiraln+Asn, Afn = AN=1chiraln−Asn. (43)

This then implies that the rational terms within and equal the ones from ,

 Agn|rational=Asn|rational,Afn|rational=−Asn|rational. (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.

Relations of this kind are rather generic and can be found for internal fermions and mixed fermion and gluon amplitudes. For analytic computations this observations were used for example in refs. [78, 24].

## 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 multi-loop 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 non-forward part of the S-matrix, , unitarity conditions , imply the nonlinear equations,

 −i(T−T†)=T†T. (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 one-loop order the unitarity relations may be written as,

 Disc(si)A1-loopn=∑states∫dΦ(ℓ1,ℓ2)Atreen1(−ℓ2,ℓ1)Atreen2(−ℓ1,ℓ2), (46)

where the state sum is over all intermediate physical states in the theory. The phase-space integral , is defined over integration contours with the intermediate momenta and on-shell; . 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 on-shell amplitudes of different loop-order in perturbation theory. For simplicity, we restrict our discussion to color ordered amplitudes.

For field-theory 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 one-loop, the discontinuities are given by phase-space integrals of multiple on-shell scattering amplitudes,

 Disc(si1,⋯sik)A1-loopn=∑states ∫d<