Multi-jet Merging with TMD Parton Branching

One of the main theoretical systematics in studies of final states with large jet multiplicities at high-energy hadron colliders is associated with the merging of QCD parton showers and hard-scattering matrix elements. We present a method to incorporate the physics of transverse momentum recoils due to initial-state shower evolution into multi-jet merging algorithms by using the concept of transverse momentum dependent (TMD) distributions and the associated parton branching. We investigate the dependence on the merging scale and illustrate the impact of the new method at the level of both exclusive and inclusive final-state observables by studying differential jet rates, transverse momentum spectra and multiplicity distributions, using vector boson + jets events at the LHC as a case study.


Introduction
Events with large multiplicities of hadronic jets in high-energy pp collisions are among the most spectacular final states observed at the Large Hadron Collider (LHC), and will continue to be explored in future collider experiments [1,2], at ever increasing jet multiplicity and with increasingly high statistics. The accurate theoretical description of multi-jet events has great practical importance, as multi-jets are used to precisely determine the properties of particles decaying to them (e.g. gauge or Higgs bosons and the top quark), to perform precision tests of strong-interaction dynamics in the Quantum Chromodynamics (QCD) sector of the Standard Model (SM), and to search for signatures of new physics beyond the Standard Model (BSM).
Theoretical predictions for multi-jet observables have relied for the last twenty years on "merging" techniques [3][4][5][6][7][8][9][10][11][12], to combine matrix-element and parton-shower event generators. The former describe the underlying hard process with bare partons providing the primary sources for widely separated jets, the latter describe the evolution of partons by radiative processes predominantly at small angles, and the two are sewn together, so as to avoid either double counting or missing events, via a "merging scheme" and merging scale.
The primary goal of the merging scheme is to ensure that the best possible approximation (fixed-order exact matrix elements or the parton shower evolution) is used at any given point of phase-space. The scheme and its parameters select the approximation to be used, and correct the event weight accordingly. Since neither of the descriptions (matrix element or parton shower) is exact, the choice of merging parameters introduces an important theoretical systematic uncertainty, discussed at leading order (LO) in [6][7][8][13][14][15][16] and next-to-leading order (NLO) in [9][10][11][12]. This systematic uncertainty reflects the mismatch between the matrix-element and the parton-shower weights assigned to a given final state. The larger the mismatch, the larger the uncertainty. The phase-space regions that are mostly affected are those describing final states for which the jet multiplicity can vary under minor changes of the merging parameters. Typically, this happens if a jet is soft or close to another hard jet. A better modeling of the emission probability for such jets by the parton-shower evolution would reduce the difference with the weight assigned to these events by the matrix-element description, reducing the mismatch and the relative systematic uncertainty.
Motivated by this, we explored in an earlier proposal [17] the role of parton-shower algorithms that go beyond the approximation of small-angle, collinear emissions, by taking into account transverse-momentum recoils in the initial-state shower through the use of transverse momentum dependent (TMD) [18] parton distributions. The reason for this is twofold. On one hand, multi-jets probe phase-space regions characterized by multiple momentum scales, in which TMDs control perturbative QCD resummations [19] to all orders in the strong coupling (e.g., in the Sudakov region [20] and high-energy region [21]). On the other hand, transverse momentum recoils can influence the theoretical uncertainty associated with combining matrix-element and parton-shower contributions [22,23], and thus affect the dependence of multi-jet cross sections on the merging scale. In this paper we further elaborate on the results of Ref. [17], providing a more detailed account of the algorithm and presenting additional studies to demonstrate its performance.
While advanced formalisms have been developed to describe TMD effects in inclusive observables (see e.g. the recent next-to-next-to-next-to-leading-order (N 3 LO) calculations [24,25] of TMD coefficient functions, and third-order studies [26,27] of transverse momentum Drell-Yan (DY) spectrum at the LHC), TMD implications on the exclusive structure of final states with high jet multiplicity have only just begun to be explored in Ref. [17], where we introduced a systematic merging method to analyze such effects in multi-jets. We aim at obtaining results for both inclusive and exclusive multi-jet observables, accounting simultaneously for the contribution of different jet multiplicities in the final state and for the contribution of initial-state TMD evolution. To this end, we employ the parton branching (PB) formulation of TMD evolution set out in [28]. Any of the existing merging algorithms can then be used, in principle, to combine samples of different parton multiplicity showered using the TMD parton branching. In this work, we focus however on the so-called MLM matching prescription [5][6][7][8]. First results from this investigation have been presented in [29].
The paper is organized as follows. We start in Sec. 2 by briefly recalling the main features of the PB formulation of TMD evolution. In Sec. 3 we describe the TMD multi-jet merging method and illustrate its features by presenting the results of applying the method to differential jet rates. In Sec. 4 we perform a detailed analysis of Z-boson production associated with multi-jets at the LHC, computing theoretical predictions for Z+jets observables with the TMD jet merging method and comparing them with LHC experimental data. In Sec. 5 we examine the theoretical systematic uncertainty associated with the merging parameters and the dependence of theoretical predictions on the merging scale. In Sec. 6 we study the sensitivity of our results to final-state parton showers, and present a comparison with Pythia Monte Carlo results based on collinear multi-jet merging. In Sec. 7 we study the sensitivity of our results to nonperturbative effects in TMD and collinear parton distribution functions (PDFs). Sec. 8 includes remarks on possible off-shell TMD effects in the evaluation of partonic matrix elements. We give our conclusions in Sec. 9.

Parton branching formulation of TMD evolution
In this section we summarize the basic elements of the PB formulation of TMD evolution. First we recall the evolution equations; next we discuss the physical picture of transverse momentum broadening from TMD evolution, and its implications for the production of multiple jets.

Evolution equation
The PB method [28] uses the unitarity picture of parton evolution [30], commonly used in showering algorithms [31,32], for both collinear and TMD parton distributions. Soft gluon emission and transverse momentum recoils are treated by introducing the soft-gluon resolution scale z M [33] to separate resolvable and non-resolvable branchings. Sudakov form factors ∆ j are used to express the probability for non-resolvable branching in a given evolution interval from one scale µ 0 to a higher scale µ, where j, are flavor indices, z is the longitudinal momentum fraction, α s is the strong coupling, and P (R) kj are resolvable splitting functions, computable as power series expansions in α s .
In this approach the TMD evolution equations can be written as is the TMD distribution of flavor j carrying the longitudinal momentum fraction x of the hadron's momentum and transverse momentum k at the evolution scale µ, µ 0 is the initial evolution scale, and µ = µ 2 is the momentum scale at which the branching occurs. The branching probabilities are evaluated using the form of the strong coupling according to angular ordering, with the scale in α s given by the emitted transverse momentum at the branching, q T = (1 − z)µ. It was observed in [28,33] that angular ordering [30,32,34,35] leads to well-prescribed TMDs, i.e., stable with respect to variations of the soft-gluon resolution scale z M . In this work we will follow this observation and will consider angular-ordered evolution.
Collinear limits may be obtained by integrating Eq. (2.2) over all transverse momenta. For z M → 1 and α s → α s (µ 2 ), the convergence to collinear PDFs satisfying DGLAP evolution equations [36][37][38] has been verified numerically [33] at LO and NLO against the evolution program [39] at the level of better than 1% over a range of five orders of magnitude both in x and in µ.
Besides the collinear limits, Eq. (2.2) can be used at unintegrated level for Monte Carlo event simulation of TMD effects [40]. The PB method enables the explicit calculation of the kinematics at every branching vertex once the evolution scale is specified in terms of kinematic variables. If the TMD distribution A j (x, k 2 , µ 2 ) evaluated at the scale µ 2 is known, the corresponding TMD parton shower can be generated by backward evolution, starting from the hard scattering process.
The basic elements of the initial state shower in the backward evolution scheme are described in Ref. [40]. By a method analogous to that used in the case of collinear showers [41][42][43], the Sudakov form factor for backward evolution is obtained from the solution of Eq. (2.2). This Sudakov factor depends on the ratio of two TMD distributions. The kinematical mapping between the branching variables and the physical variables is dictated by the angular ordering, q T = µ(1 − z), where µ and z are the branching scale and longitudinal momentum transfer, and q T is the final-state transverse momentum [28]. The evolution of the transverse momentum in the backward shower follows that of the TMD distribution. One of the features characterizing the TMD shower compared to collinear shower algorithms, in particular, concerns the treatment of the soft-gluon resolution scale z M . The resolution scale in the shower matches the resolution scale in the TMD evolution equation (2.2). This can be contrasted with the collinear shower case, in which parton distributions evolve via DGLAP equations, while the resolution scale affects the shower: see discussions in Refs. [44][45][46].
A TMD final-state shower is not available yet. Current applications of the PB method use standard Pythia [47,48] or Herwig [49,50] final-state showering algorithms. These are applied both to the final-state partons and to the timelike splittings of partons radiated from the initial state. A study of final-state showering in PB calculations is carried out in Ref. [51].
The TMD forward evolution (2.2) with backward shower for the initial state and the final-state shower will be used for multi-jet merging in this work, as will be described in the next section.
Numerical solutions to Eq. (2.2) have been used in [52] to obtain TMD densities at NLO from fits to precision deep inelastic scattering (DIS) HERA data [53], performed using the fitting platform xFitter [54,55] and the numerical techniques developed in [56] to treat the transverse momentum dependence in the fitting procedure. In [52] two sets of PB TMD distributions are described: Set 1, which uses the evolution scale as argument in the running coupling α s , similar to what is used in HERAPDF 2.0 NLO [53], and Set 2, which uses the transverse momentum in the evolution of α s . These PB TMDs have been combined with NLO calculations of Drell-Yan (DY) production in the MadGraph5_aMC@NLO [57] framework to compute vector-boson transverse momentum spectra at LHC energies [58] and fixed-target energies [59]. The theoretical predictions thus obtained are shown to provide a good description of DY transverse momentum measurements over a wide range in energies and masses, from the LHC [60,61] to PHENIX [62] to R209 [63] to NuSea [64,65].
By the same NLO + PB TMD method, angular correlations have been examined in di-jet [66] and DY + jet [51] final states at large transverse momenta. A good description of the LHC angular measurements [67,68] is found.
The PB approach has recently been used [69] to make a determination of the nonperturbative rapidity evolution kernel (see e.g. [70]) in the factorization formalism [20,71] for DY at low transverse momenta. The calculation is based on the event generator described in Ref. [40]. The results may be directly compared with extractions of the rapidity kernel from DY experimental data, e.g. [72], and with determinations from lattice calculations, e.g. [73].
An extension of the PB evolution equations has been proposed in [74] to take into account the transverse momentum dependence of splitting functions at each branching, as defined from high-energy factorization [75]. This extension is relevant for phenomenology at the highest energy frontier, where the small-x (Regge) phase space [21] opens up and highenergy resummations are called for. See Refs. [76][77][78][79][80][81][82][83][84] for shower Monte Carlo algorithms including the small-x region. In the present work, we will not consider this generalized evolution, and we will limit ourselves to using PB evolution in the form discussed in Refs. [28,33].

Transverse momentum broadening
The evolution equations (2.2) imply that a TMD distribution characterized by a given width in transverse momentum (k T = |k|) at the scale µ 0 will be subject to a broadening in k T as the evolution scale µ increases, as a result of multi-gluon emissions. An example, taken from the TMD parton distribution library [85,86], is shown in Fig. 1. Our primary interest in this work is to analyze the consequences of the large transverse momentum tails, induced by TMD evolution, on the structure of final states with multiple jets in high-energy collisions.
Consider a final state characterized by hard scale µ, e.g., the transverse momentum of the hardest jet in the event. To assess the contribution to the production of an extra jet with transverse momentum p T < µ from the high-k T tail of the TMD distribution evolved to scale µ, it is useful to introduce integral TMD distributions a j , obtained from A j in Eq. (2.2) by k T -integration as follows The distribution a j evaluated at k T = 0 gives the fully integrated initial-state distribution.
The fractional contribution to a j from the tail above transverse momentum k T , with k T of the order of the jet p T , is given by the ratio In Fig. 2 we illustrate the k T dependence of Eq. (2.3) by showing the integral TMD gluon distribution a g (x, k 2 , µ 2 ), normalized to k = 0, obtained from the PB TMD Set 2 [52] for x = 10 −2 and various values of µ. We note, for example, that for µ = 100 (500) GeV, there is a 30% probability that the gluon has developed a transverse momentum larger than 20 (80) GeV. The results in Fig. 2 illustrate that, while the distribution is falling off at large k T , for the jet transverse scales observed in the LHC kinematics the contribution from the region p T ∼ < k T ∼ < µ is non-negligible compared to the contribution of an extra parton perturbatively emitted via hard-scattering matrix elements. Both contributions need to be taken into account. To avoid the double counting between the extra jet emission induced by the TMD gluon, µ = 10 GeV gluon, µ = 100 GeV gluon, µ = 500 GeV initial-state evolution and that arising from the inclusion of a higher-order matrix element, a merging methodology is needed. This is the subject of the next section.

TMD multi-jet merging
In this section we describe the basic elements of the TMD multi-jet merging method. We illustrate its main features and apply it to the computation of differential jet rates (DJRs).

Basic elements
Parton showers are able to simulate multi-jet topologies. If an n-jet configuration is available at the matrix element level, the corresponding (n + 1)-jet configuration can be generated by a parton shower emission. However, the corresponding accuracy is limited to emissions in the soft and collinear phase space regions. A hard, wide-angle emission from an n-jet configuration will be better described by an (n + 1)-jet matrix element calculation. The naive sum of the n-and (n + 1)-jet calculations would not be correct due to regions of the (n + 1)-jet phase space which would be doubly populated, both by the (n + 1)-jet matrix elements and by the parton shower emissions off the n-jet configuration. Furthermore the (n + 1)-jet partial cross section would be unstable, with its value strongly depending on the phase space cut which prevents the extra parton from approaching the divergent soft or collinear regions. The main purpose of a jet merging algorithm is to enforce exclusivity of the (n + m)-jet matrix element calculations above a merging scale µ m , except for the highest available multiplicity which will remain inclusive. For instance, if the n-jet matrix element configuration is made exclusive, double counting with the (n + 1)-jet matrix element configuration will be avoided. This can be achieved by means of Sudakov form factors which will suppress emissions from the n-jet configuration. When properly applied to the (n + 1)-jet matrix element configuration, the Sudakov factors will also suppress the divergent phase space region, making the partial (n + 1)-jet cross section stable.
Using the Sudakov factor in Sec. 2, the (n + 1)-jet cross section dσ P S n+1 given by the parton shower approximation from an n-jet configuration can be expressed as where µ 2 refers to the scale of the emission, dφ rad = dµ 2 /µ 2 dz/z is the radiation phase space, µ 2 max is the maximum scale in the process (set at the value of the renormalization scale selected for a given process and final state, a parameter labeled as SCALUP in the conventional Les Houches file format [87]), R P S is the parton shower emission probability.
For simplicity the flavor indices and the z M parameter are not shown in Eq. (3.1). Similarly one can construct the (n + 1)-jet cross section dσ n+1 using the (n + 1)-jet matrix element, where R = σ n+1 (µ 2 )/σ n is the real emission probability. A merging algorithm would consist in using Eq. (3.1) for emissions below an arbitrary merging scale µ m and Eq. (3.2) for emissions above µ m . The residual discontinuity at the merging scale can be further reduced by reweighting the strong coupling in Eq. (3.2) to the corresponding shower history. The inclusion of TMD evolution effects in a merging algorithm can be performed according to the procedure that we describe next. This extends the standard MLM matching [5][6][7][8] (which we recall here for completeness). The merging procedure consists of the following items.
The generation cut µ c gives the lower threshold for the outgoing partons' transverse momentum. In the calculations that follow we take the generation cut to be µ c = 15 GeV. Parton samples are generated in any given kinematic configuration (e.g. within a certain pseudorapidity range) according to the matrix elements, with a probability proportional to the respective cross section.
2. Reweight the strong coupling in the matrix elements according to the values from the corresponding shower history, as prescribed in Ref. [3] and implemented in the MLM algorithm [8].
3. For each of the two initial state partons of a given event, extract values of k i (i = 1, 2) distributed according to the evolution Eq. (2.2), setting µ 2 = µ 2 max . If k i 2 ≥ µ 2 min for any i = 1, 2, the event is rejected, and its contribution to the sample cross section is subtracted. µ 2 min corresponds to the minimum energy scale for the event, defined by µ 2 min = min{p 2 ti , p 2 tij } where i, j = 1, .., n, p ti is the transverse momentum of parton i, and p 2 tij (i = j) measures the relative transverse momentum between partons i, j. In the case of Z+0 jets, we selected µ min = m Z . The overall kinematics of the final state is then reconstructed by including the transverse boost induced by the transverse momentum k = k 1 + k 2 .
4. Initial-state partons of the generated events are showered using the backward shower evolution that corresponds to Eq. (2.2). Final-state partons are showered using the standard Pythia [47,48] or Herwig [49,50] showers.
5. The MLM prescription [7,8] is applied by matching the kinematics of the showered event and that of the parton-level event after the k boost. This differs from the standard MLM procedure, where the final state of the parton-level event has no overall transverse momentum. We recall here the steps of the MLM matching: • Partons produced in the showered events are clustered with a jet algorithm (e.g. the anti-k t one [88]) defined by a "cone size" and a transverse energy. The resulting jets are characterized by the merging parameters that should be kept fixed for all samples: the cone size R clus , minimum transverse energy E clus (merging scale) and maximum pseudo-rapidity within which jets are recostructed, η clusmax . The systematic uncertainty of the method can be estimated by varying these parameters.
• For a given multiplicity n, start from the hardest parton and select its closest jet in R. The parton and the jet are matched if R < cR clus . We take c = 1.5.
After removing the matched jet from the list the matching procedure proceeds for the remaining jets and partons. The event is matched when each parton is matched to a corresponding jet. Events are rejected if not matched.
• Exclusivity is enforced for the samples with n < N by requiring the corresponding matched events not to have extra jets in order to be accepted. When n = N , additional jets are kept, provided their transverse momentum is smaller than that of all matched jets.
We will refer to the LO method described above as the TMD merging method. We have described the method using the MLM merging prescription. Analogous LO versions of the method can in principle be derived using other merging prescriptions, such as CKKW-L, and adding the corresponding steps outlined at the points 3 and 4 above.

From MLM to TMD merging
To start the assessment of the impact of TMD corrections to multi-jet merging, we consider in this section the parton-level differential jet rates (DJRs) in Z+jets final states. The DJRs represent the distributions of the variable d n,n+1 , the square of the energy scale at which an n-jet event is resolved as an (n + 1)-jet event, with parton-level jets reconstructed by the k t jet-clustering [89,90]. Since the DJRs provide the splitting scales in the jetclustering algorithm, they follow closely the measure used in the definition of the merging scale. Therefore, they have always been considered a powerful means to test the consistency and systematic uncertainties of multi-jet merging algorithms [6].
For this study, we use MadGraph5_aMC@NLO [57] to generate Z + 0, 1, 2, 3 jet samples at LO with generation cuts for the partons given by: p T > µ c = 15 GeV, |η| < 2.5, ∆R > 0.4. We consider pp collisions at a center-of-mass energy √ s = 8 and 13 TeV, for later comparisons against published LHC experimental results [60,[91][92][93]. The PB TMD evolution is implemented in the event generator Cascade [40]. We use this to generate the TMD backward shower. We use Pythia6.4 [48] for the final-state shower. We apply the parton distributions obtained from DIS fits in [52] with α s (M Z ) = 0.118. The nominal value for the merging scale is chosen to be µ m = 23 GeV.  Figure 3. The d 01 spectrum at parton level, where d 01 represents the energy-square scale at which a 0-jet event is resolved as a 1-jet event in the k ⊥ jet-clustering algorithm. The dotted blue curve represents the ME level of the Z+1 jet calculation, while the green curves represent the Z+1 jet result using the "naive MLM" (left) and TMD merging (right) methods for different values of the R clus parameter. The predictions are computed in inclusive (inc) mode.
We begin by computing the lowest-order DJR d 0,1 , giving the energy-square scale at which a 0-jet event is resolved as a 1-jet event. In Fig. 3 (right) we present the result for d 0,1 from the TMD merging method introduced in Subsec. 3.1 for different values of the R clus parameter. In Fig. 3 (left) we compare this result with the result one would obtain by changing item 5 in the method of Subsec. 3.1 to a naive implementation of the MLM merging approach, in which the jets that are reconstructed after the parton shower has taken place are matched to the original ME partons (produced in item 1 of the method in Subsec. 3.1) instead of the partons resulting from the PB-TMD evolution (as in item 3 of Subsec. 3.1). We stress that this is not the canonical MLM implementation, in which no k evolution and boost is applied. We introduce this "naive" version of item 5 in order to highlight the relevance of properly adapting the parton-level event to the new kinematics induced by the TMD evolution.
The blue curves in Fig. 3 represent the ME level of the prediction. The green curves in Fig. 3 represent the result of calculations using the "naive MLM" (left) and TMD merging (right) methods for different values of R clus . We observe in the results of Fig. 3 (left), based on "naive MLM", a strong dependence of the inclusive Z+1 jet contribution to the d 01 distribution on the merging procedure, specifically the R clus parameter. Moreover the dependence is also present at high values of d 01 , where instead the ME accuracy should be preserved. The reason for this loss of rate relative to the input parton-level distribution is the angular mismatch between the initial parton, prior to the k boost, and the jet reconstructed after radiation, whose direction is modified by the TMD boost. On the other hand, in the results of Fig. 3 (right) we observe that the TMD merging method of Subsec. 3.1 gives a contribution to the d 01 distribution from the Z+1 jet calculation that preserves the ME accuracy at high scales. Furthermore, the effect from varying R clus is localized around the merging scale, as with the standard matching algorithms.
1-jet ME 1-jet + TMD 1-jet + TMD, step1  represents the energy-square scale at which an (n + 1)-jet event is resolved as an n-jet event in the k ⊥ jet-clustering algorithm. The dotted blue curve represents the ME level of the Z+1 (left) and Z+2 (right) jet calculations. The dotted green and magenta curves represent the Z+1 and Z+2 jet ME computations respectively, after the PB-TMD evolution is applied. The solid green and magenta curves represent the Z+1 and Z+2 jet computations respectively, after step1 of the TMD merging method is applied, where step1 corresponds to the items 1, 2, 3 and 4 of the procedure in Subsec. 3.1.
We next focus on features encoded in items 1, 2, 3 and 4 of the TMD merging method introduced in Subsec. 3.1. We refer to the part of the algorithm consisting of items 1, 2, 3 and 4 as the first step (step1) of the TMD merging method. As mentioned earlier in Sec. 2, both the TMD (forward) evolution and backward shower are used for the multi-jet merging method. They are represented, respectively, by item 3 and 4 of the TMD merging algorithm in Subsec. 3.1. Item 3, in particular, implies the condition k i 2 ≤ µ 2 min applied on the forward evolution of each initial-state parton. Fig. 4 presents results for the d 01 and d 12 DJRs corresponding to step1, illustrating explicitly the role of the µ 2 min condition, which was introduced in point 3. of Section 3.1. The dotted blue curves show the ME level (1,2-jet ME) of the predictions. The dotted green and magenta curves show the result at the level in which the ME are evolved forward using Eq. (2.2) (1,2-jet + TMD) but without applying k i 2 ≤ µ 2 min . The solid green (left) and magenta (right) curves in Fig. 4 show the contribution to the d 01 (left) and d 12 (right) distributions, from the Z+1 and Z+2 jet calculations respectively, resulting from step1 (1,2jet + TMD, step1). In each case, the difference between the dotted and solid curves shows explicitly the µ 2 min effect. A more explicit representation of the impact of the TMD evolution and of the µ 2 min cut is given in Fig. 5. The red line shows the transverse momentum k t (Zj) of a parton-level Z+jet sample, following the TMD evolution that gives, by construction, k t (Zj) = |k|. The µ min cut, defined in item 3 of the merging algorithm, suppresses events where the transverse momentum built up by the evolution of either of the two initial-state partons exceeds the jet p T , leading to the blue curve.  We now move on to study the effect of the second step (step2) of the TMD merging algorithm, consisting of item 5 in Subsec. 3.1. In Fig. 6 we show the contribution from the Z+1 (left) and Z+2 (right) jet calculations resulting from step1 (1,2-jet + TMD, step1) and step2 (1,2-jet + TMD, step1+step2), to the d 01 and d 12 distributions respectively. In addition we show the computations at the ME level (1,2-jet ME). The dotted blue curves represent the ME level of the prediction (σ n ), the dashed curves represent the calculation resulting from step1, and the solid curves are the result after step2 is applied.  Figure 6. The d n,n+1 spectra for n = 0 (left) and n = 1 (right) at parton level, where d n,n+1 represents the energy-square scale at which an (n + 1)-jet event is resolved as an n-jet event in the k ⊥ jet-clustering algorithm. The dotted blue curves represents the ME level of the Z+1 (left) and Z+2 (right) jet calculations. The dashed green and magenta curves represent the Z+1 and Z+2 jet computations respectively after step1 of the TMD merging method is applied where step1 corresponds to the items 1, 2, 3 and 4 of the procedure. The solid green and magenta curves represent the Z+1 and Z+2 jet computations respectively when the full TMD merging method is applied (step1 + step2, where step1 corresponds to the items 1, 2, 3 and 4 of the procedure while step2 corresponds to item 5). The predictions represented by the solid lines are calculated in inclusive (inc) mode.
In Fig. 6 we observe that step2 introduces the merging scale cut at 23 GeV. Furthermore, when implemented in inclusive mode, step2 applies a small correction around the merging scale to the Z+1 and Z+2 jets calculations resulting from step1. The resulting correction (∼ 25%) is larger for the Z+2 jets computation than for the Z+1 case (∼ 5%) because the suppression includes not only the regions in phase space for which the partons are soft, but also the region for which the two partons are collinear, which is not present in the Z+1 calculation. Figures 4 and 6 show that at all steps of the TMD merging algorithm the predictions for the largest multiplicities agree with the ME level, which is in accordance with the expectation that fixed-order contributions dominate higher scales, while TMD effects become important at low and intermediate scales.
In the following we study the effect of applying step2 in exclusive mode compared to inclusive mode. In Fig. 7 we show the contribution from the Z+0,1 (left) and Z+1,2 (right) jet calculations resulting from step2, to the d 01 and d 12 distributions respectively. The largest multiplicity in each subfigure is calculated both in inclusive and exclusive mode, while the lowest multiplicity in each subfigure is calculated only in exclusive mode. Figure 7 shows that the suppression from the exclusive mode of the TMD merging 0-jet + TMD, step1+step2 (exc) 1-jet + TMD, step1+step2 (inc) 1-jet + TMD, step1+step2 (exc)  Figure 7. The d n,n+1 spectra for n = 0 (left) and n = 1 (right) at parton level, where d n,n+1 represents the energy-square scale at which an (n + 1)-jet event is resolved as an n-jet event in the k ⊥ jet-clustering algorithm. The predictions represent the computations of the Z+0,1 (left) and Z+1,2 (right) jet multiplicities when the full TMD merging method is applied (step1 + step2), where step1 corresponds to the items 1, 2, 3 and 4 of the procedure while step2 corresponds to item 5. The Z+1 contribution to d 01 (left) as well as the Z+2 contribution to d 12 (right) are calculated in both exclusive (exc) and inclusive (inc) modes.
method can be of the order of 0.5 at large scales while it does not affect the region close to the merging scale compared to the inclusive mode. The difference between the inclusive and exclusive Z+1 jet calculations in Fig. 7 (left) and between inclusive and exclusive Z+2 jet calculations in Fig. 7 (right) would be covered by the inclusion of higher multiplicities as will be seen in Fig. 8. We have observed in Figs. 4, 6, and 7 the effect of the different steps of the merging algorithm, with the predictions from the different steps approaching each other and the ME at large d n,n+1 values. This is a consequence of the Sudakov-like suppression resulting from the merging algorithm which follows the Sudakov factor as in Eqs. (3.1) and (3.2). The performance of the merging procedure can be better appreciated in Fig. 8, where a clear separation between the different jet samples is seen at the merging scale value while the resulting overall prediction behaves rather smoothly. A small glitch remains visible at values of the d n,n+1 variables around the merging scale. This is a common feature of all matching algorithms (see e.g. [6,8]), and reflects the residual systematic error in the matching between shower and matrix elements. In principle the matching parameters, which are varied in order to establish a range for the systematic uncertainty, can be tuned to smooth out further the transition. We return to this point in Section 5, and in particular Fig. 21, where we study the matching uncertainty for the d n,n+1 variables. In Fig. 8 the predictions include Z+1,2,3 jets computations, where all the multiplicities are calculated in exclusive mode except for the highest multiplicity which is calculated in inclusive mode. The dotted curves represent the n-jet sample contributions, while the solid curve corresponds to their sum.

Z-boson transverse momentum and φ * distributions
In Fig. 9 (left) we show the transverse momentum p T spectrum of DY lepton pairs from Z-boson decays, normalized to the inclusive DY NNLO cross section. In addition, in Fig. 9 (right) we show the φ * η [95,96] normalized distribution of DY lepton pairs from Z-boson decays, where φ * η is defined as The angle ∆φ is the azimuthal separation between the two leptons, and θ * η represents the  Separate contributions from the different jet samples are shown. All jet multiplicities are obtained in exclusive (exc) mode except for the highest multiplicity which is calculated in inclusive (inc) mode. We observe that the Z+0 jet sample constitutes the main contribution at low transverse momentum p T while the impact of larger jet multiplicities gradually increases with increasing p T . A similar behavior is observed for the φ * distribution. The merged prediction provides a very good description of the whole DY p T spectrum as well as φ * distribution, with deviations of at most 10% around the points corresponding to the 23 GeV merging scale.
The results of Fig. 9 may be compared with the results obtained in [58] by applying PB-TMD evolution matched with the NLO Z-production matrix element, without including higher jet multiplicities. The TMD merging prediction of Fig. 9 retains the good description of the low-p T region already obtained in [58] (see also the recent analysis [98]), and improves the behavior [58] in the high-p T region by merging TMD showers with higher multiplicities.

Inclusive and exclusive jet multiplicities
In Fig. 10 we show the results for the exclusive (left) and inclusive (right) jet multiplicities in Z+jets events in pp collisions at √ s = 13 TeV. For illustration, in the rest of this section we compare the predictions with the ATLAS measurements [91]. Analogous measurements of Z+multijet final states have also been performed by CMS [93].   TeV are compared to predictions using the TMD merging calculation. Separate contributions from the different jet samples are shown. All the jet multiplicities are obtained in exclusive (exc) mode except for the highest multiplicity which is calculated in inclusive (inc) mode.
The very good agreement of the prediction with the experimental measurements in Fig. 10 illustrates that not only the whole DY p T spectrum is described by the TMD merging calculation, as seen in Fig. 9, but also the number of jets which result into the lepton pair p T imbalance. What is particularly remarkable is that the agreement holds up to multiplicities much larger than the maximum number of jets (three) for which the exact LO matrix-element calculation is performed. This underscores the potential benefit of the TMD evolution in better describing hard and non-collinear emissions, compared to the standard collinear evolution.
We next compare in Fig. 11 the TMD merging calculation of the ratio of consecutive inclusive jet multiplicities to the measurement by ATLAS [91] at √ s = 13 TeV. This measurement is sensitive to the strong coupling and is also very well described by the TMD merging calculation.

Jet transverse momentum distributions
In Fig. 12 we calculate the exclusive (left) and inclusive (right) leading jet p T spectrum in Z+jets events. The calculations are compared with experimental measurements by ATLAS [91] at 13 TeV. The contributions from the different jet multiplicities to the final prediction are shown separately. While for the exclusive leading jet p T spectrum the main contribution comes from the Z+1 jet multiplicity, in the inclusive case the Z+1 jet multiplicity is only important at low p T , with larger multiplicities giving the main contributions at large p T of the jet. A similar very good description of the data is observed for both observables, notwithstanding the very different multiplicity contributions to the predictions.
Additionally, in Fig. 13 we calculate the leading jet p T spectrum in inclusive Z+2 (top left), Z+3 (top right), and Z+4 (bottom) jet events. The calculations are also compared with experimental measurements by ATLAS [91] at 13 TeV.
The very good agreement of the predictions in Figs. 12 and 13 with the experimental measurements illustrates that not only the DY p T spectrum and the jet multiplicity distribution are described by the TMD merging calculation (Figs. 9 and 10), but also the p T of the leading jet contribution to p T imbalance is well described, both in the single-jet events and in the multi-jet events.

Di-jet azimuthal separation in Z-boson events
In Fig. 14 we show the distribution in di-jet azimuthal separation ∆φ for Z+ ≥2 jets events. The calculation is compared with experimental measurements by ATLAS [91] at  13 TeV. The contributions from the different jet multiplicities to the final prediction are shown separately. Fig. 14 shows that the agreement of the TMD merging prediction with the experimental data is very good in the high ∆φ region while a 20% deficit is observed in the low ∆φ region. We expect the deficit in the low ∆φ region to be due to missing multi-parton interaction (MPI) contributions. Inclusively, Z+2 jets processes constitute the leading contribution to the ∆φ between the two jets, assuming a single scattering. If a double parton scattering is assumed instead, the leading MPI contribution at low ∆φ corresponds to a Z+1 jet interaction together with a QCD 2-jet second interaction. To explore this, we make an estimate of the MPI impact on the ∆φ distribution in Fig. 14, based on the Pythia8 [47] MPI parameter tunes CUETP8M1 [99] and CP5 [100] obtained by the CMS collaboration. Our estimate is not intended as a systematic study of MPI effects in Z + jets distributions, but simply as a proof of principle that MPI is at the origin of the behavior seen in Fig. 14.
For the tunes used to perform the estimate, the effective cross section assuming two independent hard scatterings was determined in [99] and [100] to be, respectively, σ eff = 27.9 mb (CUETP8M1) and σ eff = 25.3 mb (CP5). For each set of parameter tunes, we obtain the MPI correction as the difference between the Z+1 jet computation of the ∆φ distribution including the simulation of parton showers and MPI, and the analogous computation without the simulation of MPI. In Fig. 15 we add the resulting correction to the TMD merged prediction. The solid blue curve in Fig. 15 is the same result as in Fig. 14, while the dashed dark-green and light-green curves are the result of adding to this the MPI contribution obtained with the parameter tunes CUETP8M1 and CP5, respectively. We  see that, while the region of the largest ∆φ is not affected significantly by MPI, at low ∆φ the MPI correction contributes a 10 -15% increase to the prediction, which supports our starting hypothesis.

Di-jet mass distributions
In Fig. 16 we show the di-jet mass distribution for Z+ ≥ 2 jet events. The calculation is compared with experimental measurements by ATLAS [91] at 13 TeV. The contributions from the different jet multiplicities to the final prediction are shown separately. The agreement of the prediction with the experimental data is good. Also for the di-jet mass we estimate the MPI effects, as described in the case of the ∆φ distribution in the previous subsection. The results are shown in Fig. 17. The MPI contributions do not affect significantly the region of high di-jet masses, while they are non-negligible at low di-jet masses, particularly in the first bin of Fig. 17.

Scalar sum of transverse momenta in Z-boson events
In Fig. 18 we show the scalar sum H T of the transverse momenta of leptons and jets for Z+ jets events. The calculation is compared with experimental measurements by ATLAS [91] at 13 TeV. The contributions from the different jet multiplicities to the final prediction are shown separately.
We observe good agreement of the TMD merging calculation with the experimental data. While at low values of H T the Z+2 parton contribution plays the main role, the higher multiplicity becomes increasingly important as H T increases.

Jet rapidity distribution
In Fig. 19 we show the leading jet rapidity distribution for Z+jets events. The calculation is compared with experimental measurements by ATLAS [91] at 13 TeV. The contributions from the different jet multiplicities to the final prediction are shown separately. We observe good agreement with the experimental data. The Z+1 jet multiplicity constitutes the main contribution to the leading jet rapidity distribution, as expected.

Theoretical systematic uncertainty and merging scale dependence
In this section we investigate the theoretical uncertainties associated with the TMD multijet merging algorithm. The inclusion of transverse momentum recoils through TMD evolution influences the systematic uncertainty of the merging when matrix-element and partonshower contributions are combined. Here we focus in particular on the systematic effects occurring through the dependence on the merging scale.
In Tab. 1 we show the multi-jet rates computed with TMD merging for different multiplicities for a 10 GeV variation of the merging scale. The results are obtained for pp  collisions at 13 TeV, with the phase space selection and cuts of [91]. The rates shown are absolute, i.e., they are not rescaled to the NNLO total cross section. We observe that a 10 GeV variation of the merging scale results in less than 2% variation for all the jet multiplicities considered. This systematic uncertainty is significantly smaller than what was found with standard algorithms of collinear merging in Ref. [8], where the variation of the jet multiplicities was found to be about 10% for a 10 GeV change in the merging scale. Besides the effect on the total jet rates, we next examine the merging systematic uncertainty in the case of differential distributions. We use the TMD merging algorithm with the default value of the merging scale at 23 GeV, as in the previous calculations, and consider variations of the merging scale around the default value to 20 GeV and 30 GeV.
In Fig. 20 we show the transverse momentum (right) and normalized φ * (left) distributions of DY lepton pairs in pp collisions at 8 TeV. We observe that the effect of the merging scale variation is localized around the merging scale and it is lower than 10%.
In Fig. 21 Figure 19. Leading jet rapidity distribution for Z+jets events. Experimental measurements by ATLAS [91] at √ s = 13 TeV are compared to predictions using the TMD merging calculation.
Separate contributions from the different jet samples are shown. All the jet multiplicities are obtained in exclusive (exc) mode except for the highest multiplicity which is calculated in inclusive (inc) mode.
for DY lepton pair production in pp collisions at 8 TeV. As discussed in previous sections the d n,n+1 distribution is very sensitive to the merging scale choice. We observe that the effect of the merging scale variation is localized around the merging scale, where small kinks in the spectra appear. The effects are however small, lower than 20%, reducing the systematic uncertainty observed by all merging algorithms discussed in Ref. [8]. In Fig. 22 we show the transverse momentum spectra of the leading (top left), second (top right), and third (bottom) jets for DY production in association with jets in pp collisions at 8 TeV. We observe that the effect of the merging scale variation is localized around the merging scale and it is lower than 8%.
The above results for the total multi-jet rates, the DJRs and the transverse momentum spectra together build a picture indicating that the systematic uncertainties from multi-jet merging are reduced when the transverse momentum recoils in the shower evolution are treated through TMD distributions.

Comparison with collinear multi-jet merging
In this section we present a comparison of our results, based on TMD jet merging, with results from collinear jet merging.
To be specific, we compare the TMD merging calculation with a calculation in which the initial-state TMD shower evolution is replaced by collinear shower evolution and the TMD merging is replaced by the MLM merging, while the matrix element and final-state shower evolution are kept the same in the two calculations. Since the final-state parton shower in our TMD merging calculation uses the parton shower routine PYSHOW of Pythia6 [48], we compare our TMD merging results with the results which we obtain by using Madgraph+Pythia6 with MLM merging.
We have first verified that, as expected, at the matrix element (ME) level our results coincide with those of Madgraph when using the same generation cuts. Next we have verified that, when considering only the final-state parton shower (without including any PB-TMD evolution and initial-state shower), our TMD merging results agree with the results obtained with Madgraph+Pythia6 computations. We have finally compared the full TMD merging and the Madgraph + Pythia6 calculations. Results are reported in Figs. 23-27. Clear differences emerge in distributions that are most sensitive to higher-order shower emissions, in particular the overall jet multiplicity, shown in Fig. 23, and the p T spectrum of the leading jet in final states with at least 4 jets, shown in Fig. 25. The better agreement of the TMD merging calculation with data, relative to the canonical MLM-matching procedure implemented in the Madgraph+Pythia6 result, together with the observation made previously that the two approaches are equivalent when limited to the final-state showers only, reinforce the conclusion that the TMD initial-state evolution leads to a better description of higher-order, non-collinear emissions.

Studies of PDF and intrinsic k T dependence
In the previous sections we have studied multi-jet merging methods and we have concentrated in particular on the implications of TMD evolution (see Figs. 1 and 2) on multi-jet observables, examining in detail the example of Z + jets production at the LHC. One may wonder about the role of the nonperturbative TMD parameters at the initial (low) scale of the TMD evolution in the analysis of the multi-jet observables, and also about the role of PDFs, i.e., of the nonperturbative collinear parameters.
In traditional applications of TMD physics such as the low-p T region of DY vector-boson spectra, both nonperturbative TMD contributions and PDFs play an essential and intricate role, as recently discussed e.g. in [72,101]. Analogous scenarios are explored for processes   the calculation follows the one in [91], whose data are included in these and subsequent plots.  Figure 24. Predictions obtained using Madgraph+Pythia6 with MLM merging and the TMD merging framework are compared for exclusive (left) and inclusive (right) leading jet p T spectra in the production of a Z-boson in association with jets. The phase space for the calculation follows the one in [91].  to be measured by future DIS lepton-hadron experiments in [102,103]. In this section we see that the case of multi-jet physics is significantly different in this respect. That is, at the jet scales of the final states investigated in Secs. 4, 5 and 6, one is sensitive to the large-k T tails induced by TMD evolution (illustrated in Fig. 2), while having very little sensitivity to the intrinsic k T distributions at low evolution scales and to differences among PDF sets. This means that, on one hand, the large-k T tails due to multiple QCD radiation embodied by TMD evolution contribute to the improved description of high-multiplicity final states, as we have seen in Secs. 4 and 6, and to the reduction of the merging systematic uncertainty, as we have seen in Sec. 5; on the other hand, these predictions do not really depend on intrinsic k T distributions or PDFs.
In Figs. 28-32 we study the intrinsic k T dependence. We show results for multi-jet observables obtained in three different intrinsic-k T scenarios. The solid blue curves corre-       spond to the TMD merging results of Sec. 4 obtained with the TMD parton distributions of PB-TMD Set 2 extracted in [52] from fits to DIS data. This TMD set has an intrinsic k T distribution at the initial evolution scale µ 0 = 1.4 GeV described by a gaussian parameterization with width σ = 355 MeV. The dashed red and green curves are obtained by varying this intrinsic-k T value by a factor of 2 up and down. This amounts to a very significant distortion in the k T distribution at the initial scale. Similar variations are studied in the case of DY transverse momentum distributions in Refs. [58,59]. The results in Figs. [28][29][30][31][32] show that the predictions for multi-jet observables change little under such variations, with the changes being smaller than the experimental uncertainties.    We next investigate the effect of using different choices of collinear parton distributions in the calculation of the parton-level ME samples. The nominal choice of the integrated TMD parton density for the ME calculation is compared to the results obtained using the NNPDF2.3LO [104]

Remarks on (off-shell) TMD effects in matrix elements
In this work we have investigated the impact of TMD evolution and k T broadening, depicted in Fig. 2, on the structure of multi-jet production. We have not, however, studied the possible TMD effects which may arise at the level of hard-scattering matrix elements [21], as e.g. in the parton-level event generator described in Ref. [106]. It is known that off-shell TMD matrix elements become important in the high-energy (Regge) limit, where they control small-x logarithmic resummations of high-energy cross sections. Our work has not been aimed at treating the Regge region, and thus the off-shell matrix elements [21,106] have not been used. Nevertheless, in this section we provide a few remarks on this topic, and on the relationship between the k T -dependent (off-shell) matrix elements and the approach used in this work.
To this end, we recall [107] that the small-x resummation of high-energy cross sections is precisely achieved by convoluting the hard cross sections obtained from off-shell TMD matrix elements with the gluon Green's functions obtained from the solution of the BFKL equation [108][109][110]. Such matrix elements, although off-shell, can be understood gaugeinvariantly as high-energy limits of amplitudes for n+2 particle production, representing the reggeized initial states. As functions of the hard-scattering scale µ and the gluon's transverse momentum k T initiating the hard scattering, they behave as follows. For the gluon's off-shellness going to zero, k T µ, the collinear matrix elements are recovered. For finite k T /µ, corrections to the collinear matrix elements arise, which may be computed as a power series expansion in (k T /µ) n . The summation of these corrections leads to a "dynamical" cut-off of the partonic cross section at k T of the order of the hard scale, k T ∼ < µ, and to a falling transverse momentum tail for k T ∼ > µ. Explicit examples may be found e.g. in [111], figure 4, and [112], figure 2. Now let us convolute the cross sectionσ(k T ) obtained from the matrix elements described above with the TMD parton distribution. As shown in [107], this yields (by a suitable double Mellin transformation, performed both in transverse momentum and in energy rather than just energy as in the collinear case) the resummation of high-energy (Regge) logarithms if we take BFKL solutions for the TMD. On the other hand, if we approximate the k T -dependent cross section by a Θ function viaσ(k T ) ∼ Θ(µ − k T )σ(0), we obtain the collinear cross section times the integrated TMD density.
In the latter case, TMD effects do not change the total cross section with respect to the collinearly-factorized calculation but affect the structure of the associated jet final states through the changes in the kinematics of the initial state due to the k T in the TMD density. This is the scenario of the present work. In other words, the k T corrections to matrix elements [21,106,111,112] are essential if we are to resum Regge logarithms in the total cross section; but to analyze the kinematical effects of the k T generated by the initialstate TMD evolution, which leave the cross section unchanged but change the jet structure of the events, we may work withσ(0) (i.e., without modifying the weight with respect to collinear matrix elements) and take into account effects of k T in the TMD density. These are the effects which we illustrated earlier in Fig. 2. If we, as in the present work, aim at including the contributions of higher multiplicity matrix elements along with the TMD evolution illustrated in Fig. 2, we need an appropriate TMD merging method, and this is the main subject of this paper.
Thus, the approach followed in this paper and the TMD matrix element results [21,106,112] are to be regarded as complementary. The work in this paper does not use the complete form of the matrix element but, as outlined above, it can be viewed as stemming from a collinear approximation to it, which is sufficient for the purposes of the present study, away from regions strongly affected by small x values, and can be regarded as a starting point for a more complete treatment. Such a treatment, in order to address the truly Regge (high-energy) region, would have to incorporate explicitly the effects of k T not only in the kinematics and TMD distributions but also in the matrix elements. An appropriate generalization of the TMD merging approach presented in this paper will be required in this case.

Conclusions
We presented in this work a new development for the description of multi-jet final states in hadronic collisions, which complements the standard approach, relying on merging samples of different parton multiplicity showered through emissions in the collinear approximation, with the use of the TMD parton branching for the initial state evolution. We carried out an extensive study of the theoretical systematic uncertainties associated with this new approach, and applied our results to a first comparison with LHC ATLAS data on Z plus multijet production. Our key findings can be summarized as follows: (i) a reduced systematic uncertainty with respect to the merging parameters and (ii) an improved description of high-order emissions, giving a better agreement for final states with jet multiplicity larger than the largest multiplicity used in the generation of the matrix-element samples. Our study used MLM matching as a merging criterion; we expect that the technique can be applied also to other LO merging schemes, such as CKKW-L. Since the introduction of the LO merging, and thanks to the automation of NLO ME calculations, more powerful NLO merging algorithms have been introduced. We hope that our work will inspire the exploration of a possible extension of the TMD merging approach to the NLO case as well.