Dissecting the collinear structure of quark splitting at NNLL

We explore the collinear limit of final-state quark splittings at order $\alpha_s^2$. While at general NLL level, this limit is described simply by a product of leading-order $1\to 2$ DGLAP splitting functions, at the NNLL level we need to consider $1\to3$ splitting functions. Here, by performing suitable integrals of the triple-collinear splitting functions, we demonstrate how one may extract $\mathcal{B}^q_2(z)$, a differential version of the coefficient $B^q_2$ that enters the quark form factor at NNLL and governs the intensity of collinear radiation from a quark. The variable $z$ corresponds to the quark energy fraction after an initial $1 \to 2$ splitting, and our results yield effective higher-order splitting functions, which may be considered as a step towards the construction of NNLL parton showers. Further, while in the limit $z \to 1$ we recover the standard soft limit results involving the CMW coupling with scale $k_t$, the $z$ dependence we obtain also motivates the extension of the notion of a physical coupling beyond the soft limit.


Introduction
QCD has long been confirmed as the theory of strong interactions and is now decisively in the precision era. The demand for percent-level precision on QCD calculations is primarily being driven by LHC phenomenology with the large volume of experimental data accumulated at the LHC and the focus on firmly establishing the properties of the Standard Model Higgs sector as well as the search for possible subtle signs of new physics [1][2][3][4].
On the theoretical QCD side, in the context of precision, there has been particularly impressive recent progress in fixed-order perturbative calculations (see Ref. [5] for a review and further references). For a number of LHC applications however, one often needs to go beyond fixedorder calculations which are perturbative approximations involving a small number of partons, in contrast to the high-multiplicity hadronic final states seen in practice at colliders. The need for predictions beyond fixed-order becomes particularly obvious when one encounters observables sensitive to widely disparate scales and hence large logarithms in scale ratios, which require resummation.
Resummed calculations and the techniques behind them have also undergone substantial development over the past couple of decades. For a broad class of observables that have a property called recursive infrared and collinear (rIRC) safety [6] we can write resummed results for the cumulative distribution, i.e the distribution integrated up to some value v = e −L , in a standard exponentiated form [7]: Σ(α s , α s L) = exp α −1 s g 1 (α s L) + g 2 (α s L) + α s g 3 (α s L) + · · · , (1.1) where knowledge of the functions g 1 , g 2 and g 3 corresponds to achieving leading logarithmic (LL), next-to-leading logarithmic (NLL) and next-to-next-to leading (NNLL) logarithmic accuracy respectively. The state-of-the art has progressed from the NLL accuracy achievable for a limited number of observables in the early 1990's to achievement of automated NLL calculations for a wide range of global observables [6,8,9], and the development of NNLL (and in some cases essentially N 3 LL ) resummed calculations for the main classes of resummation , including progress on automated resummation frameworks [29,30,40,42]. However it remains true that the scope of analytic (or semi-numerical) resummation is still limited to a few types of observable which can be easily understood analytically in specific limits. On the other hand the range of observables that is relevant for phenomenology is ever increasing. Some notable examples come from the field of jet substructure in the boosted domain, where for the purposes of tagging and grooming it is common to simultaneously cut on a range of variables in order to achieve good discrimination between signal and background. In such cases, while resummation is essential to describe the results obtained, it is often hard to achieve analytically beyond the most basic modified LL accuracy (see Ref. [43] for a recent example from top tagging).
Given the limitations of current analytic resummation approaches, it is important to consider other available tools for handling general multiscale problems. This brings us to parton shower algorithms which are at the core of general purpose Monte Carlo event generators programs [44]. These, by their very nature, provide all-order perturbative results encapsulating the resummation of large logarithms for general observables. Parton showers have very broad applicability, but the question of understanding their logarithmic accuracy has proved elusive until recently. While showers are often stated to achieve leading logarithmic accuracy, and also include the main ingredients necessary for reaching NLL accuracy, Ref. [45] demonstrated that a widely used class of dipole showers, including the Pythia transverse-momentum ordered shower [46], which is the default shower in the Pythia8 program [47], failed to achieve general LL accuracy beyond the leading colour approximation and also failed to reach NLL accuracy for a number of common observables even at leading colour. Angular-ordered showers on the other hand have long been known to fail NLL accuracy criteria for non-global observables [48,49]. With a better understanding emerging of the limitations of various classes of showers, it has been possible to identify a set of principles that should be satisfied for a shower to be deemed NLL accurate [50]. NLL accurate showers, the PanScales showers, based on these principles have been constructed and numerically demonstrated to achieve full NLL accuracy for a wide range of global observables, including also subleading colour effects and spin correlations [50][51][52]. Other candidates for NLL accurate showers have also been proposed and analytically demonstrated, to achieve NLL accuracy for the thrust [53][54][55][56] and the jet multiplicity [53,54].
Given the recent spurt in activity in the context of NLL showers, which demonstrates that showers can be designed to achieve broad NLL accuracy, it is legitimate to think about whether further advances in resummation can be brought to bear on the construction of NNLL accurate showers. This is a substantially more challenging problem, but solutions to it would be invaluable in the context of extending theoretical precision for collider studies. A first obvious step in thinking about NNLL showers would be to consider the higher-order ingredients that shall be needed in order to extend shower accuracy. Once identified, the inclusion of these ingredients consistently with the existing shower framework would be required, which one can expect to be a highly non-trivial task.
In terms of the higher-order ingredients that may be necessary for extending shower accuracy, the criteria set out in Ref. [50] are useful to examine. In particular while NLL accurate configurations involve strong ordering between all pairs of emissions in at least one of two logarithmic variables (e.g. energy and angle), NNLL accuracy requires us to allow for a pair emissions that have comparable values for both logarithmic variables. Meeting this requirement entails the inclusion of higher-order splitting kernels that emerge in the double-soft and triple-collinear kinematic limits [57][58][59][60]. The inclusion of higher-order kernels in showers has been actively investigated in recent literature, both in the double-soft and triple-collinear limits [61][62][63][64] (see also related work in Refs. [65,66]). In the same context one should perhaps note that the inclusion of higher-order ingredients does not automatically improve the logarithmic accuracy of showers, which depends on the existing shower structure already being NLL accurate and its interplay with the higher-order ingredients. Given this, an investigation of the logarithmic accuracy of a shower following the inclusion of higher-order ingredients should be carried out, along similar lines to the detailed numerical tests presented in Ref. [50], before NNLL accuracy can reasonably be claimed.
In the current paper we shall carry out analytical studies, relevant to the derivation of specific higher-order shower ingredients, focussing purely on the triple-collinear limit. While the inclusion of the full triple-collinear splitting kernels [58,59] is eventually needed for general NNLL accuracy as mentioned above, we shall initially consider observables such as global event shapes, which are not directly sensitive to the splitting kernels themselves but receive an NNLL collinear contribution from an integral over the splitting kernels. 1 The coefficient that emerges from the integral of the triplecollinear splitting kernels is known in the resummation literature as B 2 [67][68][69][70]. This coefficient governs the intensity of hard-collinear radiation from an initial parton at order α 2 s . Its effective inclusion along with that of the resummation coefficient A 3 relevant to the soft limit [71,72], will be a key component of any future NNLL showers.
In this article we seek to better understand the connection between the B 2 coefficient and the triple-collinear splitting kernels by aiming to derive B q 2 , the result for a quark initiated jet, directly from the splitting kernels. In Ref. [73] we have recently derived the NNLL terms for the groomed jet mass distribution using the modified Mass Drop Tagger (mMDT or equivalently Soft Drop (β = 0) [74,75]), directly from the triple-collinear splitting functions. This work included reconstructing B q 2 as well as involving derivation of relevant NNLL jet clustering corrections. 2 In the current work we shall aim to extend and generalise the insight that emerged from the study of Ref. [73]. In particular, in order to facilitate eventual inclusion in parton showers, we shall derive here a version of B q 2 that is differential in the kinematics of a given emission and allows us to project to a definite phase-space point for a 1 → 2 branching. This in turn potentially allows us to include B q 2 by effectively giving an NLO weight to emissions that is correct in the hard-collinear limit and reproduces the standard B q 2 , familiar from resummation, upon integration over emission kinematic variables. As we shall demonstrate, examining higher-order collinear branchings in this manner also naturally produces a picture that points to an extension of the CMW coupling scale and scheme beyond the soft limit. Such an extension of the physical coupling ought to be of value for parton shower development but is also, we believe, of intrinsic theoretical interest. This paper is organised as follows: We start in section 2 with a brief reminder of the resummation coefficients, including B q 2 as defined and used in the literature. We also provide a reminder of the triple-collinear splitting functions and the phase-space we use for our calculations, as well as discussing the differential distributions we shall study. Section 3 is dedicated to our results. We start by briefly presenting leading-order results for the distributions in jet mass ρ or equivalently the angle θ 2 g , and energy fraction z for a given emission. In subsections 3.1, 3.2 and 3.3 we describe our order α 2 s results for the C F T R n f , C F (C F − C A /2) and pure C F C A terms respectively, all of which arise from the decay of an initially emitted massive parent gluon. The results we obtain are for the differential distributions in the invariant mass of the three parton system ρ and the energy fraction z of the parent gluon emission, as well as results for the distribution in θ 2 g and z, with θ 2 g being the angle of the parent gluon with respect to the final state quark. In subsection 3.4 we use our results to extract the contribution corresponding to B q 2 (z), i.e. a differential version of B q 2 , for each colour channel arising from the gluon decay subprocess. We discuss the structure of the results, including the relationship between the results for the ρ and θ 2 g distributions and explain how they are connected. We also obtain a direct connection between the z dependent functions we obtain from our differential distributions, and the NLO non-singlet timelike DGLAP splitting functions. In subsection 3.5 we combine our results with leading-order (order α s ) results, show how our results point towards an extension of the CMW coupling scale and scheme, and may be viewed as an effective extension of the concept of a web [77][78][79][80] beyond the soft limit. In subsection 3.6 we demonstrate how we may extract B q 2 (z) in the abelian C 2 F channel by noting that it ought to arise from the difference between a calculation performed with the triple-collinear splitting functions (including the one-loop corrections to a 1 → 2 splitting) and that performed using simply an iteration of successive 1 → 2 splittings. Finally in section 4 we draw our conclusions, and comment on prospects for future work. All other technical details relevant to our calculations can be found in the appendices.

B q , triple-collinear splitting functions and observable definitions
We begin with a brief reminder of the resummation literature in order to introduce the B q 2 coefficient, whose differential version we seek to compute in this article. This is the coefficient that controls the intensity of collinear radiation from a quark at order α 2 s , and hence is directly related to the quark Sudakov form factor. To exemplify this we report below the expression for the Sudakov form factor for the case of transverse momentum related resummation: where b is the impact parameter which is the Fourier conjugate of the transverse momentum or related variable. The A function encompasses soft emission effects, while the B function captures the effects of hard-collinear radiation. Each function admits a perturbative series The perturbative coefficients are then determined by matching the resummation formula onto a fixed-order computation for a given observable [76,81,82]. The form of B q 2 has long been known to have the following structure [42,76,81,83] where b 0 is the first coefficient of the QCD beta function and X v is a process and observable dependent constant. Finally, γ (2) q is the end-point contribution of the non-singlet next-to-leading order DGLAP kernel [84] Having provided a reminder of the B q 2 coefficient, we turn to the triple-collinear splitting functions that we shall use throughout this paper. The polarisation-averaged triple-collinear splitting functions we consider here were first derived in Refs. [58,59]. For an initial quark there are four distinct splitting functions to consider. In the notation of Ref. [59], these are: • Pq 1 q 2 q 3 corresponding to a 1 → 3 quark splitting involving non-identical quark flavours, associated with a C F T R colour factor.
• P (id) q 1 q 2 q 3 representing an additional interference contribution coming from identical quark flavours, with a C F (C F − C A /2) colour factor.
• P (nab) g 1 g 2 q 3 the non-abelian contribution leading to two final-state gluons and a quark, with a C F C A colour factor.
• P (ab) g 1 g 2 q 3 the abelian contribution also involving two final-state gluons and a quark but arising from independent gluon emission off a quark, with a C 2 F colour factor.
We shall work in terms of the energy fractions of the three collinear partons z i , satisfying i z i = 1, and shall label the angles between partons i and j as θ ij , such that in the collinear approximation, relevant to this work, θ ij 1. The splitting functions are reported in Ref. [59] as functions of z i and invariants s ij and s 123 where ij where E is the energy of the initial quark and s 123 = (p 1 + p 2 + p 3 ) 2 = s 12 + s 13 + s 23 . The triple-collinear phase-space in d = 4 − 2 dimensions may be expressed in the form [85] where the Gram determinant ∆ is defined as Additionally, in Ref. [73] we used a set of variables, which can be referred to as "web variables", to parametrize the triple-collinear phase space. The web variables will again prove essential to obtain analytic expressions for the distributions we are interested in. In appendix A.1 we recapitulate the phase space and recollect the physical meaning of the different variables. In Figures 1 -3 we illustrate the splitting sub-processes involved here, namely a gluon decay contribution showing a q → qg splitting followed by gluon splitting to a qq pair ( Figure 1) and a gg pair ( Figure 2) as well as the abelian contribution to q → qgg with independent gluon emission off a quark ( Figure 3). The splitting process shown in Figure 1 gives both the identical and non-identical fermion splitting functions. Figures 1 -3 also illustrate our change of independent variables from z 2 , z 3 to z and z p which are the variables associated to successive splittings, and are defined so that in the limit of strong angular-ordering between successive branchings the triple-collinear splitting Figure 2. The Feynman diagram representing the gluon decay C F C A channel. The mapping of the momentum fractions to a set of independent variables is explained in the text. The angle θ g is that between the parent gluon and the final-state quark. functions reduce (after azimuthal integration for the gluon decay channels) to a product of leadingorder splitting functions in z and z p respectively.
Finally we discuss the quantities we study here. These include the double differential distribution in ρ = s 123 /E 2 and z, where ρ is the normalised invariant mass of the three parton system that arises from the triple-collinear splitting of a quark jet with energy E, and the splitting variable z may be associated to an initial splitting as illustrated in Figures 1 -3. 3 For the gluon decay contributions, as should be evident from Figures 1 and 2, the variable z also corresponds to the energy fraction of the final-state quark so that 1 − z represents the energy fraction associated to the "parent" gluon. In addition to the ρ distribution we shall also study the distribution differential in z and angle θ g of the parent gluon. A comparison of the two distributions shall give further insight into the general structure of the result. For the abelian gluon emission process in Figure 3, we shall fix θ 13 = θ 1, the angle of emission 1 wrt the final quark, which shall set the collinearity, and then study the NLO structure induced by a smaller angle emission labelled 2, with angle θ 23 < θ.
We integrate the splitting functions over phase-space in d = 4 − 2 dimensions to obtain real emission contributions that contain poles in which reflect singularities that cancel when we combine with virtual corrections. The integrals we carry out are generically of the form where v denotes the quantity we hold fixed, P denotes the different 1 → 3 splitting functions mentioned above, α denotes the bare QCD coupling, and σ 0 is the Born cross section. Our results will be expressed in terms of a renormalised MS coupling α s , given by the relation where we have the standard MS factor and we choose µ R = E (the energy of the hard initiating parton). For the results that follow we define α s ≡ α s (E 2 ) in the MS scheme.

Results
We note that the leading-order collinear limit result for the ρ distribution is given by considering a single collinear gluon emission from a quark, with energy fraction z and angle θ g so that ρ = z(1 − z)θ 2 g . This result is equivalent to that for fixed θ g and z so that we have At order α 2 s we start by considering C F T R n f terms both for fixed invariant-mass of the three-parton system ρ and z, as well as for fixed θ g and z where θ g is the angle of the parent gluon with respect to the final quark. For the former case we have already obtained a previous fully analytic result as part of our study of the modified Mass Drop Tagger's jet mass distribution. This was reported in Ref. [73] and we shall analyse its structure in more detail here, together with that of the θ 2 g distribution calculated here.

Fixed invariant-mass
The order α 2 s result for the ρ distribution, in the collinear limit, was calculated in Ref. [73]. Two separate calculations were performed there. Firstly we carried out a part analytic and part numerical computation where the divergent terms were computed analytically and a finite remainder computed through numerical integration in four dimensions. Secondly we computed a fully analytical result using a parametrisation of the phase space based on web variables. The two calculations were found to be in perfect agreement and we report below the fully analytical result [73]:

Fixed parent angle
We now study the distribution in the angle of the parent gluon and z. The angle of the parent can be straightforwardly expressed in terms of our phase-space variables and is given, in a collinear approximation, by θ 2 2 12 with variables as illustrated in Figure 1. We note that in the limits z p → 0 or z p → 1, where all the energy is carried by a given offspring parton, θ 2 g reduces to the angle of the energetic offspring parton (1 or 2) wrt parton 3 i.e. the quark. In the collinear limit θ 12 → 0, θ 13 → θ 23 , θ 2 g reduces to the angle between the direction of either collinear parton and the quark.
We follow the same strategy as the one adopted in Ref. [73] for the ρ distribution and perform the calculation using two different methods: firstly a partly analytical calculation with a finite remainder evaluated numerically and secondly a fully analytical calculation using web variables.
In the partly analytic approach, integrating the triple-collinear splitting function Pq 1 q 2 q 3 at fixed z and θ g we obtain a result that has only a 1 pole which originates in the collinear divergence as the angle between gluon offspring partons θ 12 → 0. The collinear pole is multiplied as usual by an dependent factor and on performing an expansion we obtain a pure 1 pole term and associated finite terms.
Further, on subtracting the leading collinear divergent piece from the full integrand, we also obtain a finite term in addition to the collinear divergent term, that can be computed by integration in four dimensions, i.e. setting → 0. This finite term has two components: a term that behaves as 1/(1 − z) i.e. is singular in the limit of a soft parent and we can extract analytically, and a term that is regular as z → 1 which we can leave to numerical integration.
In this approach, the analytically determined component of our answer, which contains the collinear divergence, can be expressed as 3) The additional finite term that can be evaluated in four dimensions turns out to be identical to the corresponding result for the ρ distribution. This term takes the form below, confirmed by both our numerical evaluation and our separate fully analytic evaluation (see appendix A.2) Finally we need the one-loop virtual correction to the 1 → 2 q → qg splitting which reads: (3.5) Combining all terms we get A notable feature of our result is that it can be obtained from that for the ρ distribution via the replacement ρ → z(1 − z)θ 2 g which is an exact relationship between the two observables at leading order i.e. for the emission of an on-shell gluon. The reason for this is simply the fact that the observable dependence of our result comes solely from the divergent limit when θ 12 → 0, where IRC safe observables tend to their leading-order form. The finite terms that are computable in four dimensions are instead observable independent. More specifically, for some fixed value of observable v, once divergent contributions are removed we are left with genuine triple-collinear configurations with three energetic partons with comparable opening angles of order θ 2 ∼ v. The overall 1/θ 2 scaling of the triple-collinear splitting functions then implies that vdσ/σ 0 dv is independent of v.
Next we examine the colour suppressed identical fermion term also arising from gluon splitting to qq. This is a purely finite term and can be calculated in four dimensions. In Ref. [73] we numerically computed the contribution of this term to ρ dσ dρ as part of our calculations for the mMDT jet mass. However, in Ref. [73], we did not determine the z dependence of our result as we only required the integral over z. After the computational steps explained in appendix A.5, we provide a fully analytical result for the ρ distribution also differential in z. We also carry out a numerical calculation for fixed θ 2 g where we continue to label as θ g the angle θ 1+2,3 , the angle between the direction of the 1, 2 parton pair and the quark labelled 3. However we note that it is no longer possible to interpret this angle unambiguously as the angle of the parent gluon due to the identical quarks present in the final state. As expected, for pure finite terms, the results for the ρ and θ 2 g distributions coincide so that we obtain In deriving P (id.) (z) we have defined z 3 = z as usual but the same result is obtained for the distribution in z 2 , as one may readily anticipate for identical particles.
In Ref. [73] we provided a calculation for the C F C A contribution for the ρ distribution, that arises from the decay to gg of a parent gluon emitted off a quark. Our results were obtained in a form which was partly analytic, while finite terms computable in four dimensions were obtained through numerical integration. Here we report a fully analytic result with the details of the calculation left to appendix A.4. We obtain The results for the θ 2 g distribution give a corresponding function P nab. (z; θ 2 g ). As observed for the C F T R n f term the observable dependence in ρ/σ 0 dσ/dρ arises from the collinear divergent limit so that one may obtain the θ 2 g result via a simple substitution for ρ. Hence we have 3.4 B q 2 (z) and comments on general structure We now study the link between our results and B q 2 , the parameter in the quark form factor that relates to the intensity of collinear radiation from a quark at O(α 2 s ). We discuss also the connection between our results and the NLO timelike splitting functions [84]. Finally we show how to combine our results with the corresponding leading-order results, recovering the expected behaviour in the soft limit and giving a simple picture for z dependent corrections beyond the soft limit.

Extracting B q 2
Let us start with the θ 2 g distribution, whose C F T R n f term is reported in Eq. (3.6). In the soft limit, z → 1 the result reduces to a familiar one: (3.12) where we note the presence of a soft divergence as z → 1, reflecting the singular behaviour of the q → qg splitting function. In writing the second line we have introduced the transverse momentum of the emitted parent gluon k 2 t = E 2 (1 − z) 2 θ 2 g wrt the final quark direction. We have also defined b (n f ) 0 , which is the n f part of the first perturbative coefficient of the QCD beta function (see Eq. (2.3)) and K (n f ) is the n f term in the CMW constant [86] 14) The terms appearing in Eq. (3.12) are essentially related to NLL rather than NNLL structure, leading to the well known soft-limit prescriptions for the scale and physical scheme for the strong coupling. Indeed in terms of the type of decomposition into soft and collinear pieces presented in Eq. (2.1) these terms are all associated to the "A" series of coefficients owing to the presence of the soft divergence. Here we wish to obtain the pure collinear NNLL structure and hence we shall subtract off these terms. We shall also remove the remaining piece of the ln θ 2 g term ∝ −(1 + z) ln θ 2 g as this term is also an NLL contribution. Indeed, the −(1 + z) ln θ 2 g yields after z integration −3/2 ln θ 2 g which makes clear its association with the B 1 hard-collinear coefficient, rather than B 2 . 4 Hence we obtain where the subtracted terms are placed in square brackets and leave an NNLL pure collinear result denoted by B q,n f 2 (z; θ 2 g ), which is a differential version of the n f term of B q 2 . A similar result can of course be written for the ρ distribution by again removing soft enhanced and NLL in ρ terms . Indeed, as we have noted before, the result for the ρ distribution can be reached from that for the θ 2 g distribution via the one-gluon relation between observables i.e. via the replacement θ 2 g → ρ/(z(1 − z). This correspondence ultimately results in the relationship: We can now perform the integrals over z to write the results in a standard form: where γ (2,n f ) q is the n f part of the endpoint contribution to the DGLAP splitting kernels given in Eq. (2.4) and we determine the observable dependent constants (3.20) Notice that we inserted a factor of (2π/α s ) 2 so that our extraction of B 2 agrees with the standard definition, i.e. Eq. (2.3).. The above results are in line with the expected form of the B 2 coefficients (see Eq. (2.3)) and hence consistent with previous observations in the literature [81] that B 2 is always related to the endpoint terms of the DGLAP splitting kernels and additional b 0 X terms where X depends on the observable. Our result in Eq. (3.18) for the ρ distribution is in agreement with the collinear NNLL result in the literature for the plain jet mass [42,87] as we also pointed out in Ref. [73]. The result for the θ 2 g case has been computed here for the first time, to our knowledge. For the identical particle contribution there is no soft contribution to be removed nor any observable dependent b 0 X terms. The result derived in Eq. (3.7) corresponds to the pure NNLL constant and hence we have The above analytical result is consistent with our previous numerical evaluation in the context of the Soft Drop jet mass [73]. 5 Next we turn to the non-abelian contribution to q → qgg, involving the decay of the parent gluon to a pair of gluons, g → gg. An exactly analogous procedure can be followed as for the g → qq decay to remove soft and higher logarithmic order (i.e. NLL) contributions. In particular, in the soft limit we obtain the result given in Eq. (3.13) with b (n f ) 0 and K (n f ) replaced by the corresponding C A terms. We thus derive the results below: (1 + z) ln(1 − z).

(3.24)
We note that the sum of dilogarithms multiplying the LO splitting function p qq (z), on the second line of Eq. (3.23), is regular as z → 1 and has a small z limit given by − 1 2 ln 2 z + π 2 6 + O(z).
We then have the following integrated results: As we would expect the difference between the results for ρ and θ 2 g distributions comes only from the C F β 0 X observable dependent terms with values of X already identified in Eq. (3.18). Thus we have Note however that the integrated results given in Eq. (3.25) do not yet fully correspond to the standard form given in Eq. (2.3). In particular in order to fully recover the C F C A term of the DGLAP endpoint contribution, γ (2,C A ) q , we need to combine the C F C A results from Eq. (3.25) with the identical particle contribution computed in Eq. (3.22) which, by virtue of its C F (C F − C A /2) colour factor, contributes to both C 2 F and C F C A colour channels. Defining B q,(id.),C A 2 as the C F C A piece of the identical particle term in Eq. (3.22), and taking the example of the ρ distribution, we have 6 with γ (2,C A ) q being the standard C F C A piece of the DGLAP endpoint contribution.

Relationship to NLO timelike nonsinglet DGLAP splitting kernels
Having extracted the z dependent B q 2 pieces in the previous section, it is also of interest to consider the relationship between our results and the NLO timelike DGLAP splitting kernels themselves. 7 In order to establish a formal connection with the full structure of the mass-singularity factorisation formula in QCD, we would need to integrate our results over ρ or θ 2 g and examine the resulting pole structure to recover the NLO splitting kernels. Here it is not our intention to pursue the connection from this viewpoint but rather to study the functional dependence of our results on z, and its possible link to the NLO splitting kernels.
Returning to Eq. (3.6) and focussing on its z dependence, i.e. setting θ 2 g = 1 so as to remove the pure ln θ 2 g term, we obtain a function that may be expressed as (suppressing the overall αs 2π 2 6 Identical considerations of course hold for the θ 2 g distribution with Xρ replaced by X θ 2 g . 7 We note here that for the splitting kernels we intend to study in this subsection, the timelike and spacelike NLO splitting kernels have the same functional form.

factor)
Written in this form, the top line of Eq. (3.29) corresponds to the C F T R n f piece of the nonsinglet time-like splitting function P V (1) qq (z) [84]. When integrated with a + prescription on the 1/(1 − z) factor (or equivalently after removal of terms that diverge as z → 1) it gives the −γ The difference between the b 0 X type terms for θ 2 g and ρ arise, as we stressed before, purely from the factors of z and 1 − z in the leading-order relation ρ = z(1 − z)θ 2 g . Next, let us examine the identical fermion contribution. Here there are no b 0 X terms so we might expect a simple relationship of our result to the corresponding NLO timelike splitting function P V (1) qq (x). That is indeed the case, however the relevant function is not the P (id.) (z) obtained in Eq. (3.7). The latter is a function of z, which is the energy fraction z 3 of parton 3 which is one of the two identical final state quarks. It is therefore, as we shall verify below, a contribution to the splitting function P V (1) qq (x). If instead we carry out the calculation with fixed energy-fraction for the antiquark, i.e. set z 1 = x, we obtain a different function which we expect to coincide with P V (1) qq (x). To verify this we have simply performed the calculation numerically for several x values and checked that the result agrees precisely with the form [84,88] The above function, on restoring the qq (x). Moreover, the integral over x of P (id.) (x) also gives, as expected, after supplying the colour and coupling factors, B q,(id.) 2 as given in Eq. (3.22).
Finally we come to the pure C F C A non-abelian term from the q → qgg splitting. The NLO collinear limit result for the ρ distribution is explicitly written in Eq. (3.8). Following a similar strategy to that for the n f piece involves setting ρ = 1 and separating the result into two pieces, one of which yields the b (C A ) 0 X ρ term and the other that we may expect to be linked to the splitting kernels. It is simple to identify the terms that lead to b (C A ) 0 X ρ . These are the same as the corresponding terms for the n f piece written on the second line of Eq. (3.30) but with the replacement of 2 3 → − 11 6 . The equation analogous to (3.30) may be written as P NLO,nab. (z; ρ) = P NLO,nab. sub.
and where the suffix "sub." is a reminder that terms contributing to b 0 X have been subtracted off to obtain this function. We can indeed use this result for P NLO,nab. sub.
(z; ρ) to obtain the C F C A piece of the splitting kernel P V (1) qq (z). In order to do so we note that P V (1) qq (z) represents the splitting function where one obtains a final quark with momentum fraction z. In terms of all ways in which this can happen we should also consider the C F C A contribution from the C F (C F − C A /2) identical particle term which yields two identical quarks from the splitting of an initial quark. The momentum distribution of each of these final quarks is the same and is given by the function P (id.) (z) reported in Eq. (3.7). Thus we are led to the combination P NLO,nab. sub.
where the result on the RHS does not involve any dilogarithmic terms 8 and is in agreement with the C F C A term of P V (1) qq (z).

Combination with leading-order results
In this section we consider the combination of our O α 2 s results for the C F T R n f and pure C F C A non-abelian terms (i.e. those arising from gluon branching to a pair of gluons) with the leadingorder result Eq. (3.1), which will lead to the recovery of the running coupling in the physical scheme [42,86,89], in the soft limit. Beyond the soft limit the picture that emerges remains simple to interpret and motivates a z dependent extension of the soft limit result. Taking the θ 2 g distribution as an explicit example we combine the order α s and order α 2 s results Eqs. (3.1), (3.6), (3.10), (3.11) to define: (3.35) 8 The following identities are useful to simplify the expression: Li2(1 − z) = π 2 6 − ln z ln(1 − z) − Li2(z) and Li2 z−1 z + Li2(1 − z) = − 1 2 ln 2 z, with 0 < z < 1.
which gives (recall that we have fixed µ R = E, the energy of the initial quark) We have written the result above in a form which helps to emphasise some of its main features. Firstly the two terms on the top line of Eq. (3.36) combine to produce α CMW where we have also implicitly included terms of order α 3 s and beyond via our change of scale for the coupling.
This agrees with the expected result in the soft limit, z → 1, since the scale of the coupling is E 2 (1 − z) 2 θ 2 g i.e. the familiar transverse momentum squared of the soft emission wrt the emitting parton and we recover the CMW prescription. On the second line we note that there are two terms proportional to the beta function coefficient b 0 . The presence of a b 0 ln z term is suggestive of a redefinition of the scale of the coupling in the hard collinear region to the scale E 2 z(1 − z) 2 θ 2 g . On the other hand the b 0 (1 − z) term is universal and originates purely in the virtual corrections, implying that it can also be absorbed into the definition of the coupling via a z dependent extension of the CMW scheme. Such a redefinition of the coupling scale and scheme would imply that the QED-like n f piece of the result is fully incorporated into the definition of the coupling, consistent with suggestions made in the past literature [90,91]. One may also expect the step of defining a collinear-improved coupling to be useful from the viewpoint of inclusion of higher order (NNLL) ingredients, related to triple-collinear splitting kernels, in parton shower algorithms. 9 We shall postpone detailed discussions and definite proposals along these lines to forthcoming work [92].
The function R nab. (z) reads which is the remainder of the cross section solely with a C F C A term. The above function is not soft enhanced, as the combination of terms in parenthesis multiplying p qq (z), vanishes as 1 − z in the z → 1 limit. The R nab. (z) term can be viewed as a higher-order splitting function of pure collinear origin. We close this subsection by reminding the reader of a key property of the result Eq. (3.36). Although Eq. (3.36) represents the distribution in z and θ 2 g , it is straightforward to obtain the distribution in the mass ρ by integrating Eq. (3.36) over θ 2 g with the constraint δ(ρ − z(1 − z)θ 2 g ) exactly as one would do in a leading-order calculation. The same holds for the distribution in any observable v defined in terms of the parent gluon kinematical variables (z, θ 2 g and mass m 2 ). The fact that one can obtain one parent gluon kinematical distribution from the other at order α 2 s , by using a relationship valid in the limit of a massless gluon, implies that the effect of the gluon virtuality has effectively been absorbed into the structure of Eq. (3.36). This is reminiscent of the fact that in the soft limit and to NLL accuracy for global observables one can replace the emission of a massive gluon by a massless gluon, with the effect of the gluon branching included in the argument of the coupling and the CMW factor K. Therefore we may think of Eq. (3.36) as an extension of the web concept beyond the soft limit and into the hard-collinear region, via an extension of the CMW coupling and a higher-order splitting function.
3.6 C 2 F abelian piece Next we examine the pure C 2 F piece, i.e. the abelian q → qgg contribution. This channel is somewhat different from the gluon decay contributions we have studied thus far. Here there is no association of the results to the running coupling scale or scheme and no b 0 X terms that accompany the DGLAP endpoint contribution. Accordingly in this channel one can focus purely on illustrating the extraction of B q,(ab.) 2 (z), the abelian gluon emission piece. In Ref. [73] we arrived at B q 2 as part of the NNLL structure 10 for the mMDT jet mass distribution ρ σ 0 dσ dρ , by performing an order α 2 s calculation in the small ρ limit. For two real emissions, this involved considering two configurations. The first, called F pass (z, ρ) in Ref. [73], had a relatively energetic large-angle emission that passes a condition 1 − z cut > z > z cut , with a smaller-angle emission with no constraint on its energy. Together the emissions set a value ρ for the normalised jet mass squared, and taking ρ z cut ensures that both emissions are collinear though not necessarily strongly ordered in angle. This configuration falls entirely within the jurisdiction of the triple-collinear splitting functions. The second configuration, called F fail (z p , ρ), involves a situation where the larger-angle emission is relatively soft and fails the z cut condition. 11 This emission is then "groomed away" and does not contribute to the jet mass ρ which is then set by the smaller angle emission. The emission that is groomed away can be at any angle and its treatment requires going beyond collinear calculations. Only the F pass (z, ρ) term is directly related to B q,(ab.) 2 , while the F fail (z p , ρ) is a necessary part of recovering the full mMDT jet mass result including its leading-logarithmic terms.
Similarly for virtual corrections we needed to consider two distinct terms: firstly there is the standard one-loop correction to the Born process (qq production) and secondly there is the one-loop correction to a 1 → 2 collinear splitting. Again, it is only the latter contribution that is relevant to B q,(ab.) 2 , while the former is needed as part of obtaining the full NNLL result for the mMDT ρ distribution. Performing the order α 2 s calculation and subtracting the NLL result, which arises entirely from the approximation of emissions strongly-ordered in angle, we obtained the relevant contribution to B q 2 , alongside NNLL "clustering" corrections. Upon combining with the C 2 F term 10 We remind the reader that leading double logarithmic contributions are absent for the mMDT so that the logarithmic hierarchy starts with NLL single-logarithmic terms. 11 Note that F fail (zp, ρ) is a function of the splitting variable zp, shown in Fig. 3 from the C F (C F − C A /2) channel we recovered, though we did not study its z dependence.
In this article we demonstrate more directly, by studying a simpler example, how the B q,(ab.) 2 (z) contribution may be extracted as a difference between the triple-collinear and strongly-ordered in angle regimes, which more clearly exposes its physical origin. A key point of difference from the correlated emission calculations we have presented thus far is that when considering the decays of a parent gluon emission, requiring the parent emission to be collinear to the emitter is sufficient to constrain all three final partons to be collinear. 12 However in the C 2 F channel requiring, for instance, a small jet mass does not necessarily fix all partons to be collinear, in particular, allowing for soft wide-angle emissions unrelated to collinear physics, as already alluded to above when discussing the F fail (z, ρ) contribution to the mMDT jet mass.
We shall therefore consider a collinear gluon emission (labeled 1) with a fixed z 1 = 1 − z and angle θ = θ 13 1 wrt the final emitted quark (see Figure 3). This emission will set the collinearity in the same way as the parent emission for the correlated emission channels. We shall then examine the role of more collinear real emissions by integrating over a second emission such that θ 23 < θ, but without imposing any other constraints. Integrating over the smaller angle emission produces divergences and accompanying finite corrections. We will then combine this calculation with the one-loop corrections to a 1 → 2 collinear splitting. It is important to note here that the above mentioned angular restriction on the secondary emission is the sole reason prohibiting an analytic calculation. Essentially, the azimuthal integral is rendered incomplete due to the angular restriction, and thus we could not provide an analytic result. Finally, we will perform the exact same calculation but using strongly-ordered dynamics with factorised 1 → 2 splitting functions and phase space. We anticipate that the difference between the full triple-collinear limit calculation (including the oneloop correction to a 1 → 2 splitting) and the same calculation using the iterated 1 → 2 kernel and phase space will suffice to directly yield B q,(ab.) 2 (z). Integrating the triple-collinear splitting function P (ab) g 1 ,g 2 ,q 3 over the region θ 23 < θ 13 = θ we obtain the following result in 4 − 2 dimensions: where the label d-r denotes the double real contribution, p qq (z, ) = 1+z 2 1−z − (1 − z), the different pole terms are labelled according to the origin of the divergence i.e. whether it comes from when the second emission is soft and collinear, pure collinear or pure soft. The contribution H fin. (z) is a finite correction that we obtain via numerical integration in four dimensions.

40) with
It is to be noted that the above result is precisely the same as the F pass (z, ρ) term we obtained for the case of fixed jet mass ρ in Ref. [73] with the substitution ρ → z(1 − z)θ 2 . This relationship is just the single emission relation between ρ and θ and hence exact in the divergent limits i.e. when the second emission is vanishingly soft (z 2 → 0) or exactly collinear (θ 23 → 0) to the final quark. The finite contribution H fin. (z) is identical to the mMDT jet mass case. This structure for the answer again reflects the point made previously (for the gluon decay channels), that once divergent contributions have been removed from the triple-collinear splitting functions, the remaining finite contribution originating from relatively energetic partons, with commensurate emission angles, is independent of the observable. Accordingly one should also expect that to obtain the analogous result, from the region θ 23 < θ 13 , for other IRC safe observables v, which constrain all 3 emissions, one can use Eqs. (3.40), (3.41) and simply replace θ by its one gluon form in terms of the observable v and z.
In Ref. [73] we evaluated the integral of H fin. (z), 1 0 H fin. (z)dz ≈ 0.933 · · · , which is consistent with the analytical form 4ζ(3) − 31 8 . 13 The z dependence of this function, which we have only been able to obtain numerically is shown in Fig. 4. It shows that as z → 0 the result tends to a constant, while the steep z → 1 behaviour appears consistent with the presence of a ln z ln(1−z) We believe, on the basis of carrying out some analytical investigation, that indeed the ln z ln(1−z) 1−z behaviour is present in the answer, but nevertheless, there might be additional terms, e.g. pure ln(1 − z), that also contribute to the diverging behaviour as z → 1.
In order to obtain B q,(ab.) 2 (z) we will also eventually need to add to our calculation in Eqs. (3.40), 13 In Ref. [73] we noted that in order to recover the analytic result for B q 2 we needed to identify the numerically computed value with the aforementioned analytical form. With the help of a precise numerical evaluation and the PSLQ algorithm [93] it has been possible to analytically reconstruct this answer. We thank Pier Monni for his work to verify this. (3.41) the one-real, one-virtual contribution involving the one-loop correction to the collinear 1 → 2 splitting given by [94] θ 2 σ 0 We now compute the strongly-ordered version of the double-real emission term, which we shall use as a subtraction term to generate B q,(ab.) 2 (z). This is given simply by considering a factorised product of splitting functions and applying a factor of the single collinear emission phase-space for each emission in 4 − 2 dimensions [95]: (3. 43) in writing which we have expressed the result in terms of the renormalised MS coupling (with µ R = E). We observe immediately that Eq. (3.43) gives an identical structure to that of the pole functions in the full triple-collinear result, Eqs. (3.40) (3.41), except that the prefactor multiplying the double and single poles involves a factor z −2 rather than z −4 and there is no H fin. (z) term. Therefore on subtracting the strongly-ordered term from the full result we obtain: . (3.46) The functional dependence of B q,(ab.) 2 (z) on z is shown in Figure 5. In addition to the steep behaviour as z → 1, also seen for H fin. (z) we also note the presence of a small z divergence which comes from the small z approximation of Li 2 z−1 z ≈ − 1 2 ln 2 z − π 2 6 . Performing the integral over z of B q,(ab.) 2 (z), and combining with the C 2 F term of the C F (C F − C A /2) identical particle contribution in Eq. (3.22), we recover the expected result: is the C 2 F term of the identical particle contribution to B 2 and we have used the assumed analytic form, rather than the numerical value, for the integral of H fin. (z).
A few comments on the result in the abelian gluon emission channel are in order. As stated before, given the absence of any relationship to the definition of α s , the result obtained for B q,(ab.) 2 (z) can be entirely viewed as an effective higher-order splitting function. While in the gluon decay channels we could relate the z dependent functions we obtained to the timelike NLO non-singlet splitting kernels, our result in the present channel cannot be directly related to the C 2 F term of P V (1) qq (z). Indeed, in our current work, for all channels we have obtained the result differential in the kinematics of the first emission, i.e the emission setting the collinearity. For the abelian gluon emission case our z variable thus refers to the intermediate quark resulting from that first emission, rather than the final quark after all splittings which defines the argument of P V (1) qq (z). Moreover we have worked with a fixed first emission and imposed a constraint on the second emission to be more collinear, which crucially impacts the structure of the result we obtain and potentially modifies the calculation relative to the one needed for P V (1) qq (z).

Conclusions
In this paper we have studied collinear parton splittings at order α 2 s and NNLL accuracy, for splittings of a final-state quark. At our accuracy, i.e. to study NNLL collinear terms, the relevant limits of the QCD matrix-elements are given by the well-known triple-collinear splitting kernels. We have performed calculations for all relevant channels at order α 2 s , involving both the decay of a massive gluon emitted from the initiating quark (i.e. correlated emission terms), as well as the purely abelian term with successive gluon emissions off the initiating quark.
For the gluon decay contributions, involving C F T R n f , C F C A and C F (C F − C A /2) colour channels, we have computed the distribution differential in ρ and z, which are respectively the normalised invariant mass squared of the triple-collinear parton system and the energy fraction of the quark after emission of the massive gluon. We have also studied the distribution in θ 2 g and z, where θ 2 g refers to the angle made by the parent massive gluon with the final quark. In the abelian C 2 F channel we have again fixed the energy fraction and angle of an initial gluon emission and then examined the role of a more collinear real emission together with the one-loop virtual correction to the initial collinear splitting. The distributions we have derived give us several pieces of information that ought to be valuable when considering ingredients needed for an NNLL parton shower, and also of general theoretical interest. Firstly, for the gluon decay channels, in the soft limit z → 1, we noted that our results produced terms that reconstruct the correct scale of the coupling in the soft limit as well as giving the CMW coefficient K [86]. Further, after removal of the soft and higher logarithmic order terms we obtain a pure hard-collinear term B q 2 (z), a differential in z version of the coefficient B q 2 that governs the intensity of collinear radiation at order α 2 s . Upon integration of B q 2 (z) over z, we recovered the standard form for B q 2 including the δ(1 − z) endpoint contributions of the timelike NLO non-singlet splitting kernels accompanied by terms of the form b 0 X, with b 0 the first perturbative coefficient of the QCD β function. The value of X is dependent on the quantity whose distribution we are considering. We were additionally able to relate the z dependent functions that emerge from our calculations to the timelike NLO non-singlet splitting kernels, after removal of terms that give rise to the b 0 X contributions.
We also discussed how our results suggest a possible extension of the scale and scheme for the soft limit CMW coupling. In particular, from our results for the θ 2 g distribution, we identified a b 0 ln z term that, when combining with the leading-order result, can be absorbed in a change of scale of the coupling, thereby modifying it relative to the known soft limit result. A further term proportional to b 0 (1 − z) can be considered as an extension of the CMW scheme beyond the soft limit, and doing so, alongside modifying the scale as just mentioned, completely absorbs the n f dependent pieces into the definition of the coupling. Terms that are left over are devoid of an n f term and may be viewed as effective higher-order splitting functions. We shall leave a more detailed justification and concrete proposals for the coupling scale and scheme to forthcoming work [92].
For the abelian C 2 F channel we demonstrated how B q 2 (z) can be extracted by thinking of it as an NNLL correction that emerges naturally from the difference between the triple-collinear limit calculation (including also the one-loop correction to a 1 → 2 splitting) and the same calculation performed using NLL dynamics, i.e. with factorised splitting functions and phase-space. In the case of the abelian C 2 F piece our results gave rise to a function which we could compute only partially analytically, while we obtained the remaining part of the result through a numerical evaluation. Again, after integrating over z we recovered the known result for B q 2 , also in this channel. In conclusion, as stated before, we anticipate that the findings of this paper will be of value in incorporating the NNLL hard collinear corrections, corresponding to inclusion of B q 2 , in future higher accuracy parton showers. Prior to doing so it would be necessary to extend these results to take into account collinear splittings for gluon jets and similarly account for splittings of initial state partons, which we also leave to forthcoming work. Exploring how best to match the calculations performed here to existing shower frameworks, such as the PanScales showers [50][51][52], would also still require significant further thought and effort. Nevertheless we believe that the calculations performed here and the insights we have obtained are sufficiently general so as to constitute a key step relevant to any such future work. and the solid angle is that of q ⊥ , aligning k ⊥ along a reference axis in the transverse plane. We notice the nice feature of eq. (A.1) that the double-soft limit, z → 1, immediately recovers the double-soft phase space [96].
A.2 C F T R n f : θ 2 g distribution As we explained in the text, the fixed invariant-mass distribution has been presented in ref. [73]. Here, we give the details of the similar computation that yields Eq. (3.6). The triple-collinear splitting function reads Therefore, after carrying out the azimuthal average we find Notice that the integral over the invariant mass, s 12 is convergent from above and thus can be safely extanded to infinity. In the presence of an observable, such as the jet mass ρ, the upper limit would be dictated by the observable, e.g. max.{s 12 } = E 2 (1 − z)ρ. Performing the remaining integrals, and imposing charge renormalization, yields θ 13 = 0 or θ 23 = 0, as one can easily verify. This collinear structure can be understood, from physics standpoint, as a consequence of colour coherence (angular ordering). Now to the computation. The first thing to realize is that the total invariant mass is free from the azimuthal angle of the web phase space, in particular, and, therefore, the azimuthal average can be performed at the onset. The invariants s 13 and s 23 are given in terms of the squared transverse vectors, q 2 ⊥1 and q 2 ⊥2 , namely The full computation is somewhat tedious, so we will not proceed to give all details. Yet, we explain the necessary steps for the interested reader. First, the expansion of hypergeometric function can be performed on the integrand level and reads [97] A.5 C F (C F − C A /2): ρ distribution Let us recall the expression of the splitting function: As mentioned in the text, this colour channel contains neither soft or collinear poles, and thus one can set = 0 at the onset and perform the computation in 4 dimensions. The azimuthal average is performed using Eq. (A.8), setting = 0, followed by the s 12 and z p integrals. The final result is given in Eq. (3.7).