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.


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]. 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 [5] to more involved rescue solutions [6,7]. In [8] 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 [9], 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 [10]. This completion shall supply the functions I n 3 (j 1 , · · · , j 3 ), I n+2 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. Furthermore, unless 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. [11]. 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. 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 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: In this section, we will sketch the computation of I 4 3 and I 6 3 using the method of ref. [11]. These two integrals are defined 1 by: The Feynman contour prescription in the propagators is noted iλ in order to avoid any confusion with the parameter ǫ = (4 − n)/2 involved in dimensional regularization.
where I div 3 isolates the MS ultra violet pole in ǫ, and I 6 3 is the finite part which we will focus on. We may single out any index a in S = {1, 2, 3} and write The quadratic form z T S z becomes: The matrix G (a) is the 2 × 2 Gram matrix built from the four-vectors ∆ i a = q i − q a : G (a) ij = 2(∆ i a .∆ j a ). Its determinant does not depend on the choice of a, and it is also the determinant of the similar Gram matrix built with any subset of two external momenta. We note it simply det(G) without referring to a and unambiguously call it the 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. [11], c.f. their eq (5.2). We keep the same notations as those of ref. [11] for the different quantities, and we closely follow the strategy of ref. [11] 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.
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 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. 3 As seen by exhaustion, the only configurations leading to a positive Gram determinant would require that all three external four-momenta p 1 , p 2 , p 3 = −p 1 − p 2 of the three point function be spacelike. At the one-loop order which is our present concern, each of the three points, through which p 1 , p 2 and p 3 respectively flow, shall be connected to an independent tree. In order for p 1 , p 2 and p 3 to be all space-like, each of these trees should involve one leg in the initial state: this would correspond neither to a decay nor to a collision of two incoming bodies. 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.3), (2.4) the polynomial [a x 2 + b y 2 + c x y +d x+e y +f −i λ] in 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 such that the residue of the fake pole 1/A vanishes in each combination [ln(A + B) − ln(C)]/A and ln(A x min + B) − ln(C)]/A separately. The two terms in the r.h.s. of eq. (2.19) thus 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.
Similarly, for I 6 3 (j) we have: (2.25) with I 6 3 (j) of the form with D (j) , · · · , K (j) (α) given in table (2.24) above. 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. [11], 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.
For sector (1), K(α) = 1, A = 2 b z + e + c, C = c z + d + 2 a, and we get : 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 Fig. 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 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 Fig. 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.

det(G) → 0 whereas det(S) non vanishing
Let us first consider the polynomials g (j) (z) + 1/(2B) appearing in the denominators of the integrals I 4,6 3 (j) in eqs. (2.38), (2.39). Let us first consider γ ′′ (j) = 0, so that g (j) (z) + 1/(2B) is of degree two. Using the identitȳ where ∆ j has been defined in eq. (2.37), and the rescaled coefficients it is insightful to write the corresponding discriminant of g (j) (z) + 1/(2B) as 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 [12,13]. 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 (as commented earlier, det(G) ≤ 0). When det(G) → 0, both zeroes z ± j of 2Bg (j) (z)+ 1 are dragged away from [0, 1] towards + ∞ and − ∞ respectively. If γ ′′ (j) = 0, g (j) (z) + 1/(2B) 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 thus provide suitable integral representations in the case at hand. From a numerical point of view the explicit suppression of the ln after each term would have been 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 4 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 4 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.
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 (2) ) 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.

Characterization of the specific kinematics
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): where G (a) is the matrix of cofactors 5 of G (a) ; the superscript " T " refers to matrix transposition. The system det(G) = 0, det(S) = 0 is thus equivalent to the system the matrices G (a) andG (a) are simultaneously diagonalizable. When det(G) = 0 and G (a) has rank 6 (N − 2) -namely 1 in the N = 3 case at hand -the only eigendirection n (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 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 [15]. 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 7 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 8 , 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: 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: 7 This particular choice corresponds to singling out and erasing line and column 3 in the matrix S and considering the Gram matrix G (3) . 8 Indeed, at NLO, each of the external legs of the one loop three point function considered has to be connected to separate tree, and all the external legs of at least one of these three trees have to be final state legs. Therefore the external momentum flowing through the corresponding leg of the one loop three point function cannot be spacelike. Such configurations with three spacelike legs could appear only in higher loop diagrams, of which the one loop three point function would be seen as a subdiagram. 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:

Behaviour of I
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 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. They show that the quadratic forms g (1) (z) and g (2) (1 − z) become both equal to g(z) when s − = 0. The difference of the two integrals I 3 (1) and I 3 (2) in factor of b 1 and b 2 respectively, in eqs.  (2) ) remains bounded and has a finite limit when s − → 0, s 2 → 0.
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 the kinematic singularities are characterized by the so-called Landau conditions 9 [15,16] (see also [13]). 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.
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 10 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 11 and the same conclusions hold. 10 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. 11 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.

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 [18] T 0 of S(σ 3 = 0) given in the real mass case 12 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 14 w.r.t. those given by b = b 0 precisely Therefore, splitting S in real and imaginary parts S = S R − i S I , and writing S S † = S 2 R + ∆, with ∆ = i[S R , S I ] + 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 S I , as perturbative deformations of the eigenvalues σ k and eigenvectors v (k) of S R 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. 12 A similar discussion holds in the complex mass case as well with similar expressions cf. the previous paragraph. 13 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 b 3 , hence to B, coming from this component vanish, whereas the finite contribution to b 1 − b 2 from this component is weighted by the vanishing I 3 (1) − I 3 (2) in the decomposition (4.1).
14 More precisely b 0 is equal to the directional limit s 2 → 0 of b along the direction s − = 0 i.e. t = ∞. The discontinuities are meant for all other, finite t directions. 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.

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 3point 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 sect. 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" [12,14]. 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 (a) i generalizing eq. (2.8). Let us choose a = N to fix the ideas and make formulas simpler; the results obtained can be straightforwardly generalized to any a. 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 [12]. 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 : 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 15 G ij = 2 (p i · p j ), and the linear system given by 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 [15]. Let us focus on case (ii) assuming furthermore the {p j } to be linearly independent 16 . 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. 15 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. 16 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.
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.

B.2 Focus on N = 3
This appendix elaborates on the case N = 3 involved in subsubsec. 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 vanishing det(G) = s 1 s 3 − (p 1 · p 3 ) 2 ensures the compatibility of eqs. (B.3) and (B.4) in x and 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 subsubsec. 2.2.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 subsubsec. 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: 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: We further truncate 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: which requires A + σ 3 = tr(S) (C.14) The approximation corresponding to s − = s 2 = 0 in eqs. (C.14) -(C.16) replaces and the "zeroth order" approximations of σ 1,2 are given by:

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.

.2. Once normalized by
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 − .
Given eq. (D.25), (b 1 − b 2 )(σ 3 → 0) is given by: More precisely, since ((n 1 − n 2 ) T · l (3) ) = 0, the O(s − ) terms in the r.h.s. of eq. (D.27) cancel and, from eq. (C.29) and we get: As (e T · v (3) ) = O(s − ), the O(s 2 − ) correction in eq. (D.27) leads to a contribution to beyond the approximation retained. We thus keep: with (e T · l (3) ) given by eq. (D.9). This makes ) when both det(G) and det(S) tend to zero. 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. 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 subsubsec. 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.
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 [15,16] 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 [19] 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 headon and the two timelike external momenta are outgoing back-to-back in the transverse plane w.r.t. the incoming direction 21 .
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 20 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 u 1 and u 2 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. 21 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. 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.
(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.