A New Method for Calculating Differential Distributions Directly in Mellin Space

We present a new method for the calculation of differential distributions directly in Mellin space without recourse to the usual momentum-fraction (or z-) space. The method is completely general and can be applied to any process. It is based on solving the Integration-by-Parts identities when one of the powers of the propagators is an abstract number. The method retains the full dependence on the Mellin variable and can be implemented in any program for solving the IBP identities based on algebraic elimination, like Laporta. General features of the method are: 1) faster reduction, 2) smaller number of master integrals compared to the usual z-space approach and 3) the master integrals satisfy difference instead of differential equations. This approach generalizes previous results related to fully inclusive observables like the recently calculated three-loop space-like anomalous dimensions and coefficient functions in inclusive DIS to more general processes requiring separate treatment of the various physical cuts. Many possible applications of this method exist, the most notable being the direct evaluation of the three-loop time-like splitting functions in QCD.


I. INTRODUCTION
Achieving high precision in theoretical predictions is vital for the success of present and future collider experimental programs, as well as for the effective extraction of new physics from experimental data. A significant part of the theoretical work related to the experiment requires the evaluation of differential distributions, with most current research efforts focusing on the Next-to-Next-to-Leading Order (NNLO) or a higher level of precision. Examples of such distributions are the fully inclusive [1,2,3,4,5,6,7] and one-particle inclusive [8,9] DIS, the energy spectrum of hadrons in e + e − collisions [10,11,12], the total partonic cross-section [13,14,15,16] and rapidity distribution [17] for Higgs and vector boson [18] production at hadron colliders, Drell-Yan [19,20,21], transverse distribution of hadrons at hadron colliders [22,23,24,25,26] or particle spectra in the decays of muon [27] or heavy flavors [28,29,30,31,32,33]. Another important class of distributions that are universal and thus underlay the description of many physical processes includes the space-and time-like splitting functions [34,35,36,37,38], heavy flavor matching conditions [39] and the heavy quark perturbative fragmentation function [40,41,42].
The various distributions can be classified according to the number of kinematical variables they involve. Clearly, the larger the number of variables, the more complicated the evaluation of a distribution becomes. In this paper we will restrict our discussion to the case of distributions with a single kinematical variable. This class of distributions involves many important examples -some of them still significant open problems -like the three-loop time-like splitting functions in QCD. The extension of our discussion to cases with more than one variable will be rather transparent.
The choice of the most efficient approach to the evaluation of a particular single-scale distribution depends on its degree of "inclusiveness". The fully inclusive observables, like the fully inclusive coefficient functions in DIS [6], allow a simplified treatment based on the optical theorem. This is however a rare situation; most distributions of interest involve a specific final state, which requires that all contributing physical cuts of the relevant amplitudes be evaluated separately.
The purpose of this paper is to present a conceptually new calculational method of general applicability. As will become clear from the subsequent discussion, this method builds a bridge between two very important and seemingly unrelated calculational approaches as it provides a new perspective on the calculation of single-scale distributions. Moreover, during all stages of calculation this method requires no custom work and utilizes tools, techniques and programs that are publicly available and easy to implement in practice. Our method relies heavily on the Integration by Parts (IBP) identities [43]. It has the important feature of being formulated in terms of variables that are the most natural ones for the effective solving of the IBP identities.
With the above-described applications in mind, let us properly introduce the type of distributions σ(z) that we will be dealing with in this paper. Such distributions depend on a single kinematical variable z. For example, z can be the energy fraction of a parton produced in e + e − annihilation. We will assume that this variable is conveniently normalized: 0 ≤ z ≤ 1. The distribution σ is a scalar that is typically of the following form: The factor dPS (m) in Eq.(1) is related to the phase-space for the m-particle final state; it also contains the measure for the virtual integrations (if present). The precise form of this factor depends on the number of particles in the initial state. For a single-particle initial state processes with no virtual corrections it reads: In the case of processes with two particles in the initial state, dPS (m) has similar structure. It is detailed, for example, in [14]. Typically, expressions like Eq.(1) are UV and infrared divergent and in the following we assume that all divergences have been properly regulated by means of dimensional regularization. Besides z, the distribution σ(z) can depend on other parameters. Since their presence is irrelevant to our discussion, we will assume in the following that these have some fixed values and we will suppress them in our notations. The function f appearing in the argument of the δ-function in Eq.(1) is a dimensionless scalar. Its form is specific for each particular observable.
As a typical example we will consider the evaluation of the single particle inclusive cross-section for massless quark production in the decay of a colorless particle V → q + X (see also the appendix). Including the corrections up to next-to-leading order in the strong coupling and working in terms of bare quantities (i.e. no UV renormalization is performed) one has: In the example above, |M (k) | 2 denotes the terms proportional to α k S in the squared matrix element for the process V → q + X (see Fig.1). Clearly, the first line in Eq.(3) corresponds to the tree-level (Born) contribution while the second and the third lines respectively contain the contributions from the virtual and real-gluon emission corrections at order α S . On the above example, Eq.(1) stands for any one of the three lines in Eq. (3).
Perhaps the most elegant approach to date for the evaluation of distributions of the type Eq.(1) was proposed by Anastasiou and Melnikov [14] and further elaborated upon in [17,18]. Let us recall the salient features of this method. One uses the distributional identity: to formally replace all δ-functions appearing in Eq.(1) with propagators coinciding with the arguments of the δfunctions, i.e. one introduces the invertible mapping P acting only on δ-functions: with c an arbitrary function of the propagators. The utility of the mapping (4) is that it allows one to treat the object on the right-hand side of Eq.(4) with the usual IBP identities [43]. By solving these identities one reduces the initial distribution σ(z) to a combination of a small number of irreducible objects. Eventually, one performs the inverse mapping P −1 , thus expressing σ(z) as a linear combination (with simple known coefficients) of a small number of well defined master integrals. From the IBP identities it also follows that the master integrals satisfy a system of differential equations, which can be solved to obtain their z-dependence. To fully specify the solutions of the differential equations, one has to prescribe the corresponding boundary conditions; these can be extracted from an explicit evaluation of the master integrals in a particular kinematical point like z = 1.
In the phenomenological applications one also needs the Mellin transform 1 of the distribution in question: Performing the Mellin transform results in the evaluation of integrals over, typically, combinations of polylogarithms and rational functions of z. At present, and certainly for the case of massless distributions, there exists a very good understanding of the mapping between the classes of basic functions in z and n spaces [44,45,46,47,48,49]. In the following discussion we will consider the knowledge of σ(z) as equivalent to that of σ(n) and vise versa, i.e. we will tacitly assume that one can always perform the needed Mellin or inverse Mellin transforms. That is definitely true for the massless case. In more complicated situations one may have to resort to numerical methods to perform the inverse Mellin transform [50,51]. In any case, we need not bother about that point here. We will consider our problem as solved, as long as we know either σ(z) or σ(n).

II. THE METHOD
In the present paper we would like to advocate a new approach to the evaluation of the distribution in Eq. (1). It aims at the direct evaluation of σ(n) without calculating it first in z space as is done at present. Our proposal is to explore the obvious possibility that one can integrate over z before performing the phase-space and/or virtual integrations: We see that as a result of the interchange of the order of integrations the invariant f enters the integrals as a propagator raised to power −n. That power, however, should be treated as an abstract parameter that takes arbitrary and not fixed integer values.
By applying the mapping Eq.(4), and interchanging the order of integration as in Eq. (6), one can bring the original problem of calculating σ(z) to the following form: where D represents the appropriate measure originating from the real and/or virtual integrations, P i denote the propagators originating from the evaluation of the amplitude and x i are the arguments of all phase-space δ-functions (if present). The argument of the δ-function that defines the observed fraction z is denoted by f and the powers a i and s are some fixed integers. This way, we have effectively reduced the problem of the calculation of the differential cross-section σ(z) to the problem of evaluation of functions of the following general form: The function S appearing in Eq. (8) can also depend on other fixed parameters.
In principle, scalar terms like the one in Eq. (8) can be simplified to a minimal set of terms by applying the Integration by Parts (IBP) identities [43]. Until now, there has been no effective method for solving the IBP identities when one (or more) of the powers a 1 . . . a p is an abstract parameter (however see [52]). In the following, we present one very efficient approach for solving this problem.

III. METHOD FOR SOLVING THE IBP IDENTITIES IN PRESENCE OF AN ABSTRACT POWER
It is very well known (see for example the book of Smirnov [53] for detailed introduction) that when applied to the object S in Eq. (8), the IBP identities result in a system of linear homogeneous equations with rational coefficients that relate functions S with arguments shifted by ±1 relative to each other. If all a i 's were fixed integers, then by successively relating terms that differ with ±1 one can eventually express the original function S(a 1 , . . . , a p ) through a linear combination of several, say m, master integrals S 1 ({i 1 }), . . . , S m ({i m }). The masters S j are special cases of S(a 1 , . . . , a p ) with their arguments ({i j }) taking special values. At present, the most popular method for solving the IBP identities is the one of Laporta [54]. It is based on solving the systems of linear homogeneous equations directly, through Gauss elimination.
Clearly, if one of the parameters a i is not an integer this procedure cannot work, since: 1) with only integer steps one cannot relate the initial non-integer parameter to an element S(a 1 , . . . , a p ) with only integer a i 's and 2) the number of steps in the Gauss elimination cannot even be specified when one of the parameters is an abstract number.
In the following, we present a solution to this problem. In order to facilitate our discussion we shall assume that a 1 , . . . , a p−1 are integers having some specific values, while the last argument, a p , is an abstract parameter.
From Eq.(1) it is clear that σ(n) is a sum of a number of terms of the type in Eq.(8) that have different values of their indexes (a 1 , . . . , a p ). It is very important to observe, however, that the difference between any two values that the index a p can take is always an integer, i.e. a p − a ′ p ∈ N. This is a crucial observation, which one can use to modify the strategy for solving the IBP identities in the following way. First, one relaxes the requirement that the masters must have integer-valued indexes. Second, as we will explain in a moment, one can choose all masters in such a way that they all have the same value, say a p = r / ∈ N, of their last index i.e. the masters are all of the form: with the same r, and the j i 's being fixed integers specific to each master. It is indeed possible to arrange that all masters have the same value of the non-integer-valued index a p . That follows from the arbitrariness of the value of this parameter (we only assume that it is non-negative). Since there is no preferred value for that index, the IBP system has a sort of translational invariance along the index a p . One can understand this by saying that r and r + k, where k is a fixed integer, are equally arbitrary. Therefore, we can take as a reference value for the index a p the number r which we will consider abstract but having fixed value. Having done that, the "translational" invariance along the values of a p is now "broken". Clearly, the value r now plays the role of a zero reference point much like the value a p = 0 in the usual case when all indexes take integer values. Therefore all one need to do is to measure in integer units how much the value of the last index of an element S is displaced from the reference point r.
Next, we give a practical recipe of how to implement the above idea. Let us work with the functions 2 B: where as the "reference" point for the last index we take the Mellin variable n.
As follows from Eqns. (6) and (7) the distribution σ(n) takes the following form: where c a1,...,ap−1 are some known coefficients. To construct the needed algebraic reductions, one first applies the IBP identities on a generic monomial of the form: where all powers ν i are treated as arbitrary parameters. Next, one identifies each term of the form (11) appearing in the IBP equations, with the function B(ν 1 , . . . , ν p−1 , ν p + n), followed by the substitution ν p → ν p − n. After this manipulation the Mellin variable n is explicitly present as a parameter in the resulting equations. They can be solved in any approach available, including the one of Laporta.
A word of caution: one has to keep in mind that, as follows from Eq.(9), the functions B implicitly depend on n. Therefore, one should not confuse the integer value k in the last argument of the function B with the absolute power of the corresponding propagator 1/f , but should think of it as the "distance" -in integer units -from the reference power n.
Finally, one can map all integrals appearing in Eq.(10) to the masters obtained from the solving of the just-described reduction. This mapping is done in the standard way.
Next, we explain how one can extract the n dependence of the master integrals. Assume that an element B(b 1 , . . . , b p−1 , 0) is a master integral (with b 1 , . . . , b p−1 some fixed integers). One can inspect the already solved IBP reduction and read off from there the result for the element B(b 1 , . . . , b p−1 , −1). Note that this element differs from the master B(b 1 , . . . , b p−1 , 0) only by the value of the last index. If the element B(b 1 , . . . , b p−1 , −1) is not a masters itself, then it must be a linear combination of the master integrals: Here c(n) is a known, typically not very complicated function, and the term G(n) is a homogeneous linear combination of all master integrals, except for the master B(b 1 , . . . , b p−1 , 0). Eq. (12) is a first order non-homogeneous difference equation of the type F (n + 1) = c(n)F (n) + G(n) (recall Eq. (9)) for the master B(b 1 , . . . , b p−1 , 0). Clearly, repeating this procedure for each one of the master integrals found in the reduction run, one can derive a complete system of difference equations for all the masters. Typically, one observes certain hierarchy among the master integrals; the simplest ones satisfy homogeneous equations (i.e. G(n) = 0) that can be solved in terms of Γ-functions. These integrals then comprise the non-homogeneous terms for the equations of other masters, and so on.
In case the element B(b 1 , . . . , b p−1 , −1) is also a master integral, one should read off from the reduction the result for the yet higher term B(b 1 , . . . , b p−1 , −2). One should continue doing this until one reaches an element B(b 1 , . . . , b p−1 , −k) which is not a master itself but all elements B(b 1 , . . . , b p−1 , −s) with 0 ≤ s < k are masters. The result from the reduction for the element B(b 1 , . . . , b p−1 , −k) represents a k-th order difference equation for the master integral B(b 1 , . . . , b p−1 , 0).
To solve the resulting difference equations one can make use of existing techniques. Such equations were analyzed and successfully solved in the course of the evaluation of the three-loop anomalous dimensions in QCD [37,38] and of the two [5] and three-loop [6] coefficient functions in DIS. In most cases of physical interest the resulting difference equations can be solved after expansion in ǫ in terms of harmonic sums or their generalizations. In simpler cases, one can even solve these equations in closed form in terms of hypergeometric and/or Γ-functions. Difference equations as a way of calculating master integrals were also used by Laporta [54]. Upon solving the system of difference equations for the master integrals, one has achieved a complete extraction of the dependence of the masters on the Mellin variable n. The only remaining thing to do is to specify the initial conditions for the solutions of the difference equations. Typically, that would be the value of the masters for n = 0. This is an important fact. It implies that to completely specify the master, one need to only evaluate integrals that are fully integrated over the available phase-space. These are pure numbers that do not depend on the kinematical variable n. On the conceptual level, this is placing the evaluation of certain not-completely inclusive observables one step closer to the very familiar fully-inclusive case where, thanks to the optical theorem, one can significantly simplify the calculations by not considering separately all possible physical cuts.
Often, one can reduce the number of fixed-n integrals that have to be evaluated by hand. This follows from the property of the fixed-n reduction that not all integrals corresponding to initial conditions for the n-dependent masters are actually independent. To explore this fact one has to perform a separate fixed-n IBP reduction where n = 0 is taken from the very beginning. We have observed in simple one-and two-loop reductions as well as in rather complicated three-loop cases that this procedure indeed generates additional relations between the integrals corresponding to the initial conditions of the master integrals. One might wonder about the cost of such an additional run. That, however, should be of no concern since the fixed-n reduction is much simpler and faster than the general-n run one has to be able to perform anyway. The computer load pays off with the elimination of many of the integrals that otherwise have to be computed by hand. A fixed-n reduction can also be used as a cross-check of the general-n calculation. This is similar to the use of Mincer [55] in the three-loop DIS calculations in [36,37,38].

IV. PARTIAL FRACTIONING
Consider a case where during the evaluation of the amplitudes one gets a propagator that is not linearly independent from the constraint f defining z (see Eq. (7)). Clearly, that can happen in many ways and such linear dependence might even involve a group of several propagators. However, to simplify our point as much as possible, we will only consider a simple situation. Consider the function: where . . . stay for powers of other possible propagators. If we were to evaluate this integral in z-space we would first replace f everywhere with z, as is implied by the factor δ(z − f ). Then the factor 1/(1 − f ) becomes just the number 1/(1 − z) and drops out of the integral. However, when we work in Mellin space the constraint f = z cannot be used anymore. If n were some fixed integer, we could have applied partial fractioning n-times and split the linearly dependent propagators 1/f and 1/(1 − f ). For symbolic n, however, that cannot be done and one should again resort to solving difference equations. This can be done in the following way. The identity: immediately translates into a difference equation for the function appearing in Eq.(13): F (n + 1) = F (n) + G(n). This is a simple difference equation with non-homogeneous part given by: The integral (14) does not contain linearly dependent propagators and can be evaluated by using the procedures described previously. Another way of eliminating the linear dependence among the propagators is to expand the propagator 1/(1 − f ) in geometric series. That would completely eliminate this propagator and one would end up with a standard problem where the index n is replaced by n + k, k ≥ 0 (here one can apply the usual reduction since the index n + k is as arbitrary as the index n is). Finally, one would have to sum up the resulting expression over the index k.
Our experience shows that the combination of the above methods is sufficient to eliminate the appearance of linearly dependent propagators in any situation.

V. CONCLUDING REMARKS
The method presented in this paper represents a conceptually new approach for the evaluation of differential distributions. At the same time it has the advantage that it shares features with other available methods that previously were applied with impressive success. We would like to critically compare our approach with these methods and clarify its distinct applicability.
Our method has the following advantages over the direct momentum space approach [14]. First, the solving of the IBP reductions is much more efficient and fast. In practical terms that may not be an issue in simple one-or two-loop cases, but we have checked that it brings enormous improvement when applied, for example, at three loops. Second, the number of master integrals being produced in the course of solving the IBP's is smaller and, third, our method only involves computation of integrals that are pure numbers. This is a significantly easier task compared to the calculation of the z-dependent initial conditions at fixed values of z.
We also expect our method to bring new insights into the calculations of distributions with more than one scale. Our approach has a few common features with the method used in the inclusive DIS calculations at two-and three-loops [5,6,36,37,38,56]. Still our method has wider, in fact completely general, applicability and it relies for the formulation and solving of the IBP identities only on well established, multipurpose and publicly available methods and software. An important shared feature is the fact that the master integrals in both methods satisfy difference equations in terms of the Mellin variable n. In this respect our method benefits greatly from these existing developments, since many of the technical tools needed for its practical realization have been already established. That includes the methods for solving the difference equations, the current good understanding of the appropriate functional bases in n and z spaces and their tabulation. Finally, well established multipurpose software like FORM [58], Summer [59] and XSummer [60] exist and are publicly available. They are capable of effectively dealing with the necessary algebraical manipulations. In fact, our method opens up new venues for the application of the techniques from the two-and three-loop inclusive DIS.
As a first application of our method we have rederived the NLO coefficient functions in e + e − [10] (in fact to all orders in ǫ).
We also have experience with reductions at three-loops, where the advantage of our method becomes apparent. Our method ensures better performance of the IBP solving software and produces smaller number of master integrals.
The initial condition B(0, 1, 1, 1, 0)(n = 0) is defined through Eq.(A3) after setting n = 0 there. It is a trivial to compute number. Our work is not quite done with the solving of the IBP reductions and the evaluation of the master integral since there also appears the propagator P add = (q + k) 2 which does not belong to the set P 1 , . . . , P 5 . This propagator results from the right diagram on Fig.1; it is not linearly independent from the set P 1 , . . . , P 5 but it can appears downstairs together with these propagators. To resolve the situation one has to resort to the partial fractioning technique discussed in section IV.
Exploiting the constraints implied by the three δ-functions (with arguments P 2,3,4 ), one can easily establish that P add = 1 − P 5 . In this case one can apply the geometric series trick discussed in the previous section to all terms where P add is present downstairs: Next one takes the summation outside the integrals; the resulting integrand is of the type in Eq.(A2). With the help of the IBP reduction this integral can be reduced to the master integral discussed above. Inserting the explicit form of the master, one obtains very simple expression containing only Γ-functions. The summation over s of this product of Γ-functions can be easily performed and it results again in a product of the same type of functions. Thus, the result from the real emission radiation at order α S can be easily evaluated in closed form to all orders in ǫ. To use this expression in practical applications one has to decompose it in series in ǫ. This expansion can be easily automated with the help of the programs Summer [59] and XSummer [60].