Adaptive Integrand Decomposition in parallel and orthogonal space

We present the integrand decomposition of multiloop scattering amplitudes in parallel and orthogonal space-time dimensions, $d=d_\parallel+d_\perp$, being $d_\parallel$ the dimension of the parallel space spanned by the legs of the diagrams. When the number $n$ of external legs is $n\le 4$, the corresponding representation of the multiloop integrals exposes a subset of integration variables which can be easily integrated away by means of Gegenbauer polynomials orthogonality condition. By decomposing the integration momenta along parallel and orthogonal directions, the polynomial division algorithm is drastically simplified. Moreover, the orthogonality conditions of Gegenbauer polynomials can be suitably applied to integrate the decomposed integrand, yielding the systematic annihilation of spurious terms. Consequently, multiloop amplitudes are expressed in terms of integrals corresponding to irreducible scalar products of loop momenta and external momenta. We revisit the one-loop decomposition, which turns out to be controlled by the maximum-cut theorem in different dimensions, and we discuss the integrand reduction of two-loop planar and non-planar integrals up to $n=8$ legs, for arbitrary external and internal kinematics. The proposed algorithm extends to all orders in perturbation theory.

The decomposition of multiloop scattering amplitudes in terms of independent functions, together with the subsequent determination of the latter, is a viable alternative -often the only accessible one -to the direct integration, which, for non-trivial processes, may require the calculation of a prohibitively large number of complicated Feynman integrals. Understanding the properties of Feynman integrands has led to the development of novel algorithms aiming to the automated determination of partonic cross sections for highmultiplicity processes which have been successfully applied, in the last decade, to one-loop amplitudes. More generally, the use of unitarity-based methods and integrand decomposition algorithms has shown that exploiting the algebraic properties of the integrands may lead to the discovery of novel properties of the amplitudes, hidden beneath the superficial look of Feynman integrals' representation, which, if properly engineered, may turn into drastic simplifications for their evaluation.
In this paper, we elaborate on a representation of dimensionally regulated Feynman integrals where, for any given diagram, the number of space-time dimensions d (= 4 − 2 ) is split into parallel (or longitudinal) and orthogonal (or transverse) dimensions, as d = d + d ⊥ [1][2][3][4][5][6]. Accordingly, the parallel space is spanned by the independent four-dimensional external momenta of the diagram, namely d = n − 1, where n is the number of legs, whereas the transverse space is spanned by the complementary orthogonal directions. For diagrams with a number of legs n ≥ 5, the orthogonal space embeds the −2 regulating dimensions, d ⊥ = −2 , while, for diagrams with n ≤ 4, the orthogonal space is larger and it embeds, beside the regulating dimensions, also the four-dimensional complement of the parallel space, namely d ⊥ = (5 − n) − 2 . For this reason, the decomposition of the space-time dimensions in parallel and orthogonal directions can be considered as adaptive, since it depends on the number of legs of the individual diagram.
Decomposing the loop momenta q α i in terms of parallel and orthogonal vectors, q α i = q α i + λ α i , has the immediate advantage of exposing a subset of integration variables which can be trivially integrated away, hence they can be eliminated from the calculation before applying any reduction algorithm. In fact, multidimensional polar coordinates can be suitably introduced in order to parametrize the integral over the orthogonal space in terms of integrations over radial variables λ ii (= λ i · λ i ) and a generalised solid angle. This change of coordinates makes manifest that numerators and denominators of Feynman integrands do not depend on the same set of integration variables. Indeed, the quadratic Feynman denominators depend only on the parallel directions, on the radial variables λ ii and the relative orientations λ ij , i < j, of the transverse vectors, but they do not depend on their individual components, which can be mapped into a set of angular variables Θ ⊥ . Conversely, the numerators may depend on all variables. In the case of diagrams with n ≤ 4, the dependence of the integrand on transverse angles, say θ i , is polynomial in sin θ i and cos θ i , therefore, the integration over Θ ⊥ can be trivially performed. In this article, we show how it can be carried out by means of the orthogonality relation for Gegenbauer polynomials, as an alternative to the Passarino-Veltman tensor reduction used in ref. [2].
After integrating over the transverse angles Θ ⊥ , the integrand will solely depend on q i and on the λ ij variables appearing in the denominators. These variables correspond to (reducible and irreducible) scalar products between loop momenta and external momenta.

JHEP08(2016)164
The integration over orthogonal and parallel space has been used to evaluate multi-scale Feynman integrals, up to two-and three-point functions [2][3][4][5][6]. The goal of this communication is instead discussing how the decomposition of space-time into parallel and orthogonal subspaces simplifies the multiloop integrand reduction algorithm [7][8][9][10][11]. Namely, our objective is not the evaluation of Feynamn integrals, rather their decomposition in terms of independent integrals. We show that this procedure can be applied to arbitrarily complicated diagrams. In particular, we consider the decomposition up to two-loop eight-point planar and non-planar integrals and we discuss how the same procedure can be extended to higher orders. Previous studies of higher-loop integrands in four-dimensions can be found in [12][13][14].
Feynman integrals are multivariate integrals of rational integrands and they can be decomposed in terms of a set of irreducible integrals (IRIs) by multivariate polynomial division [9,10]. In fact, the partial fractioning of Feynman integrands amounts to iterative divisions (modulo Gröbner bases) between the numerator and the denominators, once they are written as polynomials in the components of the integration momenta in a given basis. The resulting integrand decomposition is a sum of integrands whose denominators are given by all the possible partitions of the initial set of denominators, and whose numerators correspond to the remainders of the division w.r.t. the set of denominators they sit on. The remainders of the division contain, by definition, terms which cannot be expressed in terms of denominators. In fact, since each component of a given integration momentum corresponds to a scalar product of that momentum with an element of the momentum basis, the remainder should contain only irreducible scalar products (ISPs). On the contrary, reducible scalar products (RSPs) can be decomposed in terms of denominators and external invariants.
The integrand decomposition is effectively a unitarity-based decomposition of the integrand, since each remainder can be considered as the residue of the cut identified by the simultaneous vanishing of the corresponding denominators. It should be observed that the integrand reduction can be applied as well to the case of integrals whose denominators are raised to powers higher than one [15]. Integrating the decomposed integrand over the loop momenta corresponds to the decomposition of the original integral in terms of IRIs. In fact, upon integration, some of the ISPs in the residues may generate vanishing integrals: these terms are called spurious, because although present in the integrand decomposition, they do not contribute to the amplitude. Instead, the non-spurious ISPs correspond to the (numerators of) IRIs appearing in the amplitude decomposition. Therefore, within the integrand decomposition algorithm, the reduction of any scattering amplitude in terms of IRIs turns into the algebraic problem of determining the coefficients of the monomials of the residues.
The basic elements of the integrand decomposition algorithm are: i) the space-time dimensions, namely the number of integration variables; ii) the momentum basis used for the decomposition of the loop momenta; iii) the structure of the numerators and the variables they depend on; iv) the form of the denominators and the variables they depend on; v) the structure of the residues; vi) the solutions of the cut equations. The integrand reduction algorithm was originally proposed for one-loop integrals in four dimensions [16,17] and later extended to d = 4−2 dimensions [18][19][20][21], to deal with dimensionally regulated amplitudes (see [22] for a review). In the one-loop case, the residues were built by using JHEP08(2016)164 two driving principles: on the one side, the knowledge of the set of IRIs which could appear in the decomposition of generic one-loop integrals [23]; and, on the other side, the Lorentz covariance of spurious terms which could additionally appear in the numerators.
The integrand reduction algorithm for one-loop integrals has been implemented in several public libraries, like Cutools [24], Samurai [25] and Ninja [26,27], which played an important role in the development of codes for the automatic evaluation of scattering amplitudes for generic scattering processes at NLO accuracy, as recently reviewed in [28]. In particular, Ninja implements an ameliorated integrand decomposition algorithm [26], which introduced the idea of the (univariate) polynomial division for the calculation of the residues.
In order to extend the integrand decomposition at higher orders [7,8], the same driving principles could not be applied. The first reason for this is that the basis of independent integrals is not known. Moreover, the interplay of more integration momenta makes the classification of the spurious terms less obvious. One additional difference w.r.t. the oneloop case, which was indeed to be expected, is the contribution of integrals corresponding to non-spurious ISPs [7]. Nevertheless, the systematic determination of the residues at higher order was systematized by means of algebraic geometry methods [9,10], namely the polynomial division modulo Gröbner basis. An implementation of such method is provided by the public package BasisDet [9]. Integrand decomposition beyond one-loop has been successfully applied to a first case of non trivial two-loop five-point helicity amplitude in [29,30].
One of the main outcomes of the multivariate polynomial division algorithm is the so called maximum-cut theorem [10], which can be applied whenever the on-shell conditions are sufficient in order to constrain all integration variables. In this case, the system of equations is zero-dimensional and the remainder of the division (of a numerator that depends on all variables constrained by the cut-conditions) can be cast as a univariate polynomial of degree n s − 1, being n s the number of solutions of the system. This theorem extends to all loops and to all dimensions the beauty of the four-dimensional quadruple-cut [31], which is known to have two solutions and whose residue is parametrized in terms of two monomials [16]. The number of integration variables depends on the dimensions of the loop momenta, hence, in order to freeze all of them, the number of denominators to be put on-shell depends on the space-time dimensions as well. In other words, maximum-cuts are realized by cutting diagrams with different number of external legs, according to the dimensionality of the integration momenta.
The use of the d = d + d ⊥ representation of Feynman integrals within the integrand reduction technique has several interesting consequences. To explore them, we propose a three-step algorithm, which we will refer to it as divide-integrate-divide.
Divide. First, by separating the physical directions from the (−2 )-dimensional yields simpler on-shell cut conditions, hence the division procedure becomes significantly simpler. In fact, the Gröbner basis trivialize, as they are linear in the variables to be reduced and quadratic in the irreducible variables which will appear in the residues (up to the choice of monomial order). In this case, the Gröbner bases are built from differences of denominators (basic S-polynomials). Moreover, the form of the cut conditions and of the Gröbner bases JHEP08(2016)164 are further simplified in the d = d + d ⊥ representation, due to the dependence of the denominators on a reduced set of variables, hence the determination of the cut-residues becomes computationally less expensive. We can properly talk of adaptive cutting, since the dimensions of the parallel space, i.e. the number of variables constrained by the on-shell conditions, depend on the number of legs.
Integrate. Second, after the integrand reduction, the integration over the orthogonal solid angle of the decomposed integrand allows the automatic detection and annihilation of the spurious integrals, which vanish because of the orthogonality condition enforced by the Gegenbauer polynomial integration. Within the proposed parametrization, the spherical symmetry of the transverse angular integrations offers an explicit geometric interpretation of the spurious integrals as being related to monomials which are odd under rotation group transformations, as observed in [32]. Alternatively, if the integration over Θ ⊥ is performed before the reduction, the corresponding residue will not contain any spurious term, therefore the number of non-vanishing coefficients to be determined through the reduction algorithm will be significantly smaller.
Divide. Finally, we notice that the integration of the residues over the transverse angles Θ ⊥ reintroduces, in general, a dependence on the variables λ ij . The denominators depend on these variables and, therefore, the integrated residues may be subject to a second polynomial division, which further simplifies them. In some cases, namely when the variables λ ij form a Gram determinant, this additional division can be shown to be equivalent to applying dimensional shifting recurrence relations [33,34] at the integrand level (the dimensions of any Feynman integral are controlled by the power of the Gram determinant, characteristic of each loop).
We now observe that after the integrand decomposition outlined above, the integrand will depend on a subset of the parallel space variables and on the transverse variables λ ij , which correspond just to irreducible scalar products (ISPs) between loop momenta and external momenta. Therefore, any scattering amplitude, at any loop order and with arbitrary kinematics, admits an algebraic decomposition in terms of a set of irreducible integrals (IRIs), corresponding to these ISPs.
It is important to stress that, although independent under polynomial division, the IRIs are not a minimal set. Indeed, they can be related through identities which belong to the general class of integration-by-parts relations (IBPs), hence their number can be further reduced. The amplitude, in this case, would be finally expressed in terms of a minimal set of Master Integrals (MIs). IBPs relation for IRIs can be suitably built by algebraic geometry methods through sygyzy equations [32,35,36]. In particular, the outcome of the proposed integrand reduction algorithm is suitable for an IBP-reduction in the parallel space along the lines of [36,37]. Progress on the improvement of system solving strategies for IBP equations are under intense development [38,39]. Moreover, should the reduction to MIs not be available for the process under consideration, the representation of the amplitudes in terms of IRIs can be employed together with the numerical integration of the latter. Promising advances on the numerical integration of Feynman integrals have recently been applied to a non-trivial two-loop case [40].

JHEP08(2016)164
The paper is organized as follows. In section 2, we discuss the d = d + d ⊥ representation of multiloop Feynman integrals and the integration over the transverse directions by means of the orthogonality relation for Gegenbauer polynomials. Besides analysing the properties of the transverse space for general topologies with n ≤ 4 external legs, we discuss further simplifications that can be obtained for factorized and ladder topologies. In fact, we show that the integration of Gegenbauer polynomials can be used in all cases where the numerator depends on more variables than the denominators. As an example of the considerably simplified form of Feynman integrands achieved by integrating out the transverse directions before applying any reduction algorithm, we discuss a four-point contribution to the helicity amplitude A(g + 1 , g − 2 , g + 3 , g − 4 ) at two loops. In section 3, we present the adaptive integrand decomposition algorithm for multiloop scattering amplitudes. We revisit the wellknow results for the one-loop integrand decomposition, by showing that, in d = d + d ⊥ , all unitarity cuts are reduced to zero-dimensional systems, and by providing an alternative representation of the residues, whiche read as complete polynomials in the transverse variables, due to the maximum-cut theorem. The novel parametrization of the residues emerging in d = d + d ⊥ yields a different, yet equivalent, decomposition of one-loop amplitudes w.r.t. to the known decomposition in d = 4 − 2 . At two loops, we provide a classification of the residues appearing in the integrand decomposition formula for planar and non-planar topologies with arbitrary kinematics, by considering the top-down reduction from the eightpoint maximum-cut topologies, down to the product of two one-point topologies, namely two tadpoles. As a concrete example of the application of the adaptive division algorithm, we provide the explicit expression of the coefficients of the residue of the double-box con- . In section 4, we give our summary and conclusions. We have collected in the appendices the detailed discussion of most of the calculations leading to the results presented in this work. In appendix A, we propose a new derivation of the parametric expression of Feynman integrals in terms of parallel-and transversespace variables and we discuss the change of coordinates to be performed in the transverse space in order to map, at any loop order, all integrations over the four-dimensional transverse directions into simple angular integrals. In appendices B-C, we collect some useful formulae for one-and two-loop integrals respectively, including a list of tensor integrals which can be reduced by integrating over the transverse angles. In appendix D we recall the main properties of Gegenbauer polynomials and, finally, in appendix E, we provide a representation in terms of spinor variables of the momentum-basis to which we refer throughout the text. The calculations presented in this paper have been performed with the help of Singular [41] and S@M [42].

Parallel and orthogonal space for multiloop Feynman integrals
In this section, we consider generic -loop Feynman integrals with n external legs in a d-dimensional Euclidean space, (2.1)

JHEP08(2016)164
where N (q i ) is an arbitrary tensorial numerator and the denominators D j (q i ) are defined as being {p 1 , . . . , p n−1 } the set of independent external momenta and α and β incidence matrices taking values in {0, ±1}. We first recall the usual parametrization of I d ( ) n obtained by formally splitting the d-dimensional space into the four-dimensional physical one, where external momenta and polarizations lie, and the corresponding orthogonal subspace, whose dimension is conventionally set to d − 4 = −2 . Later, we show that, when a Feynman integral has n ≤ 4 external legs which do not span the full physical space, I d ( ) n is more conveniently expressed in terms of vectors living in the d = n − 1 dimensional space described by the external kinematics and a set of transverse variables belonging to its orthogonal complement with dimension d ⊥ = d − n − 1. This alternative parametric representation of Feynman integrals remarkably simplifies, at any loop order, the integration over the transverse components of the loop momenta.

Feynman integrals in
When dealing with a regularization scheme where the external kinematics is kept in four dimensions, it is customary to split the d-dimensional loop momenta into a four-dimensional part and a (−2 )-dimensional one, so that, by defining µ ij = µ i · µ j , we have The vectors µ α i lie in a subspace which is completely orthogonal to the four-dimensional one, µ i · p j = 0, hence we can rewrite the denominators (2.2) as For the same reason, the numerator appearing in (2.1) can depend on q α i [4] and µ ij only. This implies that the integrals over the (−2 )-dimensional subspace can be expressed into spherical coordinates, and that we can integrate out all directions orthogonal to the relative orientations of the vectors µ α i , obtaining , (2.7) -7 -

JHEP08(2016)164
As it is explicitly shown in (2.5), because of this parametrization, the set of denominators which characterizes each integral depends, in general, on the same ( + 9)/2 variables as the numerator, corresponding to the 4 four-dimensional components of the loop momenta, which are decomposed into some basis of four-dimensional vectors {e α i }, x ji e α j , (2.8) and the ( + 1)/2 scalar products µ ij . It should be noticed that, the denominators of particular classes of multi=loop Feynman integrals, such as ladder topologies and factorized integrals, might depend on a reduced number of variables µ ij , due to the absence of propagators involving both loop momenta q α i and q α j . In the following, we first derive an integral parametrization alternative to (2.6), valid for Feynman integrals with n ≤ 4 external legs, by assuming that the denominators depend on the maximal number of loop variables, and then we show how a simplified parametrization for ladder and factorized integrals as well.

Feynman integrals in
For a number of external legs n ≤ 4, it is possible to reparametrize the integral (2.1) in such a way that the number of variables appearing in the denominators is reduced to ( + 2d + 1)/2. Since the numerator is a polynomial in the remaining (4 − d ) variables, their integration can be performed straightforwardly, by decomposing the numerator in terms of orthogonal polynomials. In fact, the choice of the four-dimensional basis {e α i } is completely arbitrary and, if d ≤ 3, one can choose 4 − d vectors of such basis to lie into the subspace orthogonal to the external kinematics, i.e.
In this way, we can rewrite the d-dimensional loop momenta as where q α i is a vector of the d -dimensional space spanned by the external momenta, and l=d +1 x li x lj + µ ij , (2.12) belong to the d ⊥ -dimensional orthogonal subspace. In this parametrization, all denominators become independent of the transverse components of the loop momenta,

JHEP08(2016)164
and they depend on a reduced set of ( + 2d + 1)/2 variables, corresponding to the d components of q α i plus the ( + 1)/2 scalar products λ ij . Once the decomposition (2.10) has been introduced, it can be shown that all transverse components x ji (j > d ) as well as the relative orientations of the vectors λ α i can be mapped into angular variables, defined through a suitable orthonormalization procedure described in appendix A. In particular, by introducing the angles we can define a polynomial transformation of the type so that the integral (2.1) can be rewritten as and The choice of the four-dimensional basis {e α i } and the consequent definition of the transverse space variables Λ and Θ ⊥ are determined by the external kinematics and do not depend on the specific set of denominators which characterizes the integral, therefore, the parametrization (2.16) can be used for both planar and non-planar topologies. Moreover, it should be observed that, in the case of two-point integrals with vanishing external JHEP08(2016)164 momentum p 2 = 0, the r.h.s. of (2.16) holds for d = 2, since we can define only two directions orthogonal to a massless vector.
The decomposition in eq. (2.10) allows us to express a subset of components of q α i and λ ij as combinations of denominators by solving linear equations. In fact, one can always build differences of denominators which are linear in the loop momenta and independent of λ ij , while the relation between λ ij and the denominators is always linear by definition, as can be deduced from eq. (2.13).
More explicitly, at one-loop, all the denominators can be taken to have the form where r is the total number of denominators in the loop integrand. Hence one can choose any denominator D and consider r − 1 differences of the form D j − D. These differences have no quadratic terms in the loop momenta and can thus be used to express r − 1 of the variables {x ji , j ≤ d } as linear combinations of denominators. By applying one more independent equation, given by the definition of any of the denominators, the variable λ 11 is written as a linear combination of the variables {x ji , j ≤ d }, as one can see from eq. (2.13).
At higher loops, one can split the r loop denominators into partitions identified by the subset of loop momenta each denominator depends on, and similarly consider differences of denominators belonging to the same partition which will again generate a set of linear relations between physical loop components and denominators. By solving these relations, one can express a subset of the variables {x ji , j ≤ d } as linear combinations of denominators. Finally one can, again, consider eq. (2.13) for a representative of each partition of denominators, completing the set of linear relations which can thus be solved for a subset of the variables λ ij . It is straightforward to see that the complete set of relations is equivalent to the definition of the loop denominators themselves.
For instance, at two loops, one can have at most three partitions P 1 , P 2 , P 3 , which respectively correspond to denominators having the following forms where P 1 is the set of denominators depending on q 1 only, P 2 , the set of denominators depending on q 2 only, and P 3 , the set of denominators depending both on q 1 and q 2 . Therefore one can choose a representative for each partition, say D i ∈ P i for i = 1, 2, 3, and observe that for any j ∈ P i the difference D j − D i is linear in the loop momenta. This allows us to write r − 3 linear equations which can be solved for a subset of the variables JHEP08(2016)164 {x ji , j ≤ d } in terms of the other physical directions and denominators. One can thus complete this set of relations with 3 more equations (or possibly less, if any of the partitions is empty) which are defined by eq. (2.13) applied to one denominator for each partition P i . In the case when none of the partitions is empty, these three equations can be solved for the variables λ 11 , λ 12 , λ 22 which are thus written as a combination of denominators and irreducible components of q α i by solving linear relations. If the denominators are independent of λ 12 , this variable cannot obviously be written in terms of denominators but it can be integrated out by means of the techniques presented later on in this paper. As we shall see in section 2.5, this is true at any number of loops, whenever the denominators are independent of one of these variables.
The observations made in this paragraph imply that by solving the d-dimensional cut constraints for the integrand decomposition is as complex as solving a linear system of equations. Indeedn, a similar procedure can also be applied to the decomposition of eq. (2.3), the main difference being that the resulting relations for µ ij will not only depend on the components of the loop momenta along the physical directions, but also on the orthogonal directions.

Angular integration over the transverse space
As we have explicitly indicated in (2.16), the denominators of Feynman integrals, are completely independent of the transverse components of the four-dimensional loop momenta, namely do not depend on any of variables Θ ⊥ , which are in one-to-one correspondence with {x ji }, j > d . In addition, since the numerator is polynomial in the transverse variables, after the change of variables (2.15), the integrand is mapped into a polynomial in (sine and cosine of) Θ ⊥ , with rational coefficients thed depend on Λ and on the physical directions {x ji }, j ≤ d . Finally we observe that all the integrals over Θ ⊥ are factorized one-dimensional integrals, each being of the type The values of the exponents α and β appearing in eq. (2.22) depend both on the angle θ ij under consideration and on the specific expression of the numerator. Nevertheless, these integrals can be computed once and for all, up to the desired rank, and then re-used in every calculation where they occur. One algorithmic way to perform these integrals consists first in expanding the numerator in terms of Gegenbauer polynomials C (α) n (cos θ), a particular class of orthogonal polynomials over the interval [−1, 1] (see appendix D), and then repeatedly make use of the orthogonality relation they obey, In this way, all integrations over Θ ⊥ , i.e. over all components of the loop momenta orthogonal to the external kinematics, are brought back to a unique integral formula which automatically sets to zero all spurious contributions to the Feynman amplitude. Alternatively, one can show that this angular integration is equivalent to a tensor decomposition of the JHEP08(2016)164 subspace orthogonal to the external legs of the diagram [2]. Consider a topology with n ≤ 4 external legs. A generic term contributing to a -loop integral of such topology has the form (e r · q t ) αr,t . (2.24) In the first factor on the right of the integration measure, we collected the dependence on the variables λ ij and on the components of the loop momenta along the directions of the external momenta, while the second one depends on the transverse components which can be integrated out. Because of the obvious relation the angular integration can also be performed via tensor decomposition, restricted to the d ⊥ -dimensional orthogonal subspace. In particular, this decomposition only depends on the d ⊥ -dimensional projection of the metric tensor and it is independent of the external legs of the diagram, which makes it significantly simpler than a full d-dimensional tensor reduction. This implies that we can easily perform the transverse integration by considering where α i = t α i,t (cfr. with eq. (2.24)) and S is the set of non-equivalent permutations of the Lorentz indexes ν i appearing on the l.h.s.. One can thus solve for the coefficients a σ in the traditional way, i.e. by contracting both sides of the equation with each term on the r.h.s. side and using the identities which allow us to replace the second factor in the product of (2.24) with a combination of variables λ ij . Notice that this combination only depends on the number n of external legs and on the powers of loop momenta appearing in the product of the transverse component, while it is completely independent of the expression of the other factors appearing in the integrand. As well as the explicit angular integration discussed above, this decomposition can be performed for the occurring rank once and for all, and it is independent of the internal details of the topology and of the particular process under consideration.
In the following, we use the integral representation (2.16) and apply the integration procedure described above in the case four-point integrals up to three loops. We refer the reader to appendix A, for the derivation of eq. (2.16) as well as of the explicit expression of the change of variables (2.15). General results for one-and two-loop integrals in all kinematic configurations, including a list of integrals over the transverse directions, are collected in appendices B-C.

Four-point examples
As an example, we consider the four-point topologies depicted in figure 1. Due to momentum conservation, the external momenta {p 1 , p 2 , p 3 , p 4 } span a subspace with dimension d = 3 and, as a consequence, we can build a four-dimensional basis {e α i } containing one single transverse direction e α 4 , Thus, in all the three cases, we can decompose the d-dimensional loop momenta according and λ α i are vectors in the d ⊥ = d − 3 dimensional orthogonal space, . (2.31) In this case, the set of transformations (2.15) is reduced to

JHEP08(2016)164
which expresses the transverse component x 41 of the loop momentum in terms of the single angular variable θ 11 . The numerator of a general Feynman integral corresponding to the box topology can have at most a polynomial dependence on x 41 (and hence on cos θ 11 ), so that the angular integration can always be reduced to the orthogonality relation (2.23). In particular, for the case of a scalar integral with trivial numerator, we obtain . (2.33) Moreover, as we recall in appendix D, odd powers of x 41 can be expressed in terms of (products of) Gegenbauer polynomials with different indices and vanish by orthogonality, so that only even powers of the transverse variable produce non-zero contributions.
As an example, useful for later use, let us consider the integrals . (2.34b) After writing the powers of cos θ 11 in terms of Gegenbauer polynomials, (2.35b) we can use the orthogonality relations (2.23), and obtain In the second equality, we have identified additional powers of λ 11 in the numerator, produced by the integration over the transverse component, with higher-dimensional scalar integrals, as it can be easily checked from the explicit expression of the ddimensional integral (2.33). Results for higher rank numerators can be found in appendix B.
(b) At two loops, the transverse space of the topology shown in figure 1b is described by

JHEP08(2016)164
the variables Λ = {λ 11 , λ 22 , θ 12 } and Θ ⊥ = {θ 11 , θ 22 } and we have . (2.37) In this case, (2.15) reads so that, after the change of variables, any term in the numerator depending on x 41 and x 42 is mapped into a polynomial in (sine and cosine of) Θ ⊥ , with coefficients depending on Λ, which can be easily integrated through the expansion in terms of Gegenbauer polynomials. In this way we find, for the scalar integral, , (2.39) whereas the first non-spurious monomial in x 41 and x 42 amounts to Results for higher rank numerators can be found in appendix C.

Factorized integrals and ladder topologies
The d = d + d ⊥ parametrization (2.16) applies to all Feynman integrals with n ≤ 4, nevetheless, there are special classes of multiloop integrals, associated to factorized and ladder topologies, which allow further simplifications. These integrals are characterized by a set of denominators which are independent of a certain number of transverse orientations λ ij , i.e. on a subset of the angular variables Θ Λ . This implies that, as it can be immediately understood from the properties of the change of variables (2.15) and the integration measure (2.17b), the integration via expansion in Gegenbauer polynomials can be applied, besides to all Θ ⊥ angles, also the angles Θ Λ which do not appear in the denominators. In the following, we discuss the d = d + d ⊥ parametrization in two explicit cases.

Factorized integrals
When the loop corresponding to q α i is factorized, no denominator depends on q i · q j , with j = i. In general, whether a factorized integral originates as a Feynman diagrams or from the algebraic semplification of denominators, the integrand is not necessarily factorized, since the numerator can still depend on the (d−4)-dimensional part of q i ·q j , corresponding to µ ij . Nevertheless, it can be shown that, after integrating over µ ij , by means of the orthogonality relation (2.23), the d = d + d ⊥ parametrization of a factorized integral is given by the product of the d = d + d ⊥ parametrizations of the integrals corresponding to the subtopologies. Remarkably, the transverse space of the factorized graphs can have a different dimension, since, for each of them, it depends on the number of the legs.
As an example, let us consider a bow tie integral of the type shown in figure 2a, which, according to the d = 4 − 2 parametrization (2.6), reads .
Any tensor numerator admits the generic decomposition, so that, if we introduce the change of variable cos φ ≡ µ 12 / √ µ 11 µ 22 , the integral over µ 12 can be reduced to an integral of the type (2.22), which can be evaluated by means of the by-now usual orthogonality relation (2.23), After inserting this result in (2.46), the integral over each loop momentum is completely factorized and, by comparison with the d = 4 − 2 parametrization of one-loop integrals, we can identify, for the non-trivial case α = 2n, .
(2.48) Each term in the brackets admits a d = d +d ⊥ parametrization (2.16). The decomposition can be implemented by working with two different momentum bases, each one containing two vectors orthogonal to the external legs connected to the corresponding loop. In this case, the factorized graph is obtained from the product of two identical subtopologies. In general, if the subtopologies have a different number of extenal legs, then their transverse space is not identical.

Ladder integrals
Starting from a number of loops ≥ 3, ladder topologies correspond to integrals whose denominators depend on a limited number variables λ ij . In these cases, the d = d + d ⊥ parametrization (2.16) reads exactly as in the general case (2.16), but the integration in terms of Gegenbauer polynomials can be extended to the subsets of angles Θ Λ corresponding to the λ ij which do not appear in the denominators. As an example, we consider the three-loop ladder box shown in figure 2b, for which we introduce the same set of transverse variables as for the three-loop diagram of figure 1c, and parametrize the integral exactly as in (2.41). This integral has no propagator depending on both q α 1 and q α 3 , i.e. the denominators are independent of λ 13 and hence of θ 23 , as it can be seen from (2.42). Therefore, the integral over θ 23 is reduced to the form (2.22), and it can be evaluated in the usual way, as In (2.50) the indices α and β are determined by the specific form of the numerator. In the scalar case (α = β = 0), this additional integration returns .

Simplified integrand form
The d = d + d ⊥ parametrization of Feynman integrals and the angular integration over transverse directions can be used in order to decompose scattering amplitudes in terms of a reduced number of scalar integrals without explicitly performing any tensor reduction. In fact, transverse integration can be used ab initio in order to obtain a simplified form of the integrand, free of spurious contributions, which can be more easily reduced to a combination of a minimal set of master integrals, by means of relations like integration-by-parts identities. In particular, as we show in the following example, this procedure is suited for application to helicity amplitudes which, in general, may enjoy better properties than the form factors defined in the usual tensor decomposition. Alternatively, transverse integration can be applied in tandem with algebraic methods, such as integrand decomposition, in order to achieve a step-by-step simplification of the reduction algorithm. The interplay of transverse integration and integrand decomposition will be the object of the next section.
As an example, we consider the double-box contribution to a four-gluon color-ordered helicity amplitude. For this case, we show that the integration over the transverse variables can lead to a simplified representation of the integrand, before considering the application of any reduction algorithm. The topology, in d = 4 − 2 dimensions, is defined by the seven denominators and the four irreducible scalar products where v ⊥ is orthogonal to the external momenta and it can be chosen as Notice that, conversely to the transverse vector e 4 introduced in the general discussion of section 2.4, with this definition v ⊥ is not normalized since, when dealing with realistic processes, it is convenient to use a representation which can be easily expressed in terms of spinor variables without introducing spurious square roots. It is worth observing that, while the first two scalar products in eq. (2.53) live in the physical space defined by the external momenta, the last two lie along the orthogonal direction and will be integrated out using the technique previously discussed. Finally, as it is implied by the definition of the denominators (2.52), all the external momenta p i are taken as outgoing.

JHEP08(2016)164
We consider the helicity amplitude . The double-box contribution to the amplitude is given, in a pure Yang-Mills theory, by the sum of the four diagrams shown in figures 3a-3d, namely a diagram involving only gluons and three diagrams with ghosts circulating in the loop. The calculation can be easily carried out e.g. in Feynman gauge and with an explicit choice of polarization vectors in terms of spinor variables such as We remark, however, that the final result for the on-shell residue, which we will discuss in section 3.5.1, is gauge invariant and thus independent of the previous choices. After inserting the Feynman rules for each diagram, and decomposing the loop momenta as the numerator becomes a function of the coordinates x ij appearing in eq. (2.56) and of the (−2 )-dimensional scalar products µ ij . According to eq. (2.30), the transverse vectors λ i can be identified with The d = d + d ⊥ parametrization of the integrand is simply obtained through the change of variables which, as we have already observed, makes the denominators independent of the transverse components x 41 and x 42 , where Unpon the change of variables, the numerator is given by a sum of 2025 distinct terms in the integration variables (2.59) The terms proportional to the transverse variables x 4i can be integrated out using the results listed in appendix C. Nevertheless, we want to remark that, when using a nontrivial normalization of v ⊥ , the right hand sides of the formulas for I After integrating out the transverse directions, the numerator is reduced to a sum of 773 terms in the variables which, as we will explain in more detailed in the next sections, can be easily expressed in terms of denominators and physical scalar products.

Integrand recurrence relation
In the framework of the integrand reduction method [7-9, 15, 16, 18], the decomposition of dimensionally regulated -loop integrals is rephrased as partial fractioning of the integrand, which is finally written as a sum of residues, i.e. irreducible numerators which cannot be expressed in terms of denominators D i k , sitting over all possible subsets of denominators, For an integral with an arbitrary number n of external legs, the integrand decomposition formula (3.2) can be obtained by observing that both numerator and denominators are polynomials in the components of the loop momenta with respect to some basis, which we collectively label as z = {z 1 , . . . , z ( +9) 2 }. Thus, we can choose a monomial ordering, and build a Gröebner basis G i 1 ···ir (z) of the ideal J i 1 ···ir generated by the set of denominators, being P [z] the ring of polynomials in z. By performing the polynomial division of N i 1 ···ir (z) modulo G i 1 ···ir (z), we obtain the recurrence relation whose iterative application to the integrands corresponding to subtopologies with fewer loop propagators yields to the complete decomposition (3.2). The properties of the ideal J i 1 ···ir [43][44][45][46] allow us to derive an important result concerning the parametric form of the residues corresponding to maximum-cuts. We define maximum-cut a zero-dimensional system of equations which completely constraints the loop variables z. If a maximum-cut admits a finite number of solutions n s , each with multiplicity one, it satisfies the following [15].

JHEP08(2016)164
Theorem 1 (Maximum cut) The residue of a maximum-cut is a polynomial parametrized by n s coefficients, which admits an univariate representation of degree (n s −1).
Depending on the choice of variables z and of the the monomial order, the picture presented in this section can significantly simplify. A particular convenient choice of variables turns out to be the one presented in section 2.2. Indeed, as we observed at the end of that section, we can always express a subset of the components of q α i and λ ij as a combination of denominators by solving linear relations. This set of relations is in turn equivalent to the definition of the denominators themselves. This implies that if we choose the lexicographic monomial order with λ ij ≺ x kl for k ≤ d , the polynomials in the Gröbner bases are linear in the λ ij and in the reducible components of q α i . Thus, the polynomial division can be equivalently performed by applying the aforementioned set of linear relations without explicitly computing the corresponding Gröbner basis.

Divide, integrate and divide
As we have seen in section 2, when dealing with an integral with n ≤ 4 external legs, we can introduce the d = d + d ⊥ parametrization which removes the dependence of the denominators on the transverse components of the loop momenta. Therefore, if we indicate with z the full set of ( + 9)/2 variables where are the components of the loop momenta parallel to the external kinematics, and the four-dimensional orthogonal ones, the denominators are reduced to polynomials in the subset of variables so that the general r denominators integrand has the form For convenience, we make explicit the dependence of the integrand on the component of the loop momenta τ , x ⊥ , and highlight that the denominators do not depend on x ⊥ . This observation suggests an adaptive version of the integrand decomposition algorithm, where the polynomial division is simplified because it involves a reduced set of variables τ , and where the expansion of the residues in terms of Gegenbauer polynomials allows the systematic identification of spurious terms. The novel algorithm is organized in three steps: JHEP08(2016)164 1) Divide: we adopt the lexicographic order λ ij ≺ x for the τ variables, and we divide the numerator N i 1 ...ir (τ , x ⊥ ) modulo the Gröebner basis G i 1 ···ir (τ ) of the ideal J i 1 ···ir (τ ) generated by the denominators, As a consequence of the specific monomial order, the residue ∆ i 1 ...ir depend on the transverse components x ⊥ i , which are left untouched by the polynomial division, as well as on x i , but not on the λ ij that are expressed in terms of denominators and irreducible physical scalar products. The quotient, instead, depends on the full set of loop variables. As mentioned at the end of section 3.1, the Gröbner bases do not need to be explicitly computed, since, with the choice of variables and of the ordering described here, the result of the division can be simply obtained by merely applying the set of linear relations described at the end of section 2.2.
2 In this way, we can integrate over the transverse directions through the expansion of ∆ i 1 ...ir in terms of Gegenbauer polynomials, which sets to zero spurious terms and reduce all non-vanishing contributions to monomials in λ ij , It should be noted that, due to the angular prefactors produced by the integration of the transverse directions, the integrated residue is, in general, a polynomial in τ whose coefficients depend explicitly on the space-time dimension d. The full set of ∆ int i 1 ...ir (τ ), obtained by iterating the polynomial division and the integration over the transverse space on the numerator of each subdiagram, provides a representation of the integrals of (3.2) free of spurious terms.
3) Divide: however, since ∆ int i 1 ...ir (τ ) depends on the same variables as the denominators D i k (τ ), we can perform a further division modulo the Gröebner basis G i 1 ···ir (τ ), and get

Integrate and divide
In the three-step algorithm divide-integrate-divide, outlined in the previous section, the integration over the transverse angles is performed after the integrand reduction, namely after determining the residues. This first option follows the standard integrand reduction procedure, where the spurious monomials are present in the decomposed integrand, although they do not contribute to the integrated amplitude. Alternatively, if the dependence of the numerators on the loop momenta is known, then the integration over the orthogonal angles can be carried out before the reduction. Within this second option, which we can refer to as integrate-divide, after eliminating the orthogonal angles from the integrands, the residues resulting from the polynomial divisions contain only non-spurious monomials.
In order to integrate before the reduction, the dependence of the numerator on the loop momenta should be either known analytically or reconstructed semi-analytically [47,48]. Such situation may indeed occur when the integrands to be reduced are built by automatic generators or they emerge as quotients of the subsequent divisions.

One-loop adaptive integrand decomposition
We hereby apply the adaptive integrand decomposition algorithm in order to determine an alternative parametrization of one-loop residues. As an exceptional property of one-loop integrands, we find that by working with τ variables, all unitarity cuts are reduced to zero-dimensional systems. Moreover, we show that the last step of the algorithm, i.e. the further polynomial division after angular integration over the transverse space, provides an implementation of the dimensional recurrence relations at the integrand level.
One-loop residues in d = 4 − 2 . A general one-loop integral with n external legs is characterized by a set of n denominators,

JHEP08(2016)164
The integrand depends on five variables which, in the standard d = 4 − 2 parametrization, are identified with where x i are the components of the four-dimensional part of the loop momentum with respect to a basis {e α i } of massless vectors [16,25,27,49] (one particular definition of such basis can be found in appendix E.1). The denominators D i k (z) are quadratic in z, whereas, for any renormalizable theory, the most general numerator is a polynomial of the type Higher rank numerators, such as the one appearing in effective theories can be treated in a similar way, along the lines of [15]. The polynomial division of the numerators N i 1 ···in (z) modulo the Gröbner bases G i 1 ...in (z) = {g 1 (z), g 2 (z), . . . g n (z)} returns the universal parametrization of the residues [10,16,18], (3.20) Many of the terms appearing in (3.20) are spurious, i.e. they vanish upon integration. Therefore, we can write the decomposition of an arbitrary one-loop amplitude with n external legs as a linear combinations of master integrals (IRIs), corresponding to the non-spurious terms of the integrand, where, in the second equality, we have identified µ 2 numerators with higher-dimensional integrals [50].
The particular simple form of the five-point residue appearing eq. (3.20) is due to the fact that the quintuple cut D i (z) = · · · = D m (z) = 0 is a maximum-cut which admits a unique solution (n s = 1). Hence, ∆ ijklm is parametrized by a single coefficient, which is conventionally chosen as the one corresponding to the µ 2 numerator.
One-loop residues in d = d + d ⊥ . In d = 4 − 2 dimensions the maximum-cut theorem can only fix the parametric form of the residue of the quintuple cut, since the cut conditions for all lower-point integrands form an underdetermined system for the variables z. However, all these integrands have n ≤ 4 external legs and we can introduce the d = d + d ⊥ parametrization in terms of the variables where x and x ⊥ are defined according to the four-dimensional basis given in appendix E.2. In this way, the n denominators depend on the subset of variables containing exactly n elements. Since, as explained at the end of section 2.2, these variables can be written as combinations of denominators via linear relations, and because the cut D i 1 (τ ) = · · · = D in (τ ) = 0 with n ≤ 4 is a maximum-cut, the corresponding set of cut equations is equivalent to a determined linear system and therefore has a single solution (n s = 1), which constrains all τ variables. This means that we are in the Shape Lemma position and, as implied by the discussion at the end of section 3.1, a Gröbener basis Hence, according to the maximum-cut theorem, the residues of all cuts with n ≤ 4 are constant in τ .
Although it is independent of τ , ∆ i 1 ···in can still depend on the 4 − d four-dimensional transverse variables, which are left unconstrained by the cut conditions. However, the JHEP08(2016)164 parametrization of the residue is terms of x ⊥ is completely fixed by the renormalizability condition, . . . , j 4−d ) : j 1 + j 2 + · · · + j 4−d ≤ n}. Accordingly, from the polynomial division of the numerators modulo G i 1 ...in (τ ), with lexicographic ordering λ 2 ≺ x ⊥ , we find a parametric expression of the residues alternative to (3.20), (3.26) We observe that the two-point integrand with massless external momentum p 2 = 0, whose residue is indicated as ∆ ij | p 2 =0 , is the only exception to (3.25), since the residue depends on the components x 1 parallel to p. In fact, due to the reduced dimension of the transverse space, the denominators depend on three variables τ = {x 1 , x 2 , λ 2 } so that, in this degenerate kinematic configuration, the double cut is not maximum any more.
The residues (3.26) can now be integrated over the transverse directions by means of the orthogonality relation (2.23) for Gegenbauer polynomials, which removes spurious terms and reduce non-vanishing contributions to powers of λ 2 . Hence, by making use of the results collected in appendix B, we obtain the decomposition of a generic one-loop amplitude in terms of IRIs, where, in the second equality, we have identified monomials in λ 2 in the numerator with higher-dimensional integrals. The number of IRIs in which the amplitude is decomposed can be further reduced by observing that, due to the choice of lexicographic ordering (which allows us to express λ 2 as a function of x i ), λ 2 is reducible, i.e. it can be rewritten, modulo a constant remainder, in terms of denominators. Therefore, all higher-dimensional integrals appearing in (3.27b) are reduced to a combination of the corresponding scalar integral in d-dimensions and integrals with fewer denominators. As a consequence, this additional polynomial division can be interpreted an implementation of dimensional recurrence relations at the integrand level. The final decomposition of the amplitude in terms of the minimal set IRIs reads It should be remarked that, although we have used similar a labelling, the coefficients the master integrals appearing (3.28) are different from the ones in (3.27a). Moreover, in (3.28), these coefficients can present an explicit dependence on the space-time dimension, due to the angular prefactors produced by the integration over the transverse variables. We give a summary of the results obtained from the application of the adaptive integrand reduction algorithm at one loop in table 1.

Two-loop adaptive integrand decomposition
In this section, we use the adaptive integrand decomposition algorithm in order to determine the universal parametrization of the residues appearing in the integrand representation (3.2) of the three eight-point topologies shown in figure 4a-4c. The results hereby presented are valid for arbitrary (internal and external) kinematic configuration. At two-loops, we generally deal with r denominators Feynman integrals of the type  Table 1. Residue parametrization for irreducible one-loop topologies. In the first column, τ labels the variables the denominators depend on. ∆ i1 ··· in indicates the residue obtained after the polynomial division of an arbitrary n-rank numerator and ∆ int i1 ··· in the result of its integral over transverse directions. ∆ i1 ··· in corresponds to the minimal residue obtained from a further division of ∆ int i1 ··· in . In the figures, wavy lines indicate massless particles, whereas solid ones stands for arbitrary masses.

JHEP08(2016)164
with r 1 denominators depending on q α 1 only, r 2 denominators depending on q α 2 only and r 12 = r − r 1 − r 2 depending on both loop momenta. The general numerator of a nonfactorized integrand (r 12 = 0) of a renormalizable theory is given by a polynomial of the type where z = {z 1 , . . . , z 11 } labels the full set of loop variables (which will be later specified according to the number of external legs) and J 11 (s 1 , s 2 , s tot ) denotes the vectors of integers  = (j 1 , . . . , j 11 ) satisfying the renormalizability constraints j i + 2(j 9 + j 10 + j 11 ) ≤ s 12 , with s tot = r 1 + r 2 + r 12 − 1. (3.31c) As we have already observed in the one-loop case, the present discussion can be easily extended to the case of higher rank numerators. Depending on the number of external legs, we determine the residue parametrization in two different ways: • For the parent topologies I P 1···11 , I NP1 1···11 and I NP2 1···11 , as well as for all subtopologies with n > 4 four external legs, we use the d = 4 − 2 parametrization. Accordingly, we define z as z = {x 11 , . . . x 41 , x 12 , . . . , x 42 , µ 11 , µ 22 , µ 12 }, (3.32) where {x 1 i , . . . , x 4 i } are the components of the four-dimensional vector q α [4] i with respect to the basis defined in appendix E.1. In these cases, the denominators depend on the full set of variables z and the parametric form of the residues is determined through a single polynomial division of N i 1 ···ir (z) modulo a Gröebner basis G i 1 ...ir (z) of the ideal generated by the denominators. The results obtained for the eightseven-six-and five-point integrands are summarized in tables 2-3. We observe that, according to the maximum-cut theorem, the residues of the master topologies I P 1···11 , I NP1 1···11 and I NP2 1···11 contain one single coefficient, since the zero-dimensional systems D 1 (z) = · · · = D 11 (z) = 0 admit only one solution.
• For any subdiagram with n ≤ 4 external legs, we introduce the d = d + d ⊥ parametrization and define z as definition of the basis). In these cases, the denominators depend on the reduced set of variables τ = {x 1 , x 2 , λ 11 , λ 22 , λ 12 }, τ ⊂ z, (3.34) and we can go through the full adaptive integrand decomposition algorithm described in section 3.2. We refer the reader to the appendix C for the most relevant formulae regarding the d = d + d ⊥ parametrization of two-loop integrals and the integration over transverse variables. It should be noted that, differentlt from the one-loop case, the n-ple cut D i 1 (τ ) = · · · = D ir (τ ) = 0 is, generally, non-maximum, since it does not constrain all variables τ . However, the choice of lexicographic ordering λ ij ≺ x ⊥ i guarantees that the final residues ∆ i 1 ···ir (x i ) appearing in the integrand decomposition formula (3.2) depend on the components of loop momenta parallel to the external kinematics only. All results are summarized in tables 4-7.
Finally, in the case of an integrand factorized into two one-loop diagrams, respectively with n 1 and n 2 independent external legs, we can assume, as discussed in section 2.5.1, the most general numerator to have the form α  1 ,  2 z j 11 11 . . . z j 51 51 z j 12 12 . . . z j 52 52 , (3.35) where z i = {z 1i , . . . , z 5i } labels the set of variables parametrizing q i and J 5 (n i ) is defined by (3.19). In this way, we can introduce the d = d + d ⊥ parametrization independently in any of the two loops, and then proceed with the adaptive integrand decomposition algorithm. As expected, the resulting residues, which are shown in table 8, are simply given by the product of the corresponding one-loop residues collected in table 1.
We would like to mention that the residues ∆ i 1 ···ir (x ) of non planar topologies, which are written in terms of a minimal set of physical components, produce an apparent violation of one of the renormalizability conditions (3.31a)-(3.31b) satisfied by the original numerators. This effect is due to the fact that, when the cut conditions are imposed, the presence of a number r 12 > 1 of denominators depending both on q 1 and q 2 implies the existence of linear relations between the physical components of the two loop momenta. This means that, up to subdiagrams contributions, the residues can always be rewritten in terms of a larger number of variables, in such a way to satisfy all renormalizabilty constraints (3.31a)-(3.31c).

Example: the four-point residue for
As an explicit application of the adaptive integrand decomposition, we go back to the helicity amplitude A 2−loop (p + 1 , p − 2 , p + 3 , p − 4 ) discussed in section 2.6.1 and we compute the residue of the double-box topology. We start from the full numerator, which contains 2025 terms up to rank four with respect to each loop momentum and rank six in total and we determine the residue three steps: 1) Divide: an easy way to perform the first division step of the procedure consists in observing that, since denominators are independent of x 4j , all coordinates x ij with i ≤ 3 JHEP08(2016)164  can be written in terms of differences of denominators and irreducible scalar products by solving a linear system of equation. Moreover, the variables λ ij can also be easily written as combinations of denominators and scalar products by solving simple linear relations. After this manipulations, the numerator on the cut (i.e. imposing D i = 0) is given by a sum of 70 vanishing terms in the components of the four dimensional loop momenta x = {x , x ⊥ }. The expression of the integrand on the cut found agreement with the results of [8].
2) Integrate: after integration over the transverse directions, the numerator acquires again a dependence on λ ij and it is expressed a sum of 39 non-vanishing terms in the variables x = {x , λ ij }, whose coefficient now also depend on the dimensional regulator d.   Table 6. Residue parametrization for irreducible two-point two-loop topologies. Denominators depend on the variables τ = {x 11 , x 12 , λ 11 , λ 22 , λ 12 } in the case of massive external momenta and on τ = {x 11 , x 21 , x 12 , x 22 , λ 11 , λ 22 , λ 12 } in the case of massless one. For every step of the reduction algorithm, we list the number of monomials of each residues and the set of variables appearing in it.

JHEP08(2016)164
In the figures, wavy lines indicate massless particles, whereas solid ones stands for arbitrary masses. Table 7. Residue parametrization for the irreducible one-point two-loop topology. Denominators depend on the variables τ = {λ 11 , λ 22 , λ 12 }. For every step of the reduction algorithm, we list the number of monomials of the residue and the set of variables appearing in it.
3) Divide: using the same relations as in the first step, the λ ij can be expressed in terms of denominators, completing the final division step, after which the numerator of the integrand on the cut is expressed as linear combination of 15 terms depending on the physical directions x left unconstrained by the cut conditions, i.e. on the two irreducible scalar products (q 1 · p 4 ) and (q 2 · p 1 ).
Putting everything together, after factoring out a contribution proportional to the treelevel result by means of some spinor algebra, the gauge invariant decomposition of this cut (i.e. ignoring contributions proportional to denominators) can be written as   Table 8. Residue parametrization for factorized two-loop topologies. For every step of the reduction algorithm, we list the number of monomials of each residues and the set of variables appearing in it. In the figures, wavy lines indicate massless particles, whereas solid ones stands for arbitrary masses.

JHEP08(2016)164
(the dependence on s 12 is unambiguously determined by dimensional analysis) they read

JHEP08(2016)164
where d s is the number of dimensions of the internal gluons, i.e. d s = d in the 't Hooft-Veltman scheme and d s = 4 in the four-dimensional helicity scheme. The integrals appearing on the r.h.s. of eq. (3.36), which only depend on the momenta defined by the external kinematic, can then be reduced to a minimal set of master integrals by means of traditional methods such as integration by parts. It is worth noticing that, while the original integrand had terms up to total rank six in the loop momenta, after the reduction the maximum rank is reduced down to four.

Conclusions
We presented the integrand reduction of dimensionally regulated integrals in the parallel and orthogonal space, where the number of space-time dimensions d = 4−2 is decomposed as d = d + d ⊥ . According to the external topology of each diagram, characterized by the number n of legs, the parallel space is spanned by the external momenta, d = n − 1, while the orthogonal space is spanned by the complementary directions. For diagrams with a number of legs n > 4, the orthogonal space is generated by the regulating directions, d ⊥ = −2 , while for n ≤ 4, it embeds also the four-dimensional complement to the parallel space, namely d ⊥ = 5 − n − 2 .
Owing to the representation of Feynman integrals in parallel and orthogonal space, numerators and denominators of integrands with n ≤ 4 appear to depend on different sets of variables, since the former can depend on transverse angles which are absent from the latter. Therefore, the integration over the subset of transverse variables which do not appear in the denominators can be carried out, before any reduction, simply by employing the orthogonality relation of Gegenbauer polynomials.
Because of the reduced number of variables appearing in the denominators of the diagrams with n ≤ 4 legs, the integrand reduction algorithm, which is based on the multivariate polynomial division, is simplified. In particular, the Gröbner bases generated by the denominators are linear in the variables reduced by the division algorithm, and the multivariate division is reduced to a mere substitution of the solutions of a set of linear equations, which is a consequence of the separation of the physical directions from the extra-dimensional ones. Moreover, the residues, namely the remainders of the polynomial divisions, present a novel, simpler structure. If the integration over the orthogonal directions is performed before the reduction, then the residues contain only monomials that correspond to non-vanishing integrals. On the contrary, if the polynomial division is applied to the complete numerator, the residues will contain also spurious monomials. In the latter case, the integration of the decomposed integrand over the transverse directions by means of Gegenbauer polynomials automatically detects and annihilates the spurious terms.
The outcome of the proposed algorithm is the decomposition of multiloop amplitudes in terms of a set of integrals which, beside the scalar ones, contains tensor structures corresponding to irreducible scalar products between loop momenta and external momenta. These integrals depend on the parallel directions and on the lengths of the transverse JHEP08(2016)164 vectors only. We have shown that the integration over the transverse angles, which can be systematically implemented by using Gegenbauer polynomials, plays an important role in eliminating the superfluous degrees of freedom of multiloop integrals any time that a certain subset of integration variables do not appear in the denominators. We have discussed how, in the case of factorized diagrams and ladder topologies, such integration can be applied, besides to the transverse angles, to a larger number of variables. In addition, we have shown that the integration over the transverse directions leads to integrals which can be subject to additional polynomial divisions, which in some cases correspond to dimension-shifting recurrence relations implemented at the integrand level.
We have revisited the one-loop integrand decomposition and we have shown that it is completely determined by the maximum-cut theorem in different dimensions. We have also considered the complete reduction of two-loop planar and non-planar integrals for arbitrary kinematics, classifying the corresponding residues and identifying the set of integrals contributing to the amplitude. We have discussed how the whole algorithm can be simply extended to higher loops, by giving explicit examples of four-point integrals at three loops.
The dependence of the denominators, hence of the cut-conditions, on a subset of variables determined by the number of legs suggested us to introduce the concept of adaptive cuts. We believe that the idea presented in this article of cutting diagrams in different space-time dimensions, according to the topology under consideration, can be applied, in general, to any unitarity-based algorithm.
It is known that the number of integrals emerging from the integrand reduction is not minimal. In fact, because of the properties of dimensional regularization, the number of integrals appearing in the amplitude decomposition can be further reduced by applying integration-by-parts identities and ensuing relations. We believe that novel reduction algorithms, explicitly built for decomposing integrals that depend on parallel directions and on the lengths of the transverse vectors, may lead to simplified integration-by-parts solving strategies.
As it stands, the proposed variant of a simplified integrand reduction algorithm can be used in tandem, on the one side, with automatic diagram generators and, on the other side, with codes dedicated to the automatic integrals evaluation by means of numerical or semi-analytical routines.

A Spherical coordinates for multiloop integrals
In this appendix we give a derivation the d = d + d ⊥ representation (2.16) of multiloop Feynman integrals, presented in section 2. We provide explicit formulae up to three loops and we show how these results can be extended to higher orders. We start by studying the properties of a set of auxiliary integrals that we will later identify with the integrals over the transverse space.
• For one-loop calculations, it is useful to consider integrals of the type where λ 1 is a vector of an Euclidean space, whose dimension m is first assumed to be an integer and the analytically continued to complex values. We suppose λ 1 to be decomposed with respect to an orthonormal basis {v i } as Regardless of the symmetries of the integrand, we can reparametrize I 1 in terms of spherical coordinates in m dimensions which, being {v i } orthonormal, are defined by the well-known change of variables where √ λ 11 ∈ [0, ∞) and all angles range over the interval [0, π], except for θ m−1 1 ∈ [0, 2π]. Hence, by introducing the differential solid angle in M dimensions such that we can write (A.2) as

JHEP08(2016)164
If the integrand is rotational invariant, i.e. it depends on λ 11 = λ 1 · λ 1 only, we can integrate over all angular variables in such a way to obtain, by specifying (A.6) for M = m However, in general one-loop applications, the integrand can show an explicit dependence on a subset of κ < m − 1 components of λ 1 which, with a suitable definition of the reference frame, can always be chosen to correspond to {a 11 , . . . , a κ1 }. In this way, according to (A.4), the integrand will depend only on Λ = {λ 11 } and Θ ⊥ = {θ 11 , . . . , θ κ1 } while all angles θ i1 with i > κ can be still integrated out by using (A.6) with M = m − κ, • In two-loop computations, we encounter multiple integrals of type where we suppose the two vectors λ i to be decomposed in terms of the same orthonormal basis {v i }, Analogously to the one-loop case, we would like to map all integrals associated to a subset of κ components of each vectors λ i into angular integrals. For I 1 , due to the choice of an orthonormal basis, this mapping was immediately achieved by parametrizing the integral in terms of spherical coordinates. In this case, there is an additional direction, corresponding to λ 12 = λ 1 · λ 2 , we need to trace back after the change of coordinates is performed, since the integrand will generally depend on it. The simultaneous factorization of the integral over the relative orientation λ 12 and over all relevant components of the two vectors can be obtained by expressing λ 2 into a new orthonormal basis {e i }, containing the vector e 1 ∝ λ 1 . From (A.12a) we see that, indeed, the set of vectors

JHEP08(2016)164
is a basis, although it is not an orthogonal one. Nevertheless, we can apply the Gram-Schimdt algorithm to pass from the arbitrary basis {v i } to an orthonormal one {e i }, given by (v k · e j )e j , k = 1. (A.14) By construction, the first vector of the new basis exactly corresponds to the direction of λ 1 . Applying the change of basis to (A.12b), we get where the coefficients {b i2 } are related to the components of both λ 1 and λ 2 with respect to {v i } by Since both λ 1 and λ 2 are now decomposed in two different but still orthonormal basis, we can introduce the change of variables .17) and express the integral I 2 into spherical coordinates as

JHEP08(2016)164
In addition, with some more algebra, we can express back the components of λ 2 with respect to {v i } in terms of the angular variables, cos θ 12 cos θ 11 + cos θ 22 sin θ 11 sin θ 12 (A.20) In this way, the integral over each component a i1 / ∈ {a m−1 1 , a m 1 } of λ 1 is mapped into the integral over the angular variable θ i1 whereas each component a i2 / ∈ {a m−1 2 , a m 2 } of λ 2 can be expressed in terms of the angles θ j1 with j ≤ i and θ j2 with j ≤ i + 1. Therefore, if we are dealing with and integrand depending on κ < m − 1 components of both vectors, which we can always choose to correspond to {a 11 , · · · , a κ1 } and {a 12 , · · · , a κ2 }, we can integrate out all angular variables θ i1 , j > κ and θ i2 , j > κ+1. Hence, if we define we can rewrite I 2 as • For three-loop applications, we consider integrals of the type and, as usual, we assume the vectors λ i to be initially decomposed in terms of the same orthonormal basis {v i }, When moving to spherical coordinates, we want to keep trace of the three relative orientations together with the usual subset of κ components of each λ i . The proper change of variables can be reach in two steps:

JHEP08(2016)164
1) First we express the vectors λ 2 and λ 3 in terms of the basis {e i } defined by eq. (A.14), which contains the vector e 1 ∝ λ 1 , where, similarly to (A.16), {b i2 } and {b i3 } are defined in terms of the components with respect to the basis {v i } as 2) Then we use the fact that the vectors form a (non-orthogonal) basis which can be orthogonalized by applying the Gram-Schmidt algorithm in such a way to obtain an orthonormal basis {f i }, whose first element is f 1 ∝ λ 2 , and decompose λ 3 as

JHEP08(2016)164
Eqs. (A.25), (A.27a) and (A.31) give us a decomposition of the three vectors λ i in terms of three different but still orthonormal basis. Hence, we can introduce spherical coordinates (A.33) and rewrite I 3 as In particular, one can verify that, as in all previous cases, each integral over a i1 / ∈ {a m−1 1 , a m 1 } is mapped into the integral over the angular variable θ i1 and, as we have seen for I 2 , each component a i2 / ∈ {a m−1 2 , a m 2 } can be expressed in terms of the angles θ j1 with j ≤ i and θ j2 with j ≤ i + 1. Moreover, each a i3 / ∈ {a m−1 3 , a m 3 } turns out to be function of the angles θ j1 with j ≤ i, θ j2 with j ≤ i + 1 and θ j3 with j ≤ i + 2. Therefore, if we are dealing with and integrand depending on κ < m − 1 components of all λ i , which can be always chosen to correspond to the first κ ones, we can integrate out all angular variables θ i1 , j > κ, θ i2 , j > κ + 1 and θ i3 , j > κ + 2 and obtain dcos θ 12 dcos θ 13 dcos θ 23 (sin θ 12 ) n−3 (sin θ 12 ) m−3 (sin θ 23 ) n−4 ,

JHEP08(2016)164
• It is now clear how integrals I involving a number of vectors λ i , 38) can be treated in order to define a change of variable which maps a subset of κ components of each vector as well as their ( − 1) relative directions λ ij into angular variables. Starting from the decomposition of all vectors in terms of a single orthonormal basis, one can define, by recursively applying the Gram-Schimdt algorithm, − 1 auxiliary orthonormal basis carrying information both on κ ≤ m − 1 directions of the original basis and on the relative orientations λ ij . After all vectors have been decomposed into the proper orthonormal basis, we can introduce m-dimensional polar coordinates and, by inverting the nested chain of transformations, we can obtain the expression of the components of all λ i with respect to {v i } in terms of the angular variables. The final transformation has the form where Θ Λ and Θ ⊥ label the sets of angular variables Therefore, if the integrand I only depends on κ components of each λ i , all angles θ ij , i ≥ j + κ can be integrated out, producing We can now go back to an arbitrary loop integral with n ≤ 4 external legs and, after introducing the q α i = q α i + λ α i parametrization of the loop momenta, we can rewrite eq. (A.1) as

JHEP08(2016)164
where we have explicitly indicated that the denominators depend on the d -dimensional momenta q i and on the scalar products λ ij between the transverse vectors living in d ⊥ dimensions. We now observe that the numerator can additionally depend only on the four-dimensional components of each λ α i , Therefore, the integral over the transverse vectors λ α corresponds to a d ⊥ -dimensional integral of the type I κ with κ = 4 − d so that, by substituting (A.41) in (A.43), we obtain , which reproduces (2.16).
(C.29) Moreover, in general we have

D Gegenbauer polynomials
In this appendix we recall the most relevant properties of Gegenbauer polynomials. Gegenbauer polynomials C

E Four-dimensional basis
In this appendix we provide the explicit definitions of the four-dimensional basis {e α i } used throughout the text to decompose the four-dimensional part of the loop momenta q [4] i , q α [4] i = p α 0 i + x 1i e α 1 + x 2i e α 2 + x 3i e α 3 + x 4i e α 4 . (E.1) In the following, for any pair of massless vectors q α 1 and q α 2 , we denote by ε α q 1 ,q 2 the spinor product In the d = 4 − 2 parametrization of Feynman integrals we choose, independently from the number of external legs, a basis of massless vectors {e α i } defined in terms of two adjacent external momenta p 1 and p 2 by e α 1 = 1 1 − r 1 r 2 (p α 1 − r 1 p α 2 ), e α 2 = 1 1 − r 1 r 2 (p α 2 − r 2 p α 1 ), e α 3 = ε α e 1 ,e 2 , e α 4 = ε α e 2 ,e 1 , (E.3) where In the case of two-point integrals, p 1 corresponds to the external momentum and p 2 is an arbitrary massless vector. In the case of one-point integrals, both p 1 and p 2 are chosen to be arbitrary massless vectors.