Master integrals for splitting functions from differential equations in QCD

A method for calculating phase-space master integrals for the decay process $1 \to n$ massless partons in QCD using integration-by-parts and differential equations techniques is discussed. The method is based on the appropriate choice of the basis for master integrals which leads to significant simplification of differential equations. We describe an algorithm how to construct the desirable basis, so that the resulting system of differential equations can be recursively solved in terms of (G)HPLs as a series in the dimensional regulator $\epsilon$ to any order. We demonstrate its power by calculating master integrals for the NLO time-like splitting functions and discuss future applications of the proposed method at the NNLO precision.


Introduction
After the recent success of the LHC Run I experiment, and discovery of the Higgs boson in particular, the LHC Run II is pushing the limits of higher-order calculations in QCD even further than ever. In that context, analytical calculations play a crucial role as a background for the numerical methods and phenomenological analyses in QCD.
In this paper, we focus on the analytical calculation of phase-space master integrals for 1 → n decay processes with massless particles in the final state. This type of decays within the electron-positron annihilation reactions have given us much information about the properties of quarks and gluons and the nature of their interactions as described by QCD. Moreover, they will play an outstanding role for the further precision studies of QCD at upcoming e + e − colliders at even higher energies. The classical example is jet production in e + e − annihilation, which can be used to extract values of the strong coupling constant α s from the three-jet rate and related event shape observables. In the past decade, nextto-next-to-leading order (NNLO), i.e., O(α 3 s ), contributions to the three-jet rate from the process γ * → 3 partons in e + e − annihilation were calculated [1][2][3][4]. Further improvements to this calculations at N 3 LO inevitably require analytical expression for the integrals we consider in this work. For example, three-loop splitting functions are a must-have piece for numerical calculations of N 3 LO contributions to the three-jet rate from the γ * → 6 partons process. Splitting functions for the initial-state radiation, i.e., space-like, are known exactly at NNLO [5,6]. In contrast, for the final-state radiation, i.e., time-like, they are known at NNLO only approximately [7][8][9]. Despite the fact that those uncertainties are numerically irrelevant for phenomenological applications, e.g. for the evolution of fragmentation functions [10]. The exact result are still needed, as mentioned before, for performing numerical integration in various subtraction schemes when the need to integrate local counter terms arise [11][12][13].
At the same time, a huge progress has been made in development of tools and methods for higher-order calculations in the field theory and perturbative QCD in particular. Integration-by-parts (IBP) reduction of Feynman integrals [14,15] together with differential equations for master integrals [16] proved, by state-of-the-art calculations, e.g. [17][18][19][20][21], to form a powerful framework for calculating high-order Feynman diagrams. Despite that it is usually applied for virtual integrals at the level of amplitudes, this approach can be used with the same success for analytical calculation of real phase-space integrals at the level of matrix elements, where standard approaches are usually applied, i.e., to parametrize a phase space explicitly and proceed accordingly with Feynman parameters integration and similar methods [22,23], or alternatively to work in the Mellin space with recursion relations [24] or a system of difference equations [25].
During the past few years, the method of differential equations became very popular due to the fact that a good choice of the basis for master integrals leads to significant simplifications of the differential equations [26]. Although, in general, finding an appropriate basis is not easy, the approach based on the Moser algorithm [27] was discussed in [28], which allows to reduce the system in one singular point, but not globally. A global extension of the Moser algorithm, which shows how to adjust the transformations in such a way that they do not spoil behavior in any other point was presented in [29]. It allows for systematic simplification of differential equations. Unfortunately, there is still no computer implementation of these methods available for the public use, which is very desirable to automatize the process.
In this work we propose an alternative method for calculating phase-space master integrals from differential equations and show how to fix boundary conditions. The algorithm is self-consistent, in sense that if all the prerequisites are fulfilled the proof is reduced simply to verifying that proposed solutions satisfy initial equations. The main advantage of the proposed approach is that it is relatively simple, can be easily implemented as computer code, and at the same time gives a complete solution for masters to any power in ǫ.
Although it may be not as general as other methods, it can be successfully applied for calculating splitting functions, but is not limited only to that case. As an example of practical use, we perform a detailed calculation of the master integrals for the NLO contribution to time-like splitting functions and discuss possible extensions to NNLO accuracy.
The paper is organized as follows: in Section 2 we introduce the notation and show how to calculate splitting functions from the e + e − annihilation process. In Section 3 we formulate a solution for the system of differential equations for phase-space master integrals of the topology 1 → n derived from IBP reduction rules in x-space. In Section 4 we calculate master integrals for the NLO splitting functions of 1 → 3 and 1 → 4 topologies. Finally, we discuss properties of the solutions obtained with our approach and its possible extensions to higher orders.

Splitting functions in QCD
Let us briefly review the main facts on splitting functions in the collinear factorization formalism of QCD, mainly for notation consistency. For a more detailed review we refer the reader to [25,30].
Splitting functions govern the collinear evolution in hard scattering processes with hadrons in the initial (space-like) or final (time-like) state. For processes with identified hadrons in the final state the parton-to-hadron transition is described by the parton fragmentation distributions D h f (x, q 2 ), where x represents the fractional momentum of the final-state parton f transferred to the outgoing hadron h and q 2 ≥ 0 is a time-like hard scale. The scale dependence of the fragmentation distributions is controlled by the so-called time-like splitting functions P T ba (x) 1 , and is given by where the summation runs over the number n f of effectively massless quark flavors and the gluon, b = q i ,q i , g for i = 1, . . . , n f . The splitting functions P ba can be computed in perturbation theory in powers of the strong coupling α s , where we normalize the expansion parameter as a s = α s /(4π).
As discussed at length in [30], splitting functions can be extracted using the mass factorization formalism from the electron-positron annihilation processes with photon (γ) exchange and Higgs (φ) boson exchange in the effective theory and a tagged parton p = q,q, g with momentum k 0 . For the photon-exchange process (2.3), following the notation in [31], the unpolarized differential cross-section in m = 4 − 2ǫ dimensions is given by

5)
1 Further in the text we omit the superscript T and assume all splitting functions to be time-like.
where θ denotes an angle between the beam and parton momentum k 0 . The scaling variable x is defined as For the demonstration of the method for calculating master integrals described in detail in Section 3, let us consider the time-like q → g splitting function at NLO. It can be written as where P (nr×nv) qg denotes contribution from the diagram with n r real and n v virtual legs, as illustrated in figure 1. In particular, we are interested in 1/ǫ contribution to the transverse fragmentation function, as discussed in [30], where the hadronic tensor W (2) µν (x, ǫ) for the real-virtual and real-real cases becomes and where M (3) and M (4) are amplitudes for the processes depicted in figures 1a and 1b respectively, l is a loop momentum, and dPS(n) denotes a n-particle phase-space integral In the Section 3 we provide a method to calculate this kind of integrals with some detailed examples. Below, however, we would like to present a general plan of the calculation: 1. Generate amplitudes in figure 1 with QGRAF [32] and construct from them fragmentation function F T (x, ǫ) using FORM [33].
3. Find master integrals solving differential equations in x-space as described in the next section.

Master integrals from differential equations
We consider a homogeneous system of differential equations, which in the most general case takes the form where the coefficients a ij (x, ǫ) (or the n×n matrix A(x, ǫ)) are known, f i (x, ǫ) are unknown functions, and ǫ is an infinitesimally small parameter (playing the role of a dimensional regulator in m = 4 − 2ǫ dimensions).
Assuming that the coefficients a ij (x, ǫ) are rational functions of ǫ, without loss of generality they can be written in the form where r ǫ is an integer (likely negative), which we use to denote an ǫ-rank of the matrix A(x, ǫ).
On the other hand, we restrict the matrix A(x, ǫ) to have the form where i runs over some finite set, p i is said to be the Poincare rank of A i (x, ǫ) at a singular point x i , and A i (x, ǫ) is a regular matrix at x = x i , i.e., polynomial. Such a form is imposed exclusively for a practical reason since calculations of the splitting functions are bound to the case of x i ∈ {−1, 0, 1}, which is exactly an alphabet for Harmonic Polylogarithms (HPLs) [36]. In the case of a more complex structure of denominators in the expansion eq. (3.3) the same arguments could be extended to the more general case of Generalized Harmonic Polylogarithms (GHPLs) introduced in [37], which maintain the structure and properties of HPLs [38,39]. Keeping all the above considerations in mind we proceed with providing a solution for eq. (3.1) as an ǫ-series. Taking into account a recursive definition of (G)HPLs we show that such a series can be found to any order in ǫ at a low computational price.

Solutions for ǫ-rank > 0
We are looking for the solution of the system eq. (3.1) in the form Keeping in mind the expansion eq. (3.2) it is easy to show that expansion coefficients calculated by the recursive formula lead to the desired solution, where c (k) i are integration constants determined from boundary conditions as described in Section 3.4.

Solutions for ǫ-rank = 0
There is no general solution for the system with ǫ-rank = 0, however for some special cases it is possible to write down such a solution. In this paper we consider weakly coupled systems, these are systems for which a In such a case it is possible to choose a new basis so that a new system has ǫ-rank > 0 and can be solved using the method of Section 3.1. In the remaining part of this section we provide a procedure how to accomplish that task, which consists of finding such new bases that: ij (x) = 0 for i = j; and ij (x) = 0 for i > j.

Zero-diagonal form
It is easy to verify that a system of differential equations for a new basis defined as contains zeroes as diagonal elements and has a new form

Zero-triangular form
Next, following the same strategy, we find a new basis h i (x, ǫ) which leads to the zerotriangular form of the differential equations: A complete form of the new system is rather complex and it is of no practical use to write it down here. However it can be easily obtained from eq. (3.10) after coefficients of eq. (3.11) are explicitly calculated. Let us show that such a choice indeed provides the desired zero-triangular system of equations. Taking the derivative of eq. (3.10), keeping in mind eq. (3.9) and eq. (3.11), and neglecting higher-order term in ǫ we obtain It is easy to check, by carefully switching summation variables in one of the nested sums, that right-hand side of eq. (3.12) becomes zero. At first sight, it may look that nested integrals in eq. (3.11) are way too complicated for practical calculations. In fact, they are very easy to compute taking into account the recursive nature of (G)HPLs, as was discussed earlier in this section. For our examples, discussed in the next section, we have used the HPL package [40].

Solutions for ǫ-rank < 0
As a rule, when one chooses a basis of master integrals as provided directly by the IBP rules generator, like FIRE [41], Reduze [42], or LiteRed [34,35], the system eq. (3.1) has a negative ǫ-rank. In this situation we can not proceed with the procedure described before in this section. To overcome this issue it is usually enough to adjust ǫ n factors in the masters, for example see eq. (4.3) in Appendix A.
To get a hint on how to choose n we analyze Mellin moments for the corresponding masters that leads to several possibilities:

In the presence of factors
ǫ is an ǫ-rank of the i th Mellin moment. The reason is that the logarithmic singularity in x is canceled by the Mellin moment while the second one in 1 − x introduces an additional ǫ pole. For the illustration see masters V 6 , R 7 , R 8 in the Appendices.

In the presence of factors
3. Otherwise we choose n = r (0) ǫ .

Boundary conditions
The final step of the method is to find integration constants c (k) i emerging in eq. (3.5). On the one hand, in the case of phase-space integrals we can do that by calculating Mellin moments of the solution eq. (3.5). On the other hand, the same moments can be taken from the literature or directly calculated by performing integration over the entire n-particle phase-space, i.e., As in the case of the phase-space integrals with x-space projection eq. (2.11), we can generate IBP rules for the inclusive integrals as well. That allows us to reduce the set of inclusive masters which should be calculated explicitly. Another simplification is related to the Mellin moments, which can be extracted from the difference equations. These equations in turn can be derived from the differential equations eq. (3.1), hence only one Mellin moment needs to be computed for each inclusive master.

Master integrals for NLO splitting functions
Finally, we demonstrate the practical application of the method described in the previous section. We choose to calculate two-loop contributions to the time-like splitting function P (2) qg (x) since its three-loop contribution is still not known exactly, however it will be possible to obtain them by future extension of this example to NNLO.
We follow the plan described at the end of Section 2. After IBP reduction, done with the help of LiteRed, we obtain 6 real-virtual (figure 2) and 8 real-real (figure 3) masters for contributions depicted in figure 1.
Step 1. In order to obtain a system of differential equations with non-negative ǫ-rank we choose ǫ n factors as described in Section 3.3 (see eq. (4.3)). The resulting system is given by eq. (4.9).
Step 2. We change the basis to obtain a zero-diagonal system, as described in Section 3.2.1: Step 3. We make the last change of the basis in order to obtain a zero-triangular system as described in Section 3.2.2: Step 4. We solve the resulting equations with the help of eq. (3.5) as described in Section 3.1 and return to the initial basis.
Step 5. We find the final result by fixing boundary conditions using Mellin moments given in Appendix A.

Real-real contribution
By analogy with the real-virtual case we proceed with the real-real contribution with final results given in Appendix A. We define master integrals depicted in figure 3 as where denominators D j are defined in eq. (4.6).
We tune ǫ n factors in order to obtain a non-negative ǫ-rank matrix which is given by eq. (4.10). Corresponding powers of ǫ can be seen in eq. (4.7). Next, a zero-diagonal basis reads Finally, a zero-triangular basis reads In summary, with the help of the method described in Section 3 we have found the master integrals of figures 2 and 3. The solutions are presented in Appendix A as a partlyexpanded series in the dimensional regulator ǫ with at least 3 leading terms, i.e., HPLs of weight 2. That is enough for our purpose, i.e., to extract slitting functions as discussed in Section 2. Furthermore, the higher-order ǫ-terms of the presented solutions are easy to obtain provided the corresponding ǫ-terms of Mellin moments, required to fix boundary conditions, are know.

Conclusions
In this paper we proposed a method for calculating phase-space integrals for the decay process 1 → n massless partons in QCD using integration-by-parts and differential equations techniques. The key idea of our approach is based on choosing a basis of master integrals which leads to significant simplification of differential equations. As a main result of this work, we describe an algorithm how to construct such a basis and find a solution of the resulting differential equations. The advantage of our approach comparing to available techniques is that it is relatively simple to automate for execution on a computer without loss of generality of the final solution, which is obtained to any order in the dimensional regulator ǫ in terms of (generalized) harmonic polylogarithms. That requires however to know at least one Mellin moment for every master integral in order to determine boundary conditions for the final solution.
In order to demonstrate how our method works in practice, we calculate master integrals for the decay processes 1 → 4 and 1 → 3 with a projection to x-space, needed to extract NLO time-like splitting functions from e + e − annihilation process. Analyzing this example we notice that another asset of the proposed method is that resulting master integrals are explicitly regulated in the singular points with the help of the dimensional regulator ǫ, manifested by overall factors x −1+aǫ and (1 − x) −1+bǫ in the final result.
The generalization of the results to NNLO topologies with loop insertions, needed to obtain missing n 2 f pieces of the off-diagonal time-like splitting functions, is particularly straight-forward due to the factorizability of the phase-space. In addition, master integrals with various types of projectors, not only to x-space as in the case of splitting functions, can be obtained with the described method as well.

Real-virtual case
Mellin moments of the real-virtual masters used to fix boundary conditions read [25]  Real-virtual master integrals read