On the Threshold Resummation in Forward pA Collisions

In this paper, using the Higgs production in forward rapidity region in proton-nucleus collisions as an example, we demonstrate that we can construct a systematic formalism for the threshold resummation for forward rapidity particle productions in the saturation formalism. The forward threshold jet function, which satisfies the corresponding renormalization group equation, is introduced into the new factorization formula. This calculation can be easily generalized to other processes, such as single forward hadron productions at forward rapidity region, and have important phenomenological implications.

In this paper, using the Higgs production in forward rapidity region in proton-nucleus collisions as an example, we demonstrate that we can construct a systematic formalism for the threshold resummation for forward rapidity particle productions in the saturation formalism. The forward threshold jet function, which satisfies the corresponding renormalization group equation, is introduced into the new factorization formula. This calculation can be easily generalized to other processes, such as single forward hadron productions at forward rapidity region, and have important phenomenological implications.

I. INTRODUCTION
Single inclusive particle productions in the forward rapidity region in proton-nucleus collisions, p + A → H(y, p ⊥ ) + X, is of particular importance in the search for the onset of the gluon saturation phenomenon, which occurs in high gluon density heavy nucleus target at very small x region. In this process, we measure the transverse momentum distribution of particles produced in the forward rapidity region y > 0 in the proton beam direction. It is straightforward to find that the active partons with longitudinal momentum fraction x p ∼ m ⊥ √ s e y in the proton projectile lie in the large x region, while the active partons (mostly gluons) with longitudinal momentum fraction x A ∼ m ⊥ √ s e −y in the nucleus target are from deeply small-x region when the rapidity y is sufficiently large. Here m ⊥ = m 2 + p 2 ⊥ is defined as the transverse mass of produced particle while s is the center of mass energy. Physically, partons inside the proton projectile, which act as dilute probes, can pick up sizeable amount of transverse momentum, which is of the order of the saturation momentum, after traversing the dense gluonic medium in the heavy nucleus target. Therefore, it is of great interest to study the transverse momentum distribution of particles especially in the low transverse momentum region in order to study the gluon saturation phenomenon. Normally, the so-called dilute-dense factorization, which uses collinear parton distribution functions (PDFs) for partons from the proton side and small-x gluon distributions for the low-x gluon originated from the target nucleus side, is widely adopted to formulate particle productions in the forward region. In this approach, by measuring produced particles in the forward rapidity region, one can take advantage of the extremely asymmetric kinematics (namely, x p → 1 and x A 1) to maximize the gluon saturation effects. On the other hand, x p → 1 means that this process is close to the kinematical boundary of the phase space, where the soft gluon radiations become important. The resummation of soft gluon radiations near the kinematical threshold is known as the threshold resummation. The objective of this paper is to understand the threshold resummation in the forward region in the dilute-dense factorization framework.
There have been great efforts of theoretical and phenomenological studies [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] on forward hadron productions in pA collisions using the dilute-dense factorization. In particular, the next-to-leading order (NLO) single hadron production cross section in pA collisions in forward rapidity becomes negative in the large transverse momentum region [11]. There has been some speculation [21] that the threshold resummation could in principle help to mitigate the negativity problem by systematically include higher loop contributions. Furthermore, authors in Ref. [17] specifically demonstrated that one important source of the negative cross section comes from a term which is proportional to α s ln(1 − x p ) with x p → 1. Analytically, one can trace the origin of the α s ln(1 − x p ) term and find that it is due to the so-called plus distribution in the NLO correction. Sometime, threshold resummation is viewed as a resummation of the 'plus' distributions, for example, 1 τ dξ (1−ξ)+ ∼ ln(1−τ ) in the limit τ → 1, although constant terms associated with soft gluon emissions are also resummed. This implies that the so-called threshold logarithm can start to appear and become important in the dilute-dense factorization approach. In addition, the ultra-forward rapidity hadron spectrum has been measured recently by LHCf [22], which is sensitive to the very large x region of the projectile hadron and extremely small x region of the target nucleus [23]. In this case, we believe that the threshold resummation becomes indispensable. However, the theoretical framework which can systematically resum such type of threshold logarithms in the dilute-dense factorization approach is not yet available. In the following, we will investigate the the threshold resummation in forward rapidity productions of particles in pA collisions and build such a theoretical framework.

arXiv:1806.03522v1 [hep-ph] 9 Jun 2018
In general, threshold resummation, which was originally formulated by Sterman [24], and by Catani and Trentadue [25] for the Drell-Yan process, has been a very important topic in high energy QCD studies for the last thirty years. The resummation technique developed in these two papers has been generalized and applied to many other QCD processes, and has been proved to be very useful in QCD phenomenology. For example, de Florian and Vogelsang [26] applied the threshold resummation to high p ⊥ π 0 productions in pp collisions with rapidity integrated, and found that the resummed results significantly improve the agreement between theoretical predictions and data. Similar technique was also applied to high p ⊥ Higgs productions [27]. In the context of Higgs boson production, the threshold logarithms have also been included in the low transverse momentum resummation, see, for example, Refs. [28,29]. Furthermore, with the advent of the soft-collinear effective theory (SCET) [30], simple derivation of the factorization formula for the DIS structure function in the threshold limit x Bj → 1 can be achieved as shown in Ref. [31][32][33][34]. In addition, there have been interesting studies on the joint resummation [35][36][37] of transverse momentum logarithms and threshold logarithms, which resemble some similar physical idea as the joint resummation that we are presenting below. The major difference is that we are working in the dilute-dense factorization approach with fixed forward rapidity.
Let us take the example of the Higgs boson production in gluon-gluon fusion process with an effective Lagrangian [38][39][40] to demonstrate the formulation and factorization. As the main result of this paper, in the limit which the Higgs mass M is much greater than the measured transverse momentum k ⊥ , the factorization formula for forward rapidity Higgs production in pA collisions can be written as (1) is a very interesting and elegant formula, which encodes collinear resummation in g(x, µ), small-x evolution as well as multiple interactions in S W W x A (x ⊥ , x ⊥ ) [41], transverse momentum resummation in the conventional Sudakov factor S Sud (M 2 , µ 2 b ) [42][43][44][45] and the threshold resummation in the new factor, namely the forward threshold jet function ∆(µ 2 , µ 2 b , ln x τ ). 1 C(α s ) represents the hard coefficient expanded in terms of α s without any large logarithms. Essentially, Eq. (1) also resums the threshold logarithms of type α s ln(1 − τ ) with τ ≡ M √ s e y → 1 in the forward rapidity (y) Higgs production. Interestingly, we find that the new forward theshold jet function ∆ satisfies a renormalization group equation (RGE) similar to the RGE first derived by Becher and Neubert [33,46] for jet functions in SCET, which can also be derived from the analytical solution of the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation at the end point [34].
When the anomalous dimension γ µ,b ⊥ ≡ N c z γ µ,b ⊥ −1 at the leading logarithmic level with z = ln x τ . Eq. (1) is insensitive to the choice of factorization scale (also the same as renormalization by choice) µ. If µ 2 is set to be µ 2 b , it is straightforward to show that ∆ becomes δ(ln x τ ), which reduces Eq. (1) to previous results obtained in Ref. [43]. By setting µ 2 = µ 2 b in the collinear gluon distribution g(x, µ), we implicitly encode the threshold logarithms in the collinear gluon distribution by changing its scale with b ⊥ . This requires the information of collinear PDF g(x, µ) over a large range of µ at given x. Since PDFs have large uncertainties in the large x and large µ 2 region, it is better to choose µ 2 to be a constant, and perform the threshold resummation explicitly as in Eq. (1), which is presumably more stable and accurate in the threshold limit. Furthermore, when b ⊥ gets large, µ b can become smaller than the lowest scale that collinear PDFs are defined. Without a cutoff prescription in the large b ⊥ region, this can also be a problem. In addition, in the NLO single hadron production, normally the factorization scale is chosen to be a constant in practice in Ref. [14,16] in order to avoid some technical difficulties in the numerical calculation. In addition, more close connections to the SCET formalism can be established in forward hadron(jet) productions. This will be presented in a separate work with numerical results.
The threshold resummation with respect to the double differential (rapidity and transverse momentum) cross section considered here, as shown in the LHS of Eq. (1), has been studied previously [47][48][49] in different physical framework and kinematical region. As pointed out in Ref. [49], the threshold correction becomes quite large in forward rapidity region y ∼ 4 due to small-x contributions. In our approach, Eq. (1) resums both small-x and threshold logarithms, which is in principle more suitable framework to describe forward rapidity particle productions. It would also be very interesting to see the arise of small-x logarithms and corresponding resummations in the forward rapidity production of jets following the framework developed in Ref. [49].
In the following section, we will derive ∆(µ 2 , µ 2 b , ln x τ ) in Eq. (1) and comment on the application of the threshold resummation in other processes in Sec.II. Before we conclude in Sec. IV, several comments regarding the forward threshold resummations are made in Sec. III.

II. FORWARD THRESHOLD RESUMMATION IN HIGGS PRODUCTION IN PA COLLISIONS
For the sake of simplicity, let us use the forward Higgs production as an instructive example to demonstrate the threshold resummation in saturation formalism. This process is the simplest one since there is no final state gluon radiation. In addition, one can see that the leading power contribution comes from a few diagrams, while the rest of graphs are suppressed by factors of k ⊥ M 1. We follow closely the calculation shown in Ref. [43], which shows that the LO contribution reads where resums the multiple gluon exchanges between the active gluon and the target nucleus. At LO, the transverse momentum k ⊥ of the produced Higgs solely comes from the small-x Weizsäcker-Williams (WW) gluon distribution represented by the Fourier transform of S W W x A [7]. At one-loop order, working in the leading power k ⊥ M 1 limit, we find the following corrections 2 where the terms in the first line of the above equation are conventional Sudakov type logarithms (O(α s ln i M 2 ) associated with the transverse momentum resummation [42], while the terms in the second line can give rise to the threshold logarithms in the τ → 1 limit. It is also interesting to note that the latter contribution depends on the splitting function P gg (ξ), since it is related to the renormalization of the collinear singularity in the gluon PDF. Despite the fact that the Sudakov resummation and the threshold resummation are both resummations with respect to soft gluon contribution, we will distinguish them in our work, since they have different physical interpretation.
Due to the presence of the plus-function and δ function in the gluon splitting function P gg (ξ), the end point contributions in the above one-loop corrections, in particular the threshold logarithms, become important when τ → 1.
The factorization in the threshold resummation can be illustrated and achieved in the Mellin space. Let us define the Mellin transform and inverse Mellin transform as follows where C represents the properly chosen contour which puts all the integrable poles to its left side. For sufficiently large N , the Mellin transform integral is dominated by the end point where x ∼ 1 − 1 N . In the Mellin space, the LO cross-section reads where g N ≡ 1 0 dτ τ N −1 g(τ ). The first line of the 1-loop result in Eq. (3) can be transformed similarly. The second line can be transformed as follows where one can find where ψ(N ) is the digamma function. In the large N limit, ψ(N ) has the following asymptotic expansion where B 2n is the 2n-th Bernoulli number. To proceed, we write P gg (N ) − ln N − γ E + β 0 where higher order terms O 1 N have been neglected. Because the resummation of threshold logarithms and Sudakov logarithms are with respect to soft gluons radiations which always factorize, it is expected that P gg (N ) exponentiates when arbitrary number of soft gluon emissions is resummed [50,51]. In the Appendix A, an alternative derivation of our results following the RGE idea in Ref. [34] is also provided. Therefore, we can arrive at the following resummed formula in the Mellin space Similar results have been obtained previously in the context of collinear calculations where both transverse momentum and threshold resummations are considered [28,29]. In the following, we provide a derivation in terms of the inverse Mellin transform for the logarithms associated with the threshold resummation in the above equation. For example, by performing the inverse Mellin transform, the final resummation formula can be cast into It is interesting to note that the inverse Mellin transform in Eq. (11) can be performed analytically by applying the following identity [52] , which is derived by integrating above and below its branch cut. Analogous to the analytical continuation of the gamma function, the above identity can be extended to the full complex plane. With this identity, the resummation result can be written as The resummation result as in Eq. (13) becomes singular when γ µ,b ⊥ ≤ 0 at x = τ . Nevertheless, we can use the following trick in terms of the analytical continuation to extend to the region −1 < γ µ,b ⊥ ≤ 0 In addition, occasionally, it is also useful to evolve PDFs backwards which means µ 2 µ 2 b . Then we need to further analytically continue to the region where γ µ,b ⊥ > −2 as follows In fact, by utilizing the same subtraction method, we can extend of the region of validity for γ µ,b ⊥ to any negative value, and even to −∞. At last, in the case of running coupling, the final resummed result can be written as is defined in the spirit similar to the so-called star-distribution [34,53]. The star-distribution should be understood as the last line of Eq. (14) in the region of −1 < γ µ,b ⊥ ≤ 0 and as in Eq. (15) in the region γ µ,b ⊥ > −2, which is normally sufficient for the purpose of numerical evaluations. 3 The final resummation result can be written into a compact form by introducing the threshold resummed gluon distribution g t (τ, µ b ) at the scale µ b where the forward threshold jet function ∆(µ 2 , µ 2 b , ln x τ ) is introduced for the purpose of threshold resummation. Plugging above expression into Eq. (16), we obtain the final result for Higgs boson production in forward pA collisions as described in Eq. (1) with C(α s ) = 1 + αs π π 2 2 N c + O (α s (1 − τ )) . Since the large N approximation is used in reaching above results, there are corrections of order α s (1 − τ ) which are neglected. The terms which are neglected are explicitly shown in Appendix. B.
We would like to emphasize that g t (τ, µ b ) is independent of the renormalization scale µ when τ is sufficiently close to 1. For quark distributions, similar resummation can be achieved once we replace N c and β 0 by C F and 3 4 , respectively. The off-diagonal channels are suppressed, simply because there are no plus-function or δ-function in the off-diagonal splitting functions. The consequence of the above formula is that it resums important contributions and restores predictive power in the threshold limit. It is not coincidence that the µ dependence in the collinear PDF g(x, µ) is offset by that in ∆(µ 2 , µ 2 b , ln x τ ). In fact, the convolution in Eq. (17) is determined by the end point limit of the DGLAP evolution equation [34]. Eq. (17) is a useful formula since it can provide us PDFs at any scale in the threshold limit. Detailed discussion regarding this issue and comparison with previous results are provided in the Appendix A. As to the threshold logarithms associated with fragmentation functions, an equivalent corresponding equation can be written with respect to fragmentation functions as well. In fact, as far as we know, a simplified version of this formula for valence quarks first appeared in a review paper [55] in the beginning of QCD.

III. COMMENTS ON THE FORWARD THRESHOLD JET FUNCTION
Before we conclude, several comments with respect to the forward threshold jet function ∆(µ 2 , µ 2 b , z = ln x τ ) and Eq. (16), which is our main result, are in order.
• With the identity regarding the digamma function ψ(γ) = −γ E + 1 0 du 1−u γ−1 1−u , it is straightforward to check that the forward threshold jet function introduced above 3 Analogous to the plus distribution, the star distribution is defined as follows which can be easily shown to be equivalent to the expression used in this paper for the region −1 < η ≤ 0 if one simply sets p 2 = µ 2 ln x τ and Q 2 = µ 2 ln 1 τ . Here f (p 2 ) and g(x) can be any smooth test functions.
is the solution to the following non-local RGE It is very interesting to note that the above integro-differential equation almost coincides with the RGE [33,46] developed for jet functions in SCET, once we remove the Sudakov type logarithmic terms there and identify Γ cusp = αsNc π for gluons at one-loop level. To make the connection more manifest, the scale hierarchy Q ∼ M µ i ∼ µ Λ QCD vital to the usual threshold resummation also appears in our calculation once we identify the intermediate scale µ i as µ b , which always shows up in a coordinate space formulation. The form of the RGE in Eq. (19) is specifically related to the one-loop correction of this particular process. The solution of RGE automatically contains the corresponding resummation. This interesting link with SCET can be useful for us to perform threshold resummations for other processes and go beyond leading logarithmic level with the help of the RGE technique, since Γ cusp and β functions have been calculated as high as four [56] and five loops [57], respectively.
• Again when γ i ≡ N c > 0, motivated by the discussion in Ref. [55], we can easily prove that ∆(µ 2 , µ 2 b , ln x τ ) has the following interesting propagation property This can be interpreted as that the evolution from 1 to τ represented by ∆(µ 2 1 , µ 2 3 , ln 1 τ ) can be written as the convolution of two step evolutions after summing over intermediate states.
• As a matter of fact, through mathematical induction, one can prove that the formula involving the inverse Mellin transform is valid and well-defined for any value of γ µ,b ⊥ in the complex plane with . When γ µ,b ⊥ is small, the above series converges quite fast so that the sum of the first two terms is already close to the exact result. However, when |γ µ,b ⊥ | becomes large, the complete summation should be taken into account. As expected, for γ µ,b ⊥ > 0, the above summation over k can be easily performed which gives However, when γ µ,b ⊥ becomes negative, we can no longer resum all the terms in Eq. (22). Instead, the summation over k starts from the first integer value with k > −γ µ,b ⊥ . This naturally explains the subtraction method NLO τ = 0.85 , µ 2 = 1.5GeV 2 NLO τ = 0.85 , µ 2 = 25GeV 2 NLO+0.5, τ = 0.7 , µ 2 = 1.5GeV 2 NLO+0.5, τ = 0.7 , µ 2 = 25GeV 2 NLO+1, τ = 0.5 , µ 2 = 1.5GeV 2 NLO+1, τ = 0.5 , µ 2 = 25GeV 2 for different values of τ and µ 2 calculated at NLO. All the ratios are rather close to unity. In order to separate the curves for τ = 0.5 and τ = 0.7 from the τ = 0.85 curves, we add 1 and 0.5 to the ratio for τ = 0.5 and τ = 0.7 curves, respectively. employed above and the origin of the star-distributions. For example, from Eq. (22), it is now straightforward to find that the following formula gives the star distribution in the region γ µ,b ⊥ > −3 where ∆g (3) (x, τ, µ) = g(x) − g(τ ) − g (1) (τ ) ln x τ − g (2) (τ ) ln 2 x τ . Finally, we close this section by showing some numeric results for the threshold jet functions at various value of x (τ ). In Fig. 1, we plot the ratio between g t (τ, µ b ) computed from Eq. (17) and g(τ, µ b ) at the same scale µ b for τ = 0.5, 0.7, 0.85, respectively. We use the MSTW gluon PDF [58] g(x, µ) in Eq. (17) as an input to obtain g t (τ, µ b ), and then we use again the MSTW gluon PDF g(τ, µ b ) as the denominator. The ratio is expected to be flat and close to unity when τ → 1, since g(τ, µ b ) is the solution to the exact DGLAP equation and g t (τ, µ b ) is derived from the endpoint limit of the DGLAP equation. Here, in order to get better numerical agreement, the NLO threshold resummed curves are computed with NLO running coupling and NLO DGLAP equation as given in the end of the Appendix A. Indeed, we observe that, in the range from µ 2 b = 1GeV 2 to 1000GeV 2 , g t (τ, µ b ) yields agreement with the MSTW gluon PDF at the same scale with roughly less than 10% deviations. Also, we vary the scale µ 2 from 1GeV 2 to 25GeV 2 , and find that g t (τ, µ b ) is insensitive to µ 2 choices as expected. (It seems that the scale dependence gets stronger when τ gets large. This is due to the fact that the scale evolution gets much more rapid when τ approaches 1. Nevertheless, we do see that the curves get more and more flat when τ → 1.) More interestingly, as shown in the curves with µ 2 = 25GeV 2 , we obtain excellent agreement with the PDF in the region µ 2 b < µ 2 as well, where γ µ,b ⊥ becomes negative and Eq. (15) has to be employed. This proves that the analytical continuation technique works as expected for the so-called backward evolution. Since the MSTW PDFs is provided above 1GeV 2 , we decide to make a cut at 1GeV 2 , although we believe Eq. (17) can be used to evolve g(x, µ) from µ down to the scale µ b , which can be less than 1GeV 2 .
We would like to emphasize that the above comparisons help to establish the validity of the threshold jet function, which can be applied to other forward scattering processes. In particular, in forward hadron production in pA collisions, the differential cross section can be written in terms of the parton distribution at the scale of µ b , i.e., f (a) (x, µ b ) where a represents a quark or gluon from the incoming nucleon [9]. The threshold resummation can be carried out, by applying the similar technique of this paper. Again, the final result can be cast into f (a) t (x, µ b ) multiplied by the resummation result for the hard coefficients, which will be different from the current case. In the current case, the hard part, which only depends on a term proportional to the splitting function and an α s correction with no threshold logarithms, is quite simple. However, in forward hadron production case, there exist other large threshold logarithms in the hard part, whose resummation will be also important for a precision calculation. We will leave that for a separate publication.

IV. CONCLUSION AND OUTLOOK
In this paper, we have demonstrated the resummation of threshold logarithms in the dilute dense factorization which is widely used studying small-x effects in high energy collisions by using complex analysis and RGE methods. The framework discussed above resembles a lot of similarities to the threshold resummation in SCET and traditional resummation in the Mellin space. The advantage of this approach is that final results can be expressed in momentum space analytically with decent numerical accuracy. To obtain more precise results, we believe that we need to adopt the approach derived in Ref. [28] which resums the full DGLAP splitting functions including the off-diagonal channels in the Mellin space and performs the inverse Mellin transform numerically. This also pave the way for future applications of this resummation technique in pA collisions can make phenomenological calculations in dilute dense factorizations more reliable and systematic in forward rapidity particle productions.
where P ξ→1 gg (ξ) ≡ 1 (1−ξ)+ + β 0 δ(1 − ξ) which is equivalent to τ → 1 and N 1 limits. In the Mellin space, the above equation becomes The exact solution can be written as where γ µ,µ f ≡ N c . Although we have not been able to find an analytical form for the inverse Mellin transform in Eq. (A3), we can evaluate it numerically with any positive λ when γ µ,µ f > 0. Furthermore, if one approximates ψ(N ) as ln N in the large N limit by using Eq. (9), one can find Eq. (A3) becomes the results in Eq. (17). It is also important to note that the same level of approximation has been made along the way when use the DGLAP equation at the end point. In addition, it is interesting to note that the exact solution in Eq. (A3) has the following bound which has been checked numerically for sufficiently large τ when γ µ,µ f > 0. When γ µ,µ f ≤ 0, analytical continuation has to be applied again. In the threshold limit, the exact solution is tightly bound by Eq. (A4). We have also numerically checked that the above bound (especially the upper bound) usually gives excellent numerical estimate of the exact solution. Again, similar result can be obtained for quark distributions once we change to the quark splitting function accordingly. Let us compare our result in Eq. (17) to Eq. (3.31) in Ref. [34], which can be rewritten as (converted into gluon channel in our notation) First of all, our numerical evaluation of Eq. (A3) indicates a finite difference from Eq. (A5) for τ < 1. Second, we find that the above equation is equivalent to Eq. (17) in the limit τ → 1 with corrections of order (1 − τ ). Therefore, our result, which provides a slightly different analytical formulation, is complimentary to Eq. (3.31) in Ref. [34]. It is straightforward to generalize the above calculation up to NLO with the NLO DGLAP splitting function. In the ξ → 1 limit, the gg channel NLO DGLAP splitting function reads [59,60] The NLO results for Eq. (17) can be written as Appendix B: Discussion on corrections at one-loop order Since only the dominant part of the one-loop contributions are resummed in this calculation, the difference between the exact one-loop contribution the resummed part, which vanishes in the limit τ → 1, can be computed as follows.
For the forward Higgs production in the pA collisions, the leading power one-loop contribution is proportional to where the term αs π π 2 2 N c is excluded, since it is always kept in C(α s ). For the threshold resummation formula in Eq. (17), the corresponding one-loop contribution can be obtain easily by expanding it up to first order in α s (namely, γ µ,µ b ). Alternatively, it is instructive to start from Eq. (11) after expanding exp −γ µ,b ⊥ ln N e γ E e β 0 1 − γ µ,b ⊥ ln N e γ E e β 0 . Therefore, the one-loop contribution which is resummed in Eq. (11) can be explicitly computed as follows Using the identity ln N = lim →0 The above identity can be derived in terms of expansion for any test function f (x) Therefore, it is straightforward to find that the final result is finite and it reads Then, the contributions which are not resummed in this approach can be cast into Due to the cancellation between the plus-function and the star-distribution at the end point where ξ = 1 or x = τ , we can find that the difference is no longer singular and the corresponding remaining contribution vanishes in the limit τ → 1 for fixed γ µ,b ⊥ . As a comparison, let us consider the resummation formula in Eq. (A5) derived in Ref. [34]. Similarly, one can find that the corresponding one-loop contribution and finite differences read The above difference also vanishes when τ → 1. In addition, one also needs to note that there is off-diagonal contribution from quark to gluon splittings, which is again of order α s (1 − τ ).