Constructing d-log integrands and computing master integrals for three-loop four-particle scattering

We compute all master integrals for massless three-loop four-particle scattering amplitudes required for processes like di-jet or di-photon production at the LHC. We present our result in terms of a Laurent expansion of the integrals in the dimensional regulator up to 8th power, with coefficients expressed in terms of harmonic polylogarithms. As a basis of master integrals we choose integrals with integrands that only have logarithmic poles — called dlog forms. This choice greatly facilitates the subsequent computation via the method of differential equations. We detail how this basis is obtained via an improved algorithm originally developed by one of the authors. We provide a public implementation of this algorithm. We explain how the algorithm is naturally applied in the context of unitarity. In addition, we classify our dlog forms according to their soft and collinear properties.


Introduction
Perturbative quantum field theory allows us to derive predictions for physical observable from our in-principle understanding of the fundamental interactions of nature. Experiments like the Large Hadron Collider (LHC) allow us to measure such observables and test our current conceptions of the world. One particular observable that allows us to probe the strong interactions is the production cross section of sprays of collimated hadronsso-called jets. This observable is measured with astounding precision at the LHC. Consequently, in order to maximally benefit from this measurement the precision of the theoretical prediction for this observable must at least match the experimental one. In order JHEP04(2020)167 to achieve this goal it is necessary to compute sufficiently many orders in the perturbative expansion of the cross section for the desired observable.
When physical quantities in quantum field theory are expanded perturbatively in the coupling constant, corrections beyond the leading order involve Feynman loop integrals. Examples are correlation functions depending on positions of operators, or scattering amplitudes depending on on-shell particle momenta. Feynman integrals typically evaluate to multi-valued functions, such as logarithms, dilogarithms, and generalizations thereof.
It is of great physical but also mathematical interest to understand better the connection between the Feynman integrals and the special functions that arise. In recent years, such insights allow us to predict the type of special functions, and their 'fine structure', that arise from carrying out the loop integrations, simply by analyzing properties of the Feynman integrand. These insights have already had numerous applications and streamlined many complicated calculations.
An important class of special functions is that of multiple polylogarithms [1,2]. They are iterated integrals having the same integration kernels as logarithms. The number of integrations of multiple polylogarithms is called the (transcendental) weight. For example, logarithm and dilogarithm have weight one and two, respectively. Functions with more general integration kernels may also arise in Feynman integrals, but are beyond the scope of the present paper and are not discussed here.
A heuristic observation is that L-loop integrals in four dimensions give rise to functions of weight lower or equal to 2 L. For example, at one loop in four dimensions, the maximal weight is two, which means that the space of functions is given by algebraic functions, (products of) logarithms, and dilogarithms. A special role is played by the functions of maximal weight 2L. Many examples of such functions were encountered in N = 4 supersymmetric Yang-Mills theory. It appears that many quantities in that theory are naturally expressed in terms of functions of uniform and maximal transcendental weight, see e.g. [3][4][5][6][7].
There is a conjectured connection between uniform weight integrals and properties of their integrands: the singularities of the integrand are locally of logarithmic type. This conjecture has been tested for many cases, originally in the context of planar, finite integrals in N = 4 supersymmetric Yang-Mills theory. However, this notion generalizes in a number of ways. First of all, the dual conformal symmetry of the theory (which implies a certain power counting) is not essential: for example, at one loop both box and triangle integrals give rise to uniform weight functions. Moreover, generalizations include non-planar integrals, integrals involving massive particles, for example. An important further generalization concerns dimensional regularization, where integrals are computed in D = 4 − 2 dimensions, in a Laurent expansion for small . Observing that poles such as 1/ in the dimensional regulator would correspond to log Λ for some cutoff Λ, it is natural to assign a transcendental weight −1 to . This seemingly simple concept has important repercussions. What does it mean for a function f ( , x) to have uniform weight? Writing

JHEP04(2020)167
it means that f (k) (x) has weight k, for any order in the expansion! This is a rather strong condition. In practice, the fact that properties of the loop integrand may predict which integrals evaluate to uniform weight functions is extremely helpful. The classification of integrands having d log forms can be done at the integrand level, i.e. prior to integration. This connection is well-known and has been investigated and used in a number of papers, e.g. [8][9][10][11][12]. An algorithm to do this was implemented in [13]. It is based on a suitable parametrization of the loop integrand, and analyzes the resulting rational function by taking residues iteratively. This approach is complementary to the algorithm implemented in [14] that uses methods from computational algebraic geometry to compute multivariate residues. Algorithms that can be applied to Feynman integrals (in contrast to integrands) in conjunction with the methods of differential equations [15][16][17][18][19] in order to find uniform weight integrals were discussed in refs. [20][21][22][23][24][25] In the present paper, we discuss a refined version of the algorithm of ref. [13] to find d log forms. The improvements mainly concern the following two points. Firstly, at some stage of taking residues, one may encounter integrands with denominators that are quadratic in the integration variables. We introduce a method that allows the algorithm to proceed in those cases. Secondly, the analysis performed to find integrands having d log forms is closely related to taking (generalized) cuts of integrands, and in particular to leading singularities. The latter correspond to taking the loop integrand, and performing contour integrals to take multiple residues, thereby completely localizing the integration. Obviously, doing so is much simpler than carrying out the loop integration over Minkowski space-time. We use this connection to organize the analysis of loop integrands according to different cuts, thereby simplifying each individual calculation.
It is worth pointing out that generalized cuts and leading singularities are also important methods for computing loop integrands that bypass Feynman diagrams. Given the way it is defined, the uniform weight integrands we construct are very natural building blocks for such integrand constructions, and we expect our results to be useful in this area. For recent references in this direction, see e.g. [26][27][28][29][30][31].
There is a further application of d log integrands, namely an improved control over the singularities of Feynman integrals after integration. On the one hand, it turns out that d log integrals in four dimensions are ultraviolet finite. This can be shown by a power counting argument which we explain below. On the other hand, on-shell amplitudes may have infrared (soft and collinear) divergences. For a given loop integrand, it is easy to analyze the soft and collinear behaviour responsible for divergences. By doing so one may select a basis of loop integrands/integrals with improved convergence properties. While examples of this are well known at one loop, this was first discussed systematically at higher loops in [8], with the aim of introducing finite loop integrands that are relevant for infraredfinite parts of scattering processes. The improved understanding of infrared properties of loop integrands was also used to determine the latter via bootstrap methods [32,33]. See ref. [34] for a recent application of the classification of d log intergrands according to divergence structure to four-loop form factors. It is possible to algorithmically find finite but not necessarily uniform transcendental Feynman integrals, see for example refs. [35,36]. Figure 1. The nine integral families needed to describe all master integrals for three-loop massless four-particle scattering. The external legs are associated with the momenta p 1 , p 3 , p 4 and p 2 in clockwise order starting with the top left corner.
Let us now return to the question of the evaluation of the loop integrals. As was already mentioned, knowing (conjecturally) that a given loop integrand integrates to a pure uniform weight function provides a lot of information. In fact, it is easy to see that a pure function satisfies simple differential equations. Moreover, any Feynman integral satisfies some n-th order differential equation. Equivalently, one may transform this into an n × n system of first-order differential equations for the Feynman integral and other functions (e.g. derivatives). Combining this with the information about the form of differential equations for pure functions one may conclude that one may always reach a canonical form of the differential equations [18]. The latter are very useful for computing Feynman integrals, as they are in a form where the solution in terms of special functions can directly be read off.
In this paper we apply these methods to all three-loop integrals needed for two-to-two scattering. The integrals can be arranged into nine integral families shown in figure 1. The first analytical result for three-loop ladder boxes was obtained by one of the present authors in ref. [37]. The two planar families, (a) and (e) were computed previously in ref. [38]. Some of the non-planar integrals were computed in ref. [39]. In the present paper we report for the first time on the full set of integrals.

JHEP04(2020)167
The paper is organized as follows. In section 2, we introduce our notations and conventions. Then, in section 3, we present an improved version of the algorithm of ref. [13] to find d log integrands. In section 4, we explain how this can be combined with ideas from generalized unitarity, and point out differences. In section 5, we discuss practical aspects of the application of the algorithm, and comment on the scope of applications with the current implementation. In section 6, we discuss the results of the application of the algorithm to three loops. We also classify the resulting integrands according to their soft and collinear properties. In section 7, we discuss the reduction to master integrals and the computation of the latter using differential equations. We explain how we fix the boundary conditions from physical consistency relations. Moreover, we discuss relations between integrals from different integral families, and present a minimal set of master integrals.

Conventions, notation for integrands
In this section we introduce the notation and set-up for our computation of Feynman integrals contributing to four-particle scattering. We denote the momenta of the four particles by p 1 . . . p 4 and consider all of them to be in-going such that the momentum conservation identity is satisfied. The external particles we consider are massless and on-shell such that p 2 i = 0. Furthermore, we define the Lorentz invariant scalar products Due to the specific kinematic scenario the following identity is satisfied: We always choose to eliminate the momentum p 4 using momentum conservation in our Feynman integrals. This in conjunction with the above equation allows us to express all our integrals in terms of only two variables s and t. We define If we are describing a scattering process where particles with momenta p 1 and p 2 scatter and produce particles with momenta p 3 and p 4 then both s and x are positive and x ∈ [0, 1]. The Feynman integrals under consideration in this article are plagued by ultraviolet and infrared divergences which we regulate by working in the framework of dimensional regularization and using the generalized spacetime dimension (2.5) Above, D 0 is a generic even integer and can be specified to be D 0 = 4 in order to achieve physical results. Throughout this article we will denote Feynman integrals by the letter J and differential forms that are integrated by the letters I. With this we may write

JHEP04(2020)167
In the above equation we introduce furthermore a convenient normalization factor that depends on the number of loops in the Feynman integral L.
where γ E is the Euler-Mascheroni constant.

Computing dlog forms algorithmically
Feynman integrands with integrands that can be written as d log forms are important, as they (conjecturally) evaluate to uniform weight functions after integration. In this section, we discuss a systematic way of finding Feynman integrands with this property. We introduce the necessary concepts, and illustrate the individual steps of the algorithm by examples.

dlog forms and leading singularities
We are interested in (Feynman) integrands that have the property that they can be written as a differential form that behaves as dx/x in each variable near singularities. More precisely, given a set of integration variables x i , for i = 1, . . . , n (typically, the components of the loop momentum), and external variables y j (such as Mandelstam invariants and masses, for example),we define the differential Then an integrand admitting a d log form can be written as Here and in the following the wedge corresponds to the usual definition of a differential form giving rise to an oriented volume after integration, such that e.g. dx 1 ∧ dx 2 = −dx 2 ∧ dx 1 . We see that for each term in the sum, one could change variables from x i to the set g (k) i =: τ i . The corresponding term would then look like Consequently, it is evident that all singularities of (3.2) locally (in an appropriate set of variables) behave as dx/x. The following comments are in order: • Often, one is interested in loop integrals in D 0 − 2 dimensions, for some integer D 0 . Below, we mostly consider properties of the D 0 -dimensional part of the integrand. This turns out to be sufficient for our purposes here. See [40] for a refined analysis that allows to discriminate between integrands that vanish at D = 4. • When analyzing integrands one may change variables from the loop momentum to some other convenient variable. For the question about a d log form of the integrand to be well defined it is important to allow only algebraic changes of variables.

JHEP04(2020)167
• Two integrands may lead to the same integrated function, but differ for example by a total derivative that integrates to zero. For example, we will see that the triangle integral of figure 2(c) has a d log integrand, while the bubble integral of figure 2(a) does not, although the two integrals are equivalent after integration.
This d log property of the integrand is sometimes referred to as integrands having only logarithmic singularities, as opposed to double poles. We emphasize that for this terminology to be meaningful, it is important to distinguish between integrands and integrals.
The coefficients c k can be computed, in principle, by taking multiple residues, for example by evaluating the integrand along the contour encircling the poles at τ i = 0. The coefficients c k are called leading singularities. (In some abuse of notation, sometimes the locations τ i = 0 are also called leading singularities.) Let us illustrate the d log property with some examples, following [12]. The fourdimensional integrand of the triangle integral is given by It is convenient to parametrize the loop momentum using spinor variables, p i = λ iλi , The two complex vectors multiplying α 3 and α 4 are orthogonal to p 1 and p 2 . Their normalization was chosen such that they have zero helicity weights. This implies that scalar products with other vectors can always be rewritten in terms of the standard Lorentz invariants s and t. This change of variables leads to d 4 k ∼ s 2 dα 1 ∧ dα 2 ∧ dα 3 ∧ dα 4 . (Here and in the following we tacitly drop numerical multiplicative factors.) Plugging this into equation (3.4), we obtain .

JHEP04(2020)167
One may verify (by differentiation) that this can be rewritten in the following way This is of the form of eq. (3.2). Remarkably, only a single term is needed. We also see that the leading singularity of this diagram is 1/s. 1 Of course, in this simple case this can also be seen by dimensional analysis. One may make a further interesting observation. Written in momentum variables eq. (3.7) takes the form Here k * + = βλ 2λ1 , for arbitrary β. (Obviously, (3.8) is independent of β.) A similar formula holds with k * − = βλ 1λ2 . When the triangle integrand is written in the form (3.8) we can see a close relationship between generalized unitarity and leading singularities. It might appear surprising at first sight that one may take a four-fold residue for an integral having only three propagators. To see this, it is important to realize that k * ± correspond to the two solutions of the maximal cut of the triangle integral. The leading singularity 1/s can be computed by first taking the maximal cut, which corresponds to taking the residue at k 2 = 0, (k +p 1 ) 2 = 0, (k −p 2 ) 2 = 0. Upon taking this maximal cut, a Jacobian factor is produced. For example, for one of the two possible cut solutions, this factor is 1/(sα 3 ). So we have This form has new poles at α 3 = 0 and α 3 = ∞, which were not manifest in the original integrand (3.6). The leading singularity ±1/s is then obtained by taking a further residue at either of these poles. Leading singularities involving such poles are called composite. We remark that whenever a d log representation of the form (3.2) is known, verifying it is relatively straightforward. On the other hand, determining whether such a form exists for a given integrand, and computing it, is more complicated. In the following we will present a method to derive d log forms in an automated way.
Not all Feynman integrands admit the representation (3.2). Whenever the integrand has a double or higher pole, it is impossible to rewrite it in the form dx/x (restricting ourself to algebraic changes of variables). For example, dα α 2 does not admit a d log form. Similarly, if the integrand goes to a constant in some variable, this means that there is a double pole at infinity.
Note that double poles are not always obvious and sometimes are revealed after computing residues. As an example, consider the bubble integral of figure 2

(b). Its integrand is
(3.10) 1 Note that leading singularities are only defined up to a numerical factor, since we can always rewrite dlog forms like d logA = 1 2 d logA 2 .

JHEP04(2020)167
Using again the parametrization in eq. (3.5) we have . (3.11) Taking residues at α 4 = (α 1 + 1)α 2 /α 3 , then at α 3 = 0, and finally at α 2 = −α 1 , we find We denoted the resulting form as a 'cut' integrand (in analogy with generalized unitarity). We see that the form in eq. (3.12) has a double pole at infinity, and hence I (4) 2 does not admit a d log form. Note that this also implies that any multi-loop Feynman integrand with a bubble sub-loop cannot be written as a d log form.

Partial fractioning method
In this section we show how partial fractioning can be used to systematically derive d log forms and thereby also compute the leading singularities for a given integrand. The idea is very simple: we start with one integration variable (in principle, any), and partial fraction. We then write each fraction in that variable as the differential of a logarithm. Then, we proceed with the next integration variable, and so on, until no further integration variables are left.
The question whether this algorithm terminates is closely related to the question whether the denominator is linearly reducible [41]. Making this property obvious may depends on a good parametrization of the given integrand. For on-shell integrals, the type of spinor parametrization (3.5) turns out to be very useful.
Let us illustrate the method by reconsidering the massless triangle of the previous section.
After partial fractioning I (4) 3 in equation (3.6) with respect to α 1 , and writing the corresponding terms as differentials of logarithms, we have Iterating this for the other integration variables we find the full integrand written as a sum of d log forms: This is the direct output of the algorithm, and could be simplified. In particular, although it is not obvious, this representation is equivalent to eq. (3.7). This illustrates the fact that d log representations are not unique for given integrands.

Power counting constraints on numerators
In this section we show how excluding double poles at infinity leads to certain power counting constraints. This has an important application. It will allow us to write down, for a given (Feynman) denominator, a general numerator with a finite number of free parameters. The latter can then be fixed to find all possible d log integrands for a given denominator.
As an example for the general idea, consider the following integrand where s is an external variable. We wish to construct the most general ansatz for a polynomial numerator that covers all possible d log integrands for the given denominator.
One immediate observation we can make is that if the polynomial degree (in a given variable) of the numerator is equal or higher than that of the denominator, there will be a double pole at infinity. Using this constraint we conclude that the following ansatz is sufficient to cover all possible d log forms for this denominator.
In principle, we could apply this simple power counting constraint directly to Feynman integrands, e.g. when written in the parametrization spinor variables (see eq. (3.6)). However, for integrands built with propagators we can find even stronger constraints, as we explain presently. Let us consider a general one-loop n-point integrand with loop momentum k in integer dimension D 0 , (3.17) Here, we assume the numerator N m (k) to be a monomial of factors such as k 2 or k · q i , with the total degree being m. Here q i being an arbitrary constant vector (e.g. an external momentum). It turns out to be useful to perform a conformal inversion of the loop momentum [42], k =k/k 2 , which implies This transformation reveals a double pole ink 2 for n − m < D 0 − 1. Hence we find the constraint Note that this is not the usual loop momentum power counting, since linear factors such as k · q and quadratic factors (k + q) 2 count the same. For the triangle we then find that the only d log numerator is a constant. We also find that the four-dimensional bubble integrand of eq. (3.11) does not fulfill the power counting, which is consistent with having found a double pole.

JHEP04(2020)167
Note also that the discussion so far was for a single term in the numerator. More generally, one can show that if N is expanded in a basis of the monomials k 2 and k · p i , with i = 1, . . . , n−1, the same power counting (3.19) also applies to this situation, provided that the basis terms are independent.
There is a subtlety related to the last point that we wish to address. Since we are performing the analysis in an integer dimension D 0 , it is possible to write down linear combinations of terms that are equal to zero, but in a non-trivial way. For example, consider the Gram determinant G(k, p 1 , p 2 , p 3 , p 4 ), with the loop momentum k and four independent external momenta p 1 , . . . , p 4 . It vanishes if the loop momentum is considered D 0 -dimensional, but is non-zero for (D 0 − 2 )-dimensional loop momentum. Such linear combinations may contain terms that do not fulfil the power counting constraint in equation (3.19). On the other hand, being zero, they are trivially d log forms and therefore they seem to be a counterexample to the power counting criterium. Of course, this is not so, as the requirement of independent basis terms was not met.
In practice, it is desirable to control such evanescent Gram determinants in the numerator ansatz. One may use the refined D-dimensional analysis of [40], where the integrand is written in a D-dimensional parametrization. Using again the conformal transformation one can show that linear combinations which vanish in D 0 dimensions but violate the power counting constraint have double poles also in the D-dimensional analysis.
The same power counting constraint can also be used for multi-loop integrands by applying the constraint loop by loop.
Empirically, we also found a more restrictive criterium at higher loops, namely While we do not necessarily expect this to be satisfied in general, we found it useful as a restriction of the numerator ansatz at three loops. This point will be discussed further when presenting the results.

Dealing with non-linear denominator factors
In the previous section we discussed how we are computing leading singularities by partialfractioning denominators and subsequently by taking residues. This procedure may be obstructed by denominators that are not linearly reducible. In this subsection we describe how to proceed nevertheless in certain cases. First, we discuss the case where at least one integration variable is at most quadratic in all denominator factors. Next, we discuss how to proceed in more general cases. So we start with an integrand with a denominator that is at most quadratic in all factors for some integration variable that we call x. In a first step we make a partial fraction decomposition with respect to x such that all terms are either linear or quadratic in the denominators. For the terms with linear numerators we can proceed in the standard way. Terms with quadratic numerators have to be treated differently and have the following general form: where a, b, c, u and v may depend on other integration variables.

JHEP04(2020)167
There are two residues in x, which we denote by r 1 and r 2 . In other words, the integrand can be written as where s 1 and s 2 are the two zeros of the quadratic denominator of equation (3.21). Instead of processing with the computation with the residues of r 1 and r 2 in this form, we can first simplify the expression. We do so by rewriting the last equation as where Since r 1 + r 2 is rational, for this term the computation again can be continued with our standard methods. The term r 1 − r 2 has a square root in the denominator, so we have to find a way to deal with such a term.
In case the radicand is at most quadratic in one integration variable y and all other denominator factors are linear in y, we can proceed with the help of the following formulas: To apply these formulas we first have to do a partial fraction decomposition with respect to y while treating the square root factor as a constant and possibly do a constant shift in y to get expressions of the form (3.26) and (3.27). Note that the residue in (3.27) is in general again a square root of the remaining integration variables. So we can continue using the same formulas for the next residue in case a suitable integration variable exists. It also may happen that the residue is proportional to the square root of a perfect square. In this case the square root cancels and we may choose either sign of the square root. Let us consider now two slightly more general cases. Assume we have the following integrand dy ∧ dzN (y, z) where P (y, z) is a polynomial of degree at most two in y and of degree higher than two in z. Then neither y nor z fulfil the criteria for equation (3.27) to be applied. In this special case, however, we can make the following variable transformation in z: We see that the polynomial factorizes into two linear polynomials in y. After this transformation the integrand has only linear factors in y in the denominator and the degree of y in P (y, z) does not change. This means that we can do a partial fraction decomposition in y and then apply equation (3.27). The second special case is an integrand with only one integration variable: where P (y) is a polynomial of degree two or less in y. In this case we can force a factorization of ay 2 + by + c by also allowing square root terms of the external variables. After taking the residues we get a nested square root factor in the denominator, which does not cause a problem, because we do not take further residues. Often the radicand of the square root can be written as a perfect square and hence the nested square root can be simplified.
Finally we want to discuss the case of an integrand with a square root factor in the denominator, where the radicand polynomial is at least cubic in all integration variables. In this case none of the methods discussed so far can be applied. Here we try to proceed by performing a variable transformation depending on free parameters, and then fix the latter in order to reduce the power degree for any of the integration variable in the radicand polynomial.
As an example for such a transformation, consider the following integrand dx ∧ dy (x + y) x 3 + 3x 2 y + 3x 2 + 3xy 2 + 2xy + y 3 . (3.32) The polynomial of the square root is cubic in both variables x and y, so none of the methods discussed so far can be applied. So we make a parametrized variable transformation. For this example we consider the very simple type of transformation x → x + ηy. (3.33) We find that for η = −1 the integrand simplifies to dx ∧ dy such that the radicand is now quadratic in y and we can now take the residue in y using (3.26). For a general integrand with integration variables z 1 to z n , we make the following transformations with i = 1, . . . , n and ν ∈ {0, 1}. Hereẑ i means that this variable is left out and Q is a quadratic polynomial in all integration variables except z i . We put a free coefficient η j JHEP04(2020)167 before each term of the polynomial. Since we transform only one variable at a time and Q is independent of z i we do not change the integration measure with this transformation.
After applying a transformation we check for each variable z h , where h = 1, . . . , n, if we can choose the free parameters η j such that all cubic and higher power terms of the radicand vanish. If a transformation of this type is found we apply equations (3.26) or (3.27) if the requirements to the rest of the denominator are fulfilled. If we do not find a transformation the integrand remains unsolved.

Algorithmic implementation
The input is a denominator of an integrand, and the set of integration and external variables it depends on. The denominator is required to be polynomial in the integration variables (one overall square root factor is also allowed). The algorithm makes an ansatz for polynomial numerators. It finds all numerators that have the property that the integrand can be written as a d log form with constant leading singularities.
The algorithm, with all its steps, is visualized in figure 3. Let us go through them one by one, using the example of (3.15).
Step #1 consists in finding the most general numerator ansatz subject to power counting constraints, as discussed in subsection 3.3. In our example, the result of this step is given by equation (3.16).
Step #2 consists in eliminating double poles. Note that despite the initial constraints on the numerator, there might be further double poles in the integrand which get revealed by computing the leading singularities. We will see an example of this below, in eq. (3.37).
Step #3 : we choose an integration variable that appears linearly in all denominator factors. (If this is not possible, we continue with the method described in subsection 3.4.) In the example, this is the case for both variables a and b. Let us choose b.
Next, in step #5 we find a linearly independent subset of the residues. The residues are the factors multiplying the d log factors. Choosing an independent set makes the subsequent calculation much more efficient. In our example this step is trivial because there are only two residues that are obviously linearly independent. In more complicated cases the list of residues is significantly longer. The linear relations between the different residues can be found conveniently using numerical methods. For example, one may replace all external and internal variables by random integer numbers multiple times (at least as many times as the number of residues) and then solve a system of linear equations.
Having found the relations, we express all residues in terms of an independent basis, and collect together all d log terms having the same residue as a prefactor. In practice, this step typically halves the number of terms (which is typically of the same order as the

JHEP04(2020)167
the number of parameters in the numerator ansatz). If we are interested in computing the leading singularities only, we may just keep the independent residues, dropping the d log factors.
In step #6 , we check whether integration variables are left. If so, we continue with step #2. So, in our example we again check for double poles. Indeed, at this stage there are factors (a+s) 2 in the denominator of both summands, indicating that the integrand has no d log representation for generic n i . We find the minimal constraint on the free parameters n i , such that the double pole vanishes. In other words, we demand the remainder of the polynomial division of the numerator and (a + s) to vanish. This leads to the constraints s n 2 = n 1 , n 4 = 0 . (3.37) Solving the constraints in eq. (3.37) for the n i , and proceeding with the next steps we obtain At this stage, there are no further integration variables left, so we proceed with step #7 . Here we identify the set of linearly independent leading singularities. In our example, there are two of them, n 1 /s and −n 1 /s + n 3 .
Finally, in step #8 , we find all solutions for the leading singularities to be constant numbers. We do this in the following way. We take the list of m linearly independent leading singularities and solve the system of equations where one leading singularity is one and all others are zero. In this way we obtain m independent solutions. In other words, this last step is just the inversion of a linear system of equations. Note that if this system has no solution it means that the numerator ansatz was incomplete. On the other hand, if the solution depends on a parameter, this means that the numerator terms were not independent.
In our example, this is achieved e.g. by (n 1 , n 2 ) = (s, 0) and (n 1 , n 2 ) = (0, 1). In other words, we find the following numerator solutions This means we found a basis of all d log forms with constant leading singularities for the given denominator. Let us summarize the main steps. For a given denominator we write down a numerator ansatz that includes all possible d log integrands, making use of power constraints. By repeatedly taking residues we reveal double poles that we exclude by constraining the parameters in the ansatz. After repeatedly taking residues, we eventually obtain a list of linearly independent leading singularities. We then find all solutions to the remaining JHEP04(2020)167 parameters such that all leading singularities are constant numbers. In this way we construct, for the given denominator, a basis of integrands with a d log form and constant leading singularities.

Similarities and differences to spanning set of cuts in unitarity approach
Computing residues of Feynman integrands is obviously closely related to (generalized unitarity) cuts [43,44]. This is also very natural in the context of integration-by-parts (IBP) [45,46] relations and differential equations, as the matrices can be organized according to integral sectors defined by cuts. In particular, it is possible to organize the calculation into different parts by considering a so-called spanning set of cuts [47]. This has enormous potential, as it splits the calculation into smaller parts (parallelization), and moreover each part is much simpler compared to the full calculation, and may be optimized further.
The spanning set of cuts in the context of IBP's corresponds to the maximal cuts of the master integrals that have no subsectors with further master integrals. For the leading singularities we construct the spanning set of cuts in a very similar way where instead of master integrals we consider all integrands that fulfil the power counting criterium defined in section 3.3. This leads to a different notion of spanning cuts for computing all leading singularities to that in the context of IBP relations. This can also be understood with the difference between four-dimensional integrands and integrands in D dimensions. As a consequence, for computing leading singularities in four dimensions one can in general take cuts with more propagators compared to the cuts that are used for IBP relations.
For example, in the context of D-dimensional IBP relations, the one-loop triangle integral of figure 2(a) is equivalent to the bubble integral of figure 2(b), and hence to detect it one may cut the two propagators of the bubble only. On the other hand, in four dimensions there is no such relation, and the two integrands are separate. In fact, the bubble integral is excluded by power counting. As a consequence, in this context it is sufficient to consider cuts with at least three propagators at one loop.
To compute a cut of propagators P 1 , . . . , P n , we solve the equations P 1 = P 2 = . . . = P n = 0 for some integration variables a 1 , . . . , a n and then replace these propagators by the Jacobian J = det( ∂P i ∂a j ) −1 . In this way we obtain an integrand where n variables are already integrated out. We then apply the d log algorithm to the remaining integration variables In this way we obtain the leading singularities and reveal double poles of the integrals on the cut. Leading singularities on a cut are always a subset of the leading singularities of the whole integrand. So the strategy is to combine all results of the different cuts until we have the complete list of leading singularities.
As an example let us consider the following three cuts where the indices with value 1 correspond to propagators that are cut. The only integrand we have to consider for cut c 7 is J 1,1,0,0,0,0,1,1,1 . Setting the five propagators of c 7 to zero and solving the equations with respect to five of the eight integration variables we find four solutions. The latter can be understood as the four different helicity configurations that can be chosen when all five propagators are on-shell (see [48] for a review on this topic).  We proceed to compute the leading singularities for these four integrands and find that they are all proportional to 1/(s + t). So we can normalize the integrand by (s + t) to make the leading singularities constant on the cut.

JHEP04(2020)167
Similarly we compute the leading singularities on the other cuts for the integrals in the corresponding sectors. For the three examples we find the following integrals with constant leading singularities on the corresponding cuts: c 7 : (s + t)J 1,1,0,0,0,0,1,1,1 .   Since c 7 has no subcut in the spanning cuts, we know that (s + t)J 1,1,0,0,0,0,1,1,1 is already a d log integrand with constant leading singularities. For the cuts c 13 and c 17 we have to take into account that there might be additional leading singularities or double poles on subcuts. To compensate the additional leading singularities and cancel out possible double poles on the subcuts we might have to add integrands from the corresponding subsectors. Let us consider the second integral of (4.6). Computing its leading singularities on the subcut c 7 we find an additional leading singularity. So we make an ansatz, where we add a linear combination of all d log integrals from the corresponding subsector. In this case there is just one d log integrand, which means that we have the following ansatz: sJ 1,1,0,0,1,−1,1,1,1 + n 1 (s + t)J 1,1,0,0,0,0,1,1,1 , (4.8) Computing the leading singularities on c 7 we find that after setting n 1 = − s s+t all leading singularities are constant on c 7 . Analyzing other subcuts we do not find further leading singularities, so that (4.8) is the complete d log integral for n 1 = − s s+t . The result agrees with the corresponding d log integral in table 1.

Practical comments and scope of applications
The basic version of the d log algorithm was discussed in section 3. In addition, we use as an important further improvement the cut-based organization of the calculation discussed in the last section. Here we give further practical hints on the application to specific integrals, and comment of the scope of the applications of the algorithm.

Practical hints and comments
In order to use the algorithm in a concrete application usually some preparatory steps need to be done. We discuss these, as well as some hints for its efficient use.
• Parametrization of integration variables: we find that for Feynman integrals with massless propagators, a spinor helicity parametrization such as eq. (3.5) is quite efficient. As the latter involves the choice of two special on-shell momenta, naturally, one may try different choices, as some may be better adapted to a given diagram.
(This is even more so when considering cuts.) Let us mention also that a variant of the spinor helicity parametrization can also be used in the case of massive external kinematics, by decomposing a massive momentum in terms of two (arbitrary) light-like momenta. Finally, we want to mention that another promising choice of paramterization is the 'improved Baikov' representation, see section 3.2 of [49].
• Parametrization tailored to each cut: choosing convenient parametrizations (of internal and external variables), and of the integration order, can be or practical importance. There is further potential for refinements in this direction in the cut-based approach: there, it may be natural to choose a different parametrization tailored to each cut.
• Order of integration variables: the algorithm analyzes a given integrand one integration variable at a time. After completing the analysis in one variable, it may in principle proceed with any variable that fulfils that criteria explained above. This gives a lot of possible orderings, and it may happen that the algorithm terminates for some ordering, and not for other orderings. This is closely related to the question of linear reducibility [41]. Therefore running the algorithm with different variable orderings may resolve some cases. This can naturally be parallelized.
• Dealing with square roots in the external kinematics: usually the external kinematics is expressed with a set of Mandelstam invariants and masses. In these variables, frequently square root factors appear in leading singularities . Sometimes it is possible to rationalize (some of) the square roots by changing the parametrization. See

JHEP04(2020)167
e.g. [50] for an algorithmic implementation. This can improve the performance of the algorithm, as it tends to minimize the number of square root terms encountered in intermediate steps.
• Special kinematics for problems with many variables: having many external variables may be another source of complications, as this can make intermediate expressions grow easily to such an extent that the computation is extremely slow or even not feasible. In some cases, we already have a candidate integrand, and wish to test whether it is a d log form with constant leading singularities. This can be particularly interesting with the method [24,25] that requires only a single UT integral to determine the complete UT basis. In this case, we may e.g. replace all but one external variable by numerical constants and this way prove for each variable individually that the leading singularity is independent of it.
• Integrands beyond integer dimensions: the computation of leading singularities is usually done for integrands with integer dimensions. It turns out in most cases, integrands found from an analysis in integer dimensions can be straightforwardly upgraded to integrals with full dimensional dependence without losing the uniform transcendental weight property. Whenever this is not sufficient, a refined analysis is possible, as discussed in [40]. We find that for the integrals in the current paper this is not necessary.
• Simplified dlog forms: the output of the algorithm is a d log form that can in principle be simplified further. Sometimes one can find representations with only a few or even a single term. While this can be useful conceptually, and practically for direct integration [51], this goes beyond the scope of this paper.

Scope of applications
The package provided with this paper was successfully applied to integrals with 1) up to four loops, 2) up to five external variables, 3) integrals with massive propagators. There are many examples where the computation can be done completely automatic using preimplemented routines of the package only. In these cases we apply the following (standard) procedure: • Define kinematic setup.
• Use IntegrandAnsatz to determine the set of integrands fulfilling the power counting constrains (see section 3.3).
• For up to four external momenta from which one may be off-shell, we can directly use the routine SpinorHelicityParametrization to define a parametrization and then use Parametrize to parametrize the whole integrand ansatz. For other kinematic setups the parametrization must be set individually (see section 5).
• Then use LeadingSingularities to obtain all leading singularities and double pole constraints.
• Finally use GenerataeDlogbasis to obtain the list of dlog integrands.

JHEP04(2020)167
We will now discuss the scope of application considering different integral families and discuss in which cases we used improvements to the standard procedure described above.
• Three-loop four-point : we computed d log bases of the three-loop four-point integral families (see figure 1). For all families except family (h) the computation can be done using the standard procedure. Computing on a single kernel the computation time is between a few minutes for the simplest family (a) and 9 hours for family (i) with up to 14 GB memory. For the more complicated families (c), (f), (g), and (i) we used the package together with Macaulay 2 [52] to speed up the factorization of polynomials. The d log basis of Family (h) was obtained using the cut-based approach of section 4.
• Four loops: an example for the successful application of the package to higher loop order are the four-loop form factor integral families contributing to the quartic Casimir terms of the light-like cusp-anomalous dimension in QCD [34]. Here again for the most complicated family (C) the cut-based approach is applied, while for all other families the computation takes less than 16 hours on a single core each using up to 2.2 GB memory with the standard procedure.
• Massive propagators: as a non-trivial example of integrals with massive propagators we apply the algorithm to two-loop integrals that appear e.g. in gg → gg for a massive top quark in the loop (see also [53]). In this case it is necessary to use a particular order of the integration variables. Finding a suitable order can be done by applying the algorithm for different random variable orders (each run takes approximately two minutes) until the computation is successful. The computation for this d log basis is included in an example file.
• Five-point two-loop integrals: as an example for integrals with many scales we discuss the construction of d log integrands for five-point two-loop integral families [40,54].
Here additional steps to the routines implemented in the package are needed. While for the lower sectors the package can be used with the standard procedure the d log integrals of the higher sectors are more difficult to construct. Due to Gram determinants in the integrand ansatz that vanish in the spinor helicity parametrization that was used throughout this paper, the leading singularities obtained this way are incomplete. Hence, parts of the computations have to be performed for example in Baikov parametrization where these Gram determinants do not vanish. Due to the many scales the d log integrands are constructed using the cut-based approach and the external kinematic is chosen such that all leading singularities are rational functions.
6 Results for dlog bases at three loops

Description of method and results
Here we apply our algorithm to compute d log bases for the 9 integral families shown in figure 1 (The labelling A to I follows [12].) The classification of the planar families A and E was already done in [13] and here we present the results for the non-planar families. For the different integral families we again start with constructing a numerator ansatz, subject to the power counting constraint discussed in section 3.3, including the heuristic one of eq. (3.20). In this way, it turns out that the complication of Gram determinants is avoided, as the latter would violate this condition.

JHEP04(2020)167
We did the following consistency checks: 1) for the two planar families A and E, we checked that relaxing this constraint does not lead to additional d log solutions. 2) For all families, we checked that the ansätze are closed, in the following sense: the number of independent leading singularities equals the number of free parameters. A simple counterexample is the following list of leading singularities: Clearly, no choice of n 1 (except the trivial one) renders both leading singularities constant. A complete ansatz is always closed. Therefore the fact that this does not happen supports the hypothesis that our ansatz did not miss d log terms. We found that for all families it was possible to chose a subset of d log integrals as a basis of master integrals. In some cases it is necessary to consider an integral family together with the same graph turned 90 degrees to get a complete d log master integral basis.
The second column of table 2 shows the size of the ansatz we used. The third column shows the number of d log solutions for each integral family, where also symmetric equivalent solutions are counted. The fourth column counts the number of independent d log integrals after applying integration by parts identities. The fifth column gives the number of master integrals of the corresponding family.

Classification of the dlog integrals according to their infrared properties
It turns out that all d log integrals considered in this paper are ultraviolet finite, thanks to the power counting constraints. So the only possible divergences after integration are of the soft/collinear type. The latter are encoded into properties of the integrand, and are JHEP04(2020)167 Figure 5. To reveal the −6 pole we take the following consecutive limits: 1) k 1 collinear to p 1 , 2) k 2 collinear to p 1 , 3) k 2 collinear to p 3 , 4) k 1 collinear to p 3 , 5) k 3 soft, 6) k 3 collinear to p 2 . The soft limit in k 3 contributes a pole only if applied after the series of four collinear limits in k 1 and k 2 . especially easy to study for d log integrals. It is therefore natural to classify the integrals in our basis according to their soft/collinear behavior, following [34] (cf. [8,32,33] for earlier related work.) To construct d log integrals that are finite we take the linear combination of all d log integrands and fix the coefficients such that the integrands vanish in all soft/collinear regions. We investigate these regions by parametrizing loop momenta k i with a variable x and consider the limit x → 0. For example we use k i = xk i for a soft limit and k i = αp 1 + x 2 p 2 + xk ⊥ i for the limit where k i is collinear to p 1 . Applied to an integrand, we then have to cancel out factors such as x −1−a dx, for some a, which would result in a pole in after taking the integral near x = 0. Some infrared regions require multiple loop momenta being soft or collinear simultaneously. Moreover, to reveal all poles in we find that in some cases it is necessary to consider consecutive p j and p k collinear limits of the same loop momentum k i as in eq. (6) of [34]. For an example, see figure 5. We do the analysis for all possible momentum routings where L propagators are written as 1/k 2 1 , . . . , 1/k 2 L . It is important to note that our construction corresponds to making the integrals finite locally. This is different from integrals being finite due to some cancellation of 1/ k poles after integration, which is a weaker condition. For example, classifying the 23 d log integrands of the planar double box, we find two finite integrals, in agreement with [55].
In general we expect a given L-loop integral to have a pole of order −2L at most. To find integrals that are at most of order O( −n ) we construct linear combinations of d log integrals that vanish for any valid combination of n + 1 infrared regions. For an Lloop integrand a valid combination can involve infrared limits of L independent momenta. For each loop momentum we consider any pair of two external momenta and check for all soft/collinear regions. In doing so, we find it useful to employ a spinor parametrization based on those momenta. In this way we can have poles of order −n at most.
Since the number of different combinations of infrared regions quickly becomes very large (O(10 5 ) at three loops) we first apply them to the parent integrand of the integral family to sort out the combinations that do not contribute. The remaining infrared limits can then be applied to the general linear combination of all dlog integands in a parallelized JHEP04(2020)167  computation. Note that in order to have the complete list of combined limits we also have to consider infrared regions that contribute a pole only after a certain combination of previous limits was applied. Figure 5 shows an example where the soft limit of k 3 contributes a pole only after a series of four collinear limits in the other loop momenta.
In this way, we classified all infrared poles of the integrals at three loops. The results are shown in table 3. We provide the infrared ordered dlog integrals in the supplementary material.

All three-loop master integrals from differential equations
In this section we discuss the analytic solutions of all 3-loop 4-point master integrals. First we define a set of 9 integral families that are sufficient to contain all required scalar Feynman integrals. We label an integral of a family Λ by The factor φ (D,3) was defined in eq. (2.7). We name the families by the first 9 letters in the alphabet, such that Λ ∈ {A, . . . , I}. The factors D Λ,i correspond to integer linear combination of Lorentz invariant scalar products of external momenta and the loop momenta p 5 , p 6 and p 7 . For example, We define the 9 × 15 factors D Λ,i in the supplementary material. These factors are raised to generalized powers ν i ∈ Z. However, the set of master integrals we are interested in satisfies ν i ≤ 0 for i > 10. The nine integrals J Λ 1,1,1,1,1,1,1,1,1,1,0,0,0,0,0 corresponding to the nine integral families are represented graphically in figure 1.
Next, we define a set of canonical master integrals M Λ for each integral family. Any Feynman integral expressible in terms of the definition of eq. (7.1) and with ν i ≤ 0 for i > 10 can be related to our set of master integrals via IBP relations. A master integral is a linear combination of Feynman integrals as defined in eq. (7.1) with rational numbers JHEP04(2020)167 and ratios of polynomials of Mandelstam invariants as pre-factors. Additionally, we include a normalization factor (D − 4) 6 for each master integral. We find this canonical master integrals by applying the algorithm outlined in the previous sections. In fact we find a complete basis for all families, except for one integral in family A and 3 integrals in family C using the algorithmic approach. For the missing four master integrals we select canonical integrals that have squared Feynman propagators. Such canonical integrals cannot be found by the algorithm in the form outlined above due to the power counting constraint. The number of required master integrals per family is presented in table 2. For example we choose, We give the definition of all chosen master integrals in the form of electronically readable files in the supplementary material.
In order to obtain a solution for our master integrals we apply the method of differential equations [15][16][17][18][19] in conjunction with IBP identities [45,46]. This allows us to write the total differential of our canonical master integrals in the form a, b and c are matrices with rational rational entries. We emphasize that the canonical form of the differential equation (7.4) is obtained automatically since we are using dlog integrals as master integrals. Next, we derive differential equations in the variable x by applying the momentum conservation constraint of eq. (2.3). Finally, we solve the resulting differential equations in an expansion in in terms of harmonic polylogarithms [56] of argument x up to 8 th order in the dimensional regulator. We determine the required boundary conditions from a simple physical requirement on how the solutions behave near singular points. The matrices a, b and c have integer eigenvalues. We demand that the vector of our solutions evaluated at a singular point of the differential equations is in the kernel of the space spanned by the eigen-vectors correspoding to strictly positive eigen-values of the associated matrix a, b or c. The physical explanation of this constraint can be understood as follows. The solution of our differential equations to all orders in the dimensional regulator close to the point s 12 = 0 behaves as lim Here, M Λ,s 12 =0 represents a vector of boundary constants. The matrix exponential s a 12 involves terms of the type s a i 12 , where a i are the eigenvalues of a (in general positive and negative). UV divergences are associated with infinitely small, but positive . With the analysis of Feynman integrals we carried out in previous sections we demonstrated that it is possible to choose a basis of ultra-violet finite master integrals for any kinematical point, and in fact we did. On the other hand, our solution to the differential equations would exhibit logarithmic UV divergences for positive a i at the point s 12 = 0 for a generic JHEP04(2020)167 boundary condition. To remedy this contradiction, we can choose the boundary vector M Λ,s 12 =0 such that no such divergences are present in our solution. The same has to be true for the other two singular points of our systems of differential equations, s 13 = 0 and s 23 = 0. By expanding our general solution to the differential equations around all singular points and demanding this conditions have to be satisfied we constrain all boundary conditions except for the overall normalization. We determine the latter by computing a trivial propagator type integral.
We include our solutions to the differential equation as well as the systems of differential equations in the supplementary material. Our solution is valid in the scattering region outlined above, i.e. for s 12 > 0 and x ∈ [0, 1]. Analytic continuation into another scattering region may be performed for example by following the steps discussed in refs. [57,58]. Similarly, permutations of external legs of our master integrals can be obtained using the methods detailed in refs. [57,58]. We checked that the permutations of our master integrals satisfy the permuted systems of differential equations. Many master integrals that appear in one particular family also are contained within another and thus related to the master integrals of the other family. We provide a complete set of master integrals for each family such that there are redundancies among our master integrals. In order to remove these redundancies we also provide in the supplementary material relations among master integrals across different families. These relations relate the total of 533 integrals as listed in table 2 to 221 master integrals.
While not all Feynman integrals required for massless four-point scattering amplitudes at three loops can be expressed in terms of integrals in our families, we expect that all required master integrals can. This expectation is based on the observation that all Feynman integrals that are not expressible in terms of our families contain sub-diagrams where at least one of the loops is in the form of a triangle integral. Such integrals are always reducible via IBP identities and the resulting master integrals can be included within our nine integral families.

Conclusion and future directions
Above we outlined an algorithm to find Feynman integrals with a dlog integrand within a given integral family. A preliminary version of the algorithm to find dlog forms was presented in ref. [13] and has already found multiple applications to cutting-edge problems. These include four-loop non-planar form factor integrals [34], as well as two-loop integrals with many scales [40,54]. We discussed improvements of this algorithm and provide, for the first time, a public version.
Efficient methods for obtaining d log integrals are of great value in the process of computing Feynman amplitudes. Such integrals allow to greatly facilitate the computation of Feynman integrals via the method of differential equations. Furthermore, we discussed connections of dlogintegrals and potential application within the framework for generalized unitarity. We also outlined how such integrals may be used to find Feynman integrals that are free of infrared and ultraviolet divergences or at least have lower degree of divergence.

JHEP04(2020)167
Finally, we applied the algorithm to determine a basis of master integrals required to express any amplitude for the scattering of four massless particles at three loops with only massless virtual particles. We then computed the master integrals using the method of differential equations. Our solution takes the form of a Laurent series in with coefficients that are expressed in terms of harmonic polylogarithms. In the supplementary material we provide a definition of these integrals, their explicit solution up to O( 8 ) as well as the associated systems of differential equations. With this, all integrals needed for virtual corrections to processes like di-jet or di-photon production at the LHC at next-to-next-tonext-to-leading order are known.

JHEP04(2020)167
helicity basis (see equation (3.5)). The first argument is the list of internal (and possibly external) momenta that should be parametrized. The second argument must have the same length as the first and defines the variable names for the parametrization variables. The last argument is a list of either two or three massless external momenta. The first two momenta define the basis for the spinor helicity parametrization. The third momentum is optional and defines a normalization factor to the mixed spinor vectors. Output is a list with three elements. The first is the list of integration variables, the second is the set of equations to define the scalar products and the third is the jacobian factor for transforming the differential. The output can be directly used as an input for the SetParametrization function. For this integral all possible numerators are constructet, which fulfill the dlog power counting defined with equation (3.20). An optional second argument, which has to be an integer number, specifies the dimension and its default value is 4. The following example is the massless one-loop box family: • GenerateDlogbasis[ansatz,lsing,n]: converts a given integrand ansatz and list of leading singularities into a list of dlog integrands with constant leading singularities. The input are three arguments: the first argument is the integrand ansatz. The second argument is the pair of leading singularities and double pole constraints. The third argument is the variable name n that defines the free parameters n [1],n [2],... of the leading singularities. If not all free parameters are fixed a warning is displayed. The output is a list of dlog integrands with constant leading singularities: • UseMacaulay2[True/False]: enables or disables the usage of Macaulay2 for a faster factorization of polynomials. This function requires an installed version of Macaulay2. Furthermore the path to Macaulay2 must be assigned to the variable Macaulay2Path and the variable DataPath has to be set to a directory to save temporary files: