Collinear and TMD Quark and Gluon Densities from Parton Branching Solution of QCD Evolution Equations

We study parton-branching solutions of QCD evolution equations and present a method to construct both collinear and transverse momentum dependent (TMD) parton densities from this approach. We work with next-to-leading-order (NLO) accuracy in the strong coupling. Using the unitarity picture in terms of resolvable and non-resolvable branchings, we analyze the role of the soft-gluon resolution scale in the evolution equations. For longitudinal momentum distributions, we find agreement of our numerical calculations with existing evolution programs at the level of better than 1 percent over a range of five orders of magnitude both in evolution scale and in longitudinal momentum fraction. We make predictions for the evolution of transverse momentum distributions. We perform fits to the high-precision deep inelastic scattering (DIS) structure function measurements, and we present a set of NLO TMD distributions based on the parton branching approach.


Introduction
Realistic comparisons of experimental data for hard processes at high-energy hadron colliders with theoretical predictions based on QCD factorization formulas [1][2][3][4] require Monte Carlo simulation via parton shower event generators. While great progress has been achieved in the last decade on matching and merging methods [5][6][7] to combine parton showers with perturbative calculations through next-to-leading order (NLO), important open questions still remain, both conceptual and technical, on the appropriate use of parton distribution functions in parton showers (see e.g. [8][9][10][11]) and on the treatment of the shower's transverse momentum kinematics (see e.g. [12][13][14]). The relevance of these effects is known to increase with energy [15][16][17], and they thus constitute an important theme for physics at the Large Hadron Collider (LHC) and at colliders of the next generation [18].
Candidate approaches to tackle such questions in complex, multi-scale collider processes generally include theoretical constructs designed to extend the concept of collinear parton density and decay functions, as in the case of Soft Collinear Effective Theory (SCET) [19][20][21] or of transverse momentum dependent (TMD) formalisms [22][23][24]. For instance, in Ref. [25] the TMD gluon density is determined, based on the high-energy factorization [26] and CCFM evolution equation [27], from fits to the high-precision deep inelastic scattering (DIS) data [28,29], and used in a parton shower calculation [17] to make predictions for W -boson + jets hadro-production, which can be compared with LHC experimental measurements [30,31].
Although the analysis in [17,25] proves to be successful in achieving a meaningful TMD description of both DIS and Drell-Yan measurements, it is based on a TMD form of factorization valid at high energy, and requires a matching method (provided by the CCFM equation) to include nonasymptotic contributions at collider energies. Because of the high-energy expansion, the method is predominantly sensitive to the gluon density, and quark contributions enter systematically at subleading orders.
In Ref. [32] we have proposed a different approach, based on solving coupled quark and gluon DGLAP [33] evolution equations using parton branching methods, and determining from this both collinear (integrated over transverse momenta, iTMD) and TMD parton densities. Rather than starting from high-energy resummed equations, this work relies on renormalization group evolution equations. The approach uses the unitarity formulation of these evolution equations which forms the basis of parton showering Monte Carlo simulation [2,34]. In this sense, it is close in spirit to the works in Refs. [35][36][37][38][39][40][41][42][43], in Refs. [44][45][46][47][48][49][50], and in Refs. [51,52]. Ref. [32] shows that the evolution of parton distribution functions can be calculated, including the transverse momentum dependence, from a parton branching approach, provided infrared contributions to evolution are treated by a method which takes into account consistently soft gluon emissions, near the endpoint for lightcone momentum fractions z → 1, not just at inclusive level but at exclusive level. In particular, Ref. [32] shows that this can be done by using a finite soft-gluon resolution scale in the evolution equations.
In this paper we provide details of the approach set out in [32], we show that it can be applied to higher accuracy order-by-order in the strong coupling α s , and we present numerical results at the next-to-leading order (NLO). Further we present the results of fits, based on this approach, to the high-precision DIS data [53]. First results from this work have appeared in [54].
The paper is organized as follows. In Sec. 2 we describe the main elements of the partonbranching formulation of the coupled QCD evolution equations. In Sec. 3 we present the numerical Monte Carlo solution of the coupled quark and gluon evolution equations at NLO. We compare collinear parton density functions obtained by our parton-branching solution with results obtained via the evolution package Qcdnum [55][56][57]. In Sec. 4 we illustrate an application of our method by performing a fit to the high-precision DIS data [53]. For the fit we use an updated version of the program [58] within the xFitter open-source QCD fit platform [59]. By the method of the present paper we are able to extend the fit [25] to precision DIS data significantly toward higher x and higher Q 2 . In Sec. 5 we turn to TMD parton density functions [22,60]. We discuss the identification of the transverse momentum in the initial-state parton distribution in terms of the shower's kinematics and evolution variable. We present a new set of quark and gluon TMDs including NLO evolution kernels. We give conclusions in Sec. 6.

Unitarity approach to QCD evolution equations
In this section we give the main elements of the parton-branching approach to the evolution equations. We introduce a soft-gluon resolution scale into the renormalization group evolution equations, and describe resolvable and non-resolvable emissions. We discuss the relationship of our results with the angular-ordered, coherent branching [61][62][63][64] and the behavior of the endpoint z → 1 region in transverse momentum distributions [65,66]. We construct an iterative Monte Carlo solution of the evolution equations, and apply it to the case of collinear and TMD parton densities.

The renormalization group evolution
The renormalization group evolution of parton distribution functions can be written in terms of parton splitting processes as follows where f a (x, µ 2 ) are parton distributions for a = 1, . . . , 2N f + 1 species of partons (with N f the number of quark flavors) as functions of longitudinal momentum fraction x and evolution mass scale µ, and P ab (α s , z) are the DGLAP splitting functions, depending on the running coupling α s and the splitting variable z, and computable as a perturbation series expansion We will work with the momentum-weighted parton distribution functions f a , for which the evolution equations read In the physical picture of Eq. (4), a finite resolution scale in the transverse distance between emitted partons implies, by energy-momentum conservation, that one cannot resolve partons radiated with longitudinal momentum fractions closer to z = 1 than a certain cut-off value, z > z M with 1−z M ∼ O(Λ QCD /µ), where µ is of the order of the hard-scattering scale and Λ QCD ≈ 1 fm −1 is the natural scale of the strong interactions. Removing non-resolvable radiative contributions from the evolution, on the other hand, leads to a violation of unitarity. The key idea of the parton branching method is to restore unitarity by recasting the evolution equations in terms of no-branching probabilities (Sudakov form factors) and real-emission branching probabilities. We will introduce the resolution scale parameter z M formally into the evolution equations in Sec. 2.3, and describe the unitary branching method in the subsequent sections.
To set up our formalism, we decompose the splitting functions P ab (α s , z) as where the plus-distribution 1/(1 − z) + is defined for any test function ϕ as Eq. (5) provides a classification of the singular behavior of the splitting functions P ab (α s , z) in the non-resolvable radiation region z → 1. It decomposes the splitting functions into the δ(1 − z) distribution, the 1/(1 − z) + distribution, and the function R(α s , z) which contains logarithmic terms in ln(1 − z) and analytic terms for z → 1. The δ(1 − z) and 1/(1 − z) + contributions to splitting functions are diagonal in flavor, (no summation over repeated indices). The functions D ab and K ab , or equivalently d a and k a , and the functions R ab in Eq. (5) have the perturbation series expansions The treatment which we develop in this section only relies on the decomposition in Eq. (5) and is valid at any order in α s . In practical applications one takes a given truncation of the expansions in Eqs. (8), (9). The numerical results in Secs. 3, 4 and 5 are based on the expansion to NLO (i.e., n = 2 in Eqs. (8), (9)).
Charge conjugation and SU (N f ) flavor symmetries imply that the splitting functions P ab obey the following relations to all orders, where the superscripts N S and S stand respectively for non-singlet and singlet. Therefore, P ab has three independent quark-gluon or gluon-gluon components (P qg , P gq and P gg ) and four independent quark-quark components (the N S components P N S qq , P N S qq and the S components P S qq , P S qq ). † In the next section we give explicit expressions at one-loop and two-loop orders for the D ab , K ab and R ab terms in Eq. (5).

Expansion in powers of α s
At one-loop order the coefficients of the perturbative expansions (8),(9) for d a , k a and R ab can be read from [33]. At this order, one has P N S qq = P S qq = P S qq = 0, so that all quark-quark components are degenerate. The one-loop expressions for d a , k a and R ab are given by and where the SU (N c ) color factors (with N c = 3 the number of colors) are given by At two-loop order the perturbative coefficients for d a , k a and R ab start to depend on the renormalization scheme. In the MS scheme the results can be read from [68,69]. At the level of two loops also P N S qq , P S qq , P S qq are nonvanishing so that the degeneracy of the quark splitting functions is lifted. However a residual degeneracy remains because P S qq = P S qq at this order. ‡ The two-loop coefficients for d a and k a are given by † The independent quark-quark components can alternatively be taken [1,67,68] to be the three which correspond to the three linear combinations diagonalizing the evolution of non-singlet distributions, plus the one which controls the evolution of the singlet quark distribution coupled to gluons. ‡ The degeneracy is fully lifted starting at three-loop order [70,71].
where ζ is the Riemann zeta function, and The expressions for the two-loop coefficients for the functions R ab are lengthier, and are given in Appendix A.
We will use these expansions through two loops for the numerical calculations in Secs. 3, 4, 5.

Resolvable and non-resolvable emissions
We now introduce the soft-gluon resolution parameter z M into the evolution equations (4), by splitting the integration range on the right hand side into the resolvable (z < z M ) and non-resolvable In each region, we use the decomposition (5) in the evolution equations. We include terms through Consider first the endpoint z → 1 contribution from the K ab term in Eq. (5). Using Eq. (6), we rewrite this as In the region 1 > z > z M we expand the momentum-weighted parton density as Then we see that the contribution to Eq. (17) from the non-resolvable region is of order Next, we consider the contributions to the evolution equations (4) from the other two terms, D ab and R ab , in Eq. (5). The R ab contribution can be combined with the first term on the right hand side of Eq. (19) to yield a contribution to the evolution proportional to f b (x/z, µ 2 ). The D ab contribution can be combined with the second term on the right hand side of Eq. (19), using the δ(1 − z), to yield a contribution to the evolution proportional to f b (x, µ 2 ). Further, we use that R ab has no power divergences (1 − z) −n and is at most logarithmic for z → 1, so that the integration over R ab for z > z M gives O(1 − z M ). Thus, we can write The first line in Eq. (20) contains contributions to evolution from real parton emission, while the second line contains contributions from virtual corrections. It is convenient to define the kernels in the bracket of the first line as the real-emission branching probabilities P That is, the real-emission branching probabilities P (R) ab (α s , z) are obtained from the splitting functions P ab (α s , z) in Eq. (5) by subtracting the δ(1 − z) terms and replacing the plus-distribution The virtual terms in the second line of Eq. (22) can be dealt with by using the momentum sum rule, as we see next.

Momentum sum rule
We will now use the momentum sum rule to systematically eliminate the D-terms in Eq. (5) from the evolution equations in favor of the K-terms and R-terms. To this end, we insert the momentum sum rule into the evolution equations, by subtracting the momentum sum integral in the curly bracket in the second line of Eq. (22). Recall from Eq. (7) that the D ab and K ab terms in this equation are diagonal in flavor. Therefore, by interchanging indices, we obtain from Eq.
Let us now use again the decomposition (5) for P ca (α s (µ 2 ), z) in the last line of Eq. (24). We observe that the D ca term in P ca (α s (µ 2 ), z) cancels against the first term in the curly bracket in Eq. (24), while the R ca term in P ca (α s (µ 2 ), z) may be restricted to the region z < z M , up to order O(1 − z M ). Finally, the K ca term in P ca (α s (µ 2 ), z) may be combined with the second term in the curly bracket in Eq. (24). Putting pieces together, we get We thus recognize, by using Eq. (21), that the evolution equations (22) can be written as

Sudakov form factor
Eq. (26) recasts the evolution of each parton a in terms of the real-emission probabilities P (R) ab and P (R) ba and of the resolution parameter z M . It can be rewritten in a form which has the advantage of being solvable by an iterative Monte Carlo procedure if we introduce the Sudakov form factor, defined as The Sudakov form factor ∆ a (z M , µ 2 , µ 2 0 ) has the interpretation of probability for parton a to undergo no branching between evolution scale µ 0 and evolution scale µ, where the branchings are understood to be classified according to the given resolution z M .
Noting that we obtain from Eq. (26) (removing z M and µ 2 0 from the argument list for better readability) This evolution equation can be written in a form similar to Eq. (4), but now in terms of realemission probabilities P (R) ab and Sudakov form factors: Integrating this equation we obtain, with ∆ a (µ 2 0 ) = 1, We recognize that introducing the Sudakov form factor has led to an equation which is an integral equation of Fredholm type, This can be solved by iteration as a series [32] f where We observe that while at LO in α s the splitting functions are positive definite, this is no longer the case at NLO. However, although the integrands can be negative, the integrals over the splitting functions which appear in the evolution kernels and Sudakov form factors remain positive also at NLO. We will exploit this in the next section to apply a Monte Carlo method for solving the evolution equations.

Solution of the evolution equation applying a Monte Carlo method
The solution of the evolution equation can be obtained by applying a Monte Carlo method. By this method the problem is reduced to that of generating the splitting variable z and the evolution scale µ.  1 depicts the parton evolution: we start with a parton a and evolve from scale µ i to scale µ either without any branching, or having one branching at scale µ i+1 , or having a second branching at scale µ i+2 , and so on. The probability to evolve from µ i to µ i+1 without any resolvable branching is provided by the Sudakov form factor ∆ a (z M , µ 2 i+1 , µ 2 i ). Figure 1: Illustration of the evolution process by iteration: a parton can evolve from scale µ i to scale µ without any branching (left), having one branching (middle), two branchings (right) and so on. The relevant variables are indicated.
By introducing a random number R 0 in [0, 1], we generate the value µ i+1 by solving Eq. (27) for µ i+1 at a given µ i , that is, where the upper bound may be taken to be µ max → ∞, leading to the simple expression The splitting variables z are generated from where R 1 is a random number in [0, 1], z M is the resolution parameter, and z min is the lowest kinematically allowed value.
Generating a pair of z i , µ i values many times, we obtain a true and unbiased estimate of the integrals, and a solution of the evolution equations.
We have implemented the Monte Carlo method to solve the evolution equations in a numerical program. The program is a development of the code [58] which was earlier employed by some of us for studies of the CCFM equations [25]. It is worth observing that the application to the case of the evolution equations studied in this paper presents different features with respect to the case of the CCFM equations. The differences involve especially the flavor structure of the two equations, and the behavior of the kernels at small longitudinal momentum fractions. While CCFM equations are dominated by the gluon channel, Eq. (31) has fully coupled flavor structure. The small-x behavior of CCFM kernels is controlled by the non-Sudakov form factor [27]. In the case of Eq. (31) it is essential to work with momentum-weighted distributions to improve the convergence of the numerical integration over the region of small x. In Secs. 3, 4, 5 we will employ this program to compute numerical results.

Transverse momentum distributions and ordering variables
The use of the Sudakov form factors and the iterative method of the previous subsection allows one to generate step by step each resolvable branching. In addition to providing the solution of the evolution equation, this approach has the advantage of keeping track of detailed information about each individual branching. For example, the kinematics in each branching can be calculated, similarly to what is done in a parton shower process. In particular parton distributions can be obtained, not only depending on x and µ (as in f a (x, µ 2 )), but also depending on the transverse momentum k ⊥ of the propagating parton (as in TMD parton distributions A a (x, k ⊥ , µ)).
x a p + , k t,a q t,c → µ b Figure 2: Branching process b → a + c.
Consider the splitting process b → a + c in Fig. 2. Using the notation in the figure, with plus light-cone momenta p + a = zp + b , p + c = (1 − z)p + b , we have, by applying conservation of minus light-cone momentum, where q c is the (euclidean) transverse momentum vector of particle c. For a space-like branch-ing [73][74][75] with µ 2 = −p 2 a and p 2 b = p 2 c = 0, taking z → 0 in the high-energy limit gives Eq. (39) is referred to as transverse momentum ordering. If, on the other hand, the evolution variable µ is associated with the angle Θ of the momentum of particle c with respect to the beam direction, we obtain the angular ordering relation [34,64,74] The transverse momentum of the propagating parton is calculated as The method thus enables one to determine the corresponding transverse momentum dependent (TMD) parton distribution A a (x, k, µ 2 ), in addition to the inclusive distribution f a (x, µ 2 ), integrated over k, It has been pointed out in [32] that the transverse momentum generated radiatively by the recoils in the evolution cascade depends strongly on the treatment of the non-resolvable region z → 1. Unlike the integrated distribution f a (x, µ 2 ), the TMD distribution A a (x, k, µ 2 ) is infrared-sensitive. The origin of this behavior lies with singularities, present at fixed k, which arise from branching processes in Fig. 2 with gluons at large negative rapidities, y ∼ ln q + /q − → −∞ [66]. While in the integrated distribution such singularities cancel between real and virtual non-resolvable emissions, this is not, in general, the case for the TMD distribution. As a result, supplementary conditions are needed to define the TMD distribution consistently [66].
In the framework of the parton branching solution of evolution equations discussed in the present paper, one such set of conditions is provided by the angular ordering in Eq. (40). By using the angular-ordered branching, consistent TMD distributions are defined, which are independent of the soft-gluon resolution parameter z M , for sufficiently large values of z M . In contrast, the transverse momentum ordered branching, based on Eq. (39), while entirely suitable as long as one is working at the level of integrated parton distributions, does not allow one to define TMD distributions consistently, as the transverse momentum at any given evolution scale would depend on the choice of the resolution parameter z M . In [32] the above observation is made by working at LO in the strong coupling α s . In Sec. 5 of the present paper, we confirm and extend these findings to NLO.
Using Eq. (31) and the angular ordering (40), we write the branching equation for the evolution of TMD distributions as where A is the momentum weighted distribution A ≡ xA. By applying the method in Sec. 2.6, we solve this iteratively as where and so forth.
From the solution of the branching equations (31) and (43) we will obtain collinear and TMD parton distributions. The behavior at small transverse momenta, in particular, is controlled in this formulation by the nonperturbative distributions at scale µ 0 and by Sudakov form factors, which embody the perturbative resummation and are defined as functions of the soft-gluon resolution scale separating resolvable and non-resolvable branchings. According to Eqs. (27), (31), (43) this is expressed in terms of integrals over the branching probabilities P (R) (α s , z) in Eq. (21). It depends on the residues K ab at the poles z = 1 in Eq. (21), given at one-loop and two-loop orders by the coefficients in Eqs. (12) and (16). The relationship of the z = 1 behavior with small transverse momenta is important to construct reliable theoretical predictions for transverse momentum q ⊥ spectra at low q ⊥ in the production of massive states in hadronic collisions [4,22,[76][77][78][79][80].
We next move on to numerical results based on the methods described in this section.

Numerical parton-branching solution at NLO
In this section we present numerical results from the parton-branching solution at NLO. We compare the answer thus obtained for the collinear parton density functions with the answer from the evolution package Qcdnum [55]. The corresponding comparison at LO has been shown in [32].
For the purpose of this comparison we use as input parton distributions the distributions which are the default set given in Qcdnum. These distributions are parameterized at the starting scale µ 2 0 = 2 GeV 2 as 6 , xs(x) = 0.2 (xd(x) + xū(x)), xs(x) = xs(x), xg(x) = 1.7 Three light active flavors are assumed at the starting scale, while charm and bottom quarks are produced during the evolution, for evolution scales µ > m c = 1.73 GeV and µ > m b = 5.0 GeV respectively. The running coupling α s is used at two loops, with α s (m 2 Z ) = 0.118.
In Fig. 3 we show the momentum-weighted parton densities (integrated over transverse momenta, iTMD) obtained from the parton branching solution of the evolution equations, compared to the predictions obtained by Qcdnum, starting from µ 2 0 = 2 GeV 2 for different scales µ 2 = 10, 10 3 , 10 5 GeV 2 . The parton branching solution is obtained for a fixed value of z M = 1 − 10 −5 . ¶ The overall agreement between the parton-branching and Qcdnum calculations is better than 1% throughout the whole range in x and µ 2 .    observed, confirming the findings of [32], and extending them to NLO in accord with the formalism developed in Sec. 2.

Fit to precision DIS data
The initial parton density distributions have to be determined from fits to experimental data. A general tool to perform such fits to collider measurements is the xFitter package [59]. As an application of our formalism, in this section we describe the method and results of a fit to the precision DIS data [53].
To start with, we implement the parton branching solution of the evolution equation in the xFitter package [59]. However, using a full Monte Carlo solution of the evolution equation for every new set of initial parameters would be too time-consuming to be efficient. Instead, we employ the approach developed in [25,58]. Following this approach, first a kernel A b a (x , µ 2 ) is determined from the Monte Carlo solution of the evolution equation for any initial parton of flavor b evolving to a final parton of flavor a ; then, this is folded with the non-perturbative starting distribution A 0,b (x): The kernel A b a includes the full parton evolution as in Eq. (31), with Sudakov form factors and splitting functions, and is determined with the parton branching method described earlier. The kernel A can be determined as a function of x, µ for the k ⊥ -integrated iTMD distributions, or depending on x, k ⊥ , µ for the transverse momentum dependent (TMD) distributions. The kernel is then folded with the initial condition be then used within the xFitter package to calculate the cross sections and to determine the parameters of the starting distributions A 0,b (x). We use precision measurements, in neutralcurrent and charged-current interactions at various beam energies from HERA 1+2 [53], of the reduced cross section in the range 3.5 < Q 2 < 30000 GeV 2 .
In practice, since the initial state partons can be only light quarks or gluons, it is enough to determine the kernel A b a only for one initial state quark and a gluon.
In Fig. 5 the calculated cross section σ red , obtained from the fit using xFitter, is compared with the precision measurements from HERA [53] for different values of Q 2 , showing very good agreement from low to high values of Q 2 .
Comparing this result with the fit [25] to precision DIS data based on CCFM evolution equations, note that in the case [25] the constraint x < 0.005 is applied on the data set, while no x constraint is applied in the present case. By the approach of this paper the description of precision DIS measurements can be significantly extended toward higher x and thus higher Q 2 , without on the other hand introducing any extra constraint cutting out lower x data.
We plan to analyze fits to data further in a future work.

TMD densities
By applying the parton branching method of this paper, we are able to construct TMD parton densities as described in Sec. 2.7. While large transverse momenta are generated by perturbative evolution, the nonperturbative region of small k ⊥ cannot be predicted in our approach but is parameterized by nonperturbative distributions which are to be determined from experimental measurements. For the calculations of this section we use the parameterizations given in Sec. 3 and take simple gaussian distributions exp(−|k 2 ⊥ |/σ 2 ). The widths σ are in general flavor-dependent. For the numerical illustrations that follow we take the same width for all parton species, σ 2 = q 2 0 /2 for all flavors with q 0 = 0.5 GeV. We present numerical results for the evolution of TMD parton densities, using the parton branching solution of the evolution equations with NLO evolution kernels.
It has been shown in [32] that the TMD distributions, unlike the collinear distributions, are strongly influenced by the ordering variable in the branching. In particular, the cases of the transverse momentum ordering (39) and angular ordering (40) have been examined in [32] by an explicit calculation, working at LO in the strong coupling α s . We here confirm and extend this analysis, working at NLO. We illustrate that the same behavior found at LO applies at NLO as well.
In Figs. 6 and 7 we apply the NLO numerical solution of Sec. 3 and the method of Sec. 2.7 to study the longitudinal and transverse momentum dependence of the gluon distribution, and its behavior with the soft-gluon resolution parameter z M . Fig. 6 shows the TMD gluon distribution versus the longitudinal momentum fraction x for different values of the resolution parameter, 1 − z M = 10 −3 , 10 −5 , 10 −8 . The curves are plotted for a fixed value of transverse momentum k t ≡ |k| = 10 GeV, and two values of evolution scale, µ = 100 GeV (top panels) and µ = 1000 GeV (bottom panels). * * On the right are the results for transverse-momentum ordering; on the left are the results for angular ordering. The transverse-momentum ordering, due to the effect of the non-resolvable region, does not lead to results independent of z M . On the other hand, the angular ordering correctly takes into account the cancellation of non-resolvable emissions due to soft-gluon coherence [64], and no dependence is left on the resolution parameter z M . Fig. 7 shows the TMD gluon distribution versus the transverse momentum k t , at fixed x = 10 −2 , * * The plots in Figs. 6 and 7 are produced using the plotting tool TMDplotter [60,81]. for different values of the resolution parameter, 1 − z M = 10 −3 , 10 −5 , 10 −8 . As in Fig. 6, on the right are the results for transverse-momentum ordering, and on the left are the results for angular ordering. The plots in Fig. 7 illustrate as a function of k t the same effect of the ordering and behavior in the resolution parameter z M which we have seen in the previous figure.
Analogous behavior to that in Figs. 6 and 7 was observed at LO in [32]. Figs. 6 and 7 show that the z M dependence in the transverse momentum ordering case cannot be avoided or reduced by inclusion of NLO evolution. It means that the different orderings in Eqs. (39), (40) should not be thought of as different factorization schemes, and the results in the two cases will not be related by a change in the factorization scheme.
It is worth noting that the transverse momentum ordering (39) is widely used in a variety of contexts, e.g. in low-x physics studies, as it results from approximating the exact parton-splitting kinematics in the region of strongly ordered momentum fractions z 1. The results in Figs. 6 and 7 imply that this approximation will not be valid for observables sensitive to the detailed structure in transverse momentum of the initial state. In particular, they emphasize the role of soft-gluon coherence effects leading to angular ordering in constructing well-defined TMD distributions even at low x [15,16,82,83].
In Fig. 8 we illustrate the flavor decomposition, at TMD level, resulting from perturbative evolution. We plot the TMD distributions obtained from the parton branching method for different flavors, by applying the evolution with appropriate angular-ordering condition.
To summarize, in this section we have shown that a consistent set of TMD parton distributions, valid over a large range in x, k ⊥ and µ can be determined from a parton branching solution of QCD evolution equations, as long as the soft gluon region is treated appropriately, e.g. by applying angular ordering conditions. We have shown that under these conditions the dependence on the resolution scale parameter z M drops out also for the k ⊥ -distribution, provided z M is large enough, resulting in stable predictions for the evolution of TMD parton densities.

Conclusions
Motivated by both conceptual and technical questions on the treatment of initial-state kinematics and distribution functions in QCD parton-shower calculations, we have investigated partonbranching solutions to QCD evolution equations. We have presented results of constructing collinear and TMD parton densities from this approach at NLO.
By separating resolvable and non-resolvable branchings, and analyzing the role of the soft-gluon resolution scale in the evolution, we have proposed a method to take into account simultaneously soft-gluon emission in the region z → 1 and transverse momentum q ⊥ recoils in the parton branchings along the QCD cascade.
This approach is potentially relevant for calculations both in collinear factorization and in transverse momentum dependent factorization. The starting point of the approach are the DGLAP equations, which are evolution equations valid for fully inclusive distributions. The method developed in this paper provides the branching equations which apply at exclusive level. Unlike DGLAP equations, these are necessarily sensitive to soft-gluon emission in the infrared region. We have presented the evolution equations as a function of the soft-gluon resolution scale and the ordering condition.
The branching equations for TMD densities obtained in this paper can be compared with existing TMD evolution equations: the CSS equations [4,80], which apply in the low transverse momentum region q ⊥ Q (where Q is the high mass scale of the hard scattering) and can be used in the case of low-q ⊥ TMD factorization [79]; and the CCFM equations [27], which apply in the high-energy region √ s Q, and can be used in the case of high-energy TMD factorization [26]. We have pointed to a few of the main differences and similarities in the physics described by these different approaches.
The CSS and CCFM approaches are designed to achieve high logarithmic accuracy in the resummation of higher-order logarithmic contributions in the restricted phase space regions which identify their domains of validity, respectively q ⊥ → 0 and √ s → ∞. In such approaches matching methods are required to go beyond these restricted domains and arrive at predictions valid more generally (e.g., the Y -term matching for high q ⊥ in CSS, and the large-x terms in CCFM splitting functions). On the other hand, the spirit of the approach proposed in this paper is to provide TMD distributions which can be applied over a broad kinematic range from low to high energies, and from low to high q ⊥ . We incorporate consistently renormalization group evolution, soft-gluon coherence and parton branching kinematics. The approach is general and, although in this paper we focus on longitudinal splitting functions, we believe it can accommodate further dynamical effects such as transverse splitting functions [23,66].
The formalism of this work implies that the soft-gluon resolution scale z M depends on the evolution variable µ. In the numerical examples of this paper we have limited ourselves to considering fixed values of z M . One of the main directions of development of this approach will concern the µ dependence of z M .
Furthermore, we observe that, while power-suppressed contributions of order O(1 − z M ) n ∼ O(Λ QCD /µ) n are beyond the scope of the treatment given in this paper, logarithmically enhanced contributions in ln n (1 − z M ) could be taken into account, and could be related [64] with threshold logarithms in production cross sections coupled to the parton distributions. We regard this as a further potential advantage of the formalism of this work.
Given the results presented in this paper for the parton evolution including the full flavor structure, we expect this approach to have a wide range of applications both at low energies and at high energies.