Three-loop master integrals for ladder-box diagrams with one massive leg

The three-loop master integrals for ladder-box diagrams with one massive leg are computed from an eighty-five by eighty-five system of differential equations, solved by means of Magnus exponential. The results of the considered box-type integrals, as well as of the tower of vertex- and bubble-type master integrals associated to subtopologies, are given as a Taylor series expansion in the dimensional regulator parameter epsilon = (4-d)/2. The coefficients of the series are expressed in terms of uniform weight combinations of multiple polylogarithms and transcendental constants up to weight six. The considered integrals enter the next-to-next-to-next-to-leading order virtual corrections to scattering processes like the three-jet production mediated by vector boson decay, V* ->jjj, as well as the Higgs plus one-jet production in gluon fusion, pp ->Hj.


Introduction
Feynman integrals are the building blocks for describing particle interactions beyond the tree approximation in perturbation theory. When their direct integration becomes prohibitive, because of either the number of loops or of external legs, the evaluation of Feynman integrals can be addressed in two phases, namely the decomposition in terms of basic integrals, followed by the evaluation of the latter, to be considered as an easier problem.
Within the continuous dimensional regularization scheme, Feynman integrals fulfill identities that fall in the category of the general class of integration-by-parts relations [1][2][3]. Such relations can be exploited in order to identify a set of independent integrals, dubbed master integrals (MI's), that can be used as a basis of functions for the virtual contributions to scattering amplitudes. Therefore, for any given process, it is possible to identify a minimal set of MI's, such that all Feynman integrals contributing to it admit a representation in terms of the chosen basis of MI's.
Scattering amplitudes are combinations of Feynman diagrams. Therefore, it is not unnatural that global properties of scattering amplitudes can be exploited to achieve the decomposition of the whole amplitude in terms of MI's. In fact, unitarity and factorization, both at the integral and integrand level, become suitable tools for the simultaneous decomposition of group of diagrams in terms of MI's [4][5][6] (see [7,8] for reviews). Decomposing the whole amplitude, rather than the individual diagrams, has the advantage that at each step of the computation the coefficients of the decomposition reflect the symmetries of the amplitude. Whereas, the decomposition of each integral separately may not carry such information, with the drawback of introducing spurious terms at the intermediate steps of the calculation, which necessarily disappear from the total result.
The advantage of decomposing amplitudes, rather than individual integrals, in terms of MI's clearly emerged in the last decade, as it has been the driving principle determining the breakthrough in the evaluation of next-to-leading order (NLO) virtual corrections to multi-leg processes (see for instance [9] and references therein). Activities are ongoing to extend the underlying ideas to higher orders [10][11][12][13].
The evaluation of MI's usually proceeds one-by-one in a bottom up approach, starting from the less complicated integrals and systematically enriching their structure by increasing the number of internal and external lines. MI's are functions of the kinematic invariants built with the external momenta and of the masses of the involved particles. Remarkably, the integration-by-parts (IBP) relations imply that the MI's obey linear systems of firstorder differential equations (DE's) in the kinematic invariants, which can be used for the determination of their actual analytic expressions [14][15][16].
The simplicity of the whole amplitude with respect to the individual integrals is reflected also in the fact that the analytic functions present in each MI do not necessarily appear in the final result, obtained after combining coefficients and MI's altogether. Remarkably, Hopf algebra structure (of multiple polylogarithms) can be used to simplify complicated expressions for multi-loop amplitudes [17], or global symmetries can be used to constrain the minimal set of variables an amplitude may depend on [18]. In other words, integrals depend on the kinematic invariants enforced by their external topology and masses, but amplitudes may depend on functions of such invariants, which resolve the symmetries that are not manifest at the diagrammatic level.
In this spirit, the role of differential equations is twofold, because they are not only a tool for evaluating individual integrals, but also the location where investigating how the exposed analytic properties of MI's could be transferred to amplitudes.
For any given scattering process the set of MI's is not unique [19]. Rather than considering it a limitation, this freedom can be exploited in order to select the basis of MI's that can be easier to evaluate. In the context of the differential equations method [14][15][16], reviewed in [20,21], and further developed in [22,23], simplifying the evaluation of MI's means achieving a block-triangular form of the system where: i) the dimensions of the blocks are minimal; ii) all analytic properties are exposed; iii) the solution can be determined through an algebraic procedure.
According to the study in Ref. [22], generic systems of DE's for MI's whose associated matrix depends on rational functions of the kinematics and could be brought in a canonical form where the -dependence is factorized from the kinematic. The integration of a canonical system of DE's can be carried out with algebraic patterns, where the analytic properties of its solution are manifestly inherited from the matrix associated to the sys- tem, which becomes the kernel of the representation of the solution in terms of repeated integrations [24]. These novel ideas for evaluating MI's have recently stimulated several applications [25][26][27][28][29][30][31][32][33][34][35]. In Ref. [29], we proposed to address the evaluation of MI's by means of a unitarylike formalism, where Magnus exponential [36] is responsible for the kinematic evolution of the MI's by acting on the boundary conditions. In particular, we observed that an initial set of MI's obeying a systems of DE's which depends linearly on can be brought in canonical form after applying a change of basis by means of a transformation that absorbs the constant term. The change from the initial set to the canonical basis of MI's can be implemented with Magnus exponential [36,37]. The solution admits an -expansion in terms of Dyson or equivalently Magnus series, whose coefficients are written as nested integrations [24,[38][39][40][41][42].
In this article, we apply Magnus exponential method to solve the system of differential equations fulfilled by 85 MI's required for the determination of the three-loop ladder box integrals with one massive leg, shown in Fig. 1b, and of the tower of integrals associated to their subtopologies, depicted in Figs. 3-5. All propagators are massless. The solution of the system is finally determined after fixing the values of the otherwise arbitrary constants that naturally arise from solving differential equations. In the considered case, boundary conditions are obtained by imposing the regularity of the MI's around unphysical singularities, ruling out the divergent behavior of the general solution of the systems.
Among the evaluated 85 MI's, there are 16 vertex-like integrals with two off-shell legs not yet considered in the literature, and 10 one-scale vertex-and bubble-like integrals, which had already been computed by means of Mellin-Barnes integral representation [43][44][45][46]. In general, homogeneous differential equations for single scale integrals carry only information on the scaling behaviour of the solution. The determination of the boundary constants for such differential equations amounts to the evaluation of the integrals themselves by other 2 Differential equations and Magnus Exponential All Feynman integrals belonging to the topologies which we consider in the present work can be defined in terms of the following sets of denominators D n (k i are the loop momenta). For the two-loop ladder diagram in Fig. 1a, one has supplemented by the auxiliary denominators For the three-loop ladder diagram in Fig. 1b one has instead supplemented by the auxiliary denominators D 11 = (k 1 + p 2 ) 2 , D 12 = (k 2 + p 2 ) 2 , D 13 = (k 2 + p 1 + p 2 + p 3 ) 2 , In the following we consider -loop Feynman integrals built out of p of the above denominators, each raised to some integer power, of the form where the integration measure is defined as The two-and three-loop MI's, respectively depicted in Fig. 2 and Figs. [3][4][5], are functions of the kinematic variables s = (p 2 + p 3 ) 2 , t = (p 1 + p 3 ) 2 , u = (p 1 + p 2 ) 2 , m 2 = (p 1 + p 2 + p 3 ) 2 , (2.8) with p 2 i = 0 (i = 1, 2, 3) and p 2 4 = (p 1 + p 2 + p 3 ) 2 , fulfilling s + t + u = m 2 . For convenience, we define the dimensionless ratios For planar topologies, like the ones we consider, one can always choose the Euclidean kinematic region m 2 , s, t, u < 0 such that the MI's are real. For definiteness, we work in the region 0 < y < 1, 0 < x < 1 − y (or equivalently 0 < x < 1, 0 < y < 1 − x). The analytic continuation of our results to regions of physical interests can be performed by generalizing to higher weights the procedure of Ref. [63]. The sets of MI's we choose to work with (see Sections 3 and 4) obey -linear systems of first order differential equations in the kinematic variables (σ = x, y), Given that m 2 is the only dimensionful variable, the differential equation in m 2 is related to the scaling equation. A factor (−m 2 ) − is already included in integration measure (2.7), therefore for each f i we simply have where n i is the dimension of the integral f i in units of a squared mass.

Magnus exponential
We can use Magnus exponential [29,36] to define a matrix B that implements a change of basis f → g, such that the new basis of MI's fulfills a canonical system of differential equations, ∂ σ g( , x, y) = Â σ (x, y) g( , x, y), (2.14) where the dependence on is factorized from that on the kinematic variables [22]. Accordingly, g is called canonical basis of MI's. Following Ref. [29], in order to build the matrix B of (2.13), let us introduce Magnus exponential matrix [36,37]. For a generic matrix M = M (t), this is defined as being the solution of the following matrix differential equation, The exponent Ω[M (t)] is given as a series of matrices obtained by repeated integrations, using M (t) as a kernel, where the first three terms of the series expansion read, ,

Canonical transformation
The basis change B, which brings the systems in the canonical form (2.14), can in general be obtained in both the two-and three-loop cases as the product of five exponential matrices, each implementing a change of basis according to (2.13), applied in successions. The dependence of each matrix on x and y is understood.
As a result of this transformation, A [0] m 2 vanishes and the MI's do not depend anymore on m 2 .

Then we decompose A
[0] x,0 in a diagonal term D x,0 and a off-diagonal one N x,0 , x,0 , (2.20) and use only the diagonal part for the next basis change, 3. We repeat as before and split A [1] y,0 into its diagonal and off-diagonal parts, A [1] y,0 = D [1] y,0 + N [1] y,0 , and we build the Magnus exponential again using the diagonal part 4. Now the matrix A [2] x,0 , has no diagonal term left, x,0 = N [2] x,0 , (2.24) therefore we build the Magnus exponential of the off-diagonal part [2] x,0 ] . (2.25) 5. The matrix A [3] y,0 has no diagonal term as well, A [3] y,0 = N [3] y,0 , (2.26) so we can define the last basis change After the last transformation we observe that A [4] x,0 = 0 = A [4] y,0 .
(2. 28) This means that the basis change of (2.13), with the matrix B given by

29)
absorbs the constant terms of A x and A y in the -linear systems in (2.10) and brings them to the canonical form (2.14):

Canonical system
The two partial-derivative systems satisfied by the new set of MI's, g = B −1 f , can be conveniently combined in an exact differential form, whereÂ(x, y) is logarithmic in the variables x and y, The matrices M i are n × n sparse matrices with purely rational entries, where n is the number of MI's. The value of n depends on the number of loops, and it amounts to n = 4, 18, 85, respectively at one-, two-and three-loop. The arguments of the logarithms are defined as letters, and the set of letters constitutes the alphabet. We observe that this alphabet is common both to the two-and the three loop cases. At one-loop [16], although not shown here, M 5 is absent, while at two-loop M 5 has only one non-vanishing entry. The solution of (2.14, 2.31) can be expressed as a Dyson series in , where the (matrix) coefficients of the series can be written as the iterated line integral, where dÂ i ≡ dÂ(x i , y i ). Equivalently, the solution admits a representation in terms of the Magnus exponential g(x, y, ) = e Ω[ dÂ](x,y) g 0 ( ), (2.36) where the vector g 0 ( ) ≡ g(x 0 , y 0 , ) corresponds to the boundary values of the MI's. This form is very suggestive, as Magnus exponential can be considered as an evolution operator, like in the unitary formalism, that brings the MI's g from their initial, boundary values to the considered point in the (x, y)-plane 1 . For definiteness, we integrate the exact differential form (2.31) from an arbitrary point (x 0 , y 0 ) to (x, y), along the broken path composed of the two segments in which one of the variables is kept constant. The integration is performed order by order in , up to a multiplicative vector of unknown constants. The latter are fixed by requiring the regularity of g(x, y, ) at the pseudothresholds (see Sec. 5).
The solution of the canonical system in (2.32) with the coefficient matrix given by (2.31) can be naturally expressed in terms of Goncharov's multiple polylogarithms (G-polylogarithms, for short) [24,[39][40][41], with w n being a vector of n arguments. The number n is referred to as the weight of G( w n ; x) and amounts to the number of iterated integrations needed to define it. Equivalently one has G-polylogarithms fulfill shuffle algebra relations of the type where shuffle product m n denotes all possible merges of m and n preserving their respective orderings. Because of shuffle relations, for a given alphabet and a given weight one can identify a minimal basis of G-polylogarithms. For the calculation of the considered 1 In this case, the evolution has to be understood like the variation w.r.t. the kinematic invariants that are the variables of the system of differential equations obeyed by MI's, rather than the quantum-mechanical time-evolution. Moreover, the matrices representing the systems of differential equations for MI's are not unitary.
two-and three-loop MI's the alphabet in (2.33) corresponds to the set of weights W = {0, 1, −x, 1 − x}. The number of basis elements depending on the weight n and the alphabet size α is given by the Witt formula: where µ denotes the Möbius function and the sum is done over all divisors of n. For reference we give values of N (n, α) relevant for this calculation in Tab. 1.

Two-Loop Master Integrals
A set of 18 two-loop MI's for the planar massless box integrals with one massive leg of Fig. 1a was first computed in Refs. [16,57]. We present the calculation of an alternative, yet compatible set of MI's. We begin by choosing the following basis where the integrals T i are depicted in Fig. 2 18 . The solid lines stand for massless particles; the dashed line represents a massive particle; dots indicate squared propagators; numerators may appear as indicated (p ij ≡ p i + p j ).

Three-Loop Master Integrals
The calculation of the three-loop ladder box integrals with one massive leg, represented in Fig. 1b, is the main, new result of this paper. After performing the automatic reduction

Boundary conditions
The generic solutions (2.34) of the canonical systems at two-and three-loop are written in terms of G-polylogarithms and constants to be fixed by boundary conditions. The alphabet (2.33) determines the thresholds which appear in the final result. Out of six thresholds only two are physical, since they correspond to the production of massless particles in s-and t-channels at x = 0 and y = 0. Imposing the regularity of the generic solutions at the unphysical thresholds, namely x = 1, y = 1, y = −x, y = 1 − x, amounts to ruling out the terms that give rise to divergent behaviours, hence enforcing conditions that unequivocally fix the arbitrary constants.

Relations for one-scale integrals
In general, homogeneous differential equations for single scale integrals carry only information on the scaling behaviour of the solution. Boundary constants for such differential equa-tions require the evaluation of the integrals themselves by independent methods. Within a multi-scale problem, such as the one we are considering, integrals may depend on more than one external invariant, and single-scale integrals participate in the regularity conditions of the multi-scale ones. Therefore, these relations can be exploited to determine the arbitrary constants of the single-scale integrals. Alternatively, they can reduce the number of independent single-scale integrals that needs to be independently provided. Therefore, solving multi-scale systems of differential equations yields the simultaneous determination of single-and multi-scale MI's, which are finally expressed in terms of a few single-scale MI's, to be considered as external input.
Let us discuss, as a pedagogical example, the systems of DE's for the two-loop master integrals g 2 , g 5 , g 8 and g 11 , represented by the corresponding topologies in Fig. 2, which read, and The zeroth order term of the -expansion is independent of x and y. Therefore, we can combine the equations above in order to find a relation between (the constant terms of) two one-scale integrals, At higher order in , these equations acquire a richer structure, because constants coming from the limiting values of the G-polylogarithms appear. For the considered example, the relation at the first order in is unaltered (the only constant which could appear being imaginary, hence not allowed in the Euclidean region), , (5.13) but at the second order in the relation becomes, . (5.14) Similar relations are systematically established, so that all MI's can be finally determined by providing few simple integrals as external inputs. At two-loop, the only external input is g 3 in (3.2), which can be independently computed and is given by, while g 3 and g 9 in (4.2) are the input integrals for the three-loop MI's, amounting to where we omit the common normalization factors (2.7).
We would like to observe that the relations between single-scale integrals, coming from the regularity conditions of multi-scale ones, seem not to belong to the set of IBP identities needed to derive the considered systems of differential equations. In the future, it is worth to investigate whether such relations are truly independent from IBP identities, or if they would arise when considering larger sets of identities for increasing powers of denominators and irreducible scalar products.
The sets of 18 two-loop MI's and of 85 three-loop MI's are respectively collected in two ancillary files, <2loopMIs.m> and <3loopMIs.m>. In Appendix C, we present only the analytic expression of the three loop g 83 -integral, see Fig. 5, that we consider the representative diagram of this work.

Conclusions
In this article we presented the analytic expressions of the 85 master integrals (MI's) of the three-loop ladder-box topology with one massive leg. Their calculation was performed with the method of differential equations, namely by solving a system of first order differential equations fulfilled by the MI's. The generic solution of the system was obtained in a purely algebraic way, by means of Magnus exponential method, and cast in terms of repeated integrations, according to Dyson series expansion, as recently proposed in Ref. [29]. The boundary conditions were provided by the regularity of the solutions at pseudothresholds.
The results of the considered four-leg integrals, as well as of the tower of three-and two-leg master integrals associated to subtopologies (including previously unknown twoscale vertex diagrams), were written as a Taylor expansion in the dimensional regulator parameter = (4 − d)/2. The coefficients of the series are expressed in terms of uniform weight combinations of multiple polylogarithms and transcendental constants up to weight six.
The considered integrals contribute to the N 3 LO virtual corrections to scattering processes like the three-jet production mediated by vector boson decay, V * → jjj, as well as the Higgs plus one-jet production in gluon fusion, pp → Hj, and to the three-loop one-particle splitting amplitudes.

A Canonical Matrices at Two-Loop
Here we present the sparse matrices M i (i = 1, . . . , 6) appearing in the canonical system defined in (2.31) and (2.32) obeyed by the MI's (3.2): (M 5 ) 13,13 = 4 is the only non-vanishing entry of M 5 .

B Canonical Matrices at Three-Loop
Here we present the sparse matrices M i (i = 1, . . . , 6) appearing in the canonical system defined in (2.31) and (2.32) obeyed by the MI's (4.2). We write each matrix in block form, where 0 43×42 is the null 43 × 42 matrix.
The matrix M 5 is mainly composed of zeroes and we only list its non-vanishing entries: (