Stable one-dimensional integral representations of one-loop N-point functions in the general massive case. I — Three point functions

In this article we provide representations for the one-loop three point functions in 4 and 6 dimensions in the general case with complex masses. The latter are part of the GOLEM library used for the computation of one-loop multileg amplitudes. These representations are one-dimensional integrals designed to be free of instabilites induced by inverse powers of Gram determinants, therefore suitable for stable numerical implementations.


JHEP11(2013)154 1 Introduction
The Golem project [1] initially aimed at automatically computing one loop corrections to QCD processes using Feynman diagrams techniques whereby 1) each diagram was written as form factors times Lorentz structures 2) each form factor was decomposed on a particular redundant set of basic integrals. Indeed when the form factors are reduced down to a basis of scalar integrals only, negative powers of Gram determinants, generically noted det(G) below, show up in separate coefficients of the decomposition. These det(G), albeit spurious, are sources of troublesome numerical instabilities whenever they become small. The set of basic integrals used in the Golem approach is such that all coefficients of the decomposition of any form factor on this set are free of negative powers of det(G). Let aside trivial one-and two-point functions, the Golem library of basic functions is instead made of a redundant set involving the functions I n 3 (j 1 , · · · , j 3 ), I n+2 3 (j 1 ), I n+2 4 (j 1 , · · · , j 3 ) and I n+4 4 (j 1 ). Here the lower indices indicate the number of external legs, the upper indices stand for the dimension of space-time, and the arguments j 1 , · · · , j i labels i Feynman parameters in the numerator of the corresponding integrand. The strategy is the following. In the phase space regions where det(G) are not troublesome, the extra elements of the Golem set are decomposed on a scalar basis and computed "analytically" in terms of logarithms and dilogarithms. In the phase space region where det(G) vanishes these extra Golem elements are instead used as irreducible building blocks explicitly free of Gram determinant and provided as one-dimensional integral representations computed "numerically".
Much faster and more efficient methods than those relying on Feynman diagrams techniques have been developed, e.g. based on unitarity cuts of transition amplitudes and not individual Feynman diagrams, and/or processing the decompositions at the level of the integrands [2][3][4][5][6][7][8]. Yet these methods still amount to a decomposition onto a set of basic integrals. In this respect the stand-alone relevance of the Golem library of basic functions, initially developed as a part of the Golem approach, remains. Furthermore the decompositions obtained by these new methods project onto a basis of scalar integrals and thus are still submitted to numerical instabilities caused by det(G). The issue of numerical instability is then addressed in various ways ranging from smoothing numerical interpolations over the regions of instabilities [9] to more involved rescue solutions [10,11]. In [12] the solution adopted is to provide a rescue alternative relying on the Golem decomposition to compute the amplitude in the troublesome kinematic configurations. The Golem library [13], initially designed for QCD, did not include basic functions with internal masses yet provided a convenient way of handling infrared and collinear singularities inherent in the massless case. Its completion with the cases involving internal masses, possibly complex, extends its range of use [14]. This completion shall supply the functions I n 3 (j 1 , · · · , j 3 ), I n+2 3 (j 1 ), I n+2 4 (j 1 , · · · , j 3 ) and I n+4 4 (j 1 ) in the massive cases in a numerically stable with respect to det(G) issues.
To handle det(G) issues, we advocate the use of one-dimensional integral representations rather than relying on Taylor expansions in powers of det(G). The latter may be thought a priori better both in terms of CPU time and accuracy, however the order up to which the expansion shall be pushed may happen to be rather large. a fixed large number of terms, hopefully large enough in all practical cases, be computed, it is not easy to assess a priori the optimal order required to reach a given accuracy. Actually this assessment would demand a quantitative estimate of the remainder as a function of the order of truncation, which, as with the Taylor expansion with Laplace remainder, namely requires the computation of an integral! Originally, we proposed the antipodal option of computing numerically the two-or three-dimensional Feynman integral defining respectively the three-and four point functions, more precisely hypercontour deformations thereof [1] that would be numerically more stable. Yet the computation of these multiple integrals was both slow and not very precise. It is far more efficient both in terms of CPU time and accuracy to evaluate a one-dimensional integral representation, insofar as one is able to find such a representation. In the case without internal masses, we indeed found such a representation.
The issue which we address here is the extension of this approach of one-dimensional integral representations for our set of basic integrals in the most general case, i.e. with internal complex masses. In this article we treat the case of the three point function. The case of four point functions is more involved therefore it will be elaborated separately in a companion article. We follow the approach developed by t'Hooft and Veltman in ref. [15]. In a subsequent third article, we will present an alternative approach providing integral representations for both three and four point functions equivalent to the one presented here yet with a number of new features and advantages. The present article is organized as follows. Section 2 sketches the derivation of the three point function leading to our integral representation. Section 3 treats the case when det(G) vanishes whereas the determinant of the kinematic matrix S remains non vanishing. Section 4 elaborates on the more tricky case when both det(G) and the det(S) vanish. The main body of the text presents the general arguments whereas the various technical details supporting the latter are gathered in appendices, to make the reading of this article more fluent.

Outline of the derivation
A generic three point function can be represented by the diagram of figure 1.
Each internal line with momentum q i stands for the propagator of a particle of mass m i . We define the kinematic matrix S, which encodes all the information on the kinematics JHEP11(2013)154 associated to this diagram by: The squares of differences of two internal momenta can be written in terms of the internal masses m i and the external invariants s i = p 2 i so that S reads:

JHEP11(2013)154
Gram determinant associated with the kinematic matrix S. Specifying for example a = 3, I 4 3 reads: The 2 × 2 Gram matrix G (3) and the column two-vector V (3) are explicitly given by: We then define and we get: 2 (2.12) (2.14) Eq. (2.12) is the starting point of the computation of the three point integral in ref. [15], cf. their eq (5.2). We keep the same notations as those of ref. [15] for the different quantities, and we closely follow the strategy of ref. [15] for the first integration. We only sketch these stages. An alternative strategy may be proposed which leads to the sought integral representations in a faster, more straightforward and more transparent way for three point functions, and which can be elaborated for four point functions as well thereby providing a number of interesting features. This alternative will be presented in a subsequent publication. 2 The argument of the logarithm appearing in eq. (2.13) shall be understood to contain an implicit arbitrary factor 1/M 2 with dimension -2 in order to make the argument of this logarithm dimensionless. This arbitrary M 2 dependence is cancelled by the corresponding one in the ln(M 2 /µ 2 ) involved in the term I div 3 , where µ 2 is the dimension two parameter introduced by the dimensional regularization of the ultra violet divergence subtracted in I div 3 . In practice the kinematic matrix S is rescaled from the start by its entry of largest absolute value, and so is the Gram matrix G (a) , which thereby become both dimensionless. This amounts to specifying M 2 to be this normalization parameter.

JHEP11(2013)154
The integration variable y is first shifted according to y = y ′ + α x, the parameter α is being chosen such that b α 2 + c α + a = 0 (2.15) in order that the quadratic form of x, y in the integrands of eqs. (2.12), (2.13) become linear in x. Note that the discriminant ∆ α of eq. (2.15) is minus the Gram determinant det(G). For all kinematical configurations p 1 , p 2 , p 3 = −p 1 − p 2 involved in one-loop calculations of elementary processes of interest for collider physics, det(G) is non-positive. 3 The roots α ± of the polynomial (2.15) are thus real in all relevant cases. We split the integral over y ′ and reverse the order of integrations: Since the integrand seen as a function of x and y ′ in eq. (2.16) is now linear in x the integration on x is made straightforward. For I 4 3 eq. (2.16) involves two integrals of the form where A and B are functions of y ′ and x min (y ′ ) = y ′ /(1 − α) and x min (y ′ ) = y ′ /(−α) respectively. As can be traced back to eqs.
(2.12), (2.13) has a negative imaginary part, this holds true also for complex internal masses. Therefore the numerator and denominator in the argument of the logarithm in eq. (2.17) both have a negative imaginary part, thus the logarithm in eq. (2.17) can be harmlessly split in two terms: It is convenient to add and subtract a term ln(C) in the right hand side (r.h.s.) of eq. (2.18), and split the latter into a sum of two terms lead to integrals over y ′ which are individually well defined and may be safely handled on their own. A similar treatment may be done for I 6 3 adding and subtracting a term C ln(C). We note that thus we choose In this way the integration over x yields four terms. By means of an appropriate change of variable, two of them may be further recombined so that each of the integrals I 4 3 and I 6 3 can be written as the sum of three terms. We call these terms "sector integrals" labelled I (j) , j = 1, 2, 3, they may be put in the following form. For I 4 3 we get: with the sector integrals I 4 3 (j) of the form The coefficients D (j) , · · · , K (j) (α) being provided by the following table; the dependence of the K (j) (α) on α is made explicit for further convenience.

JHEP11(2013)154
The values of the integrals I 4 3 and I 6 3 do not depend on the particular root α = α ± of eq. (2.15) chosen to perform the first integration leading to eqs. (2.23), (2.26). As in ref. [15], either of the two α roots, say α + , may be used to further compute the remaining single integrals in closed form in terms of logarithms and dilogarithms. A symmetrization over α ± would generate an unnecessary doubling of dilogarithms in the closed form that would be prejudicial regarding CPU time in practice. However the discussion of the behaviours of these integrals when det(G) → 0 is made somewhat obscure once one particular choice is made, and for this purpose it is on the contrary more enlightening to symmetrize expressions (2.23) and (2.26) over α ± , especially in the perspective of providing one dimensional integral representations free of det(G) instabilities. The α dependence comes only from the factors K (j) (α)/(D (j) z + E (j) ), not from the arguments of the logarithms in numerators. Each of the sector integrals in the decomposition of I 4 3 , respectively I 6 3 , has an explicit α dependence of the type: where L stands for the α-independent numerators in the integrands of the sector integrals I 4,6 3 (j) , and we omit the superscript (j) labelling the sector for simplicity. Symmetrizing over α ± we get: Let us introduce the following quantities: Here are the explicit forms corresponding to the different sector integrals.

JHEP11(2013)154
For sector (2), similarly, K(α) = − α, A = c z + e and C = 2 a z + d, so that: z + e and C = (c + 2 a) z + d, so that: where the coefficients b j are defined by They are such that They were introduced in the GOLEM reduction algorithm [1], and the second degree polynomials g (j) (z) are given by The polynomials g (j) (z) are namely those appearing in the integral representations of the two-point functions corresponding to the three possible pinchings of one propagator in the triangle diagram of figure 1. In what follows we parametrize the g (j) (z) generically as in order to formally handle them all at once when concerned with the zeroes of g (j) (z) + 1/(2B) further below. Let us note that the discriminant ∆ j of the second degree polynomial g (j) (z), defined by

JHEP11(2013)154
turns out to be equal to minus the determinant of the reduced kinematic matrix S {j} . This reduced kinematic matrix corresponds to the pinching of the propagator j in the triangle of figure 1, and is obtained from the matrix S by suppressing line and column j.
Correlatively γ ′′ (j) can be seen as half the reduced Gram determinant associated with the reduced kinematic matrix S {j} . Equation (2.23) can thus be written: Likewise for eq. (2.26): In eqs. (2.38), (2.39), the contour prescription inherited from (− z T Sz − iλ) in eqs. (2.3), (2.4) is implicit: the logarithmic terms ln g (j) (z) in the numerators stand for ln g (j) (z) − iλ . Let us remind that the terms ln(−1/(2B)) in the numerators have been introduced in order that the zeroes z ± (j) of the denominators (2B g (j) (z) + 1) be fictitious poles in each of the sector integrals in any case i.e. the residues vanish: hence ln(−1/(2B)) stands for ln(−1/(2B) − iλ) as well, furthermore no contour prescription around the z ± (j) is needed. Equations (2.38) and (2.39) are appealing candidates for the integral representations which we seek. Let us examine them more closely when det(G) → 0. We shall distinguish two cases: the generic case when det(G) → 0 whereas det(S) remains non vanishing, and the specific case det(G) → 0 and det(S) → 0 simultaneously which deserves a dedicated treatment. Let us subsequently examine these two cases.

JHEP11(2013)154
Identity (3.1) is derived in appendix A. It is an example of the so-called Jacobi identities for determinant ratios, relating the determinant of a matrix and related cofactors i.e. determinants of reduced matrices 4 [16][17][18]. Similar identities may be met in the treatment of the four-point function. The zeroes z ± (j) of g (j) (z) + 1/(2B) are given by is only of degree one, and its unique root z 0 (j) given by is again dragged away from [0, 1] towards ∞ when det(G) → 0. In either case, as soon as det(G) becomes small enough each of the integrals is analytically well defined and numerically safe, and furthermore the following identity holds: so that the contributions ∝ ln(−1/(2B) − iλ) sum up to zero in I 4 3 as well as in I 6 3 . In this respect, let us stress that the contributions ∝ ln(−1/(2B)−iλ) are fictitious from the start. They were introduced through eq. (2.19) with the custodial concern of separately handling integrals -the sector integrals -with integrands free of poles within the integration domain namely when either of z ± j is inside [0, 1]. When z ± j are both outside [0, 1] the introduction of the ln(−1/(2B) − iλ) terms is irrelevant and indeed identity (3.6) allows to drop them explicitly from eqs. (2.38), (2.39). The following integrals separately calculated would be submitted to numerical instabilities. Besides, if case some g (j) (z) vanishes at someẑ (j) inside [0, 1], a possible numerical improvement of the integral representation consists in deforming the integration contour in the complex z plane, to skirt the vicinity of the integrable singularity atẑ (j) , so as to prevent the integrand from becoming large and avoid cancellation of large contributions, according to a one-dimensional version 5 of the multidimensional deformation described in section 7 of ref. [1].

det(G) → 0 and det(S) → 0 simultaneously
This case is more tricky and deserves further discussion. Indeed, when det(S) = 0 and det(G) = 0, eq. (2.33) defining the parameters b j as 3 k=1 S −1 jk is no longer valid as S −1 is not defined, and the parameter B = det(S)/ det(G) is an indeterminate quantity of the type 0/0, likewise the z ± (j) are indeterminate quantities not manifestly driven away from the interval [0, 1].
In this subsection we will first characterize the specific kinematics which leads to such a case. Then we will consider kinematic configurations close to the so-called specific ones above, such that det(G) and det(S) are simultaneously small but non vanishing and we will study how I 4 3 and I 6 3 behave when det(S) and det(G) both go to zero. Anticipating on the result, we rewrite both for I 4 3 and I 6 3 the corresponding sums One of the coefficients b j , say b 3 will be shown to have a finite limit whereas b 1 and b 2 diverge towards infinity in a concomitant way such that their sum b 1 + b 2 has a finite limit. Furthermore, the difference I 3 (1) − I 3 (2) will be shown to tend to zero so that the product ) has a finite limit. A well-defined expression is thus achieved in the double limit det(S) → 0, det(G) → 0 although some of the ingredients are separately ill-defined in the limit considered. We will conclude this subsection with a comment in relation with the behaviour of the GOLEM reduction formalism in this case.
4.1 Characterization of the specific kinematics det(G) = 0, det(S) = 0 The quantities det(S) and det(G) are polynomials in the kinematical invariants. Hereafter we propose a presentation which partly linearizes the resolution of the non linear system det(G) = 0, det(S) = 0. This approach, applied here to N = 3, extends to other N , e.g. N = 4. The determinant det(S) can be written (see appendix A): In broad outline, the contour deformation is contained inside the band 0 ≤ Re(z) ≤ 1. It departs from the real axis at 0 with an acute angle and likewise ends at 1 in such a way that Im(g (j) (z)) is kept negative along the deformed contour so that the latter does not cross any cut of ln(g (j) (z) − i λ). In the case at hand this type of contour never embraces any of z ± j as soon as the latter are outside [0, 1], thus no subtraction of illegitimate pole residue contribution at z ± j has to be cared about.

JHEP11(2013)154
where G (a) is the matrix of cofactors 6 of G (a) ; the superscript " T " refers to matrix transposition. The system det(G) = 0, det(S) = 0 is thus equivalent to the system det(G) = 0, the matrices G (a) andG (a) are simultaneously diagonalizable. When det(G) = 0 and G (a) has rank 7 (N − 2) -namely 1 in the N = 3 case at hand -the only eigendirectionn (a) of G (a) associated to the eigenvalue zero is concomitantly the only eigendirection ofG (a) associated to the only non vanishing eigenvalueg ofG (a) . Forn (a) properly normalized, is thus equivalent to the following linear one: Let us now consider the condition det(G) = 0. A detailed discussion is provided in appendix B, we only summarize it here for the N = 3 case at hand. A vanishing det(G) happens (i) either when the external momenta p 1,2,3 are proportional to each other (ii) or when there exists a non vanishing linear combination of the external momenta which is lightlike and orthogonal to all of them [19]. Possibility (i) corresponds to degenerate kinematic configurations irrelevant for next-to-leading order (NLO) calculations of collider processes. Let us focus on possibility (ii) further assuming any subset of two of the three external momenta to be linearly independent. To fix the ideas, let us consider 8 p 1 and p 3 . If one of them say p 1 is lightlike it is namely (proportional to) the lightlike combination sought, whereas p 3 shall be spacelike, p 2 = −p 1 − p 3 is spacelike as well and s 2 = s 3 . If neither p 1 nor p 3 are lightlike, both shall be spacelike with s 1 = s 3 , and p 2 is (proportional to) the lightlike combination of p 1 and p 3 . Actually, configurations of type (ii) with p 3 , p 1 and p 2 linearly independent and all spacelike can also lead to a vanishing det(G), yet such configurations are not relevant for collider processes at NLO, 9 we thus discard them. Let us assume p 2 lightlike and orthogonal to p 3 and p 1 both spacelike: s 1 = s 3 ≡ s + < 0, s 2 = (p 2 · p 3,1 ) = 0, so that (p 1 · p 3 ) = − s + . We single out line and column 3 of S whose corresponding G (3) reads:

JHEP11(2013)154
The normalized eigenvectorn (3) associated with the eigenvalue zero is (up to a sign): With eqs. (2.11) and (4.6), condition (4.4) imposes the following restriction on the internal masses: 4.2 Behaviour of I 4 3 and I 6 3 when det(G) → 0, det(S) → 0 Let us assume condition (4.7) and parametrize the departure from the 'critical kinematics' using s 2 , s − ≡ (s 1 − s 3 )/2 and s + ≡ (s 1 + s 3 )/2. The determinants read: where λ is the Källen symmetric function of s + , m 2 , m 2 3 given by: The region where det(G) and det(S) are concomitantly small corresponds to |s 2 |, |s − | both small compared with the other kinematical invariants, with |s 2 /s + | and (s − /s + ) 2 of the same order, so that det(G) and det(S) can be approximated by: In order to understand in more detail the origin of the diverging contributions to the coefficients b j , the matrix S may be decomposed as follows: Let us address the real mass case first; we will briefly comment at the end of this subsubsection on how the study shall be -only slightly -modified in the complex mass case. Decomposition (4.13) corresponds to the usual diagonalization of S: the v (j) and σ j , j = 1, 2, 3 are real eigenvectors orthonormalized in the euclidean sense (v T (j) · v (k) ) = δ jk and the corresponding real eigenvalues respectively. The labelling of eigenvectors and values is chosen such that σ i=3 explicitly given by 4.14)

JHEP11(2013)154
is the eigenvalue which vanishes when s 2 and s − both vanish, whereas the two others remain finite in this limit. Introducing the column vector b = S −1 · e and the quantity B = 3 More explicit algebraic expressions of the various ingredients in the relevant regime are gathered in appendices C and D for convenience.
Let us however notice that we are taking a double limit. Properly speaking, the limits of each of these three terms in eq. (4.1) which are separately well-defined are directional limits s − → 0, s 2 → 0 in the {s 2 − , s 2 } plane keeping the ratio t = s 2 /s 2 − fixed, i.e. these directional limits are functions of t. However, the limit of the sum of these three terms in eq. (4.1) is indeed independent of t. This can be easily checked numerically, this can also be proven analytically although this is somewhat cumbersome; a proof is presented in appendix F. The ground reason why this property holds is further understood as follows: were the limit of the sum a directional one, it would imply that the three-point function would be a singular i.e. non analytical function of the kinematical invariants at such configurations. However

JHEP11(2013)154
the kinematic singularities are characterized by the so-called Landau conditions 10 [19,20] (see also [17]). For one loop diagrams, these conditions require not only that det(S) = 0, but also that the eigenvectors associated with the vanishing eigenvalue of S shall have only non negative components and that their sum be strictly positive. By contrast, in the case at hand, the eigenvector v (3) in the limit where σ 3 = 0 is, cf. appendix C: The vanishing det(S) in the present case is therefore not related to a kinematic singularity: the three-point function is regular in the limit considered, in particular this limit shall be uniform i.e. not directional.

Extension to the complex mass case
The above study was stressed to hold, strictly speaking, for real masses. Actually, it can be extended to the complex mass case with only slight modifications. Indeed in the complex mass case, the symmetric matrix S albeit complex admits a decomposition formally identical eq. (4.13): which now reflects the so-called Takagi factorization S = U · Σ · U T in terms of a real non negative diagonal matrix Σ and a unitary matrix U , instead of a standard diagonalization. The diagonal elements σ ′ j of Σ are the square roots of the eigenvalues of the hermician matrix S S † , whereas the columns u (i) of U are corresponding eigenvectors 11 of S S † . The corresponding Takagi factorization of S −1 for S invertible reads S −1 = U * · Σ −1 · U † i.e. in the tensor product notation: Identity (4.23) provides the equations for b and B which modify eqs. (4.16), (4.17) in the complex mass case. A study quite similar to the real mass case then follows 12 and the same 10 For general parametric integrals the Landau conditions provide only necessary conditions to face singularities, either of pinched or end-point type. However for Feynman integrals, Coleman and Norton [21] proved these conditions to be sufficient. 11 Let us note by passing that, unlike with standard diagonalization, the phases of the vectors u (i) involved in the Takagi factorization shall be adjusted modulo π in order to fulfill the condition (u T (j) · S · u * (k) ) = σ ′ j (u T (j) · u (k) ), because the decomposition involves the transpose of U not its hermician conjugate. 12 Technically speaking, the determination of the singular values σ ′ j and corresponding vectors u (j) may seem somewhat awkward given the algebraically more complicated form of the matrix elements of S S † . Actually we are interested in a practical case where the imaginary parts of the masses -i.e widths of unstable particles in internal lines -are much smaller than the real parts. Therefore, splitting S in real and imaginary parts S = SR − i SI , and writing S S † = S 2 R + ∆, with ∆ = i[SR, SI ] + S 2 I , the square roots σ ′ j of eigenvalues of S S † and the corresponding eigenvectors u (j) can be expanded in integer powers of matrix elements of SI , as perturbative deformations of the eigenvalues σ k and eigenvectors v (k) of SR i.e. the spectral features of the real mass case, by a straightforward application of the formalism of time-independent perturbation theory in Quantum Mechanics.

A comment on the GOLEM reduction formalism when det(S) = 0
Let us end this subsection with a comment on the applicability of the GOLEM reduction formalism [1] to configurations such that det(S) = 0. The equation S · b = e with S not invertible can still be solved e.g. introducing the so-called Moore-Penrose Pseudo-Inverse [22] T 0 of S(σ 3 = 0) given in the real mass case 13 by: provided the following compatibility condition to be satisfied: x arbitrary scalar and b 0 ≡ T 0 · e, leading to B 0 = e T · T 0 · e. The GOLEM formalism thus applies also, using b = b 0 and B = B 0 , when standing precisely at the peculiar configurations. Yet slightly away from these peculiar configurations the GOLEM ingredients defined by b = S −1 · e separately show discontinuities 15 w.r.t. those given by b = b 0 precisely at the peculiar configurations; this discontinuity comes from the contribution to b coming from the (divergent) component along v (3) , which have no counterpart in b 0 . Notwithstanding, these individual discontinuities are artefacts in the sense that they cancel out in the reduction formula when put altogether, as discussed above. 13 A similar discussion holds in the complex mass case as well with similar expressions cf. the previous paragraph.
14 The arbitrary component of b along v (3) |σ 3 =0 is irrelevant for any practical purpose. Indeed, the condition (4.27) makes the contribution to b3, hence to B, coming from this component vanish, whereas the finite contribution to b1 − b2 from this component is weighted by the vanishing I 3 (1) − I 3 (2) in the decomposition (4.1). 15 More precisely b0 is equal to the directional limit s2 → 0 of b along the direction s− = 0 i.e. t = ∞. The discontinuities are meant for all other, finite t directions.

JHEP11(2013)154 5 Summary and outlook
In this article, we provided a representation of one-loop, 3-point functions in 4 and 6 dimensions in the form of one dimensional representations in the general case with complex masses. These one-dimensional integral representation have the virtue to avoid the appearance of factitious negative powers of Gram determinants, and are therefore numerically stable and remain rather/relatively fast to compute numerically. We addressed the two cases at hand separately: the generic case when det(G) becomes small whereas det(S) remains finite, and the trickier specific case when both det(G) and det(S) become concomitantly small. Here we presented the "existence proof" for scalar integrals, but the method applies to tensor integrals as well, i.e. loop integrals involving integer powers of Feynman parameters in the denominators of their integrands.
A forthcoming article will continue the present one by the similar treatment of oneloop 4-point functions. The latter proved to be quite more involved than the 3-point case, we thus preferred to split it from the present article. These two will be supplemented by a dedicated treatment of the specific mixed case involving both finite (complex) masses, and some zero masses triggering infrared issues. In the meantime we also found an alternative approach leading to a derivation of integral representation which is perhaps simpler and also makes the algebraic nature of the ingredients involved more transparently related to the GOLEM reduction algorithm both for 3-point and 4-point functions, this approach will be presented in a separate article. Last, this approach will be fully implemented in the next version of the GOLEM95 library in Fortran 95. We will provide various numerical tests of numerical stability at this occasion.

A Useful algebraic identities among determinants
Equation (3.1) used in section 3 relates various ingredients of the reduction formula involving the one loop three point function. Similar identities can be found and used in the case of four point functions. These properties can be traced back to general algebraic identities between the determinant of a square matrix and minors of this matrix, referred to as "Jacobi identities for determinant ratios" [16][17][18]. This appendix reminds these general identities and specifies them to the case useful for the present work. Beforehand we remind a few properties useful in this respect.

A.1 Preliminaries
Let us first recall a few useful properties which we state in general for arbitrary N not just N = 3. Consider the kinematic N × N matrix S associated with a given one-loop N -point diagram generalizing eq. (2.2) for any N . We single out the line and column a, and consider the corresponding (N − 1) × (N − 1) Gram matrix G (a) associated to S, generalizing (2.7) and the (N − 1)-column vector V This defines the N × N matrix S (N ) given by: as: where the matrix G (N ) is the matrix of cofactors. Thus, using eq. (A.5):

A.2 The identity (3.1) and Jacobi identities for determinants ratios
As shown below, identity (3.1) is a special case of the following general property [16]. Let A be any n × n matrix, and A {i 1 ,··· ,ir} {k 1 ,··· ,kr} the matrix obtained from A by suppressing the lines i 1 , · · · , i r and columns k 1 , · · · , k r . Then, for any i 1 < i 2 and k 1 < k 2 :

JHEP11(2013)154
Indeed, let us specify A = S (N ) in the identity (A.16) and give the explicit forms of the other quantities obtained by suppressing appropriate lines and columns. Let us take any i = N .
Furthermore, , we first notice that, for A symmetric, and where G {j} is the (N − 2) × (N − 2) Gram matrix associated to S {j} and obtained from it via a procedure similar to the one leading to eq. (A.1). q.e.d.

B Kinematics leading to a vanishing det(G)
This appendix supplements the discussion on the kinematics leading to a vanishing det(G) provided in subsection 4.1.

B.1 General considerations
Let us consider a set {p i , i = 1 · · · , N − 1} of N − 1 four-momenta in Minkowski space, their Gram matrix 16 G ij = 2 (p i · p j ), and the linear system given by A vanishing det(G) means the existence of a set of scalars {x j , j = 1, · · · , N − 1} not all vanishing and solution of the system (B.1). Multiplying eq. (B.1) for each i = 1, · · · , N − 1 by x i and summing over i leads to the condition which means that (i) l vanishes i.e. the {p i } are linearly dependent momenta, or that (ii) l is lightlike and eq. (B.1) is the orthogonality condition (l · p i ) = 0, i = 1, · · · , N − 1 [19]. Let us focus on case (ii) assuming furthermore the {p j } to be linearly independent. 17 The orthogonality condition requires that none of the p j be timelike, and p N ≡ − N −1 j=1 p i is orthogonal to l too, thus cannot be timelike either.
If one of the p j , say p 1 , is lightlike, l is proportional to p 1 , and all the p j =1 shall be spacelike. Were p N lightlike it should be ∝ p 1 , which would contradict the extra assumption of linear independence of the {p i }, i = 1 · · · , N − 1, thus p N shall be spacelike too. These requests impose a steric constraint on N w.r.t. the spacetime dimension d = 4. As seen by decomposing p j as (p 0 j /p 0 1 ) p 1 + q j for j = 2, · · · , N − 1 in a frame chosen such that p 1 = p 0 1 (1; 0 ⊥ ; ǫ, ), with ǫ = ± and q j = (0; q ⊥ j ; 0), the maximal number of possibly independent q j is d − 2 = 2 i.e. N shall be ≤ 4. Besides, for N ≥ 4, NLO calculations involve no one-loop N -point function with external momenta all spacelike but one lightlike, neither in collision nor in decay processes: alternative (ii) only occurs for N = 3 for any NLO purpose.
If none of the p j j = 1, · · · , N − 1 is lightlike, all of them shall be spacelike. The momentum p N shall be either lightlike -hence proportional to l: one is driven back to the previous case; see the N = 3 case below -or spacelike. The latter case is submitted to a similar steric constraint as above, as seen by trading one of the p j for l; no such configuration matters at NLO whatever N .
In summary, for any practical purpose at NLO, a vanishing det(G) can happen for a linearly independent kinematic configuration only in the case N = 3. Otherwise the configurations with vanishing det(G) correspond to linearly dependent four-momenta. 16 The overall factor 2 in the definition of G is actually irrelevant in the present discussion, we keep it only for notation consistency with the bulk of the article. 17 If both properties (i) and (ii) are simultaneously fulfilled, then the rank of the (N − 1) × (N − 1) matrix G is (at most) (N − 3), corresponding to quite degenerate configurations. For example for N = 3 this corresponds to all four-momenta lightlike and collinear to each other, for which G identically vanishes. For N = 4 it corresponds to two spacelike and two collinear lightlike external momenta being a linear combination of the two spacelike ones. None of these cases are involved in NLO calculation of processes relevant e.g. for collider physics.

B.2 Focus on N = 3
This appendix elaborates on the case N = 3 involved in subsubsection 4.1, with p 1 and p 3 linearly independent and spacelike. We parametrize the lightlike combination l (defined up to an overall multiplicative constant) as l = p 1 − x p 3 . The orthogonality conditions (implying that l is lightlike) read: The The condition det(G) = 0 also implies that s 2 = s 1 + 2 (p 1 · p 3 ) + s 3 can be written Therefore s 2 = 0 iff sign(p 1 · p 3 ) = + and s 1 = s 3 , in which case x = − 1 and p 2 = − l. Otherwise s 2 < 0.

C Spectral features of S for N = 3
This appendix gathers the spectral properties of S for N = 3 which are further used in appendix D.
Accounting for the condition m 2 1 = m 2 2 ≡ m 2 and the parametrization used in subsection 4.2, the kinematic matrix S reads: Let us compute the eigenvalues σ 1,2,3 of S and the corresponding eigenvectors v (1,2,3) in the regime det(G) → 0, det(S) → 0 corresponding to s − → 0, s 2 → 0. Since σ 3 → 0 whereas σ 1,2 remain nonvanishing in this regime, in order to correctly get the coefficients b j and B in eqs. (4.16) and (4.17) respectively in subsection 4.2, we shall keep the leading dependence on s − , s 2 in σ 3 and in the components of the corresponding normalized eigenvector v (3) , whereas s − and s 2 can be safely put to zero to first approximation in σ 1,2 and the corresponding normalized eigenvectors v (1,2) . This is the approximation to which we provide the algebraic results below.

C.1 Eigenvalues
The characteristic polynomial P S (s) of S is:

JHEP11(2013)154
The small eigenvalue σ 3 . The eigenvalue σ 3 vanishing as det(S) may be extracted from eq. (C.2) rewritten by an iteration generating an expansion in integer powers of det(S). The leading term of this expansion, itself truncated to keep only the leading dependencies in s 2 and s − , is given by: Using (tr(S)) 2 = 2 2 2 m 2 2 + 2 2 2 m 2 2 m 2 3 + 2 m 2 3 2 (C.5) we have: The eigenvalue σ 3 thus has the following approximate expression: in which the terms dropped are of order s 2 2 , s 2 s 2 − , s 4 − and higher.
The two non vanishing eigenvalues σ 1,2 . The two other eigenvalues σ 1,2 are obtained from the factorization of P S (s) as:

C.2 Eigenvectors
The eigenvector v 3 associated with σ 3 . The components x, y, z of v 3 are solutions of the degenerate system: Subtracting (C.20) from eq. (C.21) yields: eq. (C.23) tells that z = O(s − (x − y)): x and y being bounded, z thus vanishes at least as O(s − ) in the limit s − → 0. We shall keep the leading dependence on s − and s 2 in the components of v (3) . Up to an overall normalization, x, y and z are given by: In eqs. (C.24)-(C.26), the dependence on s 2 has been traded for s − and σ 3 up to terms neglected at the approximation retained. 18 In what follows it is not necessary to normalize the vector l (3) to 1.

JHEP11(2013)154
the unnormalized eigenvector v unnorm (3) given by eqs. (C.24)-(C.26) can be written: v unnorm The vector v (3) ∅ is the normalized eigenvector of S associated with the eigenvalue σ 3 = 0 when det(S) = 0. Let us notice that (l T The O(s 2 − ) terms are no more explicited in eq. (C.29) as they would contribute in appendix D beyond the level of approximation retained only. The departure of v (3) from v (3) ∅ in eq. (C.29) does not depend explicitly on σ 3 , it only depends on s − .

JHEP11(2013)154
where the column vector e has been defined by eq. (4.15), and where "· · · " in eq. (D.2) stand for evanescent terms in the limit considered. As the saying goes, 'a tedious but straightforward' algebraic juggling, sketched below, leads to the following expressions for the sought coefficients.

JHEP11(2013)154
Rewriting we get: The combination of eqs. (D.18) and (D.21) using eq. (D.11) leads to: As a check, let us combine eqs. (D.12) and (D.22). We get: where we recognize the identity Notice that B 0 happens to coincide with the limit t → ∞ of B seen as a function of t = s 2 /s 2 − , as given by eq. (F.2) in appendix F.
E A relation between the zero eigenmodes of S and G (N ) when det(G) = 0, det(S) = 0 We specify a = N to fix the ideas. When det(G) = 0, condition (4.4) is equivalent to the condition det(S) = 0 according to eq. (4.2) only if G (N ) has rank (N − 2). On the other hand when G (N ) has a lower rank, its cofactor matrix G (N ) vanishes identically thus det(S) = 0 again. However, as already mentioned, Gram matrices G (N ) with ranks ≤ (N − 3) for N = 3, 4 correspond to quite peculiar and degenerate kinematics irrelevant for NLO processes, thus we do not provide any more detail about this case here. We focus on the generic case for which the Gram matrix has rank (N − 2) i.e. exactly one vanishing eigenvalue. When det(G) and det(S) vanish simultaneously the eigenvectors v (N ) andn (N ) corresponding to the eigenvalues zero of S and G (N ) respectively, happen to be simply related.

JHEP11(2013)154
To see this, using eqs. (2.7), (2.8) let us write the components i = 1, · · · , N of S · v for any N -column vector v as: As argued in subsection 4.1, the eigenvectorn (N ) fulfills condition (4.4). Therefore the ansatz v (N ) j ≡n (N ) is solution of system (E.1). Furthermore it satisfies the property N j=1 v (N ) j = 0 by construction. If S has rank (N − 1), the eigendirection of S associated to the eigenvalue zero is thus identified.
Conversely, consider v such that If δ = 0, the (N − 1)-column vectorn i ≡ v i , i = 1, · · · , N − 1 is an eigenvector of G (N ) associated to the eigenvalue zero and fulfilling condition (4.4). Let us conclude these considerations with the following remarks.
1. We recall that, at the one loop order which we are concerned with, the Landau conditions [19,20] characterizing the appearance of kinematic singularities require v i ≥ 0 for all i = 1, · · · , N and δ > 0. One may hastily infer that, an eigendirection zero of S associated with a vanishing det(G) is not associated with a kinematic singularity as characterized by the Landau conditions. This does hold true if S has rank (N − 1).
In particular, v (1) and v (2) may both fulfill the Landau conditions corresponding to piled-up singularities. The so-called double parton singularity [23] is one interesting case of this kind. It occurs for the four-point function with opposite lightlike and opposite timelike legs and with internal masses all vanishing, for which det (S) ∝ det(G) 2 , when the two lightlike momenta are incoming head-on and the two timelike external momenta are outgoing back-to-back in the transverse plane w.r.t. the incoming direction. 22 3. In practice we shall however stress that such a degeneracy of the zero eigenvalue of S is very peculiar. Beside the double parton singularity, this situation happens to occur for N = 4 with three of the four internal masses equal, for very peculiar degenerate kinematics involving two spacelike momenta, and two lightlike momenta orthogonal to the spacelike ones, collinear to each other and being linear combination of the spacelike ones. . . This quite degenerate kinematics namely corresponds to a Gram matrix with rank 1 only. As already said, such an odd case is actually irrelevant regarding NLO processes at colliders, thus do not deserve any further detail here.
F The directional limit det(S) → 0, det(G) → 0 of eq. (4.1) is actually isotropic This appendix presents an analytical proof that, whereas each of the three terms involved in eq. (4.1) are separately functions of t in the directional limit s − → 0, s 2 → 0 with s 2 /s 2 − = t fixed, the limit of their sum is actually independent of t. For this purpose we compute the t-derivative of this sum in this limit and prove it to vanish identically in t. We provide an explicit proof for I 4 3 ; the I 6 3 case, albeit more cumbersome, can be handled in a completely similar way. 21 If S has rank ≤ N − 2, such eigenvectors v (1) and v (2) can always be found even if S is complex. In the latter case, consider linearly independent eigenvectors u1 and u2 of S S † associated with the eigenvalue zero: their respective complex conjugates v (1) = u * 1 and v (2) = u * 2 are linearly independent eigenvectors of S associated with the eigenvalue zero. The matrix G (N ) being symmetric real, the eigenvectorn of G (N ) built as described shall be made real by an overall phase shift. 22 The fact that such configuration leads to a vanishing det(G) does not contradict the classification of the eligible kinematics provided in appendix B. This appendix focused on the kinematical configurations corresponding to linearly independent sets of four-momenta. On the contrary the double parton scattering singularity appears for coplanar configurations of linearly dependent four-momenta.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.