Removing infrared divergences from two-loop integrals

Feynman amplitudes at higher orders in perturbation theory generically have complex singular structures. Notwithstanding the emergence of many powerful new methods, the presence of infrared divergences poses significant challenges for their evaluation. In this article, we develop a systematic method for the removal of the infrared singularities, by adding appropriate counterterms that approximate and cancel divergent limits point-by-point at the level of the integrand. We provide a proof of concept for our method by applying it to master-integrals that are found in scattering amplitudes for representative two-to-two scattering processes of massless particles. We demonstrate that, after the introduction of counterterms, the remainder is finite in four dimensions. In addition, we find in these cases that the complete singular dependence of the integrals can be obtained simply by analytically integrating the counterterms. Finally, we observe that our subtraction method can be also useful in order to extract in a simple way the asymptotic behavior of Feynman amplitudes in the limit of small mass parameters.


Introduction
The drive for precision in collider cross sections has become a major theme in contemporary high energy physics. Precision requires elements from theory, experiment and data analysis, and a major requirement in theory involves the calculation, including numerical evaluation, of multi-loop perturbative amplitudes and multi-particle integrals over restricted regions of phase space. As a practical matter, it is often the case that amplitudes require infrared regularization, and unphysical contributions that must cancel in physical quantities. In addition, the presence of massive particles in multi-loop amplitudes in combination with massless lines, while simplifying the overall infrared structure, often results in an even more difficult burden for analytic integrations. Similarly, even at fixed loops, large numbers of external massless lines leads to further complexity.
With all this in mind, it should be helpful to develop additional methods to separate systematically, and if possible algorithmically, infrared-divergent integration regions from infrared finite, treating the latter numerically and the former, analytically. Such methods have proved invaluable at next-to-leading order [1][2][3], and considerable progress has already been made at NNLO [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] and beyond [23][24][25][26][27][28][29][30]. In the following, we explore the use of methods inspired by those that have been developed to prove the all-orders factorization of amplitudes and cross sections to a number of examples. Our goal is to provide an in principle demonstration that such an approach, based on the nested local subtractions of infrared integration regions, may be flexible, practical and useful. Although we will restrict ourselves to amplitudes in this discussion, we hope the method will provide a way toward the combination of virtual and real corrections for infrared safe quantities without the need for infrared regularization.
Nested subtractions play a central role in the proof of factorization for inclusive color-singlet hard-scattering cross sections, such as Drell-Yan production [31], and have been developed for a proof of factorization for fixed-angle scattering in gauge theories [32]. The value of factorization in organizing the calculation of higherorder cross sections has been explored in Refs. [33,34]. Here, we will concentrate on amplitudes, and illustrate the feasibility of the direct evaluation of two loop diagrams for fixed angle scattering, based on our general knowledge of infrared structure.
Our approach will be somewhat distinct from methods used in the proof of factorization in gauge theory cross sections and amplitudes. We will start with purely scalar diagrams in φ 3 theory, and indeed, we will find that some of the severe infrared behavior found in even simple scalar examples can be treated with a straightforward approach.
In the following section we recall the infrared structure of fixed-angle scattering, describe the approach we will take, based on ordered subtractions, and illustrate the method as it applies to the one-loop box diagram. In Sec. 3 we show how this method operates when applied to essential two-loop examples, including the double box and crossed double box. Finally, in Sec. 4 we discuss the use of the method in generating asymptotic expansions in small masses that regulate infrared divergences, replacing them with logarithmic dependence.

Method for subtracting the infrared divergences in loop integrals
In this paper, we treat scalar amplitudes that describe fixed-angle scattering, reorganizing them into a sum of terms, each of which is an infrared-sensitive integral that can be performed analytically, multiplied by a coefficient integral that is a tree diagram or which can be computed numerically. To these terms will be added a "remainder", which is also free of infrared divergence. In particular, we can imagine evaluating the coefficients directly in four dimensions. It is worth noting that in dimensional regularization, it would be necessary to expand these coefficients in = 2 − d/2 around = 0 in d dimensions to derive the full expression for the amplitude, but we anticipate that in the calculation of infrared-safe cross sections, this will not be necessary.

Leading regions and power counting
The essential observation that makes the subtraction method possible is that the sources of infrared divergence in fixed-angle scattering amplitudes are associated with a finite list of "leading regions" in loop momentum space, which can be enumerated for an arbitrary diagram. The generic form follows from the application of the general Coleman-Norton criterion [35] for a Landau pinch surface [36,37] to fixed angle scattering [38,39]. In this case, pinch surfaces are associated with configurations in which internal lines are either part of a "jet subdiagram", of on-shell lines connected to one of the external lines, with all line momenta parallel to that line, or are in the "soft subdiagram", and have vanishing momentum. All other lines are off-shell. The general case is illustrated for two-to-two scattering in Fig. 1.
In general, lines with vanishing momentum (soft lines) may be attached either to each other in the soft subdiagram or to the jet functions. In gauge theories, pinch surfaces at which soft lines connect to the hard subdiagram are power-counting suppressed, but in pure scalar theories this is not always the case.
The behavior of integrals in the neighborhood of an arbitrary pinch surface is particularly easy to determine in scalar theories. As illustrated in the figure, a generic pinch surface, labelled γ, can be characterized by some number of collinear loop momenta, L γ C , and by L γ S soft loop momenta. Correspondingly, at pinch surface γ, there are N γ S "soft" lines, with vanishing momentum, and N C "jet" lines, each of whose momentum is a fraction x, with 0 < x < 1 of the momentum of one Figure 1. Depiction of a general pinch surface for two-to-two scattering. Shaded blobs represent jet subdiagrams, and the open circle a subdiagram of "soft" lines, whose momenta vanish at the pinch surface. In the center is a "hard" subdiagram consisting of lines off-shell by the order of the momentum transfer. Each line connecting the soft, jet and hard subdiagrams represents an arbitrary number of lines. For the purposes of this discussion, all lines are scalars, and the dashed lines simply represent soft lines attached to jet subdiagrams. Note the possibility that soft lines as well as jet lines attach to the hard subdiagram.
of the external lines. To find the behavior of the integral at pinch surface γ, we introduce a dimensionless scaling variable δ and study the behavior of the integrand and integration volume for δ → 0, where it takes all momenta to the pinch surface. To keep track of dimensions, we label by Q the typical hard-scattering momentum scale, say Q ∼ √ŝ , for 2 → n fixed angle scaling. Now, for soft lines, which vanish in all four components at the pinch surface, we take Jet line momenta, on the other hand, approach a fraction of the corresponding external momentum according to where η µ p is a lightlike vector moving opposite to p µ , with η 2 p = 0, and where p · k ⊥ = η p · k ⊥ = 0. The scalings for these jet line components are then β j ∼ δ Q , With all x i fixed, the integrand times volume element of the remaining integrals in d dimensions then behaves as δ pγ , with p γ a "degree of divergence", given in terms of the number of loops and lines by For p γ > 0, the integral is finite near the pinch surface. When p γ = 0 there is a logarithmic divergence, and when p γ is negative, a power divergence. We note that a generic pinch surface is contained in surfaces of larger dimensionality, in which some subset of loop momenta are finite distances away from these configurations, while others remain close to them. A straightforward analysis of these integrals shows that for all such surfaces, the appropriate scaling follows the same rule [37]. An essential feature is that the pinch surfaces of the integrals that result by keeping only the leading behavior for a given pinch surface itself has (lower-dimensional) pinch surfaces and leading regions that are subsets of the pinch surfaces of the original diagram [40]. It should be noted that this is a feature of fixed-angle scattering, and that in other configurations, in particular for forward scattering, additional power counting analysis is necessary. For a massless cubic scalar theory, the degree of divergence can be arbitrarily negative at high orders, depending on the diagram. This power infrared behavior is the result of the dimensional (superrenormalizable) cubic coupling. In the cases below, for two-to-two fixed angle scattering, we will in fact encounter both logarithmic and power singularities. Although it is not our intention to propose an all-orders treatment of the scalar case, our method will deal with the power as well as logarithmic singularities.
As noted, leading regions are connected in general to other leading regions, and while some are contained within others, others may be disjoint or overlapping. The situation is analogous to the classification of ultraviolet divergences, in which divergent subdiagrams may nest or may overlap with each other. Our aim here is to show how, given any L-loop diagram defined in d = 4 − 2 dimensions, we can construct an integral that is finite in four dimensions, by a suitable subtraction, where I finite ( ) has a finite → 0 limit. The essential result of perturbative ultraviolet renormalization is that sums of products of multiple, nested subtractions produce finite Green functions. It is not necessary to make sequential subtractions involving overlapping regions. A similar structure has been developed in infrared subtraction formalisms, starting as early as Ref. [41]. Following the notation of that paper, we can represent the result as where each product is over a non-empty, ordered set N of approximation operators t a associated with pinch surfaces a, which act to the right. In Refs. [31] and [32], it was shown that sums of nested subtractions, starting from the smallest, most singular regions, can be used to separate infrared singularities from short-distance structure. We shall not review the details of these arguments, but only observe that the pattern starts by making subtractions that match the behavior of the integrand in the most singular regions of momentum space, of the smallest volume, in which the largest numbers of lines approach the light cone (or more generally, the mass shell). The nested operations then act systematically on the resulting terms to remove remaining divergences, by proceeding to subtract the next largest volume, then the next, and so on. This is possible because, as noted above, for fixed-angle scattering the pinch surfaces of the integrals after the action of the approximation operators are subsets of those of the original diagram. For proofs of factorization in gauge theories, the approximations are tailored to match leading behavior, and often at the same time to provide expressions to which the Ward identities of the theory may be applied. Generally, this results in the introduction of new ultraviolet divergences in subtractions, a feature that serves as a basis of resummation [42]. In our examples below, however, we set these considerations aside, and take a pragmatic approach to the identification of subtractions. In particular, at this stage we design subtractions to avoid induced ultraviolet divergences. Here a method introduced at one loop in Ref. [43] will turn out to be useful. This will already be apparent in our first example, the one-loop box diagram, to which we turn as a warm-up exercise in the following subsection.

Subtraction for the one-loop box
As a pedagogical example, we will apply the method of nested subtractions to the massless one-loop scalar box integral, shown in Fig. 2. We write the integral as where the propagator denominators are A j ≡ k 2 j + i0, j = 1 . . . 4, and where the internal momenta are related by with k 4 ≡ k 0 in this notation. The external momenta are taken all incoming, and satisfy with s, t two independent Mandelstam variables. In the following, we will often use the shorthand notation X ijk... = X i + X j + X k + . . .. The integral of Eq. 2.8 has infrared divergences, which fall into the classes of leading regions identified in the previous section. These leading regions are conventionally represented by "reduced diagrams", in which lines that are off-shell at the pinch surface are contracted to points. The eight leading pinch surfaces of the one-loop box fall into two categories, illustrated by the examples of Fig. 3a and b.
First, the box is divergent in the four soft limits k i ∼ δ → 0, for which the leading regions are four disjoint points in loop momentum space, illustrated by Fig.  3a. In terms of the power counting of Eq. (2.4) these regions all have L S = 1, N S = 1, N C = 2, corresponding to logarithmic divergence. Near the point k 2 = 0, for example, we have (Eq. (2.1)), k µ 2 ∼ δ → 0, and the denominators scale as (2.11) We confirm that the integrand tends to  which is of course consistent with Eq. (2.4). The integral of Eq. (2.8) is also divergent in the four collinear limits scale as in Eq. (2.2), In this region, We note the familiar double and single poles in . Our aim in this preliminary discussion is to see in practice how these poles are reproduced systematically in a subtraction formalism for this simple case. The method of nested subtractions introduces counterterms to the integrand designed to remove the four soft and four collinear divergences of the one-loop scalar box. The method removes first singularities from the smallest regions of the integration domain, and proceeds successively to remove the singularities in larger volumes. The regions of the soft singularities are clearly the smallest, since they correspond to points in the integration domain (k i = 0, i = 1 . . . 4) and so they will be removed first. In a soft limit, three of the propagators of the one-loop box are on-shell and one propagator is hard. The collinear singularities extend to larger regions k i = −x i p i , 0 < x i ≤ 1 and, in the method of nested subtractions, they ought to be removed next. In the collinear limits, two propagators are on-shell and two are hard. We note as well that each soft region is an end-point of two collinear regions.
We remove the divergence of the integral in the k 2 → 0 limit by subtracting a function that approximates the singular behavior of the integrand in that limit. We will sometimes refer to this subtraction as a counterterm. We may think of any such counterterm as the result of one approximation operator in a product of subtractions, as illustrated in Eq. (2.7). Each particular approximating operation acts to produce a new integral, which approaches a singular expression like Eq. (2.12) as δ → 0. For some purposes, in particular in proofs of factorization, a choice in which we keep only the terms with leading behavior as δ → 0 is most convenient. To be specific, let us label the subtraction operator for the k 2 → 0 as t S 2 . This operator acts as Clearly, this choice improves the ultraviolet behavior of the resulting expression, by keeping k 2 2 terms in three of the denominators [43]. It also results in a better approximation in the collinear regions, as we shall see below.
In a hopefully clear notation, we label the combination of the original diagram and the particular counterterm defined by Eq. (2.20) as Box R1 (where R1 simply denotes the remainder after the first subtraction). Exhibiting the counterterm explicitly, we have This subtraction is certainly one of the possible choices that guarantees that the integral is free of the soft singularity as k 2 → 0. The counterterm in Eq. (2.21) is chosen according to the prescription of Ref. [43], in which the denominators of eikonal propagators are not linearized. The advantages of the prescription of Ref. [43] are, first, that soft counterterms do not introduce spurious UV divergences and second, that they can be integrated analytically with standard methods. The integrand of Eq (2.21) is still divergent in other regions of the integration domain as, for example, in the remaining k µ i ∼ δ → 0, i = 1, 3, 4 soft limits. We subtract these additional soft singularities sequentially, in the same manner as above. This process is particularly simple because each of the three remaining soft limits requires the denominator A 4 to vanish for a divergent contribution in four dimensions. Indeed, none of the soft subtraction terms have further soft singularities, and all remaining soft divergences are in the first term in Eq. (2.21).
The resulting integral, subtracted for each of its four soft singularities thus has four separate subtractions, and takes the form, It is easy to verify that this integral is not singular at any of the k µ i → 0 soft limits. The subtraction in (2.21) associated with the k 2 = 0 singularity, for example, is simply 1/t times a scalar triangle. When regulated dimensionally, the explicit expression for the subtraction is easily integrated, and the four subtraction terms give In Eq. (2.22), these terms reproduce and cancel all double and single poles in the one-loop box, as given in Eq. (2.17). Evidently, the soft subtractions defined as above reproduce all the collinear as well as the soft singularities for the particular case of the scalar box. Turning our attention to the collinear singular limits, as for example in Eq. (2.14), we easily confirm that no further subtractions are necessary. The straightforward application of our method, however, would remove a remaining collinear singularity term by term, by adding an additional subtraction, determined by the collinear behavior of the soft-subtracted integral (2.22). As noted above, there is some freedom in choosing the subtraction, or counterterm, as long as it matches the singular behavior of the sum of terms in Eq. (2.22), and produces no new leading pinch surface.
Consider the limit in which the loop momentum becomes collinear to external momentum p 1 . For this example, we illustrate one of the forms of collinear counterterms that we shall use below. The subtraction acts by keeping the leading finite (δ 0 ) term in the (two) denominators that are off-shell in this collinear region (A 3 and A 4 ), and the full momentum dependence of the on-shell, collinear denominators (A 1 and A 2 ), along with the leading behavior of each term in the numerator N Box , Eq. (2.22), that defines the sum of soft subtractions, evaluated at the pinch surface, k 1 = xp 1 , 0 < x < 1. Representing the action of the p 1 -collinear approximation by t C 1 , we have, in particular, (2.25) When acting on each of the terms of N Box , however, the resulting integral, which has only two full denominators, is ultraviolet divergent. Here, we shall avoid introducing such induced divergences by adopting a slight variant of the collinear subtraction introduced in Ref. [43]. To be specific, we can introduce an extra factor that approaches unity at the relevant pinch surface, but which regulates ultraviolet behavior. . (2.26) In the nested approach, we apply the same collinear subtractions to the soft subtrac-tion terms. 1 Treating the remaining collinear regions in the same fashion, the full subtraction is We expect, of course, that since the soft subtractions already cancel all singularities, any term-by-term collinear singularities must likewise cancel among themselves. This is indeed the case, because non-zero terms in N Box cancel in the collinear limit for p 1 (where A 1 = A 2 = 0),  [44][45][46][47][48][49][50].
Although in general we would expect to evaluate the remainder with nunerical methods, as an illustration in the one-loop case of Eq. (2.22), we can introduce Feynman parameters, "complete the square" in the loop-momentum and drop numerator terms in odd powers of the loop-momentum, which integrate trivially to zero. We find, We integrate out the loop-momentum, resulting in The integration can be performed by standard methods, yielding the finite result: This is indeed the correct contribution to the finite part of the integral, as is found by comparing Eqs. (2.17) and (2.24). In summary, for our introductory one-loop example, the method of nested subtractions employed here yields the same separation of finite and divergent terms as the method of Ref. [43]. In the following, we will demonstrate that nested subtractions also allow us to treat infrared divergences in non-trivial, two-loop examples.

Application to two-loop scalar integrals
As noted above, it has been shown [31,32] that we can remove the infrared singularities of multi-loop integrals for hard scattering processes with suitable nested subtractions, which we have described in the previous section. However, beyond one-loop in multi-leg amplitudes, we are not aware of a practical construction that realizes this potential in the literature. In this section, we apply for the first time our method of nested subtractions at two loops.
We will focus on two-loop integrals with four external legs which, for light-like external momenta, already have a complicated singular structure. Explicitly, we will test that we can render integrable in d = 4 dimensions the "diagonal-box", the "bubble-box", the "planar double-box" and the "crossed double-box". These integrals represent Feynman diagrams for the scattering of massless scalar particles. In addition, they are the most complicated master integrals which appear in all 2 → 2 scattering processes in massless QCD. We believe that the set of integrals that we examine here serves two purposes: giving a pedagogical introduction to our technique at two loops, and testing it thoroughly in non-trivial applications. In particular, the planar and crossed double-box integrals have poles in the dimensional regulator of the maximum power, 1/ 4 , as they possess all the infrared singularities that are anticipated at two loops. Largely due to their complicated singular structure, the analytic evaluation of the planar and crossed double-box integrals was not amenable to traditional techniques, and was only achieved for the first time when Smirnov [51] and Tausk [52] developed powerful Mellin-Barnes methods.
In this section, we will show that the analytic structure of the 1/ poles of our two-loop examples can be derived in a simple way, by integrating less complicated counterterm integrals. In addition, our counterterm subtractions will render the remainders of the integrands free of any local singularities, and therefore amenable to direct integration methods in exactly d = 4 dimensions.
We begin our discussion with two relatively simple cases, the "diagonal" and "bubble" boxes. In these cases, no more than two nested subtractions are necessary, and the pattern follows the general considerations outlined in the previous section. We will use these examples, however, to illustrate convenient choices of finite parts for collinear subtractions. We then turn to the more complex cases of the planar and nonplanar double boxes. As our first two-loop application, we choose an example with only collinear singularities. Consider the diagonal-box integral, defined as

Subtraction for the diagonal-box
The momenta k i of the propagators are depicted in Fig. 4. One can concretely identify loop momenta with the lines k 1 and k 4 , so that The kinematics of the external momenta p i are, We have taken the p 1 , p 3 momenta to be off-shell. Our study carries through unchanged, however, in the case that one or both of them go on-shell. By inspecting all pinch surfaces, using the power counting of Eq. (2.4), we find that the diagonal-box has only collinear singularities, which we sort according to increasing volumes of their regions: • Two-collinear pairs: C 1||2 C 4||4 , in which the internal momenta k 1 , k 4 are simultaneously parallel to the external momenta p 2 , p 4 respectively. In the notation of Eq. (2.2), we parameterize the loop momenta as 5) in this region, where and • Single-Collinear: where In the C 1||2 C 4||4 limit, the momenta corresponding to We can remove this singularity by introducing a counterterm that subtracts an approximation to the integrand in this singular limit. In the notation of the previous section, we denote this as where t C 2 C 4 represents the "two-collinear pairs" approximation shown in the second equality. Here and below, we employ a notation in which the function inside square brackets is evaluated at the values of momenta shown in the subscript. We note that a factor of θ(x i )θ(1 − x i ), which ensures that along the collinear surface the momentum of a lightlike external line is shared by two internal lines, moving in the same direction, will emerge after loop integrals at these collinear limits. For now, however, the x i may be taken as unconstrained. In the example at hand, we then have , (3.12) where in the second equality, we identify a function that we will encounter in explicit integrals below, The subtraction in Eq. (3.11) removes the double collinear limit, and we take this as our starting point for the removal of the remaining singularities associated with single-collinear limits. We will then turn to the treatment of induced ultraviolet poles.
We now proceed to remove the divergence due to the C 1||2 single-collinear limit in D box | R 1 . In this limit, To treat this region, we introduce an additional counterterm to subtract the behavior of the full integrand of D box | R 1 in Eq. (3.11). The resulting subtraction automatically avoids over-counting the double-collinear limit, which is already removed. In the same notation as Eq. (3.11), we represent the resulting expression as (3.14) In the second equality, we have noted that the second and fourth terms on the right of the first equality cancel, because momentum k 1 (which becomes parallel to p 2 in this limit) flows only through propagator A 5 , according to our choice of routing of the momenta in Eq. (3.3), so that propagators A 3 , A 4 are independent of k 1 . In the subtraction notation, this amounts to t C 2 t C 2 C 4 D = t C 2 C 4 D. We retain denominators A 3 , A 4 inside the square brackets, as a reminder that they appear with denominator A 5 in a one-loop integral, evaluated at fixed k 1 = −x 2 p 2 . We remove the remaining single-collinear singularity C 4||4 from D box | R 2 similarly, giving giving At this stage, we have an integral that is free of all infrared singularities. However, the counterterms that we have introduced are divergent in the ultraviolet limit, just as they were for the one-loop box treated in the previous section. For an analysis carried out purely in dimensional regularization, this would not be a problem, but since, as above, our goal is to derive integrals that can be evaluated numerically, we need them to converge in four dimensions.
As an additional step, therefore, we modify our counterterms, so that they depend on an artificial mass µ in a manner that tames the ultraviolet behavior of the integrand. These are a variant of the subtraction in Eq. (2.26) above, still in the spirit of Ref. [43]. This gives our final expression for a fully-subtracted diagonal box, now finite in four dimensions. The subtractions are in the same pattern as in (3.15), but all integrals are now UV convergent, .
In this expression, we have added to each subtraction term in (3.15) an IR finite adjustment, in which mass dependence is introduced in the denominators that become collinear. Of course, this introduces poles associated with the new denominators. As long as µ is finite, however, these poles produce no new pinches, because lines p 2 and p 4 are lightlike. The full diagonal box equals D box | fin plus the counterterms in (3.16). The integral D box | fin of Eq. (3.16) can now be evaluated numerically (at least in principle) or analytically in exactly d=4 dimensions, since it is free of all divergences.
In fact, for this case, evaluation in dimensional regularization is particularly simple, because in dimensional regularization all mass-independent counterterms include scaleless integrals that vanish. In this way, of the nine terms in Eq. (3.16), only the first, third, fifth and the final term with four massive denominators survive. Of course, by using dimensional regularization we abandon the use of point-by-point cancellation in Eq. (3.16). Nevertheless, it will enable us to confirm the finiteness of the full subtracted form. These integrals will also come in handy in our discussion of mass-dependent integrals in Sec. 4.
We thus proceed to evaluate both the finite part and the singular parts of the original two-loop integral. The latter emerge from the integration of the counterterms, which we will perform in non-integer d = 4 − 2 dimensions. Not surprisingly, the integrations for the counterterms are simpler than the integration of the original integral.
The collinear counterterms that UV-regulate the diagonal box require integrals of the generic form where p is an on-shell external momentum, l and l + p are the momenta of the propagators attached to p and x(l) = − l·ηp p·ηp is the fraction of p that is carried by l. G(x(l)) is the internal subgraph that is attached to the propagators l, l + p, as illustrated in Fig. 5.
The masses m, M in Eq. (3.17) take non-zero values only in the terms of Eq. (3.16) that have been added to regulate the ultraviolet limit of the collinear countert-erms, in which case, m = M = µ. We can treat all these integrals in the same fashion, by first introducing a subintegral that is differential in the momentum fraction x, in terms of which, The integral of Eq. (3.18) at fixed x can be computed easily by using light-cone integration variables and Cauchy's theorem. It reads Therefore, we find that Here, the explicit denominator reduces to a constant whenever the masses m and M are equal (to µ in the case at hand). For example, Eq. (3.21) allows us to evaluate the double-collinear counterterm, for which G = [1/A 5 ], for both the k 1 and k 4 integrals, which result in x 1 ≡ x and x 4 ≡ y integrations, respectively. For this term, we find in this way, , (3.22) where A(x, y) is given by Eq. (3.13). The remaining integrals can then be done in a straightforward manner analytically in terms of rank-two polylogarithmic functions. We define the kinematic variables, and the scale µ-dependent logarithmic function, In these terms, we find which will appear as the coefficient of the double-pole counterterm for the diagonal box.
The integration of the single-collinear counterterms requires one more integral, which arises from the off-shell triangle that is opposite to the counterterm. Consider, for example, the triangle integral corresponding to loop momentum l = k 1 in Fig.  4, evaluated at k 4 = x 4 p 4 , as it appears in the third subtraction term in Eq. (3.16). The relevant integral now has three propagators, In the computation of the finite part, Eq. (3.16), the "single-collinear" term from (3.26) contributes a single pole from the expansion of A −1− that can also be computed easily in closed form, in terms of rank-three polylogarithmic functions, Quite generally, for collinear limits in two-loop diagrams, the kernel G(x) corresponds to a one-loop or a tree subgraph. G(x) is therefore a rational function of x for trees and a function of polylogarithms of x for one-loop subgraphs. The method of nested subtractions leads to functions G(x), which contain no other subdivergences.
Therefore, integrands of the form of Eq. (3.21) can be expanded as a Taylor series in , whose coefficients can be integrated either analytically or, alternatively, numerically.
Having discussed the integration of the divergent counterterms, we return to the evaluation of the finite remainder of Eq. (3.16), which in this case can be performed in exactly four dimensions. We envisage that finite remainders of two-loop integrals after the application of nested subtractions are integrated numerically in momentum space, after appropriate contour deformations away from non-pinched singularities are applied. We emphasize again that the development of an efficient numerical method requires further study, a problem that we will not address here. A method that achieves this purpose for generic multi-loop integrals has been presented in Ref [48].
For the full finite part, including the original diagram, we have from the above, with A(x, y) given in Eq. (3.13). The first term in brackets on the right of the first equality is the full diagram, the second term is the result of single-collinear subtractions, and the third term is from the double-collinear subtraction. This expression is manifestly finite in d = 4 dimensions and can also be easily integrated analytically in terms of logarithms and polylogarithms (see, for example, Appendix D of Ref. [54]). The analytic result for the finite remainder of the diagonal-box integral reads (3.29) The limit of m 1 , m 2 → 0 can be taken smoothly in Eq. (3.29). We have checked that in that limit the above result, when combined with the integrated counterterms, agrees with the analytical results of Refs. [55,56] for m 1 = m 3 = 0. Finally, we would like to comment that the terms proportional to log 3 (µ) in Eq. (3.29), log 2 (µ) in Eq. (3.25) and log(µ) in Eq. (3.27) all cancel. This is in accordance with expectations, since the strongest singularity is due to two-collinear pairs capable of producing at most 1/ 2 poles and consequently at most log 2 (µ) terms in the finite part of the integral. We now consider the bubble-box two-loop integral, which has collinear and soft IR divergences, in addition to a UV-divergent subdiagram. We give fewer details in this case, because the necessary reasoning is very similar to the one-loop and diagonal box diagrams, and is easily seen to give the desired form.

Subtraction for the bubble-box
The bubble-box diagram was computed analytically in Ref [56]. It is defined as with The momenta k i of the propagators are depicted in Fig. 6. One can concretely choose The kinematics of the external momenta p i that we choose are: where p 2 3 and p 2 4 may or may not be on-shell. The bubble-box integral possesses an ultraviolet singularity due to a one-loop (bubble) subgraph. In addition to the ultraviolet divergence, we encounter one "single-soft" (S 2 , in the notation of the previous section) and two "single-collinear" (C k 1 ||p 1 , C k 3 ||p 2 ) singularities. When momentum p 4 is lightlike, the diagram also has a "two-loop collinear" pinch surface, where both loop momenta are parallel to p 4 and lines k 1 , k 4 and k 5 , share the momentum p 4 , all moving toward the final state. A similar pinch surface arises when p 3 is lightlike. It is easy to check from Eq. (2.4), however, that, because they involve two collinear loops and only three collinear lines, these pinch surfaces do not produce a singularity in the integral in four dimensions. Therefore, it is not necessary to introduce counterterms for two-loop collinear pinch surfaces of this diagram in four dimensions, although in general it would be necessary in three dimensions.
We subtract approximations for the soft and collinear singularities sequentially, which automatically eliminates double-counting. As above, this leads to a remainder that is free of infrared singularities. We implement collinear counterterms in the same manner as in the case of the diagonal box, explicitly, The above subtractions suffice to render the integral finite also in the ultraviolet limit. Indeed, the original integrand and the soft counter-term in the first line have the same behaviour in the ultraviolet and cancel each other in that limit. Given that the bubble-box integral is known analytically [56] and that we can integrate simply the soft and collinear counterterms, it is straightforward to check that indeed the remainder of Eq. (3.34) has no 1/ poles and is finite in d = 4 dimensions. Specifically, for t = −y s, with 0 ≤ y ≤ 1, we find that The logarithmic term here that depends on the scale µ originates from the integration of the collinear counterterms. The S 12 Nielsen polylogarithm is defined in the appendix A.1. The two-loop box is defined as: where the momenta of the propagators are depicted in Fig. 1.
1 Figure 7. The two-loop double-box

Subtraction for the two-loop planar double-box integral
We now consider the planar double-box two-loop integral, which was computed for the first time analytically in Ref. [51]. It is defined as The momenta k i of the propagators are depicted in Fig. 7. One can concretely choose k 1 = l, k 2 = l + p 1 , k 3 = l + p 12 , k 4 = k + p 12 , k 5 = k + p 123 , k 6 = k, k 7 = k − l. For later use, we also define a generalisation of the scalar integral, with an arbitrary numerator N (k 2 , k 5 ) The double-box with on-shell external lines has infrared singularities, which are presented below, ordered according to increasing volumes of their regions: • Double-soft: isolated points in both loop momentum spaces, where some line has a vanishing momentum for each loop. We label these as S i S j , in which the internal lines i, j become soft.
• Soft-Collinear: an isolated point in one of the two loop spaces where one line has vanishing momentum, and a line segment of the other loop collinear to one of the external momenta. These are labelled as S i C j||l , in which the internal particle i is soft and the internal particle j is parallel to the external particle l. Note that at these pinch surfaces there is always a two-fold ambiguity in the choice of collinear line j.
• Two-collinear pairs: line segments for both loops, collinear to different external momenta, labelled C i||l C j||m , in which the internal particles i, l are parallel to the external particles k, m respectively. Two-loop-collinear: segments in which the two loops are both collinear to the same external momentum. These are identified C ijk||m , in which the internal particles i, j, k are all parallel to the external particle m. Again, there is a trivial ambiguity in the choice of lines i, j, k, because four lines go on-shell.
• Single-Soft: an isolated point in one loop, and unconstrained momentum in the other, identified by S i , in which the internal particle i is soft.
• Single-Collinear: a collinear segment in one loop, and unconstrained momentum in the other, identified by C i||l in which the internal particle i is parallel to the external particle l.
Near each of these pinch surfaces, the integral is logarithmically divergent, as may be verified by applying power counting according to Eq. (2.4).
As is characteristic of integrals with both collinear and infrared divergences, we expect both double and single poles from each loop integral, and indeed, these singularities yield up to 1/ 4 poles. To put our subtractions in context, we recall first the explicit form of the leading and next-to-leading poles of the planar box. For s > 0, t < 0, we have [51] P box [1] = Γ 2 (1 + ) In contrast to the two-loop examples of the diagonal and bubble boxes above, a direct construction of all possible nested subtractions, as in Eq. (2.7) would involve nestings with as many as five approximation operators and well over a hundred terms. We will show below, however, that the IR singularities of the planar box, associated with the pinch surfaces listed above, can be subtracted with a much smaller set of counterterms. As in the previous examples, we construct these inductively, working from smaller to larger regions. For each such region, we let the diagram, including subtractions for smaller regions, determine the subtraction for the region at hand. We thus start the process with the double soft configurations (isolated points in both loops). We design subtractions for these regions to cancel the leading poles. There are actually two types of double-soft points. In the first class, illustrated by Fig. 8(a,b), exactly two lines have vanishing momenta, while in the second class, illustrated by Fig. 8(c), three lines are joined together with zero momentum. Notice that in the latter class, one of the soft lines attaches to the "hard part" of the diagram, the two lines that are off-shell by invariant s at the pinch surface. A configuration like this is never leading in a gauge theory diagram, where the three-boson vertex is associated with a momentum factor, but it has power counting zero for φ 3 in four dimensions, and so must be included in the analysis of the planar double box.
We choose to subtract first for the former class, consisting of S 2 S 5 , S 2 S 7 , S 5 S 7 , S 2 S 4 , S 2 S 6 , S 5 S 1 and S 5 S 3 , all double-soft divergences. In these limits, as reflected in the reduced diagrams of Fig. 8(a) and (b), six of the seven propagators are on-shell and one propagator is hard: The subtractions are defined to simply replace each off-shell momentum by the appropriate invariant. For example, in the S 2 S 5 subtraction, this is accomplished by a factor A 7 /t. We label the planar box subtracted for each of these singularities as P box [N 1 ], with the numerator N in equation (3.40) given by where again A 257 = A 2 + A 5 + A 7 and so on. It is easy (and important) to check that the subtraction term associated with a given S i S j is power-counting finite in the regions around the other points S k S l in the list of Eq. (3.42). The planar box P box [N 1 ] is therefore free of all the above seven double-soft singularities. We have not, however, yet dealt with the two singular points where five propagators are on shell, and two lines are off shell. We label these points as The remaining effect of these regions can be found by direct calculation. Indeed, the subtraction terms in P box [N 1 ] can be evaluated by standard methods, and combined with the results for the full diagram, given at the first two powers in 1/ by Eq. (3.41). The complete result is In this expression, leading powers are still present, so we must certainly make further subtractions. To do so, we follow the approach mentioned above, and determine the necessary counterterms by studying the behavior of the full subtracted diagram P box [N 1 ], rather than that of the original diagram, P box [1]. We find that two counterterms with factors of 1/s in P box [N 1 ] are singular in each of the regions of Eq. (3.44). The net result is that the regions of (3.44) are actually already subtracted twice in P box [N 1 ]. Counterterms must thus add these regions back, to avoid double counting, rather than to make an additional subtraction. We thus arrive at the next stage in the subtracted planar box, labelled P box [N 2 ], with a numerator This expression is free of all double-soft singularities. Upon integration including the counterterms, it gives to order −2 the explicit form As anticipated, after the subtraction of the approximations to the integrand in all double-soft singular limits, the 1/ 4 pole is cancelled. Somewhat unexpectedly, the 1/ 3 is also cancelled. Evidently, our double-soft counterterms have removed further divergences, beyond those for which they were originally defined. In any case, P box [N 2 ] will be the starting point for the next round in the construction. Next in our ordering of subtractions of "increasing volumes" comes the subtraction of soft-collinear singularities, which may be labelled as and S i C k 5 ||p 3 , i = 2, 3, S i C k 5 ||p 4 , i = 1, 2.
We find that in these eight, potentially singular, soft-collinear limits the numerator N 2 , Eq. (3.46) that we have generated after the subtraction of the double-soft divergences vanishes: Therefore, we do not need to add any further counterterms in N 2 to remove softcollinear divergences. This result is clearly a reflection of the lack of 1/ 3 poles in P box [N 2 ]. We are ready to move on to the remaining regions, which can give at most a double pole. We now proceed with the subtraction of divergences of two-collinear pairs (in which, five propagators are on-shell). We find that the numerator N 2 vanishes in the limits, but it is finite in the remaining two limits of this type, We therefore add two counterterms to N 2 , which subtract the behavior in these limits, The integral P box [N 3 ] is then free of divergences associated with two-collinear pair singularities as well. No new singularities from smaller regions are produced, because every such singularity involves the vanishing of at least one of the pair of denominators A 1 , A 4 and A 3 , A 6 . Of a similar order are singularities due to four two-loop-collinear limits, in which N 3 scales as To remove them, we add four additional (positive) counterterms to N 3 , leading to The integral P [N 4 ] is free of all singularities due to soft/collinear configurations that involve both loop momenta, although on evaluation, it still retains double poles, This is because the resulting integrand, is still divergent from pinch surfaces in which a single momentum is soft and/or collinear to an external momentum. To remove these, we will use insight from the example of the one-loop box treated in the previous section, turning our attention to the full integrand. We deal first with the singularities in the single-soft S 2 and S 5 limits. In those, the integrand is approximated by where the resulting counterterms are given by The counterterms in F (1s) P box due to single-soft singular limits are actually straightforward to integrate. For example, the single-soft counterterm describes a one-loop subgraph with the propagators A 4 , A 5 , A 6 and A 7 , which does not contain the soft propagator, and is evaluated at a fixed value k 2 = 0. This operation factorizes the counterterm into the product of two one-loop integrals, a singular one-loop triangle (containing the propagators A 1 , A 2 , A 3 ) and an one-loop box (containing the propagators A 4 , A 5 , A 6 , A 7 ) with a numerator found by setting A 1 , A 2 and A 3 to zero in Eq. (3.53), This expression is precisely the full numerator for the subtracted one-loop box, Eq.
(2.28), which renders the k 5 integral fully finite in this case. This is not an accident, but a consequence of the fact that singularities stronger than single-soft have been removed earlier in our sequence of nested subtractions.
With the single-soft singularities subtracted, the sum of F (2) P box + F (1s) P box , Eqs. (3.55) and (3.57), should only have single-collinear singularities, which yield single 1/ poles. Indeed, upon integration we find Explicitly, the integrand of F P box is singular only in the C k 2 ||p 1 , C k 2 ||p 2 , C k 5 ||p 3 , C k 5 ||p 4 single collinear limits. In these limits, we can approximate the integrand by our final counterterm, F (1c) P box defined by (3.61) in the remaining collinear limits. This counterterm is given by an expression which, although a little long, is a straightforward generalization of the UV-finite diagonalloop box subtractions in Eq. (3.16), In the above, collinear counterterms have been introduced in the same manner as for the regulated diagonal-box, Eq. (3.16), with fractional momenta defined by We have now specified our final counterterm and the integrand is free of all singularities. We have checked that as an analytic expression, the integral of Eq. (3.64) is finite in four dimensions, We have performed the analytic integration of the counterterms in Eq. (3.65) in a straightforward manner as we will explain shortly.
To be specific, for s > 0 and 0 < y ≡ − t s < 1 we find

67)
, The finiteness of this result confirms that the counterterms reproduce all singular behavior of the analytic expression of the planar double-box P box [1], as given in Ref. [51]. As anticipated, the counterterms are simpler than the full integral. As we have noted above, the single-soft counterterms reduce to the product of two one-loop integrals. The counterterms in F (2) P box , due to double-soft, soft-collinear, double-collinear pairs and two-loop-collinear singularities contain two-loop integrals with at most six propagators. We could reduce all six-propagator integrals, following the algorithm of Ref. [56], to the diagonal-box and bubble-box master integrals, which are simpler than the original double-box integral and whose evaluation we have described earlier.
We anticipate that counterterms for such singularities (double-soft, soft-collinear, two loop collinear and double-collinear-pairs) of an arbitrary two-loop integral can be expressed in terms of integrals with at most six propagators, since this is the maximum number of propagators that can become on-shell in such configurations.
The integration of the counterterms in F (1c) P box , Eq. (3.62), due to single-collinear limits is slightly more involved, since it requires a one-dimensional convolution as described in Eq. (3.21) with a kernel that is an one-loop subdiagram. We first reduce the one-loop subdiagram to the one-loop off-shell box and bubble master integrals of the appendix A.2. The remaining one-dimensional integrations yield Nielsen polylogarithms S np (y) (appendix A.1) of uniform weight n + p ≤ 4. We have calculated the required integrals algebraically by comparing a few first terms in their series expansion around y = 0, 1, ±∞ and a general ansatz of such Nielsen polylogarithms.
While here we studied a single two-loop diagram, in a calculation of a physical two-loop amplitude it is anticipated that sums of collinear limits from all diagrams will factorize in terms of splitting functions times one-loop or tree amplitudes, simplifying the convolutions into products. We now detail the construction of local counterterms for the two-loop crossed double-box, which is depicted in Fg. 9. The external momenta satisfy,

Subtraction for the two-loop crossed double-box integral
For convenience below, and as for the planar box, we introduce the integral with an arbitrary numerator N , and define The internal momenta can be chosen as: We are interested in removing the infrared singularities of X box [1], which was computed analytically for the first time in Ref. [52]. We follow the same procedure as for the planar double-box and previous examples. Namely, we remove the singularities iteratively, following the order: double-soft, soft-collinear, two-collinear pairs/twoloop-collinear, single-soft and single-collinear. Of the sixteen distinguishable double-soft regions of the crossed box, two have the property that three lines are forced to zero momentum. In the spirit of our discussion for the planar box, we can label these zero-dimensional pinch surfaces by any two of the three lines that are coupled at a three-point vertex and have vanishing momentum. We will call them S 1 S 7 and S 3 S 6 , where we understand that these two configurations imply as well that S 5 and S 4 carry vanishing momentum, respectively. The region S 1 S 7 is illustrated in Fig. 10a. At configurations of the cross-box like this, we encounter an additional complication, due to the presence of power-like (rather than logarithmic) double-soft singularities, (3.74) In order to remove these power-like singularities, we must introduce a counterterm for which the numerator scales as δ 2 in the corresponding limits. We achieve this with which covers both cases. Notice that the integral of N 1 includes the original diagram. All other double-soft singularities are logarithmic and we can proceed to subtract them in an analogous way as in the planar double-box example. However, a subtlety remains, due to the original power-like nature of the S 1 S 7 , S 3 S 6 singularities. As we have seen for the planar box, if one encounters only logarithmic singularities in the process of performing nested subtractions, the integrands of counterterms that are introduced later for higher-dimensional pinch surfaces vanish in all limits of the lower-dimensional pinch surfaces that have been treated earlier. In other words, once a logarithmic singularity is removed at a certain step it is not introduced back spuriously with the construction of a subsequent counterterm for a different singularity. This mechanism is not automatic in the presence of a power-like singularity, which can influence the construction of counterterms for divergences which at first sight seem unrelated to it.
To illustrate this issue, let us consider the S 2 S 7 log-like singularity, in which: To remove the S 2 S 7 singularity, we could add a counterterm to N 1 , such as However, this new counterterm has spoiled the cancelation of the S 1 S 7 limit which was achieved at the first step with N 1 . Indeed, Note that there is no such problem for the region of the S 3 S 6 power-like singularity, where We can, however, easily construct a suitable counterterm that removes the S 2 S 7 singularity without introducing a spurious S 1 S 7 singularity. Such a counterterm is found in the numerator: The double-soft S 2 S 4 , S 2 S 5 and S 2 S 6 regions can be treated similarly. Counterterms that cancel this set of singularities are given by N 3 − 1, where There are two more double-soft singular limits that we have not treated so far, S 4 S 7 and S 5 S 6 . These limits leave the same propagator (A 2 ) hard, and set all other propagators on-shell, Now that we only have one hard propagator at our disposal, no single fraction with A 2 in the numerator divided by a single invariant will serve to subtract two different singularities. Instead, we introduce a counterterm that interpolates between the values of A 2 at the two singular configurations, The new counterterm, N 4 vanishes in both S 4 S 7 and S 5 S 6 limits. Unfortunately, it does not vanish fast enough in the S 1 S 7 and S 3 S 6 limits where X box [1] develops power-like singularities. We have In either of the above limits, A 2 ∼ δ and only one of A 1 or A 3 tend to the Mandelstam variable s. We can therefore modify our counterterm as follows: The integral X box [N 4 ] is now free of all double-soft singularities. We also find that is free of all soft-collinear singularities, as confirmed by explicit integration. We therefore proceed with the subtraction of two-collinear pairs/two-loop-collinear types of singularities. These singular limits do not pose any special challenges and they are subtracted along the lines of our planar double-box example. We find that the integral X box [N 5 ] with numerator is free of all singularities associated with two independent loop momenta pinched in a special kinematic configuration (soft or collinear). Finally, we need to remove the singularities due to single-soft and single-collinear limits. After these final subtractions, we find that the following integrand is free of all singularities: where, following the notation of the planar double box, (3.89) and F (1c) (3.90) In the above, B i = A i − µ 2 . Upon direct analytic integration, using the integration techniques described in the previous section for the counterterms, and the analytic result of [52] for the crossed double-box integral, we verify that Specifically, for s > 0 and y ≡ −t/s ∈ [0, 1], we find (3.97) The integration of the counterterms was performed using the same techniques as in the case of the planar double-box. A notable difference occurred in the integration of the collinear counterterms. In the case of the crossed double-box, integrals which do not have a representation in terms of Nielsen polylogarithms with a simple argument S np (y) emerge 2 . However, we have observed that the linear combination which is required in the collinear counterterm can be expressed in terms of Nielsen polylogarithms in our simple basis S np (y). Specifically, we find that

Small mass expansions
In the previous section, we rendered finite integrals that were computed in dimensional regularisation. Dimensional regularisation, however, is not an essential element; our method is in principle applicable to any infrared regulator. Infrared divergences can also be regulated by a small mass parameter. With mass regularisation, the integration over the mass-divergent regions yields logarithms that become infinite in the massless limit. The mass regulator can be artificial or physical. For example, the physical mass of the bottom-quark in processes for the production of Higgs bosons acts as a regulator for some of the infrared divergences. The logarithmic dependence of the corresponding amplitudes is of a high phenomenological interest. In this section, we will use the method of nested subtractions in order to derive simply the asymptotic behavior of certain Feynman integrals in a small-mass limit. Consider a loop integral, represented schematically as which depends on a small mass-parameter m, appearing in denominators of the standard form, k 2 i − m 2 + i0. If we take the zero-mass limit, the integral develops new infrared divergences, which are not present for finite values of the mass. In general, these appear as logarithmic corrections in the mass, which result from regions where the values of denominators are actually larger than the mass: m 2 ≤ k 2 i ≤ Q, with Q some scale fixed by the invariants. For our discussion, all invariants are of the same order. For fixed-angle scattering this is the case, and logarithmic integrals can be identified by the simple power-counting rules described in Sec. 2. At the same time, the integral can receive finite contributions from regions where for one or more denominators k 2 i = O(m 2 ). We would like to find a systematic method to isolate both the logarithmic mass dependence, and contributions that are finite for small, but nonzero mass.
To this end, we follow the method of nested subtractions, and construct an approximation f approx (k i , m) of the integrand in all the limits that become singular as m → 0. As indicated above, these limits can be identified using power counting techniques. Then, we can use these approximations to construct counterterms, keeping those mass-dependent terms that dominate each denominator in the region for which the counterterm is designed to approximate the full integral. Thus, in general, our counterterms retain mass dependence, (4. 2) The leading mass singularities are now found from the integral over the approximated integrand in the first term of the right-hand side of this relation. Because we keep all leading mass dependence, finite terms associated with momenta of the order of the mass may remain in one or more of the counterterms included in f approx . The second term in Eq. (4.2) is well-defined in the m = 0 limit, and additional finite terms, including important kinematic dependence, appear in general in It is important to note that a naive Taylor expansion of the second term in Eq. (4.2) beyond the leading order will not account for terms of O(m 1 ), which vanish as a power, but may be multiplied by logarithms. Because we drop terms that are nonleading of order in O(m) in denominators, we can only ensure the cancellation of the leading power of m in the second term, and will in general miss contributions of the form m log(m). In the following, using these ideas, we will show how this approach can be used to derive the small-mass dependence of one-and two-loop integrals.

One-loop massive triangle
Consider as our first example the scalar one-loop massive triangle of Fig. 11, taken with two light-like external lines, p 1 and p 2 here. This is a rather simple integral, which hardly needs special treatment. It is divergent in the zero-mass limit, however, and is useful to illustrate nested subtractions at finite mass. We note that this analysis depends on having, as in this case, only a single nonzero infrared-scale mass. If, for example, we were to evaluate this diagram with one or more timelike external line near the mass shell, the reasoning below would require further analysis. The integral in question is straightforward to evaluate by standard methods and is given by where the B's are denominators with regulating masses, Referring to the figure, we choose k 2 as the loop momentum, so that the remaining momenta satisfy, with p 2 1 = p 2 2 = 0, This integral is finite in d = 4 dimensions but is logarithmically divergent as m → 0: To illustrate our method applied in four dimensions with regulator masses, we construct counterterms that isolate this leading m → 0 behavior, to the level of finite terms. Of course, with an answer as simple as Eq. (4.7), the counterterms will be at least as complex as the original diagram. Our interest, however, is in the counterterm integrals themselves, because they may appear in more complex diagrams, where the asymptotic behavior of the full integral may be much more complex.
We now proceed to construct an integrand that is free of mass singularities. As is well known, and easy to confirm, the vertex with two lightlike external lines has one soft and two collinear leading pinch surfaces. The singularity of lowest dimension corresponds to the the soft limit, and this is where we begin.
In the massless limit, the soft singularity is associated with vanishing k 2 . In the presence of a mass regulator, the integrand is enhanced in the region (4.8) In this limit, We thus study the effect of a counterterm for this soft limit, subtracted from the original diagram, (4.4), at the level of the integrand, (4.10) In defining the subtraction, we have used that, for the soft region defined as in Eq. (4.8), the invariant masses of the "collinear" lines, k 1 and k 3 scale as √ mQ, with Q ∼ √ s, so that k 2 i m 2 , i = 1, 3. The counterterm is thus a good approximation to the full integral in the entire soft region. In the soft region, it reproduces both logarithmic m dependence and finite remainders. Because we are calculating in a superrenormalizable theory, the behavior from infrared regions gives in fact the full leading-power behavior of the diagram.
This would be the entire story if the soft configuration were the only pinch surface of the massless triangle diagram. As we know, however, it is not, and both the full integral and the counterterm as defined in (4.10) have pinches when the loop momentum k 2 is collinear to the external momentum p 1 or p 2 . The diagram and its soft counterterm are not guaranteed to behave exactly the same in this region. We thus expect in general a leading-power (m 0 ) contribution to I R S , and indeed, by evaluating the integral that defines I R S in (4.10), we find Evidently, the soft subtraction generates the logarithmic m-dependence of the full integral, but cannot reproduce its finite part part without aid of collinear subtractions.
The collinear pinch surface for k 2 aligned with p 1 is conveniently parameterized in the light-cone form, with Following our iterative procedure, the full, soft-subtracted triangle, Eq. (4.10) is used to define the counterterm for the collinear region C k 2 ||p 1 , using the behavior of its integrand in that region, which we may represent as with x 1 defined as in Eq. (4.13). Analogously, a collinear singularity appears in I R S when k 2 = x 2 p 2 with x 2 = 2k 2 ·η 2 2p 2 ·η 2 and m → 0. Note that although approximating k 2 3 = sx 1 produces ultraviolet singularities in each of the two terms on the right-hand side of the first equality, these divergences cancel when the two terms are combined. Thus, we may stay in four dimensions for the entire calculation, without modifying our collinear subtractions.
We now remove the collinear mass singularities by introducing to the integrand of I R S two collinear counterterms, and define the counterterms also approximate the full integrand in all regions that give leadingpower dependence. We may therefore take the limit m → 0 inside the integrand in the second integral, where it manifestly vanishes, as anticipated in Eq. (4.3). Note that although the terms proportional to m 2 appear to be multiplied by powerdivergent integrals, there is no ambiguity, because each of these two terms is the sum of collinear subtractions for the full diagram with its soft subtraction, which cancel independently in the m → 0 limit. Indeed, as we have seen, the zero-mass limit may be taken before or after the integrals.
In summary, we may write for the triangle diagram, which agrees with the result of Eq. (4.7). Our subtraction procedure has thus succeeded in reproducing the infrared behavior of this diagram, including constants, using a mass regulator. We emphasize again that although in this example the subtractions provide only a roundabout derivation of a simple result, the same relatively simple subtractions can appear in diagrams with many more lines. For our next example, we consider a slightly more complex integral. This is the one-loop scalar-box integral with a single off-shell external line. To anticipate, we will find it useful in this case to go back to dimensional regularization, but keeping the masses, with where as usual The momenta in the loop satisfy the conservation relations k 2 = k 1 + p 1 , k 3 = k 1 + p 12 , k 4 = k 1 + p 123 . (4.20) The external momenta satisfy p 1234 = 0 and p 2 1 = p 2 2 = p 2 3 = 0, p 2 12 ≡ s, p 2 23 = t, p 2 123 = M 2 . This integral contributes to the one-loop amplitude for gg → Hg, and it was computed for the first time in Ref. [57]. We will see that dimensional regularization will allow us to cancel ultraviolet divergences associated with collinear subtractions, and will facilitate the calculation of the physical mass-dependence of the diagram. In the m → 0 limit we have two soft, S 2 and S 3 , and three collinear, C k 1 ||p 1 , C k 4 ||p 3 , and C k 2 ||p 2 , singular limits when m → 0. Following our standard approach, we subtract first the soft and then the collinear limits. We note that after the soft subtractions the numerator of the subtracted diagram vanishes in the limit where the loop momentum is collinear to external momentum p 2 , C k 2 ||p 2 . Additional subtractions are therefore required only for the regions C k 1 ||p 1 and C k 4 ||p 3 .
This procedure yields an integrand that is free of leading-power mass singularities, (4.22) In this expression we have chosen x 1 and x 3 as the physical fractional momenta carried in the collinear limits, Here the η i are, as usual, arbitrary light-like vectors that are not collinear to the corresponding p i . In Eq. (4.22), the collinear subtractions induce ultraviolet divergences in four dimensions, which for this discussion we take to be regularized dimensionally. At leading order in the mass expansion, we have for J R given in case of the previous subsection, using dimensional regularization to control ultraviolet divergences that appear in intermediate steps of the calculation.
The integral we study is, The momentum assignments, k i , of the propagators are depicted in Fig. 4. One can choose, for example, k 1 = l + p 1 , k 2 = l + p 12 , k 3 = k + p 123 , k 4 = k, k 5 = k − l . (4.29) The kinematics of the external momenta p i are given by Eq. (4.27) is a master integral for the production of one or two Higgs bosons at hadron colliders. We have studied the infrared singularities of the massless diagonal-box in section 3.1, which becomes divergent in the double-collinear C 1||2 C 4||4 limit and in the two single-collinear limits, C 1||2 and C 4||4 . In the massive diagonal-box that we consider here, these singularities are screened for finite values of the mass and emerge only as we take the limit m → 0. Upon integration, these limits are responsible for generating the logarithmic mass singularities log n (m), n = 1, 2. Following steps analogous to those for the massless diagonal box in section 3.1, we construct counterterms that remove the mass singularities. Up to corrections that vanish with the mass, the fully subtracted integral is equated to its dimensionally-regulated massless limit, Let us now compare the right-hand side of Eq. (4.32) and Eq. (3.16). We observe that if we set the artificial scale µ to the physical mass µ = m in (3.16), the two expressions differ by terms that integrate to zero within dimensional regularisation. We therefore arrive at the following result for the massive diagonal-box, where D box | fin (m) is given by Eq. (3.29). We have checked this result against numerical evaluations in the Euclidean region of results in Ref. [58] after setting p 2 1 = 0 in our expression. We have also checked the coefficients of the logarithmic expansion (in the p 2 1 = 0 case) against the numerical results of an asymptotic expansion that is performed using the program of Refs [59,60], based on the strategy of regions method.

Conclusions
We have used general properties of the infrared behavior of perturbative amplitudes to develop a systematic, iterative method for the evaluation of the infrared singularities in covariant perturbation theory diagrams. As presented here, the method applies directly to fixed-angle scattering and particle production of massless and massive particles. We assumed for this analysis a single hard scale. The output for each diagram is a finite remainder, plus a sum of infrared-sensitive integrals that are generally simpler than the original diagram, and which depend on fewer external momenta.
The method proceeds through the construction of counterterms that approximate and cancel divergent behavior of integrands at a local, point-by-point, level. This approach reflects the nested structure of infrared singularities, which extends to all orders. Compared to a full "forest-like" generation of subtractions, we found that an iterative procedure, extending from regions of lower to higher dimensions, produced a very manageable set of subtractions, even in diagrams with many singular limits, such as the crossed box.
We discussed applications to a number of one-and two-loop diagrams with four external lines. The method is much more general, however, and may be particularly useful for isolating infrared poles in dimensionally-regulated multi-particle loop amplitudes.
We have shown that in massless diagrams, the procedure results in a finite remainder, which can in principle be evaluated numerically in four dimensions. It can be applied as well to diagrams with a small mass on internal lines, in which case it can be used to derive analytic expressions for leading-power mass dependence, including mass-independent terms. In this context, we studied a number of known cases, and have also derived, and numerically verified, a new result for the two-loop diagonal box diagram.
The procedure outlined in this paper can be implemented algorithmically, and as such we believe it has potential for master integrals beyond two loops, and also for applications to cross sections.