Double collinear splitting amplitudes at next-to-leading order

We compute the next-to-leading order (NLO) QCD corrections to the $1 \to 2$ splitting amplitudes in different dimensional regularization (DREG) schemes. Besides recovering previously known results, we explore new DREG schemes and analyze their consistency by comparing the divergent structure with the expected behavior predicted by Catani's formula. Through the introduction of scalar-gluons, we show the relation among splittings matrices computed using different schemes. Also, we extended this analysis to cover the double collinear limit of scattering amplitudes in the context of QCD+QED.


I. INTRODUCTION
Reaching higher-orders in the context of perturbative QCD implies a great challenge, but it becomes crucial to achieve the level of accuracy required by nowadays experiments. LHC results need to be compared to high precision theoretical predictions to unveil novel high-energy physics phenomena. In order to achieve that goal, it is necessary to understand the infrared (IR) divergent structure of QCD amplitudes. Within the framework of dimensional regularization [1,2], a lot of work has been performed at one-loop, two-loop and higher-order loops [3][4][5][6][7][8][9][10][11]. Many methods rely on the properties of the collinear/soft limit to perform the analytical subtraction of IR divergences, which allows to obtain finite cross-sections at colliders (for instance, see Ref. [12]).
Centering in the collinear limit, it is known that scattering amplitudes and squared amplitudes take a rather simple form in this class of kinematical configurations. Moreover, there are proven factorization properties 1 which show that IR collinear divergences exhibit a universal processindependent behavior, although strict collinear factorization is violated in the space-like region [14][15][16].
At the squared amplitude level, this universal behaviour is captured by the Altarelli-Parisi (AP) kernels (also known as splitting functions), which were first introduced in Ref. [17] for the doublecollinear limit at tree-level, and at the amplitude level by the splitting amplitudes [18,19]. Since then, splitting amplitudes and Altarelli-Parisi kernels have been studied at one-loop [20][21][22][23][24][25] and two-loops [26,27] in the double collinear limit. A higher-order loops analysis has been performed at the amplitude level in Ref. [28]. In the multiple collinear limit, splitting functions were studied at tree-level [29][30][31][32][33][34] and there are some results for the triple-collinear limit at one-loop order [35].
Dimensional regularization (DREG) can be implemented in various ways when performing an explicit computation. This gives rise to different DREG schemes. Since theoretical results have to be compared with experiments, it is expected that they do not depend on the scheme being used.
However, since splitting functions and splitting amplitudes are not physical observables, they can exhibit some scheme dependence. For this reason, it is necessary to understand how to relate the results obtained with different schemes. In the double collinear limit, this topic has been discussed in Ref. [36]. At one-loop level, computations were performed using several schemes choices, for both splitting amplitudes and AP kernels. In particular, in Ref. [37], the scheme dependence for 2 → 2 processes was studied at one-loop level and the authors also suggested a way to relate one-loop AP kernels in some usual DREG schemes.
In this paper, we discuss the scheme dependence of splittings amplitudes at NLO. Starting from the QCD Lagrangian in (4 − 2 ) dimensions, we decompose the gluon field and define scalargluons associated with the extra-dimensional degrees of freedom 2 (see Ref. [38,39]). Using these artificial particles we establish a link among results in several schemes, besides exploring novel DREG parameters' configuration. It is also important to mention that scalar-gluons were useful to solve some theoretical issues related with factorization in QCD when working in DREG [40,41].
This article is organized as follows. In Section 2, we briefly define the different schemes in DREG and introduce scalar-gluons in QCD. We motivate the effective Feynman rules for these objects, starting from a Lagrangian-level decomposition. In Section 3 we discuss the kinematics of the collinear limit and introduce the splitting matrices. In Section 4 we present results for the q → gq splitting at NLO. We recover known expressions, compared them with Catani's formula for the IR-divergent structure and analyze the scalar-gluon contributions. In the last part of that section, we calculate the AP kernels. In Section 5, we tackle the g → qq splitting and put more emphasis in the study of the scalar contributions. The g → gg splitting amplitude is discussed in Section 6. Since photons play a crucial role in today's collider physics (Higgs boson background, study of quark-gluon plasma and jet quenching, etc.), we extend our results to cover the q → γq and γ → qq splittings in Section 7. Finally, conclusions are presented in Section 8.

II. DIMENSIONAL REGULARIZATION AND QCD
In the context of perturbative QCD, when we want to compute higher-order contributions we have to face Feynman integrals. Generally, these are ill-defined because they involve non-integrable expressions. So, it is mandatory to introduce a regularization method to give sense to the theory.
Moreover, QCD has both ultraviolet (UV) and IR singularities, so we need prescriptions to treat them. Due to the simultaneous treatment of UV-IR divergences and gauge-invariance preserving formalism, DREG is one of the most used method to regularize QCD.
Introduced in the late sixties [1,2], the main idea of DREG is to modify the space-time Depending on the DREG scheme, it is possible to keep some quantities living in a 4-dimensional 2 For this reason, some authors also call them -scalars.
space. In some sense, this is equivalent to specify the symmetries of the extended theory. We know that QCD is a quantum field theory on a 4-dimensional space-time which is invariant under the action of the 4-dimensional Poincarè group. When we extend the theory to a D ST -dimensional space-time, it is possible to force a D ST -dimensional Poincarè invariance or just retain the physical 4-dimensional invariance. In this work we will play with these options and explore their consequences in the final results.

A. DREG schemes definition
Let's start with a general four-dimensional quantum field theory (QFT). We know that any one-loop scattering amplitude can be written in the general form where A k are coefficients depending on external momenta, color configuration and the kind of particles involved in the process (both internal and external). The DREG changes the space-time dimension to D ST in order to allow for convergence of loop integrals. Usually, we take D ST = 4−2 (with a complex number) and perform the replacement Using Passarino-Veltman decomposition or any other reduction method, we can solve tensor-type integrals and write them as q q µ 1 q µ 2 . . . q µm D α 1 1 D α 2 2 . . . D αn where {D i } and {p i } denotes the set of possible denominators (with α i ≥ 0) and physical 4-vectors defined in the 4-dimensional unregularized theory, respectively. Note that we are using the D STdimensional metric tensor in this expansion. This is an important point since we can not take the limit → 0 while computing integrals, so it is not allowed to make the replacement η D ST → η 4 until we finish the calculation.
On the other hand, DREG does not impose any specific treatment of the objects that appear in the coefficients A k . Since A k depend on the Dirac's algebra dimension and the number of fermion and boson polarizations, this means that we can change them and set an specific convention for our computations: this is called a DREG scheme. The parameters used to define a DREG scheme in the context of massless QCD (or massless QCD+QED, as we discuss in the last part of this article) are • n g : number of external gluon polarizations, • h g : number of internal gluon polarizations, • n q : number of external quark polarizations, • h q : number of internal quark polarizations, and • D Dirac : dimension of the Dirac's algebra.
There is another subtlety related with the dimensionality of particle's momenta. DREG forces loop momenta to be D ST -dimensional to ensure convergence, but there is no restriction over external momenta. For that reason, we could use them in D ST or 4 dimensions in the context of different schemes. At the amplitude level, we usually consider external particles to be physical, so their momenta are naturally 4-dimensional. However, when we compute squared matrix elements and perform phase space integrals, we can face IR singularities again. This time divergences originate in some regions of phase-space where particles have soft, collinear or soft-collinear kinematics.
DREG can be used to regularize phase space integrals [36], which implies that unobserved external momenta have to be extended to a D ST -dimensional space-time. For that reason, it is also possible to consider external momenta being D ST -dimensional.
Following with DREG schemes definition, we first study how to relate the previously mentioned parameters with the way that computations are performed. First of all, note that the D STdimensional space-time metric can be written as a direct product of a 4-dimensional and a D ST −4- On the other hand, we can introduce vectors and spinors in this space 3 . We write D Dirac = 4 − 2δ , with δ = 0 or δ = 1 to work with a 4-dimensional or a D ST -dimensional algebra, respectively. Spinors are defined starting from a representation R of Dirac's algebra; that is, we have a set of objects where Id refers to the identity in the space where representation R is defined. Since fermions are described by spinors, the number of polarizations of a massless fermion is related to Tr(Id). In 3 In the context of smooth manifolds, vector fields are defined as sections to the tangent bundle and spinor fields arise as representations of a Clifford algebra induced by the metric over the tangent space. So, D Dirac refers to the dimension of the tangent space and it must be D Dirac = D ST by definition. However, in the context of DREG we can treat them independently, since we are ultimately interested in taking the limit → 0 to recover physical results. particular, we can define where Tr Ext and Tr Int denote the trace over external and internal fermionic states, respectively, since we are treating internal and external fermions in an independent way. It is interesting to appreciate that changing h q or n q only modifies contributions which involve Dirac's traces, because using Dirac algebra and cyclic-invariance of traces, it turns out that traces are always proportional to Tr(Id). We introduced the parameters β and β R to write and control the number of fermion polarizations when performing the computations. Now let's turn to the parameter h g which is related to the gluon propagator. Working in an axial gauge, we write the sum over internal gluon's polarizations as where n is a light-like vector which verifies n 2 = 0 and n · p = 0. Here we introduced α R to take into account the number of polarization associated with internal gluons. In particular, we know Using Eq. 9, we see that if α R = 0 then h g = 2 independently of the value of D Dirac , while if we choose α R = 1 then h g = D Dirac − 2. It is important to note that this result relies in the fact that n is the D ST -dimensional null-extension of a four-vector and the metric tensor is diagonal even in D ST -dimensions. Also, we can appreciate that choosing δ = 0 (i.e D Dirac = 4) removes the dependence on α R in Eq. 9.
To control the number of external gluon's polarizations we define where Q is an arbitrary null-vector which fulfills Q 2 = 0 and Q · p = 0. When performing the explicit computation we will take Q = n with the sake of simplifying the intermediate steps.
Again, using Eq. 11 where we express the result explicitly in terms of α.
At this point, it is important to recall some properties of the gluon's polarization tensors d µν (p, n) and d Ext µν (p, n). Since we are working in the light-cone gauge (LCG), these objects should fulfill the following identities: • d Ext µν (p, n)p µ = 0 (orthogonality to external momenta p), and These requirements are related with some physical restrictions. The first condition is due to the gauge choice, while the fact that external gluons are massless guarantees the validity of the second identity. The last condition is imposed in order to recover orthogonality with external momenta when the virtual particle is on-shell. Having introduced a parametrization for polarization tensors in Eqs. 9 and 11, we show that where we use strongly that n is a light-like 4-vector (and n µ η D Dirac −4 µν = 0). Something similar happens if we consider external momenta as the null-extension of a light-like 4-vector, i.e. p = p (4) ⊕ 0. In this case, we obtain which shows that all the requirement are fulfilled for external momenta, independently of the dimension of the space in which external momenta live. However in all DREG schemes, internal momenta have to be expressed as p = p (4) ⊕ p (D ST −4) (i.e. with a non trivial component in the transverse space) and when we contract with d µν (p, n) we get which shows that, for certain combinations of parameters, propagators fail to fulfill some physical consistence conditions. In the following, we discuss deeply about this fact, performing some explicit computations and showing that failing to verify this conditions could lead to some unexpected IR divergences.
Having introduced the possible parameters that we can modify in the context of DREG, let's explain how to recover some of the most frequently used schemes. In conventional dimensional regularization (CDR) [1,2,42], internal and external particles are treated in the same way. We consider that particle's momenta is D ST -dimensional, gluons have 2−2 polarizations and massless fermions only have 2 polarizations. In other words, CDR uses n g = h g = 2 − 2 and h q = n q = 2, with D Dirac = 4 − 2 . Equivalently, we can work in this scheme setting α = 1, α R = 1 and δ = 1, according to our definitions.
On the other hand, we can set external particles in four-dimensions while keeping internal ones in D ST -dimensions. This is the 't Hooft-Veltman scheme (HV), first introduced in Ref. [2].
Sometimes it is preferable to treat all the particles as 4-dimensional objects, although internal ones have their momenta extended to a D ST -dimensional space. In the four-dimensional helicity scheme (FDH) [41,43], massless fermions and gluons have 2 polarizations states, independently of being internal or external particles. It means that FDH scheme is defined by setting n g = h g = 2, n q = h q = 2 and D Dirac = 4, or δ = 0 and α = 0 in our convention. The main advantage of this configuration is the possibility of using many identities and properties derived from the helicity method, which can be used to obtain very compact expressions.
In this scheme, both external and internal particles have 4-dimensional polarization vectors, but external momenta are continued to D ST -dimensions. This forces us to include scalar-gluons in order to achive consistency. A deep discussion about dimensional reduction schemes and its implementation can be found in Ref. [41] 4 .
Beyond these three well-established schemes, changing the parameters α R , α, δ and the value of Tr(Id), we can construct more configurations. Variations of CDR, HV and FDH with Tr(Id) = 4−4 were studied in Ref. [36]. In particular, in that paper, the authors analyzed the consequences of choosing those toy-models when performing a full NLO computation with the subtraction method. Since we can easily modify the values of β and β R in our codes, we give most of our results for an arbitrary number of fermion polarizations. In the last part of this article, we also discuss the results computed in the toy-scheme (TSC) defined by n g = n q = 2−2 = h q = h g which was introduced in Ref. [36]. Specifically, we show that this scheme preserves the supersymmetric Ward identity at one-loop level.
Aside from allowing different values of h q and n q , here we also discuss hybrid-schemes that use D Dirac = 4 − 2 (δ = 1) and h g = 2 (α R = 0), with the possibility of setting n g = 2 − 2 (α = 1) or n g = 2 (α = 0). To make the discussion easier, we call them hybrid-scheme A (HSA) and hybrid-scheme B (HSB), respectively. A summary of all the schemes treated in this paper is displayed in Table I. Although these new schemes seem to a valid choice, we can anticipate that they are inconsistent unless we add some scalar-gluon contributions.

B. Scalar-gluons: Lagrangian level decomposition
An interesting fact is related to the appearance of new scalar-type particles when we use certain DREG schemes in a D-dimensional space. We can decompose D-vectors into 4-vectors plus D − 4 scalar particles [44]; this forces us to introduce new Feynman rules for these particles and, of course, new diagrams contribute to the scattering-amplitudes. Note that this decomposition suggests that non-physical degrees of freedom can be absorbed into a certain amount 5 of scalar-particles.
To get the Feynman rules for scalar-gluons, let's start with the usual 4-dimensional massless-QCD Lagrangian density, where {i, j} are color indices, µ T a ij is the covariant derivative and we are summing over the possible quark 4 In this reference, the authors group DRED and FDH into the category of dimensional reduction schemes. Moreover, they are called old and new dimensional reduction, respectively. It is important to note that these schemes are not equivalent, but they could lead to the same results. flavors. Knowing that kinetic terms are associated with propagators, let's expand the interaction component which can be written as with which are associated to the fermion-gluon, three-gluon and four-gluon vertices, respectively. Note that we take D = 4 − 2 as the space dimension.
Now, let's consider that the space-time metric η D µν is diagonal and can be expressed as a direct product. We can perform the decomposition withÂ,γ living in a 4-dimensional space andÃ,γ in the unphysical D − 4-dimensional transverse space. Also, we have Using the expressions for the interaction terms in the Lagrangian Eq. 16 we get where we must take into account that some indices live in the 4-dimensional physical space, while others stay only in the transverse space. In Fig. 1 the available vertices are drawn. We have six possible interaction vertices which involves scalar-gluons: 2fermions-scalar, 2gluons-scalar, 2scalars-gluon, 3scalars, 2scalars-2gluons and 4scalars. After identifying the Lagrangian terms that originate the possible scalar interactions, we are able to deduce the corresponding Feynman rules. Introducing the functions and the usual Feynman rules for 4-dimensional QCD reads • fermion-gluon-fermion vertex −ıg s T aγµ , • and quadruple-gluon −ıg 2 where we are projecting Lorentz index to the 4-dimensional space through the contraction with η 4 . In Fig. 2 we write explicitly the usual QCD rules, to clarify momentum sign and ordering conventions. From these rules and conventions, we get the following associated Feynman rules for scalar interactions; • fermion-scalar-fermion vertex −ıg scalar s µ T aγµ , • gluon-scalar-gluon vertex −g scalar • scalar-gluon-scalar vertex −g scalar • 2scalar-2gluon −2ı g scalar • and 4scalar −ı g scalar where we used the labeling introduced in Fig. 1 and the conventions shown in Fig. 2. The additional factor 2 in the 2scalar-2gluon vertex comes from the expansion of the Lagrangian in Eq. 27.
As a final remark, note that when scalar-gluons appear we make the replacement g s → g scalar s .
This is done because we consider scalar-gluons as a new kind of particles which are not necessarily related with vector-gluons. So, following Landau's principle, we have to write the most general Lagrangian compatible with certain reasonable requirements. But these requirements do not exclude the possibility of having different coupling constants, so we introduce g scalar s and treat it independently of g s . Although g s = g scalar s at leading order in g s , these couplings do not evolve in the same way and, in consequence, can differ at higher-orders (see [38,45]).

C. Effective Feynman rules and other considerations for scalar-gluons
Working with scalar-gluons involves taking into account some technical details. In order to be able to write scattering amplitudes that include scalar-particles, let's motivate some effective Feynman rules and explain other useful points.
Let's start with the gluon propagator, in an axial gauge where d µν (p, n) is given by Eq. 9, using a null-vector n. Here we explicitly indicate that we are using the Feynman prescription to compute propagators. However, in the following we will omit the term +ı0 in propagator denominators, although its presence is always understood. As we shown previously, the number of gluon polarizations can be modified changing α R and δ. Working the definition of the propagator (or α R = 1 and δ = 1 with our parametrization), in which case Also it is possible to decompose them in a 4-dimensional gluon plus scalar-gluons, by setting η = η 4 in D G and using the propagator for the scalar-gluon component. If we count the number of polarizations in this case, vector- polarizations. For this reason, scalar-gluons have to be included when we set the number of polarizations of internal gluons to 2 − 2 while working with a 4-dimensional Dirac's algebra. In other words, HV results could be recovered adding to the FDH computation the corresponding scalar-gluon contribution.
It is worth noting that a completely similar analysis can be performed with external gluons. If they are treated as D ST -vectors, we can decompose them as 4-dimensional gluons plus D ST − 4 scalar particles. This implies that we can also use scalar-gluons as external legs to compensate the number of degrees of freedom of the system when working with D Dirac = 4. Explicitly, Using the definitions of V 3g and the induced rules for scalar-gluons vertices at Lagrangian level, we can get some effective rules to work with these particles. To get them we will modify the expressions given before, using the fact that η D ST is diagonal (i.e. it does not mix physical and transverse contributions). So, starting with the triple vertex we have: • g scalar s µ f abc η 4 νσ (p 2 − p 3 )μ for the 2gluon-scalar vertex; • g scalar s µ f abc η νσ (p 2 − p 3 ) µ for the 2scalar-gluon vertex; • and g scalar s µ V Cin 3g (pμ 1 , pν 2 , pρ 3 , a, b, c) for the 3-scalar vertex Note that these rules agree with the usual form of Feynman rules for vector-scalar interactions.
Also, here η ρσ can be interpreted as a delta function whose value is 1 if scalar-particles have the same index and 0 otherwise.
Following the same ideas, we can simplify quadruple interactions and we get these rules: σρ η µν for the 2gluon-2scalar vertex; • and −ı g scalar where, again, we see agreement among these expressions and the ones associated with standard quadruple scalar-vector interactions.
Finally, let's make a comment about the fermion-scalar interaction. This is the only vertex which involves dealing withγ matrices. Since D Dirac = 4, these extra-gamma matrices act trivially over spinors, so we do not have to include them inside spinor chains: this leads to helicity-violation interactions. Moreover, if we have twoγ matrices inside a chain, using the fact that {γ µ ,γ ν } = 0 and {γ µ ,γ ν } = 2η µν Id we can get ride of the transverse-dimensional indices. We give an explicit example when computing the q → gq splitting amplitude at NLO.

D. Computational implementation
We implemented the computation in Mathematica and we used FeynCalc (version 8.2) [46] to handle Dirac's algebra. FeynCalc used D as the dimension of Dirac's algebra; in particular, since it was used to perform Dirac and Lorentz algebra, then D Dirac = D. To compute integrals we used the results shown in the literature (for example, see Ref. [25]) and the integration by parts method (IBP) [47], implemented through the Mathematica's package FIRE [48,49]. We set d = 4 − 2 as the space-time dimension in this package.

III. COLLINEAR LIMITS OF SCATTERING AMPLITUDES IN QCD
To study the double collinear limit of scattering amplitudes, the first step consists in identifying the relevant kinematical variables. We describe the momenta of the external particles using the vectors p µ 1 and p µ 2 , which refer to the outgoing particles. Here it is important to note that µ is a Lorentz index which runs over the space-time dimension D ST = 4 − 2 . Since they refer to external particles, we assume that components along the additional dimensions are zero, so p µ 1 and p µ 2 behave like usual four-vectors. The momenta of the incoming particle can be obtained from momentum-conservation rules, p µ 12 = (p µ 1 + p µ 2 ), and its virtuality is given by Since we are defining η D ST µν as a flat-extension of the usual four-dimensional Minkowski metric, inner product with four-vectors behaves like a projection. Then, due to the fact that we are working with massless QCD, p 2 1 = 0 = p 2 2 . To describe the collinear limit, we introduce a null-vector n µ (n 2 = 0) that satisfies n · p 12 = 0 and that is the zero-extension of a usual four-vector (we call them with the same name to reduce the notation). Choosing that vector is equivalent to introduce a preferred direction in space-time, which allows us to get rid of unphysical degrees of freedom. In other words, we use n to settle in the light-cone gauge. Working in the light-cone gauge has advantages (for example, we do not have to consider diagrams with ghosts), but it introduces an extra-denominator in loop-integrals which makes them harder to compute. However the most important benefit of choosing a physical gauge is the possibility to exploit collinear factorization properties in an easier way.
Returning to the kinematics of the double collinear limit, we can introduce a collinear nullvector,P , which satisfiesP 2 = 0, n ·P = 0 and p µ 12 →P µ when s 12 → 0. This allows us to define the momentum fraction of particle i as where we have the constraint z 1 + z 2 = 1 and therefore we can describe the collinear limit using only the scalar variable z 1 . (Strictly speaking, we need to use also s 12 and n ·P , but since they are dimensionful we can guess their scaling properties and factorize them). We can think about z i as a measure of the contribution of particle i to the longitudinal total momentum relative to n.
In other words, n is used here to parametrize the approach to the collinear limit.
On the other hand, it is necessary to take into account the transverse component of the outgoing particles relative to the longitudinal component proportional toP . To do this we define k µ ⊥ , which verifies n·k ⊥ = 0 = k ⊥ ·P . Due to the relations among n, k ⊥ andP we can use them to parametrize the momentum of the outgoing particles [34] as where z i is the momentum fraction associated with particle i and k ⊥ 2 = −z 1 (1 − z 1 )s 12 . Note that this parametrization is consistent with the fact that both outgoing particles are on-shell and massless. On the other hand, when performing the explicit computation we do not need to express external momenta in terms ofP , n and k ⊥ : this decomposition is relevant to simplify spinor chains or scalar products that appear in matrix elements. Also, we use n · p 12 = n ·P because with / p 12 u(n) = / P u(n).
After describing collinear kinematics, let's settle some conventions to write scattering amplitudes in the context of massless-QCD with photons. Due to the presence of color charges, we will express matrix elements in color+spin space [12,50]. A general n-particle matrix ele- where |M a 1 ,...,an (p 1 , p 2 , . . . , p n ) is a vector in color+spin space. We need to remark that external legs are being considered as on-shell particles (and, moreover, QCD partons are massless).
Let's consider an n-particle scattering amplitude and assume that two particles, labeled as 1 and 2, become collinear. Since we are interested in studying the most divergent part of this kinematic limit, we will only consider diagrams in which 1 and 2 come from a parent leg P , as shown in Fig. 3. It is important to note that, in order to simplify factorization properties, we have to perform the computation in a physical gauge (for example, see Ref. [51]).   where we have introduced the amputated amplitudes A and a general propagator Prop, which depends on the kinematics and the class of particles involved in the process. It is important to note two facts: we are summing over all possible flavors of particle P , and p 12 is the momentum vector associated with the intermediate particle.
Since we are working in massless-QCD, P can be a gluon or a quark. If P is a quark, a 1 and a 2 have to be a quark and a gluon. On the other hand, if P is a gluon, a 1 and a 2 can be a quark-antiquark pair or a gluon-gluon pair. It is important to notice that a quark-antiquark pair can become collinear because we are considering them as massless particles, but gluons can always be collinear or soft.
Let's analyze what happens with each possible choice of P . If P is a quark, then its propagator where we used the definition of the light-like vectorP and we keep only the most divergent contributions in the limit s 12 → 0. Here, {i, j} are color indices associated to the fundamental representation of SU(N ). SinceP is a null-vector, it is possible to considerP as the momenta of a massless quark. So, using the completeness relation of massless spinors, we are able to use the expressions / P = λ phys.pol.
with λ being a label for possible physical polarizations of intermediate quark and antiquark, respectively. These considerations leads us to rewrite Eq. 42 as Now, going back to Eq. 41, we obtain where, in the last line, we rearranged the factors to form an n − 1 matrix element associated with a process which replaces legs 1 and 2 with a unique on-shell massless particle P . This can be done because we are working in a kinematical region where 1 2, so s 12 → 0 and we put the divergent factors to the left-side of the propagator. In other words, the replacement p 12 →P is possible in (p 12 , p 3 , . . .) because it is finite in the collinear limit.
If we consider now the case in which P is a gluon, we have to write the propagator as where η is a metric tensor which depends on the number of polarizations of gluons. As done for the quark case, we can use the definition ofP and perform the expansion and together with the completeness relation leads us to the expression which is a valid approximation in the collinear limit. Applying these results to Eq. 41 we get where, again, we are able to rearrange the expression in such a way that the first factor contains all the divergent contributions and the second one is a reduced-matrix element for a n − 1-particle process.
From Eqs. 45 and 50, we can motivate the definition of splitting matrices and amplitudes.
Working in the double-collinear limit, the quark initiated splitting matrix can be written as with a 1 and a 2 being a gluon and a quark, respectively. When the parent particle is a gluon, we get Sp g→a 1 a 2 = 1 being a 1 and a 2 a quark-antiquark or a gluon pair. In both cases, |A P,a 1 ,a 2 (p 12 , p 1 , p 2 ) is the amputated scattering amplitude associated to the process P → a 1 a 2 , without being projected over the color+spin space. If we project Sp over color+spin space we get the so-called splitting amplitudes. To recover splitting functions (as defined in Ref. [25]), we just have to remove color information from splitting amplitudes.
Now it is important to note that we left p 12 as incoming momenta in the amputated amplitude, instead of usingP . This is related to the fact that the presence of divergences in the definition of splitting matrices forces us to regularize them and keep the s 12 dependence explicitly. For that reason we must consider that the incoming particle is slightly off-shell and include all possible Feynman diagrams, also those which include self-energy corrections to the parent leg. We will emphasize this fact when computing explicitly some scattering amplitudes at NLO.
Finally, we have to remark that it is possible to get the divergent contribution to the splitting matrices at NLO without performing a full computation. For the double-collinear limit, according to Ref. [14] we can write with Sp (1) H containing only the rational dependence on the momenta p 1 , p 2 andP , and which contains all the divergent contributions and non-rational functions of z 1 . Here C i are the Casimir factors associated with the parton a i (C i = C A for gluons and C i = C F for quarks), γ i depend on the flavor of a i and c Γ is the D-dimensional volume factor associated with one-loop integrals, i.e.
If N f is the number of quark flavors, then and b 0 = γ g is the first perturbative coefficient of the QCD β function, according to our normalization. Besides that, the function f ( , z) is given by and it is associated with the kinematical behavior of matrix elements in the collinear limit.

IV. THE q → qg SPLITTING MATRIX
When working in the LCG, the presence of internal gluons makes more difficult to perform an explicit computation. So, we start with the q → gq splitting and we explain the differences among schemes. At NLO we can write the corresponding splitting matrix as where the LO contribution is Even at LO, we can decompose γ µ =γ µ +γ µ when considering D Dirac = 4 − 2 . This leads to the expression which includes an helicity-violating term that contributes only in CDR or HSA schemes. However, since gluons are treated as D-dimensional vectors in CDR, it is not required to separate explicitly the helicity-violating term. The situation is going to be different in HSA scheme because the presence of both 4 and D ST -dimensional metrics leads to a non-equal mixing betweenū(p 2 )γ µ u(P ) andū(p 2 )γ µ u(P ). It is important to appreciate that we are starting from the amputated amplitude related with q(p 12 ) → g(p 1 )q(p 2 ). This explains why we must take into account self-energy corrections to the incoming particle. In other words, to calculate the NLO corrections to the splitting matrix Sp q→gq we need to include all the diagrams shown in Fig. 4. However, it is necessary to take into account other kind of contributions to explore consistently the different schemes. As we mentioned in Section 2, when we treat QCD in the context of DREG, it is possible to decompose D ST -dimensional gluons into 4-dimensional vectors and scalar particles. Due to the fact that we can make that differentiation when drawing Feynman diagrams, it is useful to introduce the following classification of diagrams: • standard QCD contributions (STD); • helicity preserving interactions mediated by scalar gluons (SCA-nHV); • and helicity-violating interactions (SCA-HV).
To compute STD contributions, we start from 4-dimensional QCD and draw the associated Let's start describing the standard NLO QCD contribution. It can be expressed as where Sp (1,i) q→gq refers to the diagram i ∈ {A, B, C}, as shown in Fig. 4. Writing each contribution we have, where we are not making any distinction between 4 and D ST -dimensional gluons. In particular, in the context of HSA scheme, we should interpret because there are 2 − 2 gluon's degrees of freedom but vector-gluons have only two polarizations (setting α = 0 and α R = 0) while the remaining polarizations are associated with scalar-gluons.
Also, it is useful to note that Sp (1,A) q→gq can be rewritten as where Σ(p 2 12 ) is the NLO correction to quark self-energy. Except for this diagram, all the others correspond to the ones that appear when computing q → gq with massless on-shell particles. with where we used the Feynman's rules previously obtained at Lagrangian level (see Section 2). Now let's simplify this expressions using some properties of Dirac matrices in D-dimensions: we want to show explicitly that it is possible to use the effective Feynman rules introduced in the end of Section 2. Focusing in Sp (1,A ) q→gq , note that this contribution depend on However, since p 12 is a physical momenta then / p 12 = (p 12 ) σγ σ . Using that {γ α ,γ σ } = 0 and γ ργ ρ = −2 Id, then we have where we used Passarino-Veltman (PV) decomposition to write the vector-type integral in the second term. Due to the fact that vector-type integrals only depend on physical vectors then we can repeat the procedure performed in the first term and we obtain which is equivalent to use the effective scalar rules discussed in the last part of Section 2. Note that we have not used the fact that D Dirac = 4, which implies this result immediately.
The situation is analogous when we move to Sp (1,C ) q→gq , but some subtleties appear when treating Sp (1,B ) q→gq . That contribution depends on where {ρ, ν} run over the non-physical dimensions. Depending on the number of external gluon polarizations, µ can live in 4 (n g = 2) or in D ST -dimensions (n g = 2−2 ). Using PV decomposition, the tensor-type integrals can be expanded as with the inclusion of a term proportional to the D ST -dimensional metric tensor. Replacing these expansions in Eq. 74 and using that p 12 and p 2 are 4-vectors, we obtain sinceγ ρ / p i = − / p iγ ρ and A ij = A ji due to symmetry properties. On the other hand, where we used that Dirac's algebra dimension is D Dirac and it is equal to the number of gamma matrices available. So, after applying these properties and the fact that we can rewrite Int B as where we can always express 4 − D ST = η ρν (−(η ) ρν ). Note that the metric tensor inside the parenthesis comes from commuting and symmetrizing the productγ ργν . Also, the contributions involvingγ µ violate helicity conservation, so they vanish when we restrict external particles to have helicity configurations compatible with standard QCD interactions.
Summarizing these observations, we conclude that the replacement −η ρνγ ρ γ α γ µ γ βγν → (−2 )γ α γ µ γ β is valid. Thus we get an effective Feynman rule for scalar-gluons interaction with fermions, which consists in considering them as scalar particles with propagator −2ı p 2 +ı0 (see Eq. 32) and remove the corresponding Dirac matrix in the vertex. On the other hand it is useful to remember that in usual DREG schemes, if scalar-gluons are introduced then we have to set D Dirac = 4. But this limit has to be taken after replacing integrals. In other words, it is possible that some new terms (i.e. not present in the expressions when using effective Feynman rules for scalar-gluon) survive when applying directly Lagrangian level Feynman rules. But this terms are always proportional to integrals which vanish in the limit D Dirac → 4. This situation occurs in q→gq because there is a term proportional to γ α γ µ γ α = −2(1 − )γ µ (see Eq. 81). So, after this discussion, we can rewrite the SCA-nHV contributions as where we used the effective Feynman rules for scalar-gluons setting D Dirac = 4.
Finally, we want to make a brief comment about SCA-HV components. When working in HSA/HSB schemes it is possible that STD contributions mix helicity-preserving and helicityviolating terms, whose origin is the contraction of 4-dimensional metric tensors (coming from the gluon propagator) with D ST -dimensional structures. We discuss this point in the next subsections, using the results for Sp q→gq to give an explicit example.

A. Amplitude level results
Following with the study of q → gq splitting amplitude, we performed the computation without specifying the polarization of the involved particles. This implies having larger spinorial structures and more complex tensor-type integrals, but this will allow us to compute contributions to the NLO Altarelli-Parisi kernel in an easier way.
Let's start with the NLO standard-QCD contribution to the splitting matrix. After writing explicitly the corresponding Feynman diagrams and replacing the involved loop-integrals, we find that T a C (STD,1) q→gqū (p 2 )/ (p 1 )u(P ) where the coefficients C (STD,i) q→gq are given by where δ controls Dirac's algebra dimension and we left α R as a free parameter. Note that there is a term that explicitly involves an helicity-violating interaction. It is proportional to 1 − α R and only contributes when we work in HSA scheme (α = 1) because external gluons must have 2 − 2 polarizations in order to allow for this kind of interactions. Also, it is worth noting that modifying α R only introduces O( 2 ) differences in coefficients C (STD,1) q→gq and C (STD,2) q→gq . However, expanding C (STD,3) q→gq we find which implies that Sp (1,STD) q→gq acquires an additional contribution to the double pole which is If we want to check our calculations, we can set α R = 1 to recover well known results in FDH (δ = 0) and CDR/HV schemes (δ = 1). In particular, when using CDR we have to assume that µ (p 1 ) is a D ST -dimensional vector, while in the remaining schemes µ (p 1 ) is associated to a 4dimensional space. It is important to appreciate that we used the properties (p 1 ) · n = 0 (related to the definition of the null-vector n) and (p 1 ) · p 1 = 0 (because the outgoing gluon is a physical massless vector particle, with transverse polarization) to simplify the expressions.
Following with the study of different contributions to the splitting amplitude, we can compute Sp (1,SCA−nHV) . After replacing integrals and performing some simplifications, it can be expressed as where we consider a 4-dimensional Dirac's algebra. Note that this expression is simpler than the STD contribution presented before. This is due to the absence of two-gamma matrices in the spinorial chain, which were replaced by a -dimensional metric, and the simplification of some gluon-propagators. Also, it is worth noting that SCA-nHV terms are finite in the limit → 0, so they can be added to the other contributions without modifying the divergent structure. This allows us to interpret the addition to the SCA-nHV terms to the splitting as a DREG scheme choice. Moreover, note that from Eqs. 85 and 90 we can recover the relation which tells us that HV results can be obtained from FDH ones by just adding SCA-nHV contributions. This is a really interesting property, because sometimes it is easier to perform the computation using 4-dimensional algebra. Moreover, this relation is still valid when we set the polarization of external particles to the possible 4-dimensional physical values. And, in that situation, we can take advantage of working in FDH scheme because we can apply a wide range of novel techniques, such as the helicity method.

B. Scheme dependence and divergent structure
Following with the analysis of our results, we can test the decomposition suggested in Eq. 53.
First, we assume that α = 0 (i.e. we neglect HSA scheme) and use only STD diagrams. If we expand in series around = 0 and rearrange divergent contributions, we find H,q→gq + I C,q→gq p 1 , p 2 ;P Sp (0) q→gq , with I C,q→gq p 1 , p 2 ;P = c Γ g 2

Sp
(1) where we left δ and α R as free parameters. The structure of I C,q→gq exactly agrees with the expected singular behavior of unrenormalized splitting amplitudes. However, some discrepancies appear in the finite contribution. According to Ref. [14], Sp H only contains rational functions of z and . This is completely true when α R = 1, since it reduces to and the finite remainder becomes But when considering α R = 0, this contribution involves a non-vanishing combination of hypergeometric functions, which can not be expressed using only rational terms. So, when we work in HSB scheme, Sp H is no longer a pure rational function. The situation becomes worse if we choose to work in HSA scheme, setting α = 1 and α R = 0.
In that case, it is not possible to cast Sp (1,STD) in the form expressed in Eq. 92 because the divergent structure verifies which involves additional poles that can not be absorbed in any term proportional to the LO splitting amplitude. This indicates that something else has to be added when performing computations inside HSA scheme (or, conversely, that the definition of HSA scheme must be different).
In fact, we need to take into account all the scalar-gluon contributions, both SCA-nHV and SCA-HV. To understand this, we remind the reader that HS schemes assume D Dirac = 4 − 2 = D ST This relation tells us that the consistent version of HSB is the HV scheme (δ = 1 and α = 0). On the other hand, in HSA scheme we must allow the presence of scalar-gluons as external particles.
Again, this is equivalent to add the same kind of diagrams but decomposing the outgoing gluon polarization vector as µ =˜ µ +ˆ µ . In other words, if we add all the contributions required to cure the inconsistencies of HSB, we just end in CDR scheme (δ = 1 and α = 1). We will emphasize this point in the following subsection, when computing Altarelli-Parisi kernels.
In summary, after analyzing the scheme dependence of our results for q → gq splitting amplitude and comparing them with Catani's formula (Eq. 92), we conclude that HSA/HSB configurations are not suitable choices for performing calculations. Instead, we will use CDR, HV and FDH schemes, with the possibility of changing the number of fermion polarizations (playing with the parameters β and β R previously defined).

C. NLO corrections to AP kernels
Having LO and NLO contributions to the splitting matrix we can obtain the NLO correction to the Altarelli-Parisi (AP) kernel q → gq. In order to do that, we use the expansion where we must consider the regulator as a complex-valued parameter. If we sum over the physical polarization states of outgoing particles, sum over colors (averaging the incoming ones) and project over the helicity-space of incoming particles, we obtain the polarized AP kernels. Also, it is possible to sum and average over the physical polarizations of the parent parton, which leads to the definition of the unpolarized AP kernels. 6 As expected, the sum over polarizations depend on the scheme being used. If we consider FDH or HV, external particles have physical 4-dimensional polarizations, but when we set in CDR, they live in a D ST -dimensional space. So, in the last scenario, a scalar-gluon can be considered as an external particle, which implies that we must also consider spin-flip contributions at amplitude level. If we compute STD contributions to Sp q→gq , we can obtain AP kernels in any scheme. It is important to note that, when considering CDR scheme, spin-flip contributions are hidden inside the definition of the D ST -dimensional polarization vector, as we saw in Eq. 60. So, we do not need to include explicitly external scalar-gluons, but we can use them to give a physical interpretation to some contributions.
After this brief discussion, let's show explicit results. Starting at LO, we get for the polarized and unpolarized kernels, respectively. Note that when summing over external fermions polarizations, we get a global factor Tr(Id) = 4 − 4 β multiplying our results, but it cancels with the average factor. So, q → gq AP kernels are independent of the number of fermion polarizations. Also, we can prove that since the kernel is diagonal in helicity space. For this reason, we only present the NLO correction to the unpolarized kernel, which is given by where α = 1 in CDR and α = 0 in FDH/HV schemes. As expected, we can appreciate that NLO corrections are independent of β R and β. On the other hand, it is important to take into account that we must consider only the real part of the r.h.s.
where φ denotes external scalar-gluons and we use the replacement suggested in Eq. 34. To obtain the NLO correction to this result, it is necessary to take into account some SCA-HV diagrams and compute the corresponding splitting matrix. Since we are decomposing only external gluons, the required contributions can be recovered from Sp (1,STD) by just making the replacement µ (p 1 ) →ˆ µ (p 1 ). So, we can write with the corresponding Feynman diagrams shown in Fig. 6. After summing over external particles polarizations and averaging, we get We can appreciate that which reflects the fact that additional gluon polarizations can be interpreted as scalar particles, and, in consequence, that it is possible to recover CDR results working with external 4-dimensional gluons and adding the remaining degrees of freedom treating them as scalar-particles. Of course, this separation has to be performed with each external gluon to be consistent, which makes a bit cumbersome to carry out this analysis in general.

V. THE g → qq SPLITTING MATRIX
In the previous section we treated in great detail the splitting amplitude q → gq. Here we focus in the process g → qq, which is closely related to the first one. However, due to the fact that it is initiated by a vector particle, there are some differences.
As a starting point, we write where the LO contribution is where p i is the physical momentum of particle i and we associate the massless vectorP to the incoming gluon in the collinear limit, in spite of having a momenta p 12 which verifies p 2 12 = s 12 . The NLO standard-QCD contribution can be expanded as where Sp (1,i) g→qq refers to the diagram i ∈ {A, B, C, D}, as shown in Fig. 7. Note that diagrams A and D expands the self-energy correction to the incoming gluon with a tiny virtuality s 12 . For that reason, we can rewrite their contribution as where Π(p 2 12 ) can be extracted from Π µν (p 12 ) after contracting with two gluon polarization vectors * µ (P ) ν (P ). (See Appendix B for further details on the computation of Π(s 12 ) and Σ(s 12 )). For this process, there are four possible SCA-nHV diagrams which contribute to the amplitude.
Following Fig. 8 and using Feynman rules at Lagrangian level, we have with

Sp
(1,C ) When we discussed the structure of the contributions to q → gq splitting amplitude, we mention the possibility of having q 2 -type integrals. Here we face the problem explicitly when analyzing g→qq . If we expand the triple gluon vertex, we find Since the scalar contribution is computed setting D Dirac = 4, we can write the involved integral as where P is a permutation of Lorentz indices {ρ, ρ 1 , µ, σ}, k i ∈ {p 12 , n} and Q is a ordering of The important fact here is that Int A only has 4-dimensional components, which implies (1,A ) g→qq = 0 when using a standard scheme for scalar-gluon contributions. The remaining terms of the splitting matrix can be written as

Sp
(1,C ) where we used the same argument presented in the previous section to make the replacementŝ Related with the scalar-gluon contributions, here we saw an important fact. Although many diagrams can be constructed by using the effective rules, some of them are going to be zero due to the presence of only q 2 -integrals. This integrals appear when a transverse index contracts with the loop-momentum q. So, to avoid them, transverse indices should form closed chains, that is which is equivalent to say that each chain is going to be proportional to the trace of the transverse metric tensor (Tr [η ] = −2 = (η ) µ µ ).

A. Amplitude level results
Before showing the explicit results for the g → qq splitting matrix, let's work out the possible spinorial structures which are going to appear. First of all, LO contribution is proportional to FIG. 8. Feynman diagrams associated with SCA-nHV contribution to g(P ) → q(p 1 )q(p 2 ) at NLO. u(p 1 )/ (P )v(p 2 ) and there will be a term proportional to this in Sp (1) g→qq . Due to the symmetry p 1 ↔ p 2 , the propertiesū and the presence of only two physical vectors (p 12 and n), we can only have one additional spinorchain with one gamma matrix inside:ū(p 1 )/ nv(p 2 )p 1 · (P ). Although there can be spinor-chains of up to five gamma-matrices, Dirac's algebra and the previous properties allow to reduce them to combinations ofū(p 1 )/ (P )v(p 2 ) andū(p 1 )/ nv(p 2 )p 1 · (P ). For these reasons, after replacing Feynman integrals in the expressions for Sp for the NLO standard contribution, where the coefficients C (STD,i) g→qq are given by where we set α R = 1 since we will not use HSA/HSB schemes here. It is interesting to note that the full NLO correction to the splitting matrix is proportional to Sp g→qq . Besides that, we can appreciate that C (STD,1) g→qq is symmetric when interchanging particles 1 and 2.
Again, when using α = 1 we have to assume that µ is a D ST -dimensional Lorentz index, while in the remaining schemes µ is associated to a 4-dimensional space. Moreover, if we rearrange the contributions to Sp (1) g→qq in the last scenario, we find that it verifies with

Sp
(1) as expected according to Eqs. 53 and 54. In contrast to the q → gq splitting, Sp g→qq depends on β R . However, this parameter seems to define a well-behaved scheme since it respects the universal divergent structure of splitting amplitudes and the finite remainder is kept composed only by rational functions.
On the other hand, for the scalar-gluon contribution we have and we can recover the relation which tells us, again, that HV results can be recovered from FDH ones by just adding SCA-nHV contributions.

B. NLO corrections to AP kernels
Finally we can compute the contributions to both polarized and unpolarized AP kernel. For the LO contribution we get where we can appreciate that the results depend explicitly on β (i.e. the number of external fermions polarizations). Due to the fact that Sp (1) g→qq is proportional to LO, NLO corrections to AP kernels can be written as where we kept only the real part of the r.h.s. We can appreciate that this expressions depends on both β and β R , and it is not possible to cancel this dependence by setting β = β R . But it is interesting to appreciate that the additional factors in Eq. 133 disappear in TSC scheme.

VI. THE g → gg SPLITTING MATRIX
Finally, we arrive to the g → gg splitting amplitude. It is worth noticing that this case involves dealing with many properties of polarization vectors, but it has the advantage of being free of spinor chains. For that reason, here we deal only with scalar products which are well-defined in DREG for every value of D.
As done with the previous configurations, the splitting matrix can be decomposed as where the LO contribution is where p i is the physical momentum of particle i and (T a (A)) bc = ıf abc are the generators of SU (3) C in the adjoint representation.
The NLO standard-QCD contribution can be expanded as FIG. 9. Feynman diagrams associated with the standard QCD contribution to g(P ) → q(p 1 )g(p 2 ) at NLO. Here the incoming gluon is off-shell and its virtuality is (p 12 ) 2 = s 12 .
being Sp (1,i) g→gg associated with diagram i ∈ {A, B, C, D, E}, as shown in Fig. 9. We have to remark that due to symmetry properties, diagrams C and D only describe the associated topology. In other words, there are two diagrams C (and D), which are obtained from the displayed graph by interchanging particles 1 and 2; Sp (1,C) and Sp (1,D) include the sum over all the possible relabellings of final-state particles associated with the process.
On the other hand, diagrams A and B expands the self-energy correction to the incoming gluon with a tiny virtuality s 12 . As we have done in the g → qq splitting, we can rewrite their contribution as with Π(p 2 12 ) given in Appendix B. When dealing with the scalar-gluon contribution, we find many possible diagrams. However, as we have seen in the previous computations (explicitly in Sp (1) g→qq ), the only non-trivial terms arise from taking the trace of transverse metrics. In other words, transverse indices have to form a closed chain and be completely contracted with metric tensors; otherwise, we will have q 2 -integrals, which are set to zero when D Dirac = 4. So, following Fig. 10 and using effective Feynman rules for scalar-gluons, the SCA-nHV contribution can be written as FIG. 10. Feynman diagrams associated with the scalar-gluon contribution to g(P ) → g(p 1 )g(p 2 ) at NLO.
We only consider diagrams which contribute non-trivially to the splitting amplitude. with It is important to note that Sp (1,E ) g→gg is zero due to color properties. In fact, we get where we have contracted the effective 2scalar-2gluon vertex with a factor f ade coming from the triple-gluon interaction.

A. Amplitude level results
As a first step, let's study the possible structure of the splitting matrix. In this process we have three physical momenta (p 1 , p 2 andP , or equivalently, n) and three physical on-shell polarizations vectors. Since external legs are massless on-shell particles we have the constraints P · (P ) = 0 = n · (P ) ⇒ p 12 · (P ) = 0 , where we have forced all the polarization vectors to vanish when contracted with the null-vector n, relying in the gauge invariance. So, we have the following non-zero scalar products: and where we are using p 1 · (P ) = −p 2 · (P ). Now we have to form all the possible structures that involve the three polarization vectors and that are compatible with the symmetry of the system when interchanging particles 1 and 2. Thus we get and notice that E + 2 is symmetric while E 1 ,E − 2 ,E 3 are antisymmetric. After replacing Feynman integrals in the expressions for Sp (1,i) g→gg and summing all the contributions, we realize that only two structures survive: E 1 + E − 2 (this is proportional to LO splitting) and E 1 − 2 s 12 E 3 . So, we can write for the NLO standard contribution, where the coefficients C (STD,i) g→gg are given by where we set α R = 1 to exclude HSA/HSB schemes.
Following Eq. 53, Sp (1) g→gg can be rewritten as H,g→gg + I C,g→gg p 1 , p 2 ;P Sp (0) g→gg , with I C,g→gg p 1 , p 2 ;P = c Γ g 2 as expected. Moving to the scalar-gluon contribution we obtain and comparing it with STD contributions in different schemes we get which agrees with the relation found for q → gq and g → qq splittings.

B. NLO corrections to AP kernel
Finally we can compute the contributions to the Altarelli-Parisi kernels. At LO we have for the polarized and unpolarized kernels, respectively. Note that when we set α = 1, then external gluons have 2 − 2 polarizations and they are treated like D ST -dimensional vectors. So, we must set δ = 1, which allows us to cancel the α dependence in the unpolarized kernel.
Moving to NLO, we obtain where α = 1 in CDR and α = 0 in FDH/HV schemes. It is worth noticing that the polarized kernel contains some terms proportional toP µ and n µ , but sinceP · (P ) = n · (P ) = 0 we neglect them to simplify the result.

VII. SPLITTINGS MATRICES INVOLVING PHOTONS
Let's consider an extension of massless QCD with the inclusion of a QED photon. This is a natural step when we want to study photon-production in the context of hadron colliders, since photons represent a very clean signal in the detector and QCD corrections can not be ignored.
This model can be described by extending the gauge group to SU(3) C × U(1) E which involves adding a new vector field A µ . The associated D-dimensional Lagrangian reads where {i, j} are color indices, g e is the electromagnetic coupling (i.e. the absolute value of electron charge), E Q is the charge of quark's flavor Q (E u,c,t = 2/3 and E d,s,b = −1/3) and F µν = ∂ µ A ν − ∂ ν A µ is the gauge-field strength tensor for the Abelian group U(1) E . From the interaction term, we can deduce that the Feynman rule for the quark-photon-quark vertex is −ıg e µ E Q γ µ and it is proportional to the identity matrix Id C in the color space. Since quarks belong to the fundamental representation of SU(3) C , then Tr(Id C ) = N C = C A which is going to be important when computing AP kernels.
In the next subsections, we show the associated splitting functions at NLO in the QCD coupling constant α s : Sp q→γq and Sp γ→qq . It is worth noticing that processes involving two photons and one gluon (i.e γ → γg or g → γγ) vanish due to color conservation, because they are proportional to Tr(T a (F )) = 0. On the other hand, there are not splittings with one photon and two gluons, because they involve a fermion loop with three vectors attached to it and, after summing all diagrams, we arrive to an expression which is again globally proportional to Tr(T a (F )) = 0.
It is worth noticing that we can check the divergent structure of splitting matrices involving photons using a formula similar to Eq. 53. For 1 → 2 processes, any splitting can be written as Sp (1) p 1 , p 2 ;P = Sp (1) H p 1 , p 2 ;P + I γ C p 1 , p 2 ;P Sp (0) p 1 , p 2 ;P , with Sp (1) H finite in the limit → 0 and containing only rational functions of p 1 , p 2 andP , and associated with the divergent contributions. Note that Eq. 164 is very similar to Eq. 54, with the exception of single pole proportional to b 0 . Explicitly, which is related to the fact that this kind of splitting only involves two colored-particles and we have to remove the single -pole coming from the renormalization of QCD coupling. Also, we have to take into account that C γ = 0 = γ γ because photons do not carry color.
Finally, since we are interested in studying the scheme dependence of splitting amplitudes, we are treating external photons and gluons in the same way. In other words, we can adapt the conventions shown in Section 2 for gluons to obtain phys.pol.
where (p) denotes the polarization vector associated to photons. The advantage of choosing this gauge is that it allows us to make a straightforward reduction from the pure QCD splittings, since this implies that (p) · n = 0 also for photons.

A. q → γq
This process can be considered as an Abelianization of q → gq, because it is not possible to have a triple-gluon vertex contribution. So, having performed a detailed study of q → gq in previous sections, we are able to extract some important results for q → γq without doing a full computation again.
First of all, the list of possible Feynman diagrams up to O (α 2 s ) is shown in Fig. 11. Note that they are essentially the same that we used for q → gq (see Figs. 4 and 5), except for the diagrams that include a triple-gluon vertex. At LO we have FIG. 11. Feynman diagrams associated with q(P ) → γ(p 1 )q(p 2 ) at LO and NLO. We include also the SCA-nHV contributions.
while the NLO standard-QCD corrections can be written as where the coefficients C (STD,i) q→γq are given by Analogously, for the NLO scalar-gluon contribution we have and it is straightforward to verify that which shows that the cancellation of scalar degrees of freedom occurs separately in Abelian and non-Abelian vertices.

Sp
(1) where we see that the divergent part (which contains -poles and branch-cuts) is isolated into I γ C , while Sp (1) H only contains rational functions and is finite in the limit → 0. Moreover, the new spinor chain which appears in the NLO computation is entirely contained in Sp (1) H . Finally, we can compute the corresponding contributions to the Altarelli-Parisi kernel. Since it is a quark initiated process, the polarized kernel verifies due to helicity conservation. So, the unpolarized kernel at LO is given by where we can appreciate that the result is independent of the number of fermion polarizations.
On the other hand, at NLO we have where α is a parameter that allows us to change between CDR (α = 1) and HV/FDH (α = 0) schemes.

B. γ → qq
Finally, we arrive to Sp γ→qq . Starting with g → qq, we have to replace the incoming leg with a photon, which forces us to eliminate self-energy correction (diagrams A and D in Fig. 7) and other term which includes a triple-gluon vertex. So, up to O (α 2 s ), we only have the diagrams shown in Fig. 12. The LO contribution reads while standard NLO correction is with On the other hand, for the NLO scalar-gluon contribution we have and, again, we find that the relation is fulfilled.
Testing the divergent structure of Sp (1) γ→qq we find that
Finally, the corresponding contributions to the Altarelli-Parisi kernels are for the LO terms and for the unpolarized NLO correction. The NLO polarized kernel can be expressed as because Sp (1) γ→qq is proportional to Sp (1) γ→qq .

VIII. CONCLUSIONS
In this work we have studied the double collinear limit and we computed the associated splitting matrices at NLO in α s for both pure QCD and QCD plus photon-quark interactions. As a first consistency check, we have verified that the divergent structure of splitting matrices agrees with the general form shown in the literature (for example, in Refs. [3,14,16]). Moreover, we found that the scheme dependence can be predicted up to O ( 0 ) using Eq. 11 in Ref. [35]. Also, we compared our results for usual DREG schemes with those available in Ref. [25], and, again, we found an agreement.
Besides the comparison of explicit results, we shown that FDH and HV schemes can be related at the amplitude level by introducing scalar-gluons. In fact, we verified that the relation is always fulfilled. Moreover, if we only consider fixed helicity configurations allowed by standard 4-dimensional QCD interactions, then we can extend the validity of Eq. 193 to include the CDR scheme. This is an important fact because it allows us to perform the same computation following two different paths, each of them having advantages in certain situations. For example, if we want to compute fixed-polarization splitting amplitudes (or matrices) in the CDR/HV scheme, it is more suitable to work with the r.h.s. of Eq. 193, because we settle D Dirac = 4 and many useful identities can be used. In particular, we can use Fierz identities to contract spinor chains and reduce them to bispinor products. The improvement in the treatment of results can be much better when more particles are involved (for instance, when studying the multiple-collinear limit).
On the other hand, if we want to compute Altarelli-Parisi kernel corrections, it is better to use the l.h.s. of Eq. 193 and work with D Dirac = 4 − 2 . The reason is that when we close spinor chains and sum over polarizations, we get rid of spinors and obtain traces which involves Dirac's matrices. Since the relations that we use to solve Dirac's traces are valid with any value of D Dirac , then we can simplify them and the final result only contains scalar products. Also, we do not have to compute each helicity configuration separately, which makes the computation straightforward.
This can be considered a great advantage, even if this procedure involves dealing with tensor type integrals which can have up to three free Lorentz indices. (In Appendix A we collect all the integrals required for the double collinear limit).
In the context of AP kernels, we also showed that it is possible to relate CDR and HV computations by just taking into account external scalar gluons. In fact, for the q → gq process, we find P CDR q→gq = P HV q→gq + P q→φq , which is a complement to Eq. 193 at the squared-amplitude level. Of course this relation can be extended to more general processes: we just have to decompose external gluons into 4-dimensional vectors plus scalar particles and compute each contribution separately.
Finally, let's make some comments about the alternative schemes studied in this article. In Section 2 we introduced some parameters that allowed us to control Dirac's algebra dimension (δ), the number of gluon polarizations (α R and α for internal and external particles, respectively) and the number of fermion polarizations (β R for internal fermions and β for external ones). By examining the behavior of Sp (1) q→gq with different parameter's values and comparing the divergent structure predicted by Eqs. 53 and 54, we conclude that the hybrid-schemes (i.e. α R = 0) are not consistent unless we include the corresponding scalar-gluon contributions. But, after adding these contributions, we get the same results provided by HV and CDR schemes. In other words, we show that the consistent version of HSA and HSB schemes are CDR and HV, respectively.
As anticipated in Section 2, FDH and TSC schemes are compatible with the supersymmetric Ward identity, even at one-loop level. In Ref. [36], it was shown that tree-level Altarelli-Parisi kernels computed in FDH and TSC schemes fulfilled this identity, i.e. P g→gg (z) + P g→qq (z) = P q→qg (1 − z) + P q→qg (z) , given that we set C A = C F = T R = N f . In this situation, if we identify quarks and gluinos then QCD is similar to N = 1 super Yang-Mills theory. From a physical point of view, this is possible because we consider the same number of bosonic and fermionic degrees of freedom. However, from Eqs. 103, 134 and 161 we can explicitly show that Eq. 195 is verified at one-loop level, for both FDH and TSC schemes. This result makes TSC an interesting choice, since it has a very symmetric and democratic way of treating all the particles involved in the computation.
It is interesting to appreciate that we performed the computations following a path that allowed us to keep track of Lorentz indices and metric tensors. In other words, we replaced integrals in Sp (1) before contracting with Sp (0) and summing over polarizations. This involved dealing with tensor-type integrals, which makes the calculation more complicated. If we were only interested in obtaining NLO corrections to AP-kernels, we could have first performed the contraction, and then replace the corresponding scalar integrals. However, scalar q 2 -integrals could appear in all schemes, with the exception of CDR. In spite of that, this approach is better suited when considering multiple-collinear splitting amplitudes, because tensor-type integrals become very lengthy and complicated when increasing the number of external physical momenta.
The natural next-step of this work is to extend the analysis to cover the multiple-collinear limit, and the possibility of computing them using recursion-relations [52], even at loop-level.
(228) APPENDIX B: PARTON SELF-ENERGIES Gluon self-energy When computing the gluon self-energy at one-loop level, we find that there are two Feynman diagrams which contribute to Π µν . They are shown in Fig. 13. Using conventional Feynman rules, we define × V Cin 3g (−p, q, p − q; µ, σ, ρ)V Cin 3g (−q, p, q − p; σ , ν, ρ ) , where we are using p as the external momenta which verifies p 2 = s. After integrating the loop-momentum we arrive to