Two-loop scalar functions with $I$ internal lines, $I \le 5$

This article displays a proof of concept of the mixed analytical/numerical method, presented in previous publications, to compute two-loop functions with up to five massive propagators in a scalar theory having three- and four-leg vertices as the Higgs sector of the Standard Model. Several amplitudes are considered with two, three and four external legs. Some of them diverge in the UV region and we demonstrate that the method works in that case. It is shown that all these classes of amplitudes can be generated by four master topologies. Results of a numerical evaluation for some kinematics are presented, they are compared to a public code and agree well within the error bars quoted by the different programs.


Introduction
The control of quantum corrections in quantum field theory requires high performance tools and techniques to calculate the physical observables as precisely as possible. This is necessary to reduce the unphysical scale dependency from the theoretical predictions and make comparison with the precise data collected at TeV colliders. These data are decisive in inspecting the Standard Model and exploring the physics beyond the Standard Model. To achieve precise predictions, it is required to evaluate scattering amplitudes, which are the quantities connecting theory to experiments, beyond leading and next-to-leading orders. However, the increasing number of produced particles in the final state (processes with many jets), and the number of required loops at the wished order (order beyond NLO) lead to a dramatic increase in the number of contributing Feynman diagrams. These diagrams can have very complicated structure due to multi-loop topologies with massive internal propagators. Therefore, the only way out of this serious challenge is the automatic computation of the scattering amplitudes. So, it is of great importance to be able to compute scattering amplitudes beyond leading and next-to-leading order accuracy in an efficient, fast and highly automatic way.
The automation of the next-to-next-to leading order calculation (NNLO) has received a growing attention in the last few years. The major challenge of automated computation of two-loop multileg processes, which is one of the fundamental ingredients of NNLO calculation, is the stable and fast numerical calculation of the scalar Feynman integrals, since the full analytic evaluation is not available, so far, in the general multi-leg massive case. There has been an enormous progress in this task but it still extremely challenging, and the available results are derived for very specific processes, see for example refs. [1,2]. There are many approaches to evaluate two-loop Feynman diagrams, each one has its own advantages and disadvantages. It seems that the numerical approach for multi-loop amplitudes calculation is more practical and useful, since the analytical way, especially the differential equation method, which is most promising in this domain [3,4,5], becomes very complicated and remains always an issue with the growing number of external legs and scales in the general massive case because of the emergence of elliptical integrals. In alternative methods, based on Mellin-Barnes techniques [6,7,8,9,10] part of the integration is performed analytically and the remaining part is integrated numerically. As far as we know, the number of integrals left depends on the topologies under consideration which might be rather costly. The sector decomposition [11,12,13,14,15], which is a fully numerical multidimensional integration based method, gives a reliable results but it has a computing cost which can be high. This method is implemented in the public codes SecDec [16,17,18] and Fiesta [19]. There are other methods, developed in the last few years, to deal with the numerical multi-loop Feynman amplitudes, we cite for example: the numerical extrapolation approaches [20,21], the numerical solution of the differential equation [22,23], the numerical unitarity inspired methods [24], etc. It would therefore be useful to perform part of the Feynman parameter integrations analytically in a systematic way to reduce the number of numerical quadratures as proposed in ref. [25] which is the main subject of this paper. This paper is a proof of concept of the applicability of the semi-numerical method of the two-loop multi-leg calculation, in the general massive case, presented in [25]. In this approach, every twoloop scalar Feynman diagram, with N external legs in n dimensions, is expressed as a sum of some 1 two dimensional integrals as follows (2) where W (ρ, ξ) are some weighting functions and the (1) I n N (ρ, ξ) are"generalised one-loop type" N -point Feynman-type integrals. The "effective number of external legs" N and the "effective dimension" n depend on the two-loop topology in particular the number of internal lines I, and on the dimension n. The generalised one-loop function (1) I n N (ρ, ξ) is calculated analytically by the help of the "Stokes-type" identity introduced in ref. [26], see also [27,28]. Once this function is calculated analytically, the two-loop function is obtained by numerical integration only over the sole two remaining variables ρ and ξ, which, compared to the direct numerical methods, may represent a significant gain on numerical and computational efficiency. In this article, we present several two-loop functions with a number of internal lines less or equal to five (I ≤ 5) in a scalar theory having three-and four-leg vertices. Such a scalar theory could be, for example, the Higgs sector in the Standard Model. The goal of this paper is not to cover all the possible topologies and all possible kinematic configurations but rather to discuss specific ones as a proof of concept for the method presented in ref. [25]. We construct two codes, one written in Fortran and the other written in Mathematica, to perform the numerical integration over the two remaining variables. We cross checked our results, obtained by these codes, with SecDec. We showed that the results obtained by the three different codes are in a excellent agreement within the error bars.
The outline of this article is the following. The section 2 concerns preliminaries. We discuss the superficial degree of UV divergences of the two-, three-and four-point two-loop scalar Feynman integrals, in theories involving three-and four-leg vertices. Then, we remind the reader of how the re-parameterisation of the Feynman parameters works. This is discussed in more details in [25] leading to an integrand which is a second order polynomial in some new parameters. The coefficients of this polynomial may factorise under certain conditions which are presented here. The different topologies are classified in four main categories to which corresponds a formula in terms of weight functions and "generalised one-loop functions". In section 3, the first master topology is presented, this is a two-loop two-point function with three internal lines which is globally UV divergent. It is the only one of its category. Section 4 is dedicated to topologies with four internal lines. There is essentialy one master topology, a two-loop two-point function also globally UV divergent. Another topology, a two-loop three-point function, is a child of the first. Some other topologies are just mentioned because they are trivial extensions of the one presented. In section 5, we discussed the topologies involving five internal lines. There are two master topologies both two-loop two-point functions, one contains a UV divergent sub-graph, the other is UV finite. Four child topologies, two-loop three-and four-point functions, are also evaluated because the kinematics related cannot be guessed from the mother topologies. Here also, other child topologies, corresponding to trivial extensions of the mother topologies, are just quoted. In section 6, we give the numerical results for the nine topologies, presented in the previous sections, for five different kinematic configurations. To facilitate the reading of this paper, we have postponed details of the calculations to many appendices. The appendix A contains the details of the computation of the master topology of section 3. In appendices B and C, we give the results, and how to derive them, for many integrals appearing in globally UV divergent topologies. The appendix D shows how to calculate the oneloop two-point functions at order O(ε 0 ) and O(ε 1 ) which appear in the ε expansion of the UV divergent topologies of section 4. In appendix E, we present the evaluation of the one-loop threepoint function whose integration domain is the square and not the usual simplex as in the case of the ordinary one-loop functions. The appendices F and G are dedicated to the computation of the oneloop three-point function at order O(ε) appearing in some topologies of section 5. In appendix H, we provide further checks of formulae found for the different scalar two-loop amplitudes.

Preliminaries
Let us discuss the different topologies presented in this article. The scalar topologies are labelled as T 2n i j k where n is the number of external legs, i is the number of vertices with three-legs, j is the number of vertices with four-legs. The knowledge of these numbers does not fix uniquely the diagram, thus another integer k is introduced labelling the different two-loop diagrams having the same number of external legs, three-leg vertices and four-leg vertices. We also provide in the table 1, the correspondence between the labelling used in this article and the Nickel index following ref. [29].

Superficial degree of UV divergence
Let us discuss in a general way the superficial degree of divergence. For that, let us note g 3 , resp. g 4 the coupling constant for the three-leg (resp. four-leg) vertices, the dimension (in terms of energy) of these two coupling constants is: where N is the number of external bosons of the diagram G, n i is the number of vertices having the coupling constant g i .

Two-point functions
For a two-point function, the superficial degree of divergence becomes If n 3 > 2 then w(G) < 0 and the Feynman integral associated to the diagram G is convergent in the UV region 1 . So a UV divergent two-loop two-point function, in this theory, has 2, 1 or zero three-leg vertex. But once the number of three-leg vertices is determined, since the number of loop and the number of external legs are fixed, the number of four-leg vertices is also known. Indeed, let us note I the number of internal lines of the diagram G and L the number of internal loops, there are two equations which give some constraints on the different number of vertices In our case, taking n 3 = 0, eqs. leading to the following solution: n 4 = 1 and I = 4. Therefore, the two-loop two-point functions which diverge globally in the UV domain have either zero three-leg vertex and two four-leg vertices with three internal lines or two three-leg vertices and one four-leg vertex with four internal lines.

Three-point functions
For a three-point function, the superficial degree of divergence becomes: If n 3 > 1 then w(G) < 0 and the Feynman integral associated to the diagram G is globally convergent in the UV region. Taking n 3 = 1, eqs. (2.3) reduce to which implies that n 4 = 2 and I = 4. Note that for n 3 = 0, the system of equations (2.3) has no solution in N. Therefore, the two-loop three-point functions which globally diverge in the UV domain have one three-leg vertex and two four-leg vertices with four internal propagators.

Four-point functions
The same study can be followed for a two-loop four-point function. In this case, the superficial degree of divergence of a two-loop four-point Feynman diagram G is given by The only possibility to have a global UV divergence is to set n 3 = 0. Let us determine I and n 4 with the help of eqs. (2.3), we have to solve which implies that I = 4 and n 4 = 3.

Skeleton for the re-parametrisation of the Feynman parameters
The goal of this section is to give details on the way the Feynman parameters are re-parameterised for a two-loop scalar function. All the common formulae will be collected here, the specific ones will be discussed for each topologies.
We mainly follow the discussion of ref. [25]. Namely, the n-momentum q i carried by each propagator i of a two-loop scalar function is parameterised in the following way q i =k i +r i , where the nmomentar i depend on the external n-momenta p i whereas thek i are some linear combinations of the loop n-momenta k 1 and k 2 which are specific for each topology 2 . Once this linear combination is fixed, the energy-momentum conservation at each vertex determines ther i except two, let us note themr 1 andr 2 3 . Each n-vectorr i can be decomposed aŝ where α i , β i = 0, 1 and Q i ({p k }) is a linear combination of the external n-momenta p k . We further introduce a 2 × 2 matrix A, two other n-momenta r 1 and r 2 and a quantity C such that the The symbol τ j , resp. m j , is the Feynman parameter, resp. the mass, associated to the internal line having a n-momenta q j . According to its definition, cf. eq. (2.11), the matrix A is equal to where the three sets S 1 , S 2 and S 3 are defined as: S 1 is the set of propagator labels whose momentum involves only k 1 not k 2 , S 2 is the set of propagator labels whose momentum involves only k 2 not k 1 and S 3 is the set of propagator labels whose momentum involves both k 1 and k 2 4 . Let us choose N − 1 independent linear combinations of the external n-momenta 5 t i i = 1 · · · N − 1. This choice is of course specific of the topology too. The n-momenta r 1 and r 2 can be expressed in terms ofr 1 , r 2 and the n-momenta t i as andB is a 2 × N − 1 matrix depending on the Feynman parameters τ i . Furthermore, the term C which is equal to can be written as where Γ is a N − 1 × N − 1 matrix depending only on the Feynman parameters τ i . Once the integration over the loop n-momenta k 1 and k 2 is performed, the two-loop amplitude is a I-dimensional integral over the Feynman parameters τ k , constrained by I k=1 τ k = 1, of an integrand proportional to det(A) I−3/2 n [F({τ k }) − i λ] n−I cf. [25]. This function F({τ k }) is given by where Cof[A] is the matrix of the cofactors of A. One can trade the n-momenta r 1 and r 2 in eq. (2.17) againstr 1 andr 2 using eq. (2.13). But this new dependence in the n-momentar 1 andr 2 cancels against the one coming from the term C (cf. eq. (2.16)) in the polynomial F({τ k }) and the latter becomes Thus the knowledge of the matricesB and Γ which are specific to each topology enables to determine F({τ k }). Note that although the polynomial F({τ k }) is unique for a given Feynman diagram, the matricesB and Γ are not because parts ofB can be reabsorbed into Γ, cf. the discussion in the Appendix A of ref. [25]. An example will be given in the next section. Finally, the amplitude related to a scalar two-loop diagram with N external legs is given by (2) The symbol κ stands for the kinematics, it is a set containing the independent invariants and internal masses squared while the symbol T indicates the topology under consideration.
The next step is to introduce three auxiliary parameters ρ i for i = 1, 2, 3 defined by Because of its definition in eq. (2.12), the matrix A is given in terms of the parameters ρ i by and thus Note that, due to their definitions, the ρ i parameters sum to 1. Then the Feynman parameters τ k whose label belong to the set S j are re-parameterised such that τ k = ρ j u k . Since the k∈S j τ k = ρ j the parameters u k sum to 1. Eq. (2.19) becomes then (2) I n N (κ; T ) = (−1) I+1 (4 π) −n Γ(I − n) In eq. (2.23), |S k | denotes the cardinal of the set S k . Lastly, to get rid of the constraint on the sum of the parameter ρ j , the following change of variable is performed ρ j 1 = ρ ξ, ρ j 2 = ρ (1 − ξ) where j 1 and j 2 are two distinct elements of the set {1, 2, 3}, this yields where (1) I n N is the "generalised" one-loop N -point function given by (1)  with N = I − 2 and n = 2 (n − 2). The weight function W (ρ, ξ), in eq. (2.24), has the following expression In eq. (2.25), the quantityF({u k }, ρ, ξ) which is a second order polynomial in the u i , is given by 6 where U is a I − 3 vector gathering the leftover parameters u i after one got rid of the constraints k∈S j u k = 1. In eqs. (2.24) and (2.25), the argument E denotes the integration domain over the u i once these latter constraints have been taken into account, this domain is not always a simplex as in the case of the genuine one-loop functions. In the cases treated in this article it can be the vector andC is a scalar. They all depend on the parameters ρ and ξ and on the kinematics. Note that, since the dependence on the variables ρ and ξ is only implicit, in the cases where these variables are fixed, we drop them from the list of the arguments of the generalised one-loop functions for the sake of simplicity. Note also that some quantities introduced in this section depend implicitly on the topology. Strictly speaking, they have to be labelled differently for each topology. To make the notations lighter, we will ignore this except for the matrixG, the vectorṼ and the scalarC to which we will associate a different letter for each topology.
The "generalised one-loop functions" in space-time of dimension n = 4 − 4 ε will be expanded around ε = 0 as (1) where the quantities I Since we limit ourselves to topologies with a number of internal lines I ≤ 5, N ≤ 3. For N = 2, the coefficients of the expansion of the right-hand side of eq. (2.28) is given by Whatever the choice of j1 and j2, it is always possible to factorise ρ from F({u k }, {ρ l (ρ, ξ)}).
8 while for N = 3, they are given by In eq. (2.25), the way to choose the parameters u i is guided by the fact that in some case det(A) can be factorised from the vectorṼ and the scalarC. In ref. [25], we gave a sufficient condition for this property to hold. If we achieved to build aB matrix homogeneous of degree 1 in the parameter u i 7 then, the first term of eq. (2.18) which is not proportional to det(A) will contribute only to U T ·G · U , and thus the factorisation of det(A) will hold for the vectorṼ and the scalarC.
Let us try to draw the condition to build aB matrix homogeneous in the u i parameters. For that, let us discuss two-loop planar diagrams in a scalar theory having three-and four-leg vertices. These vertices are either external vertices where one or two external legs can be branched, they are N of them and are represented by blobs; or internal vertices which connect only internal lines. The topologies of these diagrams fall into three categories depicted on figs. 1 and 2.
Class of topologies T 2 In figs. 1 and 2, the dots represent trees whatever they are which contain the propagator whose labels are missing on the picture. For the class of topologies T 1 , the different sets S k are equal to The key idea is that, in order to be able to build aB matrix Class of topologies T 3 which is homogeneous in the parameters u i , the sets S i i = 1, 2, 3 must contain at least one label of a n-vectorr i whose linear combination of the external momenta p k , as defined in eq. (2.10), is zero. If we chooser 1 =r k with k ∈ S 1 ,r 2 =r l with l ∈ S 2 , then the discussion comes down to know if a n-vectorr j with j ∈ S 3 which is a linear combination of onlyr 1 andr 2 can be built. It is clear that takingr 1 =r i andr 2 =r i+1 ,r N +3 (resp.r N +2 ) in the class of topologies T 1 (resp. T 2 ) is equal tō r 1 +r 2 and thus, for these classes of topologies, it exists a n-vectorr whose label belongs to the set S 3 which has this property. On the contrary, for topology T 3 , there is no way to build the n-vector r N +1 as a linear combination ofr 1 andr 2 only because on the two vertices connecting three internal lines, there is an external n-momentum connection. In this case, all the Feynman parameters τ k with k ∈ S 3 will appear in the entries of theB matrix. And thus, after the re-parameterisation in terms of the parameters u i and ρ j , some entries of this matrix will not be homogeneous of degree 1 in the parameter u i . The same conclusion can be reached by studying the class of non planar topologies. To end up the discussion, the condition for which the factorisation of det(A) does not hold for the vector V and the scalar C for a given diagram is that it has two four-leg vertices connecting three internal lines.
In the next sections, we present a set of scalar topologies having a number of internal lines I less than or equal to five 8 . For a fixed value of I, a detailed calculation is given only for the primary topologies, they are parent topologies from where child topologies can be generated. This is due to the fact that in eq. (2.23), the complexity of the integrand depends only on I and not on the number of external legs. Thus the way to generate other topologies having the same structure as an initial N -point topology is to change the three-leg vertices into four-leg one. But it appears two kinds of three-leg vertices: the one connecting three internal legs and the one connecting two internal lines and one external line. Changing the first kind of three-leg vertex into a four-leg vertex amounts to add an extra external leg to the initial topology leading to a N + 1 topology whose kinematics cannot be guessed from the one of the initial topology without performing the explicit calculation. On the contrary, changing the second kind of three-leg vertex into a four-leg vertex leads also to a new N + 1 topology but with a kinematics which is a trivial extension of the one of the initial topology. These latter topologies will just be mentioned in the article.
In the rest of the article, we will come across typical integrals as where f is a regular function as ρ → 0. Up to terms of order of O(ε 2 ), the distribution ρ −1+ε can be replaced by The (+) distributions are defined as usual

Topologies with I = 3
There is in fact only one topology depicted in fig. 3 having three internal propagators. From the discussion of sect. 2.1, this topology has zero three-leg vertex and two four-leg vertices and the superficial degree of divergence is 2. ¤ q 2 For this topology, we parameterise the internal momenta in the following way Thus the sets S k , k = 1 . . . 3 are given by which fills the matrix A, cf. eq. (2.12). Furthermore we chooser 1 =r 1 andr 2 =r 2 . T is a one dimensional vector T = [p] and the matricesB and Γ are given bȳ Expanding the right-hand side of eq. (2.18) gives There is no need to introduce the ρ i parameters in this case, obviously ρ i = τ i for i = 1, 2, 3. Then we set τ 1 = ρ ξ and τ 2 = ρ (1 − ξ), the two-loop two-point function then becomes and κ a = {p 2 , m 2 1 , m 2 2 , m 2 3 }. Eq. (3.8) represents a two dimensional integration of a "generalised one-loop" tadpole function in a space-time of dimension 4 − 4 ε, the functionF(ρ, ξ) plays the role of the mass of the effective particle which propagates through the tadpole.
The details of the computation have been put in appendix A to facilitate the reading. Summing the different parts and gathering terms proportional to ln( where ξ ± are the roots of the polynomial ξ 2 − ξ + 1 and the integrals I j j = 1, 7 are given in the appendix B. From eq. (3.10), the divergent parts are obviously symmetric under the exchange of two masses. This is also the case for the finite part even if it is difficult to read it from the last equation. This symmetry can be seen from eq. (3.7), it is the way the Feynman parameters τ i are re-parameterised which breaks this symetry. The fact that the coefficient of the ε −2 term is proportional to masses can be understood as follows. The diagram depicted in fig. 3 has divergent subdiagrams which are one-loop bubbles. Shrinking one of them to a point leads to a tadpole which is proportional to a mass.
On a practical point of view, the two double integrals on the variables ρ and ξ appearing in the eq. (3.10) are computed numerically. Nevertheless, an improvement can be made here by noticing thatF(ρ, ξ) is only quadratic in ρ, thus the ρ integration could be performed analytically leaving only one dimensional integral to be computed numerically.

Topologies with I = 4
These topologies have a superficial degree of divergence w(G) = 0. There is only one primary topology: topology T 22 2 1 1 depicted in fig. 4, the others can be built out of it as will be shown later on. The internal momenta are parameterised in the following way and about the matrix A, cf. eq. (2.12). We chooser 1 =r 1 andr 2 =r 3 . This choice leads tō and Γ = [τ 2 ] (4.4) As in sec. 3, T is a one dimensional vector whose element is p. Using eq. (2.18), we get Parameterising the Feynman parameters τ i in the following way and performing the change of variable ρ 1 = ρ ξ and ρ 3 = ρ (1−ξ), then the two-loop scalar amplitude becomes With the parameterisation (4.6), all the non zero entries of the matrixB are homogeneous in u, thus the factorisation of 1 − ρ + ρ ξ (1 − ξ) for the coefficients ofṼ b andC b in the polynomial of eq. (4.8).
After a partial fraction decomposition, eq. (4.7) can be re-written as Notice that only the first term of eq. (4.12) diverges. Indeed, there is a global factor 1 − ρ in G 1 and G 2 and thus G 2 (1, ξ) = 0 which implies that the second term is finite.
Let us start with the first term of eq. (4.12) and for that let us introduce The function G 1 (ρ, ξ) can be expanded around ε = 0 Using eq. (2.32) and the definition of the + distributions (eqs. (2.33) and (2.34)), this leads to From eq. (4.8), one can see thatF(u, 0, ξ) is a function of u onlỹ so eq. (4.19) becomes with I 1 is given in appendix (B).
For the second term in eq. (4.12), we have to study This term does not diverge and so expanding around ε = 0 and keeping the relevant terms yields where the integrals I 2 , I 3 and I 4 are given in appendix B. Then summing up the results for M 1 and M 2 and gathering terms proportional to ln F (u, ρ, ξ) − i λ leads to the scalar two-loop amplitude does not depend on ξ (cf. eq. (4.20)). In eq. (4.24), the "generalised one-loop functions" I  As already explained, the determination of the new matrixG, the new vectorṼ and the new scalar C for this topology will be done explicitly, they cannot be extracted from those of the primary topology T 22 2 1 1 . The internal momenta are parameterised as follows and thus to the matrix A through eq. (2.12). We chooser 1 =r 1 andr 2 =r 2 . This choice leads tō The vector T is taken to be Using eq. (2.18), we get The Feynman parameters τ i are parameterised in the following way Then, the following change of variables is performed ρ 1 = ρ ξ, ρ 3 = ρ (1 − ξ) leading to ρ 2 = 1 − ρ because the ρ i 's sum to 1. This yields As discussed previously, the new kinematics carried byG c ,Ṽ c andC c cannot be guessed from the kinematics of the primary topology T 22 2 1 1 . Indeed, since the topology T 23 1 2 1 has two four-leg vertices connecting three internal lines, the factorisation of (1 − ρ + ρ ξ (1 − ξ)) for the coefficients ofṼ c andC c in the polynomial of eq. (4.33) does not hold (cf. sec. 2.2). Note that the kinematics of the primary topology (eqs. (4.9), (4.10) and (4.11)) can be recovered by letting p 1 = 0 because of the choice of the labelling of the internal legs in fig. 5. Making the change of variable u = 1 − w in eq. (4.32) leads to a new polynomial in w whose coefficients arẽ Thus letting p 3 = 0 leads to the kinematics of the primary topology also but with m 2 ↔ m 3 which is not straightforward by looking at eqs. (2) I n 3 (κ c ; T 23 1 2 1 )

Other topologies with I = 4
For the sake of completeness, let us mention two other topologies generated by changing the threeleg vertex having an external leg in topologies T 22 2 1 1 and T 23 1 2 1 into a four-leg vertex, they are depicted in fig. 6 Figure 6: New topologies generated from T 22 2 1 1 .
The amplitudes for these two topologies are still given by the right-hand side of eq. (4.24). In these cases, no need to redo the explicit computation for the kinematics. The kinematics of T 23 1 2 2 is the one of topology T 22 2 1 1 with p = k 1 = −(k 2 + k 3 ) and the kinematics of the topology T 24 0 3 1 is the one of the topology T 23 1 2 1 with p 1 = k 1 , p 3 = k 3 and p 2 = k 2 + k 4 .

Topologies with I = 5
For amplitudes having I = 5, the superficial degree of divergence w(G) = −2 but some topologies can have UV divergent subdiagrams. There are two primary topologies, one having a UV divergent subdiagram T 22 4 0 1 depicted in fig. 7 and another one with a UV free topology T 22 4 0 2 depicted in fig. 13.

Topology T 22 4 0 1
The internal momenta are parameterised in the following way and about the matrix A through eq. (2.12). We chooser 1 =r 1 andr 2 =r 2 . This choice leads tō The vector T is taken to be Using eq. (2.18), we get for the polynomial F({τ k }) Then, we parameterise the Feynman parameters τ i in the following way and performing the following change of variable ρ 1 = ρ ξ, ρ 3 = ρ (1 − ξ) the two-loop scalar amplitude becomes We, then, proceed following the strategy developed for the other UV divergent diagrams. After partial fraction decomposition, eq. (5.8) can be re-written as Notice that only the first term of the eq. (5.13) diverges. Indeed, there is a global factor (1 − ρ) 2 in K 1 and K 2 and thus K 2 (1, ξ) = 0 makes the second term finite.
Let us focus on the first term of eq. (5.13) and for that let us introduce Using eq. (2.32) and keeping terms up to order ε 0 , N 1 becomes Note thatF(u 1 , u 2 , 0, ξ) does not depend on ξ, it is given bỹ For the second term of eq. (5.13), we just have to take ε = 0. Putting the two contributions together and gathering terms proportional to 1/F(u 1 , u 2 , ρ, ξ), the two-loop amplitude reads The "generalised one-loop three-point functions" appearing in eq. (5.20) are a bit special because the matrixG d is not invertible. The general case for the computation of "generalised one-loop three-and four-point functions" where det(G) = 0 has been treated in the Appendix C of ref. [26]. In this case, the one-loop three-point function boils down to a difference of one-loop two-point functions as can be seen from eq. (2.25) of ref. [26] where the coefficients b i and B are given by eqs. (C.45) and (C.46) of the same reference, namely

Topology T 23 3 1 1
This topology is built by inserting an external leg to one three-leg vertex connecting three internal lines in the topology T 22 4 0 1 . The kinematics will be computed explicitly. The internal momenta are parameterised in the following way and to the matrix A through eq. (2.12). We chooser 1 =r 1 andr 2 =r 2 . This choice leads tō The vector T is taken to be Using eq. (2.18), parameterising the Feynman parameters τ i in the following way and performing the following change of variable ρ 1 = ρ ξ, ρ 3 = ρ (1 − ξ) the two-loop scalar amplitude becomes The coefficientsG d ,Ṽ d andC d of the preceding subsection can be recover by setting p 1 = 0 and p 2 = p 3 = p.
Note that due to the choice of the n-momentar 1 andr 2 the matrixB is not explicitly homogeneous of degree 1 in the variables u 1 and u 2 . Nevertheless the factorisation of 1 − ρ + ρ ξ (1 − ξ) for the coefficientsṼ e andC e holds. Using the transformation given by eq. (A.12) of Appendix A of ref. [25], we can introduce a matrixB (0) given bȳ The new matrixB is explicitly homogeneous of degree 1 in the variables u 1 and u 2 . With the help of eq. (A.16) of ref. [25], the new matrix Γ reads This is an example on how the transformation given by eq. (A.12) of ref. [25] works.
Since this topology is a "child" of the topology T 22 4 0 1 , the two-loop amplitude is given by the right-hand side of eq. (5.20) but with different coefficients for the polynomial in u 1 and u 2 (2) I n 3 (κ e ; T 23 3 1 1 ) 3 (Σ (2) ;G e ,Ṽ e ,C e , 0, ξ) Some comments are in order. Firstly, the coefficient of the divergent term is proportional to the one-loop three-point function obtained by shrinking the bubble in fig. 8 as it should be. Secondly, the generalised one-loop three-point function I coincide with the genuine one-loop three-point function with no restriction on the kinematics, it is given in ref. [26]. But the generalised one-loop three-point function of order ε, I Its analytical formula as well as the derivation to get it are given in appendix F.

Topology T 24 2 2 1
This topology is built by inserting an external leg to the three-leg vertex connecting three internal lines in the topology T 23 3 1 1 . The kinematics will be computed explicitly. The internal momenta are parameterised in the following way and to the matrix A through eq. (2.12). We chooser 1 =r 1 andr 2 =r 2 . This choice leads tō For this topology, the vector T is chosen to be Using eq. (2.18), and parameterising the Feynman parameters τ i in the following way and performing the change of variable ρ 1 = ρ ξ and ρ 3 = ρ (1 − ξ) as for the topology T 23 3 1 1 , this leads to For this topology, the functionF(u 1 , u 2 , ρ, ξ) is given bỹ Note that this case has two four-leg vertices connecting three internal lines thus there is no factorisation of (1 − ρ + ρ ξ (1 − ξ)) in the components of the vectorṼ f or in the scalarC f . This topology is similar to the topology T 23 3 1 1 , the only differences are the coefficients of the polynomial F(u 1 , u 2 , ρ, ξ) in u 1 and u 2 . The two-loop four-point function is thus given by the right-hand side of eq. (5.20) using the kinematics given by eqs. (5.45), (5.46) and (5.47) Note that the result of sec. 5.2 can be retrieved by setting p 4 = 0 (t = p 2 1 and s = p 2 3 ) in the latter equations. In principle, setting p 1 = 0 would also lead to a kinematics of the same type as the one of sec. 5.2 but it is not easy to see that on eqs. (5.45), (5.46) and (5.47). Nevertheless, the transformation lets invariant the integration volume. At the end of this transformation, the new polynomial in the integration variables w 1 and w 2 has the following structurẽ The last equations (5.51), (5.52) and (5.53) are more handy for the limit p 1 = 0. This is due to the mirror symmetry of the diagram depicted in fig. 9, namely As in sect. 4, other UV divergent topologies can be generated by replacing the three-leg vertices connecting two internal legs with one external leg in the topologies T 23 3 1 1 and T 24 2 2 1 by fourleg vertices leading to the following topologies. They are presented in figs. 10, 11 and 12. They correspond to a trivial extension of the kinematics from the one related to topologies T 23 3 1 1 and T 24 2 2 1 .

Topology T 22 4 0 2
This is the second primary topology with I = 5 mentioned at the beginning of sec. 5, the corresponding amplitude is free of UV divergences.
Parameterising the Feynman parameters τ i in the following way and performing the following change of variable ρ 1 = ρ ξ, ρ 2 = ρ (1 − ξ) the two-loop scalar amplitude becomes It is clear from eq. (5.63) that there is no UV divergence in this case, thus we can safely take ε = 0 in this equation. Note that the choice of the parameterisation defined in eq. (5.62) makesB homogeneous of degree 1 with respect to the variables u i , thus the factorisation of 1 − ρ + ρ ξ (1 − ξ) for the coefficients ofṼ g andC g in the polynomial of eq. (5.64).
Thus the two-loop amplitude reads 3 (K (2) ;G g ,Ṽ g ,C g , ρ, ξ) (5.68) The "generalised one-loop function" I 3 (K (2) ;G g ,Ṽ g ,C g , ρ, ξ) looks like a one-loop three-point function but the Feynman parameters runs through the unit square instead of the usual simplex. Its computation is presented in appendix E.

Topology T 23 3 1 2
This topology is generated from the primary topology T 22 4 0 2 by replacing one of the three-leg vertex connecting three internal legs by a four-leg vertex. Let us compute the related kinematics.

Topology T 24 2 2 2
This topology is generated from the primary topology T 22 4 0 2 by replacing the two three-leg vertices connecting three internal legs by four-leg vertices. In this case also let us compute the related kinematics. The internal momenta are parameterised in the following way For this topology, we choose The use of eq. (2.18) and the following parameterisation of the Feynman parameters τ i define the polynomial in u 1 , u 2 in the integrand. Furthermore the following change of variable is performed ρ 1 = ρ ξ, ρ 2 = ρ (1 − ξ). Thus letting ε → 0, the two-loop amplitude reads where κ i = {p 2 1 , p 2 2 , p 2 3 , p 2 4 , s, t, m 2 1 , m 2 2 , m 2 3 , m 2 4 , m 2 5 } and s = (p 1 + p 2 ) 2 and t = (p 1 + p 4 ) 2 . The polynomial in the variables u 1 , u 2 is given bỹ The Mandelstam variable u is taken to be u = (p 1 + p 3 ) 2 . This topology, also, has two four-leg vertices connecting three internal lines, thus there is no factorisation of 1−ρ+ρ ξ (1−ξ) inṼ i andC i .
The two-loop amplitude is still given by the right-hand side of eq. (5.68) but again with a different kinematics The results of sec. 5.5 can be obtained by setting p 3 = 0, i.e. s = p 2 4 , t = p 2 2 and u = p 2 1 (p 4 for this kinematics plays the role of p 3 for the kinematics related to topology T 23 3 1 2 ). One could recover the results of sec. 5.5 by letting p 1 = 0, it is not straightforward to see that from eqs. (5.89), (5.90) and (5.91). It is better to use the transformation on the integration variables U = L − W where This transformation lets the integration volume invariant. The polynomial in u 1 and u 2 becomes in terms of the new variables w 1 and w 2 This is more handy to use eqs (5.94), (5.95) and (5.96) to take the limit p 1 = 0. This is related to the mirror symmetry of the diagram depicted in fig. 15, namely The kinematics of the sec. 5.4 can be recovered by letting p 1 = p 3 = 0, p 4 = −p 2 , s = p 2 2 , t = p 2 2 and u = 0.

Other topologies with I = 5 free of UV divergences
As in the preceding cases, some new topologies can be generated from the topologies T 22 4 0 2 , T 23 3 1 2 and T 24 2 2 2 by changing one or two three-leg vertices connecting one external leg and two internal lines into four-leg vertices. The kinematics related to these topologies are trivial extensions of those related to the starting topologies, the former will not be discussed further. For the sake of completeness, they are presented in figs. 16, 17 and 18.

35
The results are presented under the following way (2) We keep unexpanded around ε = 0 the overall factor (4 π) −4+2 ε Γ(1 + 2 ε). For each topology, two tables are presented. The first one gives five points of the phase space, they are chosen more or less randomly avoiding exceptional kinematics since the purpose of this publication is a "proof of concept" of the method. The second one gives, for each phase space point, the real and imaginary parts for the coefficients a −2 , a −1 and a 0 with the errors related. These errors are the errors reported by the numerical integration package 9 used to performed the ρ and ξ integration. Note that the numerical values for the coefficients a −2 and a −1 come from the evaluation of analytical formulae, there is no numerical integration here. This explains that the errors for these coefficients are set to zero.
To compute the numbers in the second table we built a Fortran program to perform the numerical evaluation of the different formulae presented in this article. The results of this program are cross checked by another program built in Mathematica [31]. In addition, for each phase space point, the coefficients a −2 , a −1 and a 0 are also compared with the results of the SecDec program [11]. Our results are in perfect agreement with those of SecDec within the error bars quoted by the three programs. Note that we differ by a global −1 factor with SecDec because of the use of a different measure for the loop integration.

Summary
This article can be seen as a proof of concept of the validity of the novel framework of mixed analytical/numerical approach, introduced in ref. [25], for the computation of scalar two-loop multi-leg Feynman diagrams appearing in scalar theories with three-and four-leg vertices. This method allows us to express any scalar two-loop N -point function in n dimensions as a sum of double integrals of some generalised scalar one-loop type functions multiplied by some weight functions. In the current work, we limited ourselves to topologies with five internal propagators at most, where all the internal masses are taken to be real. To show the generality of our method, we considered several N -point two-loop functions, where some have global and/or sub-leading UV divergences.
More complicated topologies will be considered in future publications. The studied topologies are classified in four categories. The topologies of each category are generated by inserting some extra external legs to the master topology. The first class (with I = 3) contain only one topology, denoted T 22 0 2 1 , it corresponds to the two-loop two-point function with three internal lines and two external legs, which is globally UV divergent. The second class (with I = 4) contains four globally UV divergent topologies generated by the master topology T 22 2 1 1 , the later one corresponds to the two-loop two-point function with four internal lines and two external legs. The third and the fourth category (with I = 5) have two master topologies, one of them, denoted T 22 4 0 1 , contains a sub-leading UV divergence and the other one, denoted T 22 4 0 2 , is finite. Both master topologies correspond to two-loop functions with five internal lines and two external legs. One can generate from the divergent topology nine other topologies, and from the finite one eight other topologies belonging to the same class. The topologies generated by the same master topology have the same formula with different coefficients for the polynomialF. We implemented the two-dimensional integral representation of nine of these topologies belonging to the four classes in two different codes, one is written in Fortran and the other one is written in Mathematica, to perform the numerical integration over the two remaining variables (ρ and ξ). We showed that our results are in good agreement with the public code SecDec, which is based on sector decomposition method. This paper is accompanied with many appendices which are dedicated to the fully analytical calculation of the divergent parts and some related finite contributions of these topologies. Furthermore, we presented the full analytical calculation of the one-loop two-point function at O(ε) (Appendix D), the one-loop three-point function integrated over a square (Appendix E) and the one-loop three point function at O(ε) (Appendix F).
This method opens a new window towards the automation of two-loop computation, since traditional reduction methods can be applied to reduce tensor integrals to scalar integrals and scalar integrals to scalar integrals with less propgators. In principle, it is competitive with the fully numerical approaches since it keeps only two numerical integrations for every two-loop amplitude in a systematic way. It worthwhile to mention that, in the case where a sum of two-loop diagrams have to be computed, the numerical integration can be easily factorised out and performed once on the sum of the integrands which could lead to an extra gain in cpu time and precision. It is hard to do in other numerical methods. Further developments can be carried out in different directions. An immediate extension of this work would be to consider complex masses. Almost all the results are there. Indeed, all the new integrals presented in the appendix D and F trivially extend to complex masses and the results for the "generalised one-loop functions" are given in ref. [27]. These formulae for complex mass cases as well as numerical applications are not shown here mainly because this article contains already a lot of information. They will be postpone to other publications. One could, inside this scalar theory, tackle amplitudes whose number of internal lines is greater than five. For I = 6, the "generalised one-loop function" associated is a " generalised one-loop four-point function", the formulae for an extended kinematics being given in ref. [26,27,28] (see also [32]). While for I = 7, we have to deal with a "generalised one-loop five-point functions" which are a sum of "generalised one-loop four-point functions" as can be seen with standard one-loop reduction methods. Another direction would be to extend the method for theories including particles with spin 1/2 or 1, gauge theories for example. In this case, at the end of the re-parameterisation the numerator will be a polynomial in the parameters u i as well as in the parameters ρ and ξ. To each monomial, will correspond a "generalised one-loop functions" with a numerator which can be reduced to "generalised scalar integrals" -we think that any reduction method works, certainly the one proposed in [30] will do. The extra dependance on ρ and ξ will be taken into account in the numerical integration. One has to care about the proliferation of terms in the integrand which may slow down the numerical integration. Lastly, a third direction would be to consider two-loop amplitudes with IR divergences, some internal scalar particles are massless, where some extra work may have to be done on the "generalised one-loop functions" in order to push further the expansion around ε = 0. Our method works in principle but it would be good to have a proof of concept in this case too. Along these lines, the IR cases with photons or gluons could be also investigated.

Acknowledgements
This article is the continuation of a work initiated by Prof. Shimizu who passed away too early, we keep in memory a physicist who pushed his projects wholeheartedly. We would like to thank P. Aurenche for his support during this project and for a careful reading of the manuscript. This appendix presents the gory details for the intermediate computation leading to eq. (3.10). By splitting the factorF(ρ, ξ) and using partial fraction decomposition, eq. (3.8) can be written as The first term in eq. (A.1) diverges at ρ = 0. Using eq. (2.32) and expanding the coefficient 1 (ρ, ξ), the ρ integration is now of the type Injecting the results of eqs. (A.7) in eq. (A.6), we end with where I 1 is given in appendix B.
The third term diverges at ρ = 1 and ξ = 0 or ξ = 1. The same care as for the second term has to be taken here. Indeed, H 3 (ρ, ξ) is not a regular function, thus we re-write the third term as where Now, H 3 (ρ, ξ) is a regular function except may be for exceptional kinematics. So the third term T 3 can be written as The first term of eq. (A.15) is finite whereas the second one diverges. Let us focus on this second term. Integrating easily over ρ, we end for this term with The integral over ξ of the second term in the square bracket can be obtained from the integral of first term by changing m 1 ↔ m 2 . The third term in the square bracket integrated over ξ leads to finite terms. Thus, let us focus on the integral of the first term in the square bracket which is of the type Using the definition of the Dirac and the "plus" distributions, eq. (A.17) becomes Playing with the fact that As already said, the integration of the second term of eq. (A.16) can be extracted from the result of eq. (A.20) by changing m 1 into m 2 and vice versa. The integration of the third term is just given by Putting everything together and after some algebra, we get for T 3 The remaining integration over ξ can be easily performed. Adding the finite part and using the result given by eq. (B.27), the third term T 3 reads The definition of the quantities ξ ± is given in the appendix B.
The fourth term does not diverge either, it is enough to make an expansion around ε = 0. Nevertheless, it is safer to re-write this term as where The expansion around ε = 0 yields where I 2 , I 3 and I 4 are given in appendix B.

B Integrals appearing in globally UV divergent topologies
In this appendix, different integrals appearing in the UV globally divergent topologies are computed. Let us introduce the two roots of the polynomial ξ 2 − ξ + 1: which obey to the following relations ξ ± = e ±i π/3 Let us also define B = 1 − ξ (1 − ξ)) throughout this appendix.
B.1 I 1 The ρ integration can be performed easily: The first integral is done in appendix C, the other ones are not difficult to perform and finally we get This integral can be re-written as: Integrating by part the second and the fourth terms in the squared brackets leads after the ρ integration to B.3 I 3 Integrating by part the second and the fourth terms in the squared brackets leads after the ρ integration to Integrating by part the second and the fourth terms in the squared brackets leads after the ρ integration to B.5 I 5 B.6 I 6 Let us introduce n a strictly positive integer (n = 2, 3) and integrating by part over the ρ variable leads to With this last result, the integral I 6 can be put under the following form Then, using where A is a complex number with a non zero imaginary part, and the properties of the roots ξ + and ξ − yields For any integer n > 1, one can show that With this result, the ρ integration can be easily performed leading to The integrals over ξ are of the same type as those for I 6 , thus using the same techniques we readily get B.8 I 8 The ρ integration is performed using eq. (B.18), yielding The ξ integration can be readily performed leading to In several places, the following integral appears This integral is difficult to perform directly, nevertheless one can study the integral Performing the ρ integration first leads to But now, the integral L t can be computed by integrating firstly on ξ. Let us note the two roots of the polynomial −ρ ξ 2 + ρ ξ + 1 − ρ bȳ these roots obey 1 −ξ + =ξ − , thus L t becomes Since the argument of the logarithm is a product of positive factors, it can be split and the ξ integration is then performed easily, yielding To integrate over ρ, we use the Euler third change of variable: ρ = 4/(t 2 + 3) which rationalises the square root and with some algebra, eq. (C.6) becomes Then, the integral over t is split into three parts Let us then compute the different parts.
For the first integral L (1) t we make the following change of variable t = 1/ √ v, this leads to A partial fraction decomposition and a splitting of the logarithm because the numerator and the denominator of its argument are positive yields All the integrals can be easily performed and we end with For the second one L (2) t , the integrals over the two terms, taken separately, in the squared brackets of eq. (C.12) diverge at t = ∞ but the integral of their sum is convergent. Thus, let us introduce a large t cut-off Λ such that The logarithm can be split into two others and for the last integral of eq. (C.15), we set t = √ u, this leads to The integrals can be easily performed, the ln(Λ) terms vanish and we end with Thus, taking the limit Λ → ∞, we get Lastly, for the third integral, we set t = 1/u, and we have to compute A partial fraction decomposition and a splitting of the logarithm leads to Then, a change of variable is performed in each integral in such way that after this change of variable the argument of the logarithm is just the new integration variable v The integration between 1 and 2 can be safely split because the pole is far away from the real axis: So putting things together, and expressing all the arguments of the dilogarithms in terms of ξ + and ξ − , the following result is obtained

D One-loop two-point function at O(ε)
In this appendix, the computation of the one-loop two-point function at order O(ε 0 ) and O(ε 1 ) is presented. These functions appear in the ε expansion of the UV divergent topologies T 22 2 1 1 and T 23 1 2 1 . Let us remind their expression 10 where G is real and V and C are either real or complex. In the latter case, the vanishing quantity λ, inherited from the Feynman prescription for the propagators, can be neglected with respect to the non zero imaginary parts of V and C. Eqs. (D.1) and (D.2) become where u + and u − are the two roots of the polynomial in argument of the logarithm in eqs. (D.1) and (D.2) Expanding the terms in square bracket yields The integrals in the left-hand side of eq. (D.6) can be readily computed leading to Given eq. (D.7), the new integrals, with respect to the one-loop two-point function, to be calculated are 10 As the ρ and ξ dependence is implicit through the argument G, V and C, contrarily to the main text we drop them from the arguments of the functions I and thus the integral I 1 can be readily computed. For the integral I 2 , after an integration by part, eq. (D.10) can be put under the following form The first integral in the square bracket appears in the computation of the one-loop three-point function for example and can be expressed in terms of dilogarithms. Using this known result leads to Putting everything together, we end with The formulae given by eqs. (D.8) and (D.13) may not be suited for numerical applications in the case of special kinematics: G → 0 or u + → u − , etc. They have to be replaced by dedicated formula which will not be presented here for the sake of simplicity.

E One-loop three-point function integrated over a square
Despite the fact that the Feynman parameter integration domain is not the usual 2-simplex, the method developed in refs. [26], [27] and [28] can be applied to compute this function analytically. This will be the purpose of this appendix.
To start with, we slightly change the notation to match the one of these articles 11 .
where the short-hand notation D(x 1 , x 2 ) means Changing the power of the denominator in the integrand of eq. (E.1) and applying once the "Stokes-like" identity, we get At this stage, one can integrate over the variable χ using appendix D of ref. [26]. Then introducing the following quantities to get a more compact formula (E.5) 11 We also drop the dependence on ρ and ξ which plays no role here.
Another way to evaluate eq. (E.3) is to apply once more the "Stokes-like" identity on the integration over x. To do that, let us note that the four terms in the square brackets of this last equation have the form where the polynomial E(x) has the following structure Then by changing the power of the denominator in eq. (E.8) and applying the identity [A.5] of ref. [26], we end with where the function L 4 3 (∆ 2 , ∆ 1 , D) is the same as the one which appears in the computation of the genuine one-loop three-point function (cf. eq. (2.35) of ref. [26]). Let us give it here for the sake of completeness and The function I The different quantities depending on the sector i are given as a function of the elements of the matrixǦ, the elements of the vectorV and the scalarČ.
For sector 1 For sector 2 det (G [2] For sector 3 For sector 4 det (G [4] Note that despite the fact that it is difficult to express the coefficients b i in terms of elements of the inverse of the kinematical matrix S as in the genuine one-loop case, the magic identities given in [26] (cf. eqs. (2.36), (2.37)) related to the property of determinants still holds, namely Thus the different roots are given again by Note that the results presented in this appendix are valid in the real mass case. In the case with complex masses, similar formulae can be derived along the same lines but using results of ref. [27] instead of ref. [26].

F One-loop three-point function at O(ε)
This appendix is dedicated to the computation of the one-loop three-point function at O(ε) which appears for the topologies T 23 3 1 1 and T 24 2 2 1 To compute I 3 , a new integral, related to the former and on which it will be more handy to apply the method developed in ref. [26], is defined The quantity J β is linked to I by the relation The computation of J β is similar to the infrared case treated in ref. [28]. Using eq. (2.45) of this reference, multiplying by −2 −1−ε /Γ(1 + ε) to take into account the different normalisation and setting ε = −β, we get and ν = 1/(1 + β).
Using appendix D of ref. [26] to trade the ρ integration for an integration between 0 and 1 and appendix A of ref. [28] to integrate over χ, we get But what is needed is the derivative of J β with respect to β taken at β = 0 which leads to Thus we have to compute the integral between 0 and 1 of the ratio of the logarithm squared of a second order polynomial in the integration variable over another second order polynomial in the same variable.
Let us describe the method used to compute analytically the kind of integrals appearing in eq. (F.8) in terms of Nielsen polylogarithms. For that, let us study the following integral along the same line as the quantity K R 0,1 (A, B, u 2 0 ) appearing in the appendix E of ref. [26] (2) K R 0,1 (A, B, where A is real, B and u 2 0 are complex with a vanishing imaginary part. u 2 R 0 denotes the real part of u 2 0 . Let us nameū a root of the polynomial A u 2 + B,ū = −B/A, eq.(F.9) can be re-written as where S u = sign(Im(−ū 2 )) and u 0 is defined by is the sign of imaginary part of u 2 0 , namely Im(u 2 0 ) = i sλ. Expanding the terms into the curly brackets of eq. (F.10) leads to where two new integrals have been defined Using parity arguments, cf. appendix E of ref. [26], some integrals, for exampleR (2) (u 0 ,ū, u 0 ) and −R (2) (−u 0 , −ū, − u 0 ) can be glued together to give an integral between −1 and 1. In addition, the u 0 appearing at the denominator can be traded safely for u 0 (cf. appendix E of [26]). Let us introduce new notations, namelȳ where the path of integration for the two integrals is along the real axis. After some algebra, eq. (F.11) becomes where the function F( u 0 ,ū) is given by eq. (E.17) of ref. [26]. At this point, it remains to computē R (2) (ū, u 0 ) andQ (ū, u 0 ).

F.1 Computation ofQ (ū, u 0 )
For this calculation,ū and u 0 are complex numbers, the imaginary part ofū never vanishes while the one of u 0 may be zero. The product ln(u −ū) ln(u +ū) can be written in the following way ln(u −ū) ln(u +ū) = ln(−ū) ln(ū) + ln(ū) ln 1 − ū u + ln(−ū) ln 1 + ū u Furthermore, the last term of eq. (F.17) can be rewritten as It is easy to show that the η(1 − u/ū, 1/(1 + u/ū)) = 0 because u is real, thus the functionQ (ū, u 0 ) can be expressed in the following waȳ Let us define and compute the following integral where the integration variable u runs along the real axis. Setting z = u − u 0 , the integral L 2 (ū, u 0 ) becomes Then, the idea is to split the integral from −1 − u 0 to 1 − u 0 in two parts But because of the structure of the argument of the logarithm squared in eq. (F.21), the z variable cannot run along the straight line [0, 1 − u 0 ] (resp. [0, −1 − u 0 ]) for the first (resp. the second) integral of the right member of eq. (F.22). Indeed, one of the cuts of the logarithm may cross these segments as shown in appendix G and thus the contour must be deformed. In order to avoid this crossing, the strategy to deform the contour is to go through one of the branch points. Instead of integrating on a segment between 0 and 1 − u 0 (resp. −1 − u 0 ), the integration contour is chosen to be a segment between 0 and one of the branch points, then a segment from this branch point and 1 − u 0 (resp. −1 − u 0 ), cf. fig. 19. Figure 19: Example of a contour deformation running through the branch point zB. The blue (resp. red) segments are the contour deformation for the integration between 0 and 1 − u0 (resp. −1 − u0) In this way, the integration contour never crosses the branch cuts of the logarithm. Of course, any of the two branch points can be chosen equivalently and even the branch point used for the contour deformation needs not to be the same for the two paths 0, 1 − u 0 and 0, −1 − u 0 .
To fix the idea, we will choose for the two contour deformations the branch point z B =ū − u 0 (cf. appendix G). Let us introduce the integral J 1 given by which depends on a real parameter β. As discussed previously, the integration path is chosen to therefore eq. (F.34) becomes L 2 (ū, u 0 ) = 2 Li 3 ū + 1 u − 1ū If the branch point z B had been chosen instead of z B for the two contour deformations, we would have got the right-hand side of eq. (F.36) withū changed in −ū, but since L 2 (ū, u 0 ) = L 2 (−ū, u 0 ) as it is obvious from its definition, cf. eq. (F.15), the results are the same as it should be.
G Study of the cut of the logarithm squared G.1 Case of L 2 (ū, u 0 ) In this appendix, we will study the following integral with w = (z + u 0 )/ū where u 0 andū are genuine complex numbers. For a generic complex number w, let us note w R = Re(w) and w I = Im(w). The real and imaginary parts of w are given by w R = (z R + u 0 R )ū R + (z I + u 0 I )ū I |ū| 2 (G.2) and the argument of the logarithm in eq. (G.1) is given in terms of w R and w I by The cut of the logarithm is given by the two conditions Im 1−w 1+w = 0 → w I = 0 and In terms of the variable z these two conditions translate into with w 0 R = (z R + u 0 R )/ū R . Thus the cuts are some parts of a straight line z I = A 2 z R + B 2 with A 2 =ū I /ū R and B 2 = (ū I u 0 R −ū R u 0 I )/ū R in the z plane and the two branch points which correspond to w 0 R = ±1, are given by G.2 Case of L 1 (ū, u 0 ) In this case , let us study the cuts of the following function L (z) = ln 2 (1 − w) (G.7) with w = (z + u 0 )/ū. The real and imaginary parts of w is given by eqs. (G.2) and (G.3) and the branch cut of the function L (z) is given by the conditions Im (1 − w) = 0 → z I =ū Ī u R (z R + u 0 R ) − u 0 I and Re (1 − w) ≤ 0 → z R >ū R − u 0 R ifū R > 0 and w 0 R ≥ 1 z R <ū R − u 0 R ifū R < 0 and w 0 R ≥ 1 In this case, there is only one branch cut and so only one branch point which correspond to z B .

H Further checks of the two-loop scalar amplitudes
In this appendix, we provide further checks of formulae found for the different scalar two-loop amplitudes. The idea behind these checks is the following. A two-loop N -point scalar amplitude in momentum space is of the type The last term in eq. (H.2) can be viewed as two propagators of a scalar particle having a mass m k connected through a three-leg vertex whose external leg has a zero momentum, in other words, it represents a zero momentum insertion in the internal line propagating the scalar particle of mass m k . We will use this idea to relate a N -point two-loop scalar amplitude to a N + 1-point two-loop scalar amplitude whose one of the external momenta is taken to be zero. This idea is far from being new, cf. [34,35] for example. Several examples of such relations are given. Notice that, in this appendix, the different functionsF will be labelled for each topology with the same letter used to label the kinematics.
H.1 Relations between topology T 22 0 2 1 and T 23 1 2 1 A quick proof of such a relation can be provided by taking the derivative with respect to m 2 3 of eq, (3.8). Using eq. (3.9), we readily get whose right-hand side has the same structure as the right-hand side of eq. (4.32). Now, concerning