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 (N$^3$LO). From this, we extract the exact LL and NLL small-$x$ contributions to the yet unknown N$^3$LO splitting functions, both in the standard $\overline{MS}$ scheme and in the $Q_0 \overline{MS}$ scheme usually considered in small-$x$ literature. We show that the impact of unknown subleading logarithmic contributions (NNLL and beyond) at N$^3$LO is significant, thus motivating future work towards their computation. Our results will be also needed in future to match NLL resummation to N$^3$LO 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 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 next-to-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 structurefunction 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 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 fixedorder result at small x. Second, because the resummation also includes subleading effects, mostly related to the running of the strong coupling, we 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 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 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 threeloop 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, 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 App. 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 App. 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 can be derived from the resummed results of Ref. [38]. Their computation is presented in App. 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 App. A) 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, 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 is (see App. A) 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 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 App. A. With those results, we can construct all the other entries using Eqs. (2.9) and (2.7). Denoting with γ 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 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] .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 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 Sect. 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 Fig. 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 Fig. 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 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 Fig. 1, P qq and P gq behave similarly to P qg and P gg at small x, due to the colour-charge 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 quark-initiated 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 Fig. 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 Fig. 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 Sect. 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 , 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 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 App. 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 Fig. 4, while the T parameter variation has no appreciable effect at this order on any of the two variants.
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 (Fig. 1) and their expansions (Fig. 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 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łodowska 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 byγ In the next, we start from LL resummation, and then move to NLL. The computation follows closely the one presented in Sect. 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] where the function χ s (α s /M ) is the dual of the LO anomalous dimension α sγ0 , where in the last equality we have used the definition Eq. (A.6), and the formulae for the derivatives 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 As already noted in Ref. [38], the fact thatγ 1 vanishes is not surprising: indeed, the LL pole of the exact NLO γ (1) + and NNLO γ (2) + are accidentally zero, so the only part which is supposed to be predicted correctly by LL resummation was indeed expected to vanish.
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.  [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 Sect. 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 The function χ s,NLO (M, α s ) is the generalization of χ s to the next order, which is obtained as the exact dual of the NLO anomalous dimension, This kernel can be expanded as  [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] 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 Sect. 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 The anticollinear χ s,NLO expands as The kernels give instead where implicit derivatives denoted with a are with respect to M , i.e. the first argument. The new functions appearing in the equations above are given by Plugging each expansion in Eq. (A.4) we find 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 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 Sect. A.3.
Putting everything together according to Eq. (A.1), we obtain γ res NLL (N, α s ) = α sγ0 (N ) + α 2 sγ1 (N ) 9 Note that this expression differs from the analogous in Ref. [38] by NNLL terms, due to the correction to the function χ corr . 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, 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 Sect. 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ñ and further expanding we get 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 Note that ultimately one wants to construct the anomalous dimensions which are dual to the BFKL kernel in asymmetric (or DIS) variables. In Ref. [38] we argued that ∆ DL-(N)LO γ (N)LL rc , which contains the RC resummation corrections to be added to the DL result, could be computed directly from the kernel in symmetric variables, as the conversion from symmetric to asymmetric amounts to adding +N/2 to the anomalous dimension, which is then subtracted again by the matching. However, the argument was superficial, and induced by the collinear approximation of the kernel used for computing the resummation of RC effects. Indeed, the actual DL BFKL kernel in symmetric variables behaves asymptotically as −2M in the collinear region, which by duality (either fixed-coupling or running-coupling) reproduces the −N/2 behaviour which is then removed by the conversion from symmetric to asymmetric variables. Therefore, using a collinear approximation to the kernel is more appropriate if the kernel is the one in asymmetric variables. Thus, the RC parameters (c, κ and M min ) should be the ones of such kernel. While the curvature and the value of the kernel at the minimum are not sensitive to the conversion, the position of the minimum is, and thus in the 3.0 version of the HELL code we use the results of Ref. [38] with

B Running coupling corrections to the DGLAP-BFKL duality
The BFKL kernel (in symmetric variables) is not fully symmetric for the exchange M ↔ 1 − M because of a number of effects, all induced by the running of the strong coupling. The origin of these effects can be traced back to (see e.g. [29]): • the scale at which α s is computed in different pieces of the fixed-order kernel; • the conversion from the unintegrated (k t -dependent) PDFs, to which the original BFKL equation refers, to the integrated (collinear) PDFs, to which the DGLAP equation refers; • the so-called running coupling corrections to the duality, namely the correct derivation of the duality between DGLAP and BFKL taking into account the scale dependence of the strong coupling.
These effects have been extensively discussed in the literature and they will not be rederived; here we limit ourselves to present them. All these effects can be implemented as corrections to the BFKL kernel, which then gives by naive (i.e., fixed-coupling) duality the correct NLL contributions to the resummed anomalous dimension. These corrections amount to the following expression to be added to the symmetric kernel, 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, . 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 Sect. 3.