Four-loop splitting functions at small x

We consider the expansion of small-x resummed DGLAP splitting functions at next-to-leading logarithmic (NLL) accuracy to four-loop order, namely next-to-next-to-next-to-leading order (N3LO). From this, we extract the exact LL and NLL small-x contributions to the yet unknown N3LO splitting functions, both in the standard MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} scheme and in the Q0MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} scheme usually considered in small-x literature. We show that the impact of unknown subleading logarithmic contributions (NNLL and beyond) at N3LO is significant, thus motivating future work towards their computation. Our results will be also needed in future to match NLL resummation to N3LO evolution. In turn, we propose an improved implementation of the small-x resummation and therefore release a new version of the resummation code (HELL 3.0) which contains these changes.


Introduction
The data thus far collected by the experiments at the CERN Large Hadron Collider (LHC) have reaffirmed the Standard Model as a remarkably successful theory of fundamental particles and their interactions. Thus, in absence of striking signature of new physics phenomena, the theoretical community is compelled to perform calculations with ever smaller uncertainties so that predictions with ever increasing accuracy and precision can be compared to data of outstanding quality, thereby exposing subtle differences and discrepancies that may reveal the presence of physics beyond the Standard Model.
In the context of strong interactions, accuracy is usually achieved by computing predictions that include an increasing number of terms of the perturbative expansion in the strong coupling α s (henceforth, the fixed-order expansion). Leading-order (LO) cross-sections in QCD can be computed for an essentially arbitrary number of external particles. Automation has been achieved in recent years also for NLO calculations and an increasing number of NNLO calculations is now available in computer programs. Moreover, for hadron-collider processes with simple topologies, recent milestone calculations have achieved N 3 LO accuracy [1,2]. This is particularly important because the main production channel of the Higgs JHEP06(2018)145 boson, i.e. gluon-gluon fusion [3][4][5], falls under this category. Furthermore, precise theoretical predictions for LHC processes also require precise and reliable parton distribution functions (PDFs). In particular, the lack of a N 3 LO determination of PDFs is an important source of uncertainty on the Higgs cross-section [6]. Although a global determination of such PDF set cannot be foreseen in the near future, several ingredients are either already available, or focus of current research. For instance, deep-inelastic scattering (DIS) coefficient functions with massless quarks have been known at three loops for a long time [7], and a lot of progress has been done in the context of heavy quarks, e.g. [8][9][10][11][12]. Another important ingredient of this rather ambitious task is the determination of the DGLAP kernels, which control the scale dependence of the PDFs, at N 3 LO. Recent progress with four-loop splitting functions [13,14] suggests that this calculation could be completed rather soon.
A complementary approach to the fixed-order expansion consists of exploiting all-order resummation. In the context of PDF determination, small-x (or high-energy) resummation is of particular relevance. Small-x resummation of DGLAP evolution is known to nextto-leading logarithmic accuracy (NLL) and it is based on the BFKL equation [15][16][17][18][19][20]. However, the proper inclusion of LL and NLL corrections is far from trivial, due to the perturbative instability of the BFKL evolution kernel. This problem has been tackled by more than one group in the 1990s, see for instance refs. [21][22][23][24], refs. [25][26][27][28][29][30] and refs. [31][32][33][34], and resulted in resummed anomalous dimensions for PDF evolution. These techniques have recently been applied in the context of PDF determination in refs. [35,36], where small-x resummation at NLL accuracy has been included in the PDF evolution and in the computation of DIS structure functions through the public code HELL [37,38].
It has been found that small-x resummation stabilises the perturbative behaviour of both evolution kernels and partonic coefficient functions, thereby improving the description of structure-function data at low x. In particular, it is well-known that potentially large logarithms at small-x are absent in NLO splitting functions due to accidental cancellations, while they start to contribute at NNLO. As a consequence, PDFs determined with NNLO theory improved by NLL small-x resummation differ rather significantly from the ones determined with NNLO alone. Furthermore, while at NNLO the most singular term in the gluon splitting function is of order α 3 s x log 1 x (the term with two logarithms being again accidentally zero), at N 3 LO the most singular term is of order α 4 s x log 3 1 x . Hence, the aforementioned instability at low x is very likely to become rather worse at N 3 LO. Small-x resummation would then be mandatory for improved precision. We note that in order to resum all small-x logarithms that appear at N 3 LO, one would have to consider NNLL resummation, which would be based on the three-loop BFKL kernel, which despite a lot of recent progress [39][40][41][42][43][44][45] is not yet fully known.
In this work we examine in some detail the fixed-order expansion of the NLL resummed splitting functions up to four loops. This exercise is interesting for several reasons. First, it enables us to predict the coefficients of the leading and next-to-leading small-x contributions to the yet-unknown N 3 LO splitting functions, thus offering either a strong check or a way of complementing the fixed-order result at small x. Second, because the resummation also includes subleading effects, mostly related to the running of the strong coupling, we JHEP06(2018)145 are able to assess the impact of unknown NNLL (or higher) contributions on the four-loop result. Third, although we predict that N 3 LO splitting functions will be unstable at small-x (much more than the NNLO ones), their inclusion will be most likely beneficial at moderate and large x, and therefore we conclude that the most reliable result in future will be obtained by using N 3 LO evolution provided it is supplemented by small-x resummation. The expansion of the resummation to O(α 4 s ) presented here is also a crucial ingredient for the N 3 LO+NLL matching procedure. Finally, by explicitly studying the behaviour of subleading contributions up to forth order in perturbation theory, we are able to identify a potential source of instability in our previous implementation of the resummation. We propose here an improved way of dealing with this class of subleading contributions and consequently we release a new version of the resummation code HELL 3.0, where these changes are implemented.

DGLAP evolution at small x
In this section we summarise small-x resummation of the DGLAP splitting functions. Small-x logarithms appear in the singlet sector and we have therefore to consider a 2 × 2 evolution matrix that couples together the quark singlet and the gluon. Currently, small-x resummation is known to NLL and we find convenient to express resummed and matched results as where P N k LO ij are the (k+1)-loop splitting functions and ∆ k+1 P NLL ij represent the resummed predictions P NLL ij minus their expansion up to order α k+1 s , namely Eq. (2.1) is valid, in principle, for any value of k. Matching of the resummation to NNLO (k = 2 in the above notation) was achieved in ref. [38] and later applied in refs. [35,36] for PDF determination. In this work we instead focus on the matching to the next perturbative order, namely N 3 LO (k = 3). We note, however, that in order to really improve the quality of the result, one should also increase the logarithmic accuracy of the resummation contribution so that no potentially large logarithm is left unresummed. Therefore, one would like to reach at least N 3 LO+NNLL: we will leave this rather ambitious goal to future work.
Small-x resummation of DGLAP evolution is usually performed in a conjugate (Mellin) space. Therefore, we define the entries of the anomalous dimension matrix in the singlet sectors as In this non-standard notation, usually adopted in the small-x resummation literature, the leading small-x logarithms of the form 1 x log k 1 x are mapped into poles in N = 0.

Brief recap of small-x resummation of DGLAP evolution
We now recall how the resummation of DGLAP splitting functions is constructed, mainly following ref. [38]. First, one considers the plus eigenvalue γ + (N, α s ) of the singlet anomalous dimension matrix eq. (2.3). This is resummed by first exploiting the duality between DGLAP evolution and BFKL evolution, and then supplementing the result by the resummation of a class of subleading contributions originating from the running of the strong coupling. Additionally, requiring the symmetry of the resummed BFKL kernel and imposing momentum conservation leads to perturbatively stable results. Since the knowledge of the BFKL kernel at N k LO allows the resummation of γ + at N k LL, at the moment we can only reach NLL accuracy, γ NLL + . Once the eigenvalue γ + is resummed, one proceeds with the resummation of γ qg . Its all-order behaviour at NLL is described by the equation [46,47] where the square-bracket notation is defined by the recursion [30,48] γ k+1 Note that in refs. [37,38] a variant of the resummation where r(N, α s ) → α s β 0 , which corresponds to a limit in which γ + is assumed to be proportional to α s , was used to infer an uncertainty on γ qg . Because only a finite number of coefficients h k are known, the implementation of the resummation, described in ref. [37], is only approximate. However, for not-too-large values of α s , the implementation is numerically stable and reliable. Simple power-counting at small-N shows that the quark anomalous dimension γ qg starts at NLL. Therefore, at this accuracy we have some freedom in how we choose the logarithmic accuracy of γ + appearing in eq. (2.4). In ref. [37] a dedicated anomalous dimension, denoted LL , was constructed specifically for this purpose. This LL anomalous dimension is essentially a LL anomalous dimension, but its singular structure, which at resummed level is encoded in the position of the rightmost pole in N space, is taken from the NLL result. The reason for using this hybrid object can be summarised as follows: • on the one hand, it is preferable to use the γ NLL + anomalous dimension in order to avoid singularity mismatches between different entries of the anomalous dimension matrix; • on the other hand, since using γ + at LL in eq. (2.4) is formally sufficient to achieve NLL accuracy in γ qg , it was convenient from a numerical point of view to use as much of the LL result as possible, because of its better numerical stability.
However, recently in ref. [38] various improvements in the construction of the resummed anomalous dimensions have been proposed and implemented in the numerical code HELL. With these developments, the computation of the full NLL anomalous dimension is faster and much more stable and reliable, and it is therefore now possible to either use the hybrid

JHEP06(2018)145
LL result or the full NLL one in eq. (2.4). In particular, the latter choice corresponds to the original approach of ref. [30]. We will explore the effects of both options in the following. All the other entries of the singlet anomalous dimension matrix can be derived from the plus eigenvalue and γ qg . In particular, the results can be written in a rather simple form if we consider the resummed contributions ∆ k γ ij , defined according to eq. (2.2): We recall that these relations are able to predict only the LL part of γ gq , while we do not have enough knowledge to predict its NLL part, which is then only approximate in the equation above. Having all the anomalous dimensions in the singlet sector, we can construct the splitting functions by Mellin inversion. 1

Perturbative expansion of the resummation
We now consider the perturbative expansion of the resummed result presented in the previous section. The goal is twofold. On the one hand, the expansion of the resummation is needed to construct the resummed contributions ∆ k γ NLL ij , ∆ k P NLL ij , eq. (2.2), namely for the matching of the resummed result to fixed order. On the other hand, the α s expansion of the resummed results provides a prediction for the small-x behaviour of the fixed-order splitting functions.
In ref. [38] we have already determined the expansion of the NLL resummed splitting functions to O(α 3 s ), which was needed to match resummation to NNLO. In that case, there was no point in using the result of the expansion to predict the NNLO behaviour at small-x, as the three-loop splitting functions are known [49,50]. In this work we push the expansion to one extra order, O(α 4 s ). These results would be needed to match resummation to N 3 LO, and specifically to construct N 3 LO+NLL resummed results. The four-loop splitting functions, however, are not yet fully known [13,14]. Therefore, at the moment our results can be used to construct approximations, valid at small x, of the unknown N 3 LO splitting functions, or simply to supplement the ongoing computation with the knowledge of the exact small-x behaviour. In future, when the four-loop splitting functions will be computed, it will also serve as a cross check.
In order to obtain the expansion of the resummed entries of the anomalous dimension matrix, we have to expand both γ NLL + and γ NLL qg to the desired accuracy; the other anomalous dimensions are recovered using eqs. (2.7). Let us first introduce a generic notation for the expansion of the plus anomalous dimension, In phenomenological applications, a damping at large x is added to make the transition from the small- x region (where the resummation is relevant) to the large-x region (where the fixed-order description is appropriate) as smooth as possible, and finally momentum conservation is reimposed. These details have been described in ref. [38] and are not repeated here.

JHEP06(2018)145
which is valid both for the NLL anomalous dimension and for the auxiliary LL one. The expansion of the qg anomalous dimension, according to eqs. (2.4), (2.5) and (2.6), is given by where the anomalous dimensions γ 0,1,2 are the coefficients eq. (2.8) of the expansion of either γ NLL + or γ LL + , depending on the choice adopted in eq. (2.4). Note that the expansion eq. (2.9) depends on a parameter T . This parameter has been introduced to account for the two variants of the resummation of running coupling contributions described above. In particular, when r(N, α s ) as given by eq. (2.6) is used then T = 2, while for the variation r(N, α s ) → α s β 0 then T = 1. The first h k coefficients appearing in eq. (2.9) are [47] Note that the knowledge of γ NLL qg to O(α 4 s ) requires the expansion of the plus anomalous dimension up to O(α 3 s ), as it does not depend on γ 3 . This is due to the fact that γ qg is a pure NLL quantity, as it is clear from the factor of α s in front of eq. (2.4).
We thus need to compute the first four orders of γ NLL + , while we just need the first three orders of γ LL + as it only possibly enters in the expansion of γ qg . The precise construction of these resummed anomalous dimensions was presented in detail in ref. [38], and we do not repeat it here (some details are given in appendix A). We just recall that due to the actual construction of the plus eigenvalue, which is based on the duality between DGLAP and BFKL evolutions, the LL and LL resummed anomalous dimensions automatically contain the fixed LO anomalous dimension, while the NLL anomalous dimension contains the NLO one. 2 However, the qg anomalous dimension, eq. (2.4), requires a purely resummed anomalous dimension, which goes to zero at large N , in order to avoid producing spurious large-N terms (see discussion in appendix B.2 of ref. [37]). Therefore, we first define the resummed contributions ∆ 1 γ LL + and ∆ 2 γ NLL + to be the resummed results at LL and NLL minus the LO and NLO anomalous dimensions, respectively (the notation is the same as in ref. [38]). Then, we construct the purely resummed LL and NLL anomalous dimensions as The functions γ LL 0 , γ NLL 0 and γ NLL 1 are not fixed by the resummation, and we thus have a degree of arbitrariness in how to define them. Instead, the expansions of the resummed contributions

JHEP06(2018)145
can be derived from the resummed results of ref. [38]. Their computation is presented in appendix A. Part of them, specifically γ LL 1 and γ NLL 2 , have been already computed and presented in ref. [38]. The next terms, γ LL 2 and γ NLL 3 , are reported here for the first time. Starting from LL resummation, the first order of the anomalous dimension was chosen in ref. [37] to include the LL and NLL contributions of the LO anomalous dimension. The NLL term, being it a constant at this order, is further multiplied by a function 1/(N + 1) to make it vanish at large N . The expression, which we adopt also here, is where a ij are defined in eq. (A.3). The next orders, as predicted by the resummation, are given by (see appendix A) where all the coefficients are defined in appendix A. Eq. (2.15) is a new result. Note that, by construction, both γ LL 1 and γ LL 2 vanish in N = 1, as the resummation is built to preserve momentum conservation.
Moving to the NLL resummation, we need to choose both the LO and NLO contributions. We decide to adopt the same strategy as in the LL , i.e. keeping only the LL and NLL contributions from the fixed orders. In particular, for the LO term we use exactly the same approximation used for the LL , and at NLO, due to the fact that the LL term is accidentally zero, we simply have a NLL term, where a 21 , the NLL coefficient of the NLO, is defined in eq. (A.3). In γ NLL 1 we have also included a subtraction term of the same form as the NLL term in γ NLL 0 , which restores momentum conservation, i.e. γ NLL 1 (1) = 0. We have decided to add this feature as the effect is formally NNLL, and it makes it more in line with the LL case where the O(α 2 s ) term vanishes in N = 1. We stress that we have played with variants of both γ NLL 0 and γ NLL 1 , adding momentum conservation to the first, relaxing it in the second, varying the way it is implemented, and so on: the effect at the level of the splitting functions is moderate. The third and fourth coefficients are instead found expanding the resummation, and their form where again the various coefficients are defined in appendix A, and we have left implicit the functionχ 1 (0, N ), eq. (A.49). Eq. (2.18) was already presented in ref. [38], 3 while eq. (2.19) is a new result of this study. Before moving further, some comments are in order. We observe that, due to accidental zeros of the LL singularity both at NLO and NNLO, the leading singularity at these two orders is the NLL one, namely 1/N and 1/N 2 , respectively. These NLL poles are predicted correctly by NLL resummation, so in particular γ NLL 2 has the exact leading singularity (γ NLL 1 has it by construction). In contrast, the leading singularity of both γ LL 1 and γ LL 2 is predicted by the resummation and thus, being the resummation just accurate at LL, is not exact. While the full LL anomalous dimension, being an all-order result, is reliable, each term of its perturbative expansion may not be. In particular, these two terms, γ LL 1 and γ LL 2 , do not contain anything of the exact result, because they are zero at LL. Since the impact of these two orders in the expansion of resummed splitting functions (in particular P qg and P qq ) may be substantial, we may expect that using LL resummation may give rise JHEP06(2018)145 to somewhat unreliable resummed contributions when matched to high orders, e.g. NNLO or N 3 LO. More precisely, we may expect NLL resummation, which has the exact leading contributions at these two orders, to lead to more reliable matched results to high orders.
These considerations suggest that the use of the NLL anomalous dimension in the construction of γ qg (as originally suggested in ref. [30]) is preferred. Thus, from now on we will adopt the NLL anomalous dimension as default ingredient for the resummation of anomalous dimension matrix, and possibly consider LL resummation as a variant to estimate the impact of subleading logarithmic contributions. Both options are now available in the new 3.0 version of the HELL code.

The four-loop splitting functions at small-x in the Q 0 MS scheme
The results of the previous section allow to construct all ∆ 4 γ NLL ij , and thus by Mellin inversion all ∆ 4 P NLL ij , needed to match NLL resummation to N 3 LO evolution. While these results will be in practice of no use until the computation of the N 3 LO splitting functions will be completed, it is interesting to extract from them the small-x terms at LL and NLL of the yet unknown four-loop splitting functions. These may be also useful to construct approximate predictions of the four-loop splitting functions while waiting for their full computation.
We remind the reader that the resummation procedure previously described lead to all-order results which are not in the traditional MS scheme, but rather in a related factorization scheme called Q 0 MS [39,47,51,52], which is particularly suitable for small-x resummation. Indeed, in the MS scheme there are some cancellations of large small-x contributions taking place between parton evolution and coefficient functions. The Q 0 variant of the scheme automatically removes these large contributions from both objects, leading to resummed predictions which are perturbatively much more stable. For this reason, here as well as in previous studies, we always use the Q 0 MS scheme for results including small-x resummation. The difference between the MS and Q 0 MS factorization schemes influences the resummation of the anomalous dimensions beyond the leading logarithmic accuracy, as well as the resummation of the coefficient functions. We first concentrate on Q 0 MS, while we present results for MS in the next section.
The small-N expansion of the expansion terms of γ NLL + is given in appendix A. With those results, we can construct all the other entries using eqs. (2.9) and (2.7). Denoting with γ

JHEP06(2018)145
The gq and qq anomalous dimensions are obtained by simply multiplying the gg and qg ones by C F /C A , even though we stress again that only the 1/N 4 pole of the resulting gq anomalous dimension is correct. Hence we have The corresponding small-x logarithms in the four-loop splitting functions are easily obtained by Mellin inversion of eqs. (2.20), according to where M −1 denotes the inverse Mellin transform. Thus, four-loop splitting functions exhibit a much stronger growth at small-x than those at a previous order, which behave like α 3 s x −1 log x. To our knowledge, the NLL contribution to P gg is explicitly presented here for the first time.

The four-loop splitting functions at small-x in the MS scheme
The effect of scheme change between Q 0 MS and MS turns out to be of relative order α 3 s , and thus all the fixed-order results considered in previous studies on small x resummation happened to be identical in either scheme. However, in this work we are considering resummation matched (or expanded) to N 3 LO, thus becoming sensitive to the scheme choice even at fixed order. It is thus important to recall how the conversion is performed. The goal of this subsection is also to provide the small-x contributions of the N 3 LO splitting functions in the MS scheme, namely in the scheme in which the full four-loop computation will likely be carried out.
A factorization scheme change is a multiplicative redefinition of the PDFs f and coefficient functions C. Focussing on MS and Q 0 MS, and considering both processes with one or two hadrons in the initial state (i.e. coefficient functions with one or two flavour indices), we have where α s = α s (Q 2 ) and we denoted with a MS label quantities in that scheme and without label quantities in the Q 0 MS scheme. Accordingly, the anomalous dimensions change as The function Λ ij is a matrix in flavour space implementing the scheme change. As far as small-x scheme changes are concerned, this matrix is trivial in its non-singlet part, so we

JHEP06(2018)145
focus only on the singlet. Up to NLL, its form is given by [30] where both R and v are LL functions, i.e. functions of α s /N to all orders. The form of the LL part of the matrix is such that the scheme change has no effect on the LL part of the anomalous dimension matrix. The three empty slots in the NLL part of the matrix have an effect only on the γ gq entry at NLL, which is however not determined by NLL resummation (as we already stressed), and are thus of no relevance. Furthermore, they also affect the resummation of partonic coefficient functions but this effect is beyond the accuracy currently achieved in the context of small-x resummation. Thus, at the currently available logarithmic accuracy, the three empty slots can be any LL function of α s /N .
The scheme-change function was calculated long ago [47] R  .27), where γ + is the resummed one (in either scheme, as only its LL part needs to be correct, and at LL the scheme change is ineffective on the anomalous dimensions). The form of v, eq. (2.29), is such that at NLL both γ qg and γ qq are the same in MS and Q 0 MS. Additionally, we have already noted that the matrix structure of Λ ij at LL is such that is has no effect on the LL anomalous dimensions. Thus, the only anomalous dimension which is sensitive at NLL to the scheme change is γ gg , and we have having used the expansion of the function R in powers of M , The scheme change contribution is entirely due to the derivative term in eq. (2.26). Of course, also the NLL part of γ gq changes by the scheme change, but as we already repeated several times NLL resummation is not able to predict it. To conclude, we report the actual JHEP06(2018)145 expansion to NLL of the gg anomalous dimension in the MS scheme: At this accuracy, all the other entries are identical to their Q 0 MS counterparts given in section 2.3. To our knowledge, the NLL contributions to P gg in the MS scheme are explicitly presented here for the first time.

Numerical results and discussion
Thus far we have presented analytical results. We now concentrate on numerics and we illustrate the difference between the two variants of the resummation discussed above. We also present approximate results for the four-loop splitting functions which are based on the expansion of the resummation and we critically assess the trustworthiness of this construction.

Resummed splitting functions at NNLO+NLL
First, we consider the four singlet splitting functions at fixed order and with resummation using the NLL anomalous dimension, which is the new default in HELL 3.0. In figure 1 we show P gg (blue) and P qg (orange) in the left plot, and P gq (blue) and P qq (orange) in the right plot, at LO (dotted), NLO (dashed), NNLO (dot-dot-dashed) and NNLO+NLL (solid). The resummed result is supplemented with an uncertainty band, which aims to estimate the impact of unknown subleading logarithmic contributions. Following ref. [38], this band is obtained by considering variations of the way RC resummation of γ + is implemented and of the way the resummation of γ qg is performed, and summing in quadrature the two effects. 5 The qualitative aspect of these results is the same of those obtained with the HELL 2.0 settings of ref. [38], i.e. using the LL anomalous dimension. 6 To better appreciate similarities and differences, we compare the two variants of the resummation in figure 2, focussing on P qg on the left and on P gg on the right. The current default, denoted with "NNLO+NLL" (solid red) is compared to the choice we made in ref. [38], which has been labelled "NNLO+NLL (LL )" (dot-dashed green). For completeness, we also show fixed-order results (gray). We see that the difference for P gg are very 5 In fact, in ref. [38] we considered only the second variation for Pqg and Pqq; we now use a more symmetric approach and use both for all the splitting functions. 6 We warn the reader that we have discovered a bug in the implementation of our NLL results. The numerical impact is not dramatic and it is discussed in detail in appendix B. All numerical results presented here, including the ones with HELL 2.0 settings (i.e., using the LL anomalous dimension), have been obtained with the corrected implementation of the resummation. small and well within the uncertainty band. The difference for P qg is not large either, even though the uncertainty bands are smaller and comparable in size with such a difference. Similar considerations hold for the other splitting functions because, as it is clear from figure 1, P qq and P gq behave similarly to P qg and P gg at small x, due to the colourcharge relations, eqs. (2.7). We have checked that the comparison between the predictions obtained using the NLL anomalous dimension versus the LL remains equivalent also at NLO+NLL accuracy.
A striking feature of the result is the small size of the uncertainty band that we obtain for P qg . This is rather counterintuitive because P qg (and P qq ) start at NLL and they are therefore known only at their leading non-vanishing logarithmic accuracy. Perhaps this signals the limitation in our way of estimated theoretical uncertainties, which is currently purely based on the variation of subleading contributions related to the running of the strong coupling. Thus, in order to better assess the impact of subleading logarithms, it would be important to push the accuracy of the resummation, at least for the quarkinitiated splitting functions, one logarithmic order higher.
We now consider the resummation matched to one order higher, in view a future combination with N 3 LO splitting functions. Also in this case, we compare the two variants of resummation, which led to very similar results when matched to NNLO. We have already argued from theoretical grounds that usage of the LL variant is less favourable because one has less control over the subleading poles that appear in the expansion of the resummation. In figure 3 we see that this worry is indeed justified. In these plots we show the resummed contributions x∆ 4 P qg (on the left) and x∆ 4 P gg (on the right), which would have to be added to the N 3 LO splitting functions according to eq. (2.1). The solid red curve denotes the default resummation in HELL 3.0 based on the NLL anomalous dimension, while in dot-dashed green we show the LL variant. Both plots show an issue of the LL , which was absent when matching at lower orders: the resummed contributions give rise to a (most likely) spurious contribution at moderate x, which is instead absent if full NLL is employed in the resummation. Indeed, for x 10 −3 ÷ 10 −2 we expect to be outside the resummation region, and the effect of resummation should be smaller compared to the fixed-order contribution, which is more reliable in this region. This behaviour is violated in ∆ 4 P ij when using the LL anomalous dimension, due a large contribution of γ LL 2 , eq. (2.15), entering in eq. (2.9), which makes ∆ 4 P ij even larger than ∆ 3 P ij in this region, despite it being of higher order in α s . Consistently, the green curve also has a rather large uncertainty band in that region, which makes it almost compatible with the red curve, which has instead a smaller uncertainty, as one would expect. We interpret this behaviour as further confirmation of the aforementioned theoretical arguments in favour of our new default. However, it should be stressed once again that both curves feature the same logarithmic accuracy, and hence this discrepancy contributes to our theoretical uncertainty.

Approximate N 3 LO
We can use the expansion of the resummed splitting functions to O(α 4 s ) to make an approximate prediction of the N 3 LO splitting functions, simply by adding it to the exact NNLO ones. This is shown in solid red in figure 4 for P qg (left plot) and P gg (right plot), and in dot-dashed green for the LL variant. According to the discussion in the previous section, the latter curve is not expected to be accurate in the region of moderately large x, where it has an unphysically large effect. These predictions are further supplemented by the same uncertainty band that appears on resummed results, which happens to be invisible for P qg when using our new default implementation 7 (and thus confirms that our uncertainty band underestimates the actual uncertainty on P qg ). Additionally, we also show in light solid black the asymptotic small-x behaviour at N 3 LO, as obtained by adding to the exact NNLO the pure NLL contributions, eq. (2.20), without any subleading effects. Note that these results are obtained in the Q 0 MS scheme. While we have not implemented the resummation (and hence its expansion) in MS, we can easily plot the asymptotic behaviour of the splitting functions in this scheme, exploiting the results of section 2.4. The quark splitting function P qg is unaffected, while the MS asymptotic result for P gg , eq. (2.32), is shown in solid cyan.
In the small-x limit, approximate predictions behave as their asymptotic expansions, by construction. The difference is due to subleading NNLL contributions, behaving as 1 x log 1 x at this order (a straight line in the plots). While these NNLL contributions are subleading at asymptotically small x, their effect is sizeable for all the x range shown in the plots, which is rather large, reaching x = 10 −9 . This is true in particular for P gg , 7 Indeed, at this order, our uncertainty band originates from the parameter T appearing in eq. (2.9), and on the potential dependence of the anomalous dimensions γ0,1,2 on the parameter T defined in appendix A. Since the NLL anomalous dimension is more precise than the LL one, none of the γ NLL 0,1,2 depends on T , while γ LL 2 does (and also γ NLL 3 , which contributes to Pgg). It is the latter (T ) dependence that generates the uncertainty bands in figure 4, while the T parameter variation has no appreciable effect at this order on any of the two variants.

JHEP06(2018)145
where the pure NLL asymptotic curve is very different from the approximate N 3 LO, so much that in order to display the asymptotic behaviour we had to plot xP gg (x) in a rather extended range. For this reason, we also added an inset which zooms in to the region 10 −5 < x < 0.1, most relevant for HERA and LHC physics, to better appreciate the perturbative behaviour of the fixed-order splitting functions. This shows once again that subleading contributions have a very important role at intermediate and moderately small values of x. Similar conclusions were reached some time ago in ref. [53] (see also ref. [54] for a similar study in the non-singlet case). This also suggests that the approximate N 3 LO prediction that we plotted has a huge uncertainty, likely larger than what we estimate with our uncertainty bands.
By comparing resummed results (figure 1) and their expansions (figure 2), we can conclude that, while the small-x contributions to P qg , as obtained from the expansion of the resummation, behave in a perturbative way, the prediction of the N 3 LO contribution to P gg , which is more directly sensitive to the (perturbatively unstable) BFKL kernel, is very different from its all-order counterpart, as it was the NNLO contribution. We must conclude that these approximate N 3 LO predictions cannot be regarded as a faithful estimate of the actual N 3 LO (especially for P gg and P gq ) due to potentially underestimated subleading contributions. However, what we can certainly conclude is that exact N 3 LO evolution (when available) will be unreliable at small x, and thus it will necessarily have to be supplemented with the all-order resummation of small-x contributions.
As we have argued before, because of the sensitivity of the N 3 LO splitting functions to subleading logarithmic contributions, it would be important to push the resummation to NNLL accuracy. However, NNLL resummation requires at least the knowledge of the NNLO BFKL kernel, which is so far only known in a collinear approximation [39]. It will be important in the future to explore the possibility of computing the BFKL kernel to NNLO [40][41][42][43][44][45], and perhaps to consider the option of using its collinear approximation.
Finally, it interesting to note that the asymptotic behaviour in MS appears to be closer to the all-order result (albeit computed in a different scheme). This suggests that a future study of the resummation in MS may reveal interesting properties in terms of the size of subleading contributions, despite the fact that in this scheme we expect stronger cancellations between coefficient functions and parton evolution.

Conclusions
In this paper we have discussed the role of higher-order corrections to the splitting functions, which govern the evolution of the parton distribution functions. In particular, we have exploited results in small-x resummation to study the behaviour of the yet-unknown four-loop splitting functions in the singlet sector. Our results stress once again the fact that small-x singularities lead to loss of perturbative stability, when higher orders are considered. This has been masked so far by accidental cancellations at NLO but it becomes apparent at NNLO, even though also there the strongest singularity is accidentally zero, thus slightly mitigating the perturbative deterioration. Instead, at the next orders there are no accidental cancellations, so the perturbative instability is no longer moderated, and

JHEP06(2018)145
we estimate it to be rather severe at N 3 LO. Thanks to this work, this potential instability can be solved by adding the resummation to the full four-loop splitting functions, when these will become available.
We have also investigated the possibility of using the expansion of the resummation to construct approximate N 3 LO splitting functions. Unfortunately, we have found that subleading corrections, which are only partially included in our approach, have a sizeable impact at moderate x, thus rendering the construction of approximate fixed-order splitting functions rather uncertain. However, our asymptotic results can be used as a check on the full four-loop calculation, or for complementing an approximate computation based on integer Mellin moments.
While performing these studies we have encountered a potential source of instability in the way the resummation of the quark anomalous dimension γ qg was implemented in HELL 2.0 in ref. [38], which was based on a hybrid resummation formula denoted LL . Therefore, we have adopted as a new default a resummation fully based on NLL and consequently released a new version of the resummation code HELL 3.0: www.ge.infn.it/∼bonvini/hell As the distinction between the two choices is beyond the accuracy of the calculation, the old option can be, and should be, still used to estimate theoretical uncertainties. Furthermore, we anticipate that an analogous situation appears in the resummation of partonic coefficient functions. This issue will be discussed in a forthcoming study [55]. Finally, these results have been recently exploited in a double-resummed calculation of the Higgs production cross section [56].

Acknowledgments
We thank Andreas Vogt for encouraging us to write up these results and Richard Ball for useful discussions. The work of MB is supported by the Marie Sk lodowska Curie grant HiPPiE@LHC.

A Perturbative expansion of resummed anomalous dimensions to N 3 LO
In this appendix we derive the expansion in powers of α s of the resummed plus eigenvalue of the singlet anomalous dimension matrix presented in ref. [38]. Specifically, we provide the detailed computation of the expansion of the NLL resummed anomalous dimensions up to O(α 4 s ), and of the LL anomalous dimension up to O(α 3 s ), as needed for matching NLL resummation in DGLAP evolution to N 3 LO. We recall that the anomalous dimensions at LL and NLL are constructed as [38] γ res LL (1, α s ) = 0, through a subleading function f mom defined in eq. (A.8). All these ingredients have been presented in ref. [38] and will be used in the following.
Before starting, we recall that the DL anomalous dimension is constructed starting from the fixed-order BFKL kernel matched to the fixed-order anomalous dimension. Thus, one of the ingredients of (N)LL resummation is the (N)LO anomalous dimension. In ref. [38] an approximate form for the input anomalous dimension was suggested to facilitate the numerical implementation and to solve a potential issue. The approximation does not represent any loss of accuracy, as the only requirements needed for the input anomalous dimension are to be accurate at NLL and to conserve momentum, both of which are satisfied in the approximation of ref. [38]. At LO and NLO they are given bŷ In the next, we start from LL resummation, and then move to NLL. The computation follows closely the one presented in section 3 and 4 of ref. [38], extending it to one extra order.

A.1 Expansion of the LL anomalous dimension
We start expanding the LL anomalous dimension up to O(α 3 s ). The first ingredient for resummation is the DL resummed anomalous dimension γ DL , which is obtained from the implicit equation The function χ Σ (M, N, α s ) is the so-called off-shell BFKL kernel [29,37]. For LL resummation it is given by [29] χ

JHEP06(2018)145
where the function χ s (α s /M ) is the dual of the LO anomalous dimension α sγ0 , is the off-shell extension of the LO BFKL kernel after subtracting double counting with χ s . The last term restores the momentum conservation constraint γ DL (1, α s ) = 0, namely by duality χ Σ (0, 1, α s ) = 1, through a function and with the coefficient The coefficient χ 01 appearing in eq. (A.7) is the first of the expansion of χ s , where in the last equality we have used the definition eq. (A.6), and the formulae for the derivatives

JHEP06(2018)145
which can be derived from the very same definition. The prime denotes a derivative with respect to the argument of the function, so χ s (1/γ 0 ) is a derivative with respect to 1/γ 0 , andγ 0 is a double derivative with respect to N . The anticollinear χ s gives instead The kernel eq. (A.5) expands as where the derivative is with respect to M , i.e. the first argument. Putting everything together eq. (A.4) brings to the expanded equality from which it immediately follows Note that the O(α 0 s ) term cancels automatically, becauseγ 0 in eq. (A.12) is the one used in the definition of χ s , eq. (A.6). The expansion terms of the off-shell kernelχ 0 (M, N ), eq. (A.7), are given bỹ which lead to the predictions As already noted in ref. [38], the fact thatγ 1 vanishes is not surprising: indeed, the LL pole of the exact NLO γ + are accidentally zero, so the only part which is supposed to be predicted correctly by LL resummation was indeed expected to vanish.

JHEP06(2018)145
Having computed the expansion of the DL part, we now move to the RC contributions. The function that resums the running-coupling effects is given by [38] γ rc (N, α s where k ν (z) is a Bateman function, with In the equations above, M min (α s ) is the position of the minimum of the BFKL kernel, 8 and c(α s ) and κ(α s ) are the value and curvature of the kernel at such minimum. In deriving eq. (A.23), the α s -dependence of the BFKL kernel has been approximated linearly, keeping the value of the kernel and its α s -derivative correct. In ref. [38] a variant of this approximation in which the kernel is assumed to be proportional to α s (i.e., as if it was a purely LO kernel) was considered to study the impact of subleading logarithmic contributions. This variant is recovered by letting c → c/α s , κ → κ/α s , and represents an equally valid alternative. The RC contribution to be added to the DL-LO result is given by where κ 0,1 and c 0,1 are the O(α 1,2 s ) terms of κ and c, and m 1 is determined from M min (α s ) = 1/2 + α s m 1 + O(α 2 s ). To cover both α s -dependence approximations, we have introduced a parameter T which equals 2 in the default approximation and equals 1 in the limit c → c/α s , κ → κ/α s . The values of κ 1 , c 1 and m 1 depend on the actual kernel used. For LL resummation, the DL-NLO is used for RC resummation, differently from pure LL resummation which uses the DL-LO kernel. Thus, the second order coefficients in eq. anomalous dimension, there is a mismatch in the singularities at small x which are cured by the matching function γ LO+LL match defined in ref. [38]. Its expansion gives where are the differences between the coefficients computed with the DL-NLO and the DL-LO BFKL kernels. Explicit results are given in section A.3.
Putting everything together, the expansion in powers of α s of the full LL anomalous dimension is given by Using eq. (A.11) and replacing the explicit form ofγ 0 (N ), eq. (A.2), we obtain the results presented in eq. (2.14) and eq. (2.15).

A.2 Expansion of the NLL anomalous dimensions
We now move to the NLL anomalous dimension. This time, we need to push our expansion to one order higher. The off-shell kernel needed for NLL resummation is  [37], and can be written as where the functioñ is regular in M = 0 and 1 − M + N = 0. The functionχ 1 is given by [37]

JHEP06(2018)145
The extra term χ corr 1 (M, N ) was corrected in ref. [38]; however, there was another issue, which we have discovered only now. This issue has an effect at NLL, namely the claimed accuracy, but it manifests itself (at NLL) only in terms of O(α 4 s ) and beyond of the anomalous dimension (for this reason, the comparison of the expansion of the anomalous dimension with the exact one at three loops was successful). The origin of the issue and our solution are discussed in detail in section B. The actual expression that we use here is given by Finally, the momentum conservation coefficient is given by We now consider the expansion of the DL-NLO anomalous dimension where implicit derivatives denoted with a are with respect to M , i.e. the first argument.

JHEP06(2018)145
The new functions appearing in the equations above are given by +γ 1χ 0 (0, N )+ 1 2γ where The N 3 LO functionγ 3 (N ) cannot be simplified significantly, so we do not manipulate it further.
We now move to the running-coupling contributions. The RC correction to the duality is implemented through the function Note that this expression differs from the analogous in ref. [38] by NNLL terms, due to the correction to the function χ corr 1 .

JHEP06(2018)145
which contributes at NLL. Further RC corrections from RC resummation, eq. (A.23), are NNLL corrections. They amount to where as before T is either 2 or 1 depending on the approximate α s dependence chosen.
To cancel the singularity mismatch between eq. (A.56) and eq. (A.55) we further need the matching function The coefficients above are the ones obtained from the NLO kernel, given in section A.3.
Putting everything together according to eq. (A.1), we obtain . Specifically, we compute the residues of all the poles in N = 0, −1, −2, and construct an approximation based only on these terms, This approximation has the advantage of describing in x space all contributions behaving as x k log j x for all non-vanishing terms with j ≥ 0 and for k = −1, 0, 1, namely all non-vanishing terms at small x plus the leading corrections to them. Of course this approximation will not be accurate at large x, but since the final results will be damped at large x the inaccuracy should be negligible. To verify the quality of our approximation, we have compared it with a simpler approximation which does not include the contributions from the poles in N = −2, i.e. those terms behaving as x log j x in x space. The difference between the two results is almost imperceptible. Before concluding, it is useful to extract from eq. (A.58) its small-x behaviour up to NLL, which provides a prediction for the yet unknown four-loop anomalous dimensions. Expanding the anomalous dimensions eq. (A.2) in N = 0,  10 We could not compute the derivatives analytically, so we have used the PSLQ algorithm to find the results of the integrals.

JHEP06(2018)145
Denoting with γ To our knowledge, the NLL term at N 3 LO was never explicitly presented in the literature. Note that the scheme is Q 0 MS. The conversion to MS was presented in section 2.4.

A.3 Computation of the coefficients describing the minimum of the BFKL kernel
In this section we compute the expansion coefficients of c(α s ) and κ(α s ) of the BFKL kernel in symmetric variables. We will assume that the minimum is not in M = 1/2, but iñ The expansion coefficients of the minimum can be found by solving perturbatively the minimum condition

JHEP06(2018)145
Note that the last term vanishes due to the fact that the LO kernel is symmetric, and hence all its odd M -derivatives are zero. Putting everything together, we then obtain where every off-shell kernel is implicitly assumed to be computed in M = 1/2, N = 0. Explicitly, we have at lowest order where the three contributions in square brackets correspond to the three items above, respectively. The last term contains a χ 0 (M ) in the denominator, which is singular in M = 1/2, and thus produces a new singularity that makes the result perturbatively unstable. It has been realised long ago [28,29] that these singular contributions can be removed (resummed) if the running coupling corrections to duality are accounted for to all orders rather than included perturbatively. This is the goal of the RC resummation, which provides an anomalous dimension where these contributions are resummed by solving the BFKL equation with running coupling (in a given approximation). Thus, it is convenient to translate the last correction to the kernel in a correction to the anomalous dimension, which is then regulated when the RC resummed result is added. This is always possible, and a simple computation reveals that one can pour the last term of eq. (B. This implies that the naive dual of the kernel (the symmetric kernel plus ∆χ, eq. (B.1), without the last term) does not coincide with the high-energy limit of the exact anomalous dimension at O(α s ). This is problematic when one wants to resum the collinear singularities of the kernel using χ s , the dual of the LO anomalous dimension, as the singularities would not cancel. One could in principle cure this problem by computing χ s from a modified fixed-order anomalous dimension to which the −α s β 0 term, eq. (B.3), has been subtracted, namelyγ 0 →γ 0 + α s β 0 . However, this becomes too artificial and not very physical. A better way to solve the issue is to pour back the lowest order term of the expansion of ∆γ into ∆χ. We would thus have (using the same names for the modified objects) In our previous implementation of the resummation, ref. [38], we used a different (wrong) form of χ corr 1 , taken from ref. [29], which differs from the correct eq. (B.6) by the last term. The singularity of the last term is the correct one, and thus leads to a proper resummation through χ s , which is the argument used by ref. [29] to introduce such a "subtraction" term. However, the missing M -dependent contributions at O(α 2 s ) produce spurious NLL contributions in the anomalous dimension. This NLL effect starts to appear at O(α 4 s ), which explains why the comparison of the three-loop anomalous dimension obtained from the resummation to the exact result was successful. For this reason, the full resummed result after the correction is not very different from the previous bugged one.
To conclude, we show in figure 5 the comparison between the P qg and P gg splitting functions obtained using the LL anomalous dimension before (dot-dashed gray) and after (dot-dashed green) the correction of the bug, and the new definition of M min , eq. (A.91). Both splitting functions appear to be harder after bug corrections, which is due to both γ NLL + and γ LL + being indeed harder. We stress however that much of this effect is induced by the change in M min , rather than to the correction in χ corr 1 . The uncertainty band on P qg is significantly reduced, again mostly due to the different M min used. Overall, however, the effect is not dramatic, and likely comparable to (if not smaller than) unknown subleading corrections, which in P qg are likely underestimated by our uncertainty band, as we already commented in section 3.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.