Matching NLO QCD with parton shower in Monte Carlo scheme - the KrkNLO method

A new method of including the complete NLO QCD corrections to hard processes in the LO parton-shower Monte Carlo (PSMC) is presented. This method, called KrkNLO, requires the use of parton distribution functions in a dedicated Monte Carlo factorization scheme, which is also discussed in this paper. In the future, it may simplify introduction of the NNLO corrections to hard processes and the NLO corrections to PSMC. Details of the method and numerical examples of its practical implementation, as well as comparisons with other calculations, such as MCFM, MC@NLO, POWHEG, for single $Z/\gamma^*$-boson production at the LHC, are presented.


Introduction
Higher-order perturbative corrections in Quantum Chromodynamics (QCD), important for the LHC data analysis, are calculated order by order in the strong coupling, α s , while some numerically important ones, related to soft and collinear singularities, can be resumed to the infinite order. The most valuable, albeit technically difficult, way of QCD resummation is in form of a Monte Carlo (MC) event generator [1]. It is widely recognized that the most promising way of getting high precision QCD calculation for hadron collider data analysis is a common framework of fixed-order calculations combined with a parton shower Monte Carlo (PSMC). The pioneering work, in which the complete first-order QCD corrections to the hard process of heavy boson production in hadron-hadron collision were combined with PSMC, was that of Ref. [2]. Shortly later another interesting variant was proposed in Ref. [3]. Presently both methods are available for many processes, see Ref. [4].
It is worth to mention that partial efforts in this direction were done earlier in the course of development of the most popular PSMCs, like for instance in Ref. [5] or [6] and in many other works. In these earlier attempts, the distributions generated in PSMC were improved using a tree-level exact matrix element (ME) of QCD, while virtual corrections were neglected or taken in the leading-logaritmic approximation [7]. It was also known for a long time, that an ad hoc approach, in which PSMC differential distributions were corrected using the exact ME and the overall normalization was corrected by hand to fixed-order next-to-leading (NLO) integrated cross section, was providing distributions in a quite good agreement with experimental data. It may be therefore a little bit surprising that it took two decades to work out a systematic method of combining the NLO-corrections to the hard process, known from the early 1980's, see for instance Ref. [8], with the leading-order (LO) PSMC, also dating from the early 1980's. Apart from the lack of interest in more precise QCD calculations due to poor data quality, main reasons for this much delayed development can be seen from problems addressed in Refs. [2,3]. Namely, any such a method requires a very good NLO-level analytic understanding/control of distributions from PSMC and, either NLO-level complete phase space coverage for the hard process or a practical methodology of correcting for the lack of it. Luckily, a new wave of developments of the LO parton showers, see Refs. [9][10][11][12], has lead to modernized PSMCs, better suited for merging/matching with the fixed-order QCD calculations, in particular with better or even complete coverage of the hard process phase space.
It is now obvious that the next challenge on the way to even higher-precision perturbative QCD calculations needed until the end of the LHC era two decades from now, is to combine the fully exclusive NNLO corrections to the hard process and the NLO parton shower. The fixed-order NNLO corrections to many processes are already well established, see for instance Refs. [13,14], but the NLO PSMC needed for such a progress is still not available, except of feasibility studies summarized in Ref. [15]. Interesting partial solution of combining NNLO-corrected hard process with the LO parton shower can be found in 1 refs. [16][17][18]. The present work is relevant for the above future developments it the sense that it presents a simplified method of correcting the hard process to the NLO level in combination with the LO parton shower (PS). In other words, it offers a simpler alternative to the MC@NLO and POWHEG methods of Refs. [2,3], which may hopefully pave the way to the NNLO hard process combined with NLO PSMC.
The new method described here, nicknamed as KrkNLO, was already proposed in Ref. [19], where its first numerical implementation was done on top of the dedicated toy model PSMC and was limited to the gluonstrahlung subset of the NLO corrections. Later on it was tested numerically in a more detail in Refs. [15,20]. In the present work the KrkNLO method is implemented within the standard PSMC Sherpa 2.0.0 [11]. A pilot study of KrkNLO implementation outside PS MC, using MC event encoded in the event records produced by Herwig++ 2.7.0 [9] and Sherpa 2.0.0 was also done. Let us stress, however, that the overall simplifications of the KrkNLO method comes not completely for free, as it requires to use parton distribution functions (PDFs) in a special Monte Carlo (MC) factorization scheme (obtained, however, easily from reprocessing the MS PDFs), and it is required that the basic LO PSMC provides for the NLO-complete coverage of the hard process phase space (this is also not a problem for all modern PSMCs). Our method is simpler to implement in the case of PSMC with an ordering based on the transverse momentum k T or a q 2 variable of Ref. [21], inspired by the classic Catani-Seymour work [22]. However, it can be also easily implemented on top of PSMC that uses the angular ordering -without the need of the so-called truncated showers required in the POWHEG method, see Refs. [20,23] for more discussion on that.
The main advantage of the KrkNLO method is a simplification of the NLO corrections due to the use of PDFs in the MC factorization scheme. The implementation of the entire NLO corrections with help of a single multiplicative simple weight on top of the LO distribution is quite unique feature of the KrkNLO method.
Numerical studies presented here will be extended in the future publications to a wider range of distributions, energies, implementation variants, including comparison with the data.
The paper is organized as follows. In Section 2 we introduce the kinematics and the phase space parametrization for the considered process. In Section 3 we describe in detail the KrkNLO method. Sections 4 and 5 contain some numerical results of the KrkNLO implementation: Section 4 from the fixed-order NLO cross-checks while Section 5 from the NLO+PSMC comparisons with the MCFM, MC@NLO and POWHEG programs for the main Z-boson observables. Section 6 summarizes the paper. In Appendix A we add some details on the first gluon emission in the backward evolution in PSMC.

Kinematics and phase space parametrization
In the present work we are going to concentrate on the Drell-Yan process, specifically production and decay of the heavy boson Z/γ * 2 in proton-proton collisions. At the leading order (LO), qq → Z is the only partonic subprocess that contributes. At the next-to-leading (NLO) level, the real correction qq → Zg and the virtual correction to qq → Z contributes -to be referred to collectively as qq channel. The remaining NLO contributions qg → Zq andqg → Zq are tree-level only -to be called the qg channel. Fig. 1 illustrates part of the notation that will be used throughout the paper. The diagram shows the real correction to the Z-boson production in the qq channel with the gluon four-momentum denoted by k 1 . The four-momenta of the incoming forward and backward partons, p 1F and p 1B are related to the four momenta of the incoming protons, P F and P B , as follows

Single emission
where The invariant masses of the incoming parton pair prior to gluon emission and that of the produced Z boson are denoted by s 0 andŝ = s 1 , respectively. Their ratio 3 can be related to light-cone variables of the emitted gluon, as seen from the kinematics of the emitted gluons expressed in terms of the light-cone Sudakov variables The subscript in k 1 and s 1 is kept to underline that, in the context of PSMC, there is more parton emissions further away from the hard process qqZ vertex. spanned on the four-momenta p 1 F and p 1 B of the incoming partons, prior to gluon emission. As seen in Eqs. (2.4), the variables α 1 , β 1 are simply the fractions of the light-cone components of the gluon four-momentum to the centre-of-mass (CM) energy of the incoming partons (prior to gluon emission). The ratio z 1 , the gluon transverse momentum k 1T and the gluon rapidity y 1 are related to the α 1 , β 1 variables as well.
In PSMC the (eikonal) phase space measure of the emitted gluon is always split one way or another into two parts belonging to the quark and antiquark emitters. A sharp division along the y 1 = 0 angular boundary was used in Ref. [19], while in modern PSMCs, see Refs. [9,11,12], a more gentle division, introduced in Refs. [24,25] and inspired by the Catani-Seymour work [22], is exploited where m F , B are the so-called soft partition functions. What is important for our methodology in the following is that the sum of the two parts attributed to two emitters/showers reproduces the original 1-particle phase space without any gaps and (or well-controlled) overlaps. Clearly the m F + m B = 1 property of the functions used in Eq. (2.6) takes care of that. Once the above (overlapping) separation of the emission phase space into two part is applied, a different evolution variable of PSMC in each of them is defined instead of a common one, like k 2 T = α 1 β 1 s 0 or rapidity. Using the above evolution variable and z 1 , the single-emission (eikonal) phase space (2.5) is easily re-parametrized The relation between the old and new variables are illustrated graphically in Fig. 2. The transformation back to the Sudakov variables is different for each part: The upper phase-space limit α 1 + β 1 ≤ 1 transforms into while the positivity conditions, α 1 > 0 and β 1 > 0, enforce the IR-boundary cut-offs for the two parts, correspondingly. Also, in most of the phase space region populated according to ∼ m F factor we may approximate q 2 in the ∼ m B phase-space sector. The above kinematical limits are also shown on the logarithmic Sudakov plane in Fig. 2 for the ∼ m F sector. NB. The IR cut-off k 2 1T > k 2 T min marked in this figure translates into a slightly stronger cut-off on 1 − z 1 , easily calculable.
The same kinematic limits for one emission are also illustrated directly in terms of 1 − z 1 and q 2 1 F variables in Fig. 3, including also second emission for the purpose of the following discussion.
The essential ingredient of the 1-particle phase space reorganisation towards PSMC is an introduction of the on-shell "effective beams" p 0 F and p 2 B such that p 2 F + p 2,B = p 1 F + p 1,B − k 1 . Their definition is not unique. For example for the ∼ m F branch one may choose (2.12) Figure 3: The illustration of the kinematic boundaries in the forward evolution (FEV) algorithm with n = 2 particles. The parton no. 1 is emitted as a first one, within the biggest shaded triangular (blue) area defined by (1 − z 1 ) 2 > q 2 1 /s 0 and q 2 1 > q 2 s . The parton no. 2 is generated as a smaller triangular area marked with dashed (red) line according to (1 − z 2 ) 2 > q 2 2 /s 1 and q 2 2 > q 2 1 . Third parton is not generated, but its would-be-allowed space is marked by the smallest triangle below s 2 =ŝ. The vertical line marking the upper phase-space boundary q 2 < Q 2 in D(Q 2 , x) of Eq. (3.24) is also marked.
The same for the ∼ m B branch modulo obvious variable interchange. All that in the rest frame of p 1 F + p 1 B . However, in practice of a typical PSMC with backward evolution, all four-momenta are reconstructed starting from q 2 1 F ,1 B and z 1 variables in the rest frame of the effective beams p 0 F + p 0 B , which are constructed in the first place.
Technical details of the construction of effective beams are not so important for our analysis. Just as an illustration let us define explicitly one possible construction, which was introduced in Ref. [21]. In the rest frame of the hard process and effective beams P 2 = p 2 F + p 2 B one may construct all four-momenta starting from the (q 1 F , z 1 , φ 1 ) set, translating it into (α 1 , β 1 ) and using (for ∼ m F ): (2.13) Finally, knowing P 1 = P 0 − k 1 one may transform all newly constructed four-momenta k 1 , p 1 F , p 1 B to their common rest frame. Altogether, the complete reorganisation of the 1-real emission phase space from a simple form based on the Sudakov variables to equivalent parametrization using the variables of PSMC, based on backward evolution algorithm (applying the Catani-Seymour softpartition factor m F ), keeping track of the kinematical limits, and definingx = x 1 = x 0 z 1 , s = s 1 = sx 1 , looks as follows: (2.14) where q 2 1 = q 2 1 F ,P (z) ≡ (1 − z)P qq (z) and P qq (z) is the DGLAP [26] splitting function It is important to see that the full phase coverage requires integration over q 2 F to extend above the effective mass squaredŝ = s 2 of the LO hard process. Since q 2 F k 2 1T , it means that the transverse momentum aboveŝ is included in the above phase space.
In Section 5 on numerical results it will be commented more on what kind of evolution variable is chosen in the parton shower generation of Sherpa and Herwig++.

Multiple emissions
Although the above 1-emission kinematics is enough for most of our prescription for the NLO-correcting of the hard process, for better understanding of the role of the underlying multi-emission PSMC, it useful to extend the 1-emission treatment of the kinematics to 2 and more emissions in the initial-state parton showers. The important technical point is the choice of the numbering (labelling) of the emitted particles. In Fig. 3 we use the numbering in the emission chain (ladder) starting at the incoming hadron and ending next to the hard process, that is the labelling of the forward evolution (FEV) algorithm. From now on we switch to backward evolution (BEV) labelling which starts next to the hard process and ends at the incoming hadron, see Fig. 4 for illustration.
Let us summarize briefly on the effective beam technique for simplicity limiting its description to a single tree F of emissions (shower). The sequence of the effective emitter beams (p i F ,p i B ) is defined starting from the hard process, with the four-momentum P i , such that and they are used to span the four-momentum of the emitted gluon introducing the Sudakov variablesα i ,β i . These Sudakov variables are related to the evolution variable q 2 i F (a factorization scale) and the light-cone variable z i of PSMC as follows:β Finally, in the recursive backward reconstruction of the four-momenta starting from the hard process, one employs the Sudakov-like decomposition of the emitted parton in terms of the emitters after the emission 4 : (2.16) All exact kinematical limits (including ordering in the evolution variable) are represented by the following inequalities: In particular, the kinematical limits for the first emission in the backward-evolution (BEV) labelling are The above kinematic limits in terms of the BEV variables q 2 i and z i , look more complicated than in FEV scenario, although they represent exactly the same phase space region, and are illustrated graphically in Fig. 5.
As already underlined, in terms of the BEV variables, the phase space gets apparently widened after each emission, for instance q 2 >ŝ is already available for the first emission and due to lowering of the IR boundary by 1/z 1 factor more phase space is available for the second emission. This phenomenon, important for the full coverage of the phase space, is illustrated graphically in Fig. 6. It is, of course, an artefact of the BEV phase-space parametrization, which in the FEV world corresponds to the phase-space reduction due to energy conservation. particles. The parton no. 1 is emitted within the semi-triangular shaded (blue) area defined by (1 − z 1 ) 2 > q 2 1 /(s 0 /z 1 ) and s > q 2 1 > q 2 s . The parton no. 2 is generated within the second semi-trapezoid area marked by the dashed (red) line according to (1 − z 2 ) 2 > q 2 2 /(s 1 /z z ) and q 2 1 > q 2 2 > q 2 s . The third parton is not generated, but its would-be-allowed space is marked by the leftmost almost-trapezoid line. Figure 6: The illustration of the kinematic of the BEV algorithm with n = 2 particles in the extreme case when the parton no. 1 is emitted with q 2 1 > s 0 =ŝ and the parton no. 2 is generated within the area not accessible for the parton no. 1 due to a higher IR boundary.
We omit here the discussion of the "kinematical cross-talk" between two parton showers, which means that for emissions with the common q 2 -ordering in two initial showers (as in any realistic PSMC), the emission in one shower reduces the available phase space in the other shower. This effect is easily incorporated in the kinematical construction (mappings) of PSMC. The only thing one has to watch out is the correctness of the soft limit in the case of two and more emissions, see for instance the discussion in Ref. [21]. This subject will be covered in a more detail in our future publications.

KrkNLO methodology
Very briefly, the essence of the KrkNLO prescription defined in Ref. [19] is that NLO corrections are introduced by a multiplicative weight on top of distributions from the LO PS, acting either inside PSMC or outside it, on the MC event record. This NLO weight is sensitive to the parton four-momentum with the highest q 2 (or maximum k 2 T in case of k 2 T -ordering), although in Ref. [19] it was demonstrated that such a NLO multiplicative weight works also in the case of angular ordering, provided that summation over all emitted partons is performed.
In MC@NLO the correcting weight is essentially additive and the NLO − LO correction to PSMC is provided from outside in form of additional MC events, with a non-positivedefinite weight. In parts of the phase space which are not covered by the LO PS, extra events provide the entire NLO distributions (positive weights), otherwise extra events, correcting the LO distributions to the NLO level, have typically (inconvenient) negative weights.
In POWHEG the entire NLO correction to LO PS is provided by an external MC generator -the highest k 2 T emission is isolated/subtracted from PSMC, following the doublelogarithmic Sudakov exponential factor, and generated according to the NLO distribution outside PSMC, while trailing emissions with lower k 2 T (suppressed by the Sudakov exponent) are left for generating within LO PSMC.
Both KrkNLO and MC@NLO require good analytical control of the LO PS distributions, while for POWHEG it is less important. In addition, KrkNLO requires that LO PS fills in the entire NLO phase space with some non-zero distribution. The Sudakov suppression is also exploited in KrkNLO, but in a different way than in POWHEG.
In the standard NLO corrections to a hard process with MS PDFs, part of the NLO corrections feature degenerate 1-dimensional longitudinal phase space with k 2 T = 0 exactly. These corrections enforce in practice certain convolution on the MS PDFs, which in the case of POWHEG and MC@NLO is done in-flight inside MC. In the KrkNLO prescription, the implementation of these corrections is moved outside MC (simplifying it). This leads to the use in KrkNLO of the modified LO PDFs, in the so-called MC factorization scheme. This reorganization is mandatory because one of the main aims of the KrkNLO method, that the NLO corrections are implemented with a well-behaved multiplicative positive weight, is not compatible with such a degenerate collinear phase-space contribution to the NLO MC distributions.
Last but not least, any scheme for correcting the hard process of LO PSMC to the NLO level requires a formal proof that the resulting distributions are indeed of the NLO class (no double-counting, no NLO leak). Such a proof in an algebraic form is not trivial, not only because the NLO total cross section has to be verified, but also any NLOclass observable (experimental event selections) has to be properly reproduced. In other words, it has to be done using functional space of all the NLO observables. In the case of KrkNLO such a proof was done in Ref. [19], both algebraically and numerically, for a toy-model PSMC with the angular ordering. Here, we shall provide such an algebraic proof starting from the NLO-corrected multiparton distributions for realistic PSMC of the kind implemented in Sherpa and Herwig++ using the BEV algorithm.
In the following, we are going to collect building blocks for the NLO weight, then we shall elaborate on the multiparton distributions of LO PSMC without and with the NLO weight of KrkNLO. We shall pay particular attention to the question of the completeness of the phase space in PSMC and to the equivalence between the backward and forward evolution algorithms in PSMC. Finally, we shall show that for an arbitrary NLO-class observable, KrkNLO gives the same result as simpler NLO calculations with the collinear PDFs, instead of PSMC, such as for instance MCFM [27].

NLO-correcting weight
Let us collect the ingredients for construction of the NLO corrections to the hard process of the Z-boson production and decay in the proton-proton collisions.
The fully differential NLO cross section of the production and decay of the Z boson in the quark-antiquark annihilation process, with the simultaneous emission of a single real gluon, can be cast in a well-known compact form, see Ref. [19]: where the Sudakov variables 5 (α, β) are spanned on momenta of the effective beams of q andq prior to the gluon emission, see Eq. (2.4). The Born differential cross section for Z-boson production and decay is well known (see for instance [23] for the exact expression). The solid angle Ω = (θ, φ) is the direction of the lepton from the decaying Z boson in its rest frame andŝ = (1 − α − β)s. The angles θ F and θ B depend on α and β as well -their precise definition was given 6 in Ref. [28]. The integration over luminosities of q andq is not yet included.
How does the above compare with the distributions from PSMC also restricted to the single gluon emission? The differential cross section for the gluon emission from the quark emitter (i.e. the ∼ m F part in Eq. (2.6)) in PSMC reads where q F is the evolution variable defined in Eq. (2.7) andθ is another effective angle in Z decay specific to LO PSMC, for instance the so-called Collins-Soper angle [29]. Adding the gluon emission fromq simply amounts to the α ↔ β symmetrization, resulting in The integration limits are not explicit, but they are the same as in Eq. (2.4).
In our discussion, we shall often use the following objects: the additive NLO correction and the NLO multiplicative weight for the qq channel The above weight is especially simple in the case of averaging over the angles in Z decay: and it can be used in approximate implementations of the NLO corrections. The analogous weight for the qg channel is (3.8) After summing up the contributions from the two emitters, d 5 σ LO qq is obtained, which is exactly the same 7 as in Ref. [19]. The important consequence of the above is that many details of the matching between the MS NLO corrections and the LO PSMC in the present KrkNLO implementation remains the same as in Ref. [19]. In particular, the virtual plus soft-real correction, when the PDFs in the MC factorization scheme are used, is the same as in Ref. [19] and reads 8 As one can see, this virtual+soft-real correction is constant, i.e. kinematics-independent.

MC factorization scheme
In this subsection we extend the definition of the MC factorization scheme, introduced in Ref. [19] for the quark-antiquark channel, to the NLO DY process with the quark-gluon initial state. For the completeness and convenience of the reader we first provide the main formulae for the qq initial state. The NLO qq-channel coefficient function for the DY process in the MS factorization scheme is given by [19] The corresponding coefficient function in the MC factorization scheme, defined in Ref. [19], reads V S is the virtual plus soft-real gluonstrahlung correction given in Eq. (3.10). From the above equations, following Ref. [19], we can obtain a qq contribution to the relation between the MC-scheme and MS-scheme quark (antiquark) PDFs: Similarly, for the NLO qg-channel contribution to the DY process we have: Using Eqs. (3.13) and (3.16), we can relate the MC-scheme quark (antiquark) PDFs to the corresponding MS-scheme PDFs in the following way (3.17) The above relation is universal, i.e. process-independent, at the NLO level. It is simply because it is defined uniquely with respect to the MS scheme.
The gluon PDF is equal between the MC and MS schemes up to the O(α s 2 ) corrections for processes with no gluons at the Born level, such as the Drell-Yan process considered in this work. Hence for the DY process we may use As one can see, in the MC factorization scheme the NLO coefficient functions in both the qq and qg channels are substantially simpler than the corresponding ones in the MS scheme; in particular they are free of logarithmic singular terms. Since the latter terms are absorbed into the MC PDFs, i.e. they are included in a resummed way, one may also expect that the higher-order QCD corrections in the MC scheme are smaller than in the standard MS scheme.

Multiparton reality of PSMC
In the above we have restricted ourselves to a single emission and hence omitted all the multi-emission reality of PSMC. Let us elaborate on that in a more detail now, because the NLO correcting weights are put not on top of single-emission distributions but on top of multiparton distributions of LO PSMC, so they have to be known and controlled for the hard process within the NLO precision, preferably in a closed algebraic form.
Let us start from LO PSMC in the forward-evolution (FEV) formulation. The equivalent backward-evolution (BEV) formulation will be presented later on. Restricting ourselves temporarily to the pure gluonstrahlung case, the FEV differential cross sections of gluons emitted from the q andq emitters 9 reads as follows 9 We adopt a convention in which . The principal evolution variable q 2 i was introduced in Section 2, Eq. (2.7), and the labelling of the emissions starts from the hard process 10 , as in Fig. 4. The emission distributions for the ladder labelled with F are the following: where the Sudakov function reads and for the ladder B they look the same, except for theα i ↔β i swap.
The important variable s ij entering d 3 ρ F , B i and S F , B for the single shower/ladder was already defined in Section 2 as s i =ŝ/Z i , withŝ = sx F x B . For two showers, in any realistic PSMC, the emissions are generated (and the four-momenta are reconstructed) simultaneously in both showers using the competition algorithm, in which a common q 2ordering in both showers is emerging in a natural way 11 . Within such a common ordering for the two showers, the variable s ij =ŝ/(Z F i Z B j ) includes all z's from the emissions in both the ladders, starting from the hard process 12 .
Strictly speaking, the above implicit "kinematical coupling" of the two showers through s ij variable prevents us from rewriting the distributions of Eq. (3.19), without any approximation, into the traditional convolution of the two LO PDFs and the hard process, as it was possible in the toy PS MC in Refs. [19,20]. However, a slight modification of the kinematical coupling (modulo N 3 LO corrections) allows us get from Eq. (3.19) the following equivalent factorized inclusive formula where dσ/dΩ is the hard cross section from Eq. (3.1) and the LO PDF D F MC is resulting from the FEV algorithm run separately for each single shower, written in form of the 10 This is unnatural for the present FEV scenario, but better suited for the BEV algorithm in the following. 11 This method leads to forward-backward symmetric distributions, contrary to generating first the emissions from q and later on all the remaining emissions fromq. 12 In the PSMC jargon this is referred to as a recoil mechanism.
following time-ordered (T.O.) exponential 13 , is defined analogously. It is important to note that the objects D F,B MC appearing in Eq. (3.22) are not just scalar functions but they have non-trivial and well defined internal structure, as explicitly seen in Eq. (3.23). In particular, the MC PDFs, D F,B MC , result from the Markovian process and thanks to kinematical mappings they respect the phase-space constraints exactly. Therefore, they are not equal to the standard DGLAP parton distributions functions, in particular they integrate emissions up the absolute phase-space limit, c.f. Eqs. (2.4) and (2.10), rather than stopping at some arbitrary scale Q 2 = µ 2 F . The reader may check, analysing one and two emissions in a detail, that the effect of the above kinematical coupling of the two showers through the variable s ij is conveniently absorbed by the construction of the four-momenta defined in Section 2.2, hence Eq. (3.22) is equivalent to Eq. (3.19) up to the N 2 LO level, i.e. neglecting the N 3 LO and higher corrections. The above equivalence could be also tested numerically, similarly as it was done in Ref. [19].
Why do we insist on the FEV representation of PSMC, knowing that any typical PSMC is built using the BEV algorithm? The important reason is that in any methodology of combining fixed-order perturbative corrections with PSMC one has to make an algebraic contact with the standard diagrammatic perturbative calculations, including factorization theorems, resummation techniques, etc., which are all defined within the standard Lorentz-invariant phase space (LIPS). The FEV parton shower works directly within the LIPS 14 , while in the BEV algorithm the relation between LIPS and the PSMC distributions is obscured due to the presence of ratios of the PDFs in the BEV distributions. The way out is to define a framework in which the BEV and FEV distributions are by construction exactly identical. This requires certain non-trivial extra care and effort. We shall follow this path in the following.
The PDF of Eq. (3.23) looks pretty standard, for instance it is directly implementable in form of the FEV Markovian MC. Its peculiarity is already quite obvious for a single emission n = 1, where the upper phase-space boundary limit is not just regular q 2 1 <ŝ, but more complicated q 2 1 < (1 − z 1 ) 2ŝ /z 1 , reflecting the realistic phase-space boundary of the hard process, while the lower boundary is just regular q 2 1 > q 2 s , i.e. the same as from the ordering in PSMC. This peculiarity influences mainly the first emissions, closest to the hard process 15 , as discussed in Section 2. The consequence of the above peculiarity is that for our aim of the rigorous correspondence of the multiparton distributions between the FEV to BEV algorithms we need an auxiliary PDF with two competing factorization scales, Q 2 andŝ: where q n+1 ≡ q s and it coincides with the PDF of Eq. (3.19) for The main difference between this new PDF and the one introduced earlier in Eq. Fig. 3, that restrict emissions only to the region below the scale Q 2 , which in principle is an arbitrary parameter, as long as Q 2 ŝ. In that sense the new PDFs, D F,B MC , are closer to the standard DGLAP parton distribution functions.
With the above definitions we may work out the BEV algorithm providing the LO distributions, which by construction are exactly the same as from the FEV algorithm 16 : where x F i = i m=1 z m includes z i from d 3 ω F i and similarly for x B j from d 3 ω B i . The singleemission distributions and form factors are defined as follows: 27) and the shower label H = F , B was temporarily omitted wherever it was unambiguous to do so. 15 Already for the second emission, the limit q 2 2 < q 2 1 is more important than q 2 2 < (1 − z 2 ) 2ŝ /(z 1 z 2 ), which can be neglected modulo NNLO. 16 Modulo the kinematical coupling mentioned above.
Although it may not look obvious, the integrated cross section and all the multiparton distributions are exactly the same as in Eqs. (3.22) and (3.23). Let us present a sketchy proof of that 17 . Directly from the BEV algorithm for the PDF of Eq. (3.23) of the showering quark we get (3.28) In order to prove that the above PDF is the same as from the FEV algorithm one exploit the fact that in the T.O. exponential representation of the double-scale PDF of Eq. (3.24) we may detach (factorize off) the first emission, obtaining the following integral equation: . (3.30) The above identity can be used many times in order to eliminate the ratios of the PDFs in the BEV formula and transform without any approximation the BEV distributions of Eq.(3.28) into the FEV distributions of Eq. (3.23).

NLO weight in parton shower and its algebraic validation
Having defined the framework in which the FEV and BEV distributions of LO PSMC are identical, we may define a MC weight correcting the hard process in the (notoriously) inefficient FEV LO PSMC, but having advantage of a straightforward connection to the perturbative expansion, and apply it within LO PSMC implemented using the efficient BEV algorithm. The NLO completeness of such a scheme can be proven analytically within FEV, without any troubles due to the absence of the 'obscure' ratios of the PDFs of BEV. The only extra cost is that in order to gain a perfect FEV ↔ BEV compatibility, we had to introduce the special LO PDFs D(ŝ|Q 2 , x) with two factorization scales, at least for the purpose of the discussion of the NLO completeness of the KrkNLO scheme, while in the practice we shall be able to obtain D from the standard MS PDFs The differential cross section in which the hard process is NLO-corrected for the qq channel in the KrkNLO method reads as follows 18 : where see eqs. (3.6) and (3.10). The structure of this weight, with the summation over the emitted gluons, is a straightforward generalization of the weight used in many MC programs with multiphoton exponentiated corrections, see for instance Ref. [32].
As it was pointed out in Ref. [20], in order to get the complete NLO corrections to the hard process, it is enough to retain in the sum in Eq. (3.31) only a single term, the one with the maximum k 2 T . Since our q 2 -variable is practically identical to k 2 T we may retain only one term for a gluon with the maximum q 2 , from one of the two showers. In the case of the BEV algorithm with competition, this gluon is just the one which was generated first.
Similarly, as we have detached the first emission from the T.O. exponential representation of the PDF in Eq. (3.29), we now factorize off explicitly the gluon emissions with the maximum q 2 F and q 2 B from the product of two PDFs in Eq. (3.22): where s 1 F =ŝ/x F and s 1 B =ŝ/x B . Moreover, we restore also the kinematical coupling of the two showers, that is we are going to replace in the following s 1 F , s 1 F → s =ŝ/(x F x B ). The above decomposition is easily exploited in calculating the NLO cross section with the MC weight of Eq. (3.31) 18 We adopt a convention that 2 n=1 d n = 0. truncated to a single term W [1] qq with the maximum q 2 : where we have introduced x F and x B before the emission, as in a typical fixed-order NLO calculation. They are related tox F andx B from Eq. For the purpose of the proof of compatibility of the above formula with the fixedorder calculation (MCFM) for the entire class of the LO and NLO observables, we have introduced in the above the jet functions, J LO (x F , x B ) and J NLO (x F , x B , z 1 , k 2 1T ), which satisfy the following properties: often referred to as an infra-red safety requirement. In the experimental practice, the function J NLO represents the most general 2-dimensional histogramming in the transverse momentum of the gluon (or the heavy boson) and the z 1 variable is related to the rapidity difference between the heavy boson and the emitted gluon (or the ratio of the effective mass of the heavy boson and the heavy boson plus gluon system). The meaning of the above limit k 2 1T → 0 is that for the first bin in k 2 1T , the one including the k 2 1T = 0 point, we are not allowed to do any binning in z 1 (we have to sum up inclusively over all z 1 , similarly as for J LO we sum up inclusively over all k 2 1T ). In order to recover the fixed-order NLO formula of MCFM, all terms O(α 2 s ) have to be carefully eliminated. It is easy to do it for ∼ d 3 ρ W [1] qq parts, which are formally O(α 1 s ). We replace in them D where we were also allowed to replace z 1 F , z 1 B → z 1 , because in terms of the Sudakov variables before the emission it is the same variable. We have also introduced s 1 = sx F x B . The remaining parts ∼ (1 + ∆ V S ) are less trivial. In order to use again the identity of Eq. (3.33) for folding in three integrals within the {. . . } braces back into a simple product of two PDFs, we have to do something with the z 1 -dependence in dσ/dΩ and the k 2 1 -dependence of the jet function J NLO in the second part. The solution is to add and subtract a term proportional dσ dΩ (s 1 ,θ)J LO (x F , x B ) in the second part, regroup and use again Eq. (3.33): (3.37) 19 The collinear singularity in q 2 → 0 is killed by W [1] qq , so dq 2 /q 2 e −S(q2...) cannot give a ∼ 1/α s contribution.
where we have profited from finiteness of the integrals, in order to omit (1 + ∆ V S ) and the Sudakov exponents wherever possible and to set the lower integration limits of q 2 to zero.
Final clean-up of O(α 2 s ) terms involves the change q 2 1 F → s and q 2 1 B → s in the PDFs in the second integral and recombining it with the third one: (3.38) We may finally eliminate the MC weight of Eq. (3.32) going back to the NLO distributions: (3.39) The above is the typical fixed-order NLO calculation formula with the explicit definition of the NLO observable, as implemented for instance in MCFM, following the Catani-Seymour work [22], albeit with the PDFs in the MC scheme and absent the strictly collinear terms. Eq. (3.39) cannot be implemented, of course, in form of an efficient MC event generator with positive weights. In this way we have proven algebraically that the KrkNLO scheme is equivalent to the fixed-order NLO calculation in the entire functional space of the NLO-class experimental observables. It seems that this kind of a rigorous, albeit tedious, mathematical proof is not available for the other methods of combining the NLO corrections with the LO parton shower.
What remains to be done is to relate D F MC (ŝ|q 2 , x) to the standard PDFs of the MS scheme and summarize on the overall logics of the KrkNLO method of correcting the hard process to the NLO level, before describing its implementation and numerical tests.

Summarizing the KrkNLO method
In the following, we summarize on the key elements of the KrkNLO method, starting from the qq channel only and then adding the second qg channel.
1. The first element is the reorganization of the NLO corrections, the same as Ref. [19].
Starting from the standard MS calculation, strictly collinear part of the NLO corrections is removed thanks to the clever redefinition of the LO PDFs from the MS scheme to the MC scheme, see eq. (3.17). In a nutshell, in the standard MS procedure, one subtracts from the diagrammatic real + virtual results a pure ∼ P (z)/ pole, while in the case of KrkNLO one subtracts the Γ MC (z, ) function, which contains extra non-pole terms. As discussed in Ref. [19], this Γ MC function is just the integral over the LO distribution of the single-gluon in the parton shower MC 20 , written in 4 + dimension, together with the Sudakov form factor, obeying typical parton shower sum rules. This element of the KrkNLO prescription remains the same as in Ref. [19], provided the LO distribution is the same as in the present study. This assumption is valid because the LO distribution of Eq. (3.3) is indeed the same in Ref. [19] and in the Catani-Seymour inspired implementations of the single-gluon emissions in Sherpa and Herwig++. What remains to be checked is whether the upper phase-space limit, α + β ≤ 1, is not spoiled in a given PSMC by the use of the backward evolution algorithm. We shall come back to this point shortly.
2. The second element of the KrkNLO prescription is the construction, implementation, and algebraic validation of the multiplicative weight introducing the NLO corrections in the multiparton environment of PSMC. The weight of Eq. (3.6) is implemented following Eq. (3.31). The summation over gluons is necessary only for the angular ordering, while for the q 2 -ordering or k T -ordering it is enough to keep just one term for the gluon next to the hard process. Validity of above method can be proven rigorously, see Sec. 3.4, in the case when the initial-state PSMC is realized using the forward Markovian evolution (FEV) algorithm, in the presence of any observable-defining function J, see Eq. (3.35). Such a FEV algorithm would be terribly inefficient for the resonant Z-boson production, but it is mathematically perfectly well defined. NB. In such a gedanken FEV scenario with the complete phase-space coverage, the output PDF from the MC at the hard process scale will be automatically in the MC scheme, following Eq. (3.17), provided that the input PDF at low q 2 s is in the MS scheme 21 .
3. The above validation proof of the NLO weight has to be extended to the (efficient) backward evolution (BEV) scenario of the typical PSMC, such as Sherpa or Her-wig++. It was shown in Section 3.3 that it is feasible, by means of formulating the twin FEV and BEV algorithms, using the same PDFs, which produce exactly the same multiparton exclusive distributions. What is highly non-trivial is that the full phase space coverage is not lost. The price for this was that we had to introduce the auxiliary PDF D(ŝ|Q 2 , x) of Eq. (3.24) with two competing factorization scales, in addition to the single-scale PDF of Eq. (3.23). Luckily, at the end this auxiliary PDF can be finally eliminated from the BEV algorithm. Its main role is to provide theoretical control for the FEV → BEV transition, desired in view of the previous point.
4. The elimination of the auxiliary parton distribution function D(ŝ|Q 2 , x) can be seen best with the analysis presented in Appendix A, where the transition FEV → BEV is analysed in a fine detail for n = 0, 1 emissions. In Eq. (A.1) for n = 0, we see that we may replace D(ŝ|s, x) → D(ŝ, x) using Eq. (3.25). For n = 1 in Eq. (A.6), one may show that, modulo O(α 2 s ) terms, the ratio of D-functions can be replaced with any kind of the LO PDFs and the exponents, exp(−∆), can be neglected. What cannot be approximated there, it is the exact implementation of the evolution kernel, keeping the correct phase-space limits due to the θ-functions. This is however not the problem, because any standard BEV implementation with the veto algorithm does that correctly 22 . We conclude that, in the BEV scenario, all D-functions go away and only the single-scale PDFs in the MC factorization scheme seen in Eq. (A.1) are left finally in the game 23 . For the purpose of the next section we denote (Q 2 =ŝ): In the above analysis we have omitted qg channel. For the Z/γ * -boson production process however, is not difficult to include it following the same steps described above 24 In our implementation of the KrkNLO method in PSMC, whenever in the first step of the BEV PS algorithm the transition from the quark (antiquark) to the gluon is generated, we associate with such an event the weight computed according to Eqs. (3.8) and (3.9). This weight corresponds to observables averaged over the Z/γ * -boson decay angles is sufficient for the current paper, as we are going to concentrate on the Z/γ * -boson transverse momentum and rapidity distributions. The discussion of the Z-decay leptonic observables, as well as more details on the implementation of the qg channel, is reserved for a separate publication.  22 More detailed analysis of this kind was done but is left beyond the scope of the present paper. 23 NB. One can also show that the effective PDF generated by the above BEV at low q 2 s will be in the MS scheme -the BEV algorithm is not only undoing the LO evolution, but also the transition from the MS to MC scheme. 24 All ingredient distributions are well known in the literature. 25 Unless stated otherwise in the text.

MCFM
Sherpa Herwig++ σ tot [pb] 936.9 ± 0.1 937.2 ± 0.2 937.0 ± 0.6 Table 1: Values of the total cross section at the Born level for the Drell-Yan process in the MS scheme. and the G µ -scheme [33] for the electroweak sector of the SM. To compute the hadronic cross section we also use the MSTW2008 LO set of parton distribution functions [34], and take the renormalization and factorization scales to be µ 2 In order to check that our settings are identical in all used programs, we began with the comparisons at the Born level. The results presented in Table. 1 show a very good agreement (within statistical errors) between different programs.

PDFs in MC scheme
In Section 3.2, we have introduced the MC factorization scheme and explained why this scheme is better suited for the matching NLO results with the parton shower than the standard MS scheme.
The MC factorization scheme comes with a new set of parton distribution functions that can be obtained from the standard MS PDFs via the relations given in Eqs. (3.17) and (3.18). Note that the convolution terms in Eq. (3.17) introduce the dependence on the renormalization scale via the strong coupling multiplying the expressions for ∆C 2q and ∆C 2g in Eqs. (3.13) and (3.16). In what follows, we shall always choose that scale to be equal to the scale Q 2 = µ 2 F of the MS PDF. For clarity of the discussion, before presenting the complete NLO corrections to the Drell-Yan process, we shall be starting in what follows from results limited to the pure qq channel. This case corresponds to setting the gluon coefficient functions of Eqs. (3.14) and (3.15) to zero. That in turn implies that, for the two schemes to be consistent at O(α s ), the last term in the relation between the quark PDFs in the MS and MC schemes from Eq. (3.17) should be omitted.
The MC PDFs for light quarks and light antiquarks, obtained from the MSTW2008 LO set [34] are shown in Fig. 7 as a ratio to the MS PDFs. The curves correspond to the choice of Q 2 = 100 GeV. The same scale was used as the argument of α s . Each panel shown in Fig. 7 contains two curves. The solid line anticipates including both the qq and qg channels, hence it was obtained exactly with the formula (3.17). On the other hand, the dotted line corresponds to the case of the pure qq channel, hence only the convolution with ∆C 2q was used, as explained earlier.
As seen in Fig. 7, the MC PDFs without the gluon, e.g. without the last term in Eq. (3.17), are very similar to the MS quarks at low and moderate x. The MC PDFs in this region are only about 2% higher. At large x the ratio increases, reaching the value 2-3 near x = 1. Note that the MC PDFs are always higher than the standard MS quark PDFs, which results from the fact that the second term on the RHS of Eq. (3.17) is always positive. Adding the gluon component to the MC quark makes a little difference at large x, whereas the small and moderate x region is affected significantly. We see that the convolution with ∆C 2g is negative and, at small x, it can reduce the MS quark PDFs even by 20%. We have found a very similar picture for other quark flavours as well for other Q 2 values.
Before moving to NLO, it is interesting to check how the LO result changes when switching from the MS to MC PDFs. From the LO point of view, both PDFs are equivalent, as the differences are formally of the O(α s 2 ) order and higher. We can see, however, some numerical differences in the LO distributions between those cases. And, indeed, as shown in Fig. 8, the differential distributions of the dilepton mass and the dilepton rapidity differ at LO depending on whether the MS or the MC PDFs are used. In the case of the MC PDFs without the gluon convolution, c.f. Eq. (3.17), the mass spectrum is ∼ 5% higher than that with the MS partons over a broad range around the peak. On the other hand, when we include the gluon convolution, the result with MC PDFs gets up to ∼ 20% below that of MS. The dilepton rapidity distribution with the MC PDFs with (without) the gluon is smaller (larger) by a similar amount for central rapidities and Here x 1 and x 2 are the usual hadron's energy fractions of the incoming partons, while y Z is the rapidity of the produced Z-boson. The large rapidities in the forward or backward direction correspond to one of the x'es being large and that is the region where the differences between the MS PDFs and the MC PDFs are substantial.

MC and MS schemes at NLO
The cross section for the Drell-Yan process at NLO in the MS scheme can be schematically written as where f i (x 1,2 ) are the standard MS parton distribution functions, whose scale dependence is understood and ⊗ denotes the convolution via integration over z variable. The first term in the square brackets is just the Born contribution, which is followed by the NLO correction in the qq channel. The second sum corresponds to the NLO contribution from the qg channel. In both cases, the functions C MS q,g have non-trivial z-dependence (suppressed here for clarity of notation) and they are related to the coefficient functions from Eqs. (3.11) and (3.14) as C MS q,g = 1 α s C MS 2q,2g (z). Hence, the dependence on α s is explicit in Eq. (4.3) and all the following formulae of this section.
When going to the MC scheme, the MS coefficient functions in Eq. (4.3) need to be replaced by their MC scheme counterparts and the PDFs need to be transferred to the MC scheme as well according to Eq. (3.17). This leads to the formula where, again C MC q,g = 1 α s C MC 2q,2g (z), with the MC coefficient functions defined in Eqs. (3.12) and (3.15) and similarly for ∆ C q,g = 1 α s ∆C 2q,2g (z), with the latter given in Eqs. (3.13) and (3.16).
By construction (c.f. section 3.2), terms proportional to a given partonic luminosity, f i f j , are identical in both factorization schemes up to the order O(α s ). The results in Eqs. (4.3) and (4.4) differ however at O(α s 2 ), which is beyond NLO and therefore such a difference is allowed.
In order to validate our transformation from the MS to the MC scheme, we have performed an explicit check of the equivalence of the two schemes up to O(α s ), as well studied numerical importance of the higher order terms. The calculations were performed using our standard setup defined at the beginning of Section 4.
Let us start from comparing the cross sections at LO. This corresponds to setting the coefficient functions C q,g = 0 in Eqs. (4.3) and (4.4). As shown in the previous section, the LO cross sections will differ between the two schemes because of the PDFs. To check the extent to which this happens, we performed an explicit computation with MCFM [27] in both factorization schemes for either the pure qq channel or both channels.
Indeed, as demonstrated in Table 2, the LO MS and MC cross sections are not identical, both in the case of pure qq channel, where the LO result with MC PDFs is ∼5% higher, and in the case where both channels are included, where the MC result is ∼20% lower. The beyond-LO, O(α s ) and O(α s 2 ), terms are also given in the table. We see that most of the difference comes from O(α s ). The O(α s 2 ) terms are in fact very small, below 1% in both cases. It is interesting to note that the difference between σ (0) DY in the two schemes, for the case with both channels, comes primarily from the term proportional to α s f q ⊗ ∆C g ⊗ f g , hence it originates from large gluon luminosity.
We turn now to a similar comparison at NLO. Table 3 gives the NLO-only results in the two schemes for the cases of qq and both channels. At NLO, we need to be careful what we take for the coefficient functions. In the MS scheme, we just use the original MCFM 6.6 [27] implementation. However, in order to carry out consistent NLO 6.70 ± 0.20 f q ∆ C q fq 25.79 ± 0.04 f q ∆ C q fq + fq∆ Cf q 355.00 ± 0.29 16.61 ± 0.02 MC 147.9 ± 0.20 Table 3: Values of the NLO contribution to the total cross section for the the Drell-Yan process in the MS and MC factorization schemes. The results were obtained with MCFM 6.6 [27] and its version adjusted to the MC scheme, MCFM * . The O(α s ) corrections in the MC results is split into contributions proportional to various terms of Eq. (4.4).
calculations in the MC scheme, we had to change the coefficient functions from C MS 2q,2g (z) to C MC 2q,2g (z). We dubbed this modified version of the program MCFM * and used it together with the MC PDFs to compute our NLO predictions in that scheme.
We see that also the sub-leading corrections to the Drell-Yan process are not the same in the MC and MS schemes (numbers in the two first lines of Table. 3). This is allowed provided that all the difference comes from terms of the order α s 2 and beyond. As shown in Table 3, by extracting only the O(α s ) terms of the MC result and summing them up, we recover the MS cross section exactly. There, we also give subleading correction introduced via the MC PDFs and we see that they contribute at most 5% to the NLO correction in the MC scheme in the case of the pure qq channel but can be quite sizeable when both channels are considered.
Finally, it would be interesting to fit the quark and gluon PDFs in the MC factorization scheme directly to experimental data, in order to minimize higher-order effects. For that we would need to define the gluon PDF in the MC factorization scheme using a process in which the gluon PDF enters already at the LO level, e.g. the Higgs-boson production process. This will be done in our future publication.
Let us conclude that the we have successfully validated the MC factorization scheme by explicitly showing that all numerical differences w.r.t. the MS results are of the order of α s 2 and α s 3 . We emphasize that the validation of the new MC scheme presented in this section provides a highly non-trivial check as various components, given schematically in Eq. (3.17), come from different parts of the calculation (PDFs, coefficient functions) and the agreement up to α s is much more sophisticated than a purely algebraic relation.

Results for NLO with parton shower
Numerical implementation of KrkNLO is presently done mainly using the Sherpa PSMC of Ref. [11] version 2.0.0 26 with the dipole organization of the parton shower distributions inspired by the Catani-Seymour work [22], see Ref. [35]. However, the evolution variable in Ref. [35] is chosen to be the transverse momentum distribution, while in the actual Sherpa 2.0.0 implementation it is the q 2 ∼ α(α + β) variable of Section 2, see Ref. [36]. The dipole shower implemented in Herwig++, see Ref. [21], is quite similar. Actually, the choice of the evolution variable is not critical in the KrkNLO method, as long as the coverage of the phase space for the emission closest to the hard process is assured, and the PSMC distribution is under the perfect control. In fact, for any choices of the dipole PS evolution variable in Sherpa or Herwig++ the same distribution of the single gluon emission is obtained, provided that the contributions for the quark and antiquark emitters are added. For the Sherpa MC program, we have checked numerically the completeness of the coverage of the phase space of the first emission from the backward evolution by examining the gluon distribution on the α, β plane.
In the following, the KrkNLO-matched results will be compared mostly with the results of the MC@NLO technique, implemented on the top of the same parton shower within Sherpa 27 . Using the same parton shower in both cases ensures that all differences in the matched results come purely from the differences in the methods themselves rather than details of the shower implementation.
Nevertheless, the difference between dipole showers currently implemented in Sherpa and Herwig++ is limited, and it amounts to using the q variable in the former and k T in the latter. These two evolution variables are however very close in practise. Hence, in our final discussion, we shall also compare the KrkNLO results with those of POWHEG, as implemented in Herwig++ 28 using an automated setup based on Matchbox framework [37] with adaptive sampling [38].
In the numerical validation of the KrkNLO approach, we shall first compare it to the fixed-order NLO calculation, as implemented in the MCFM program. Later on, the comparison with MC@NLO implementation of Sherpa will be discussed. In order to minimize the influence of various higher-order effects on the difference between KrkNLO and MCFM, we will initially limit our numerical exercises to a simplified version, in which only the qq channel is kept. Another initial limitation, with the same aim in mind, will be to profit from the option in Sherpa to stop the backward evolution in the parton shower just after the first emission (starting from the hard process). This will help us to eliminate the influence of the differences in the so-called recoil schemes [39], which start to play a role from the second emission onwards. The running α s will be also introduced gradually. Since we study the differences in the matching methods, we have switched off in Sherpa non-perturbative effects including the intrinsic k T , multiparton interactions and hadronization. However, we would like to stress that there is no problem with switching them on for the KrkNLO method.
In all the following results for the KrkNLO method the PDFs in the MC factorization scheme, discussed in Sec. 4.1, will be used. The factorization scale will be set to µ F = M Z . As for the choice of the renormalization scale in the matched NLO + PS results, we note that there is always some level of arbitrariness. In the pure fixed-order NLO calculation, in the case of DY, one usually sets µ R = µ F = M Z . On the other hand, the LO parton shower uses µ 2 R = q 2 , with the latter being closely related to the evolution variable of the shower. The parton shower is unitary, hence the choice of the argument of α s (µ 2 R ) does not influence the total cross section, it changes, however, shapes of some distributions, e.g. of p T,Z . The differences between α s (q 2 ) and α s (M 2 Z ) are formally of the higher order, that is beyond NLO and, from that perspective, they are equivalent in the context of NLO + PS matching. They can however be numerically relevant.
In what follows, for KrkNLO method we adopt the following procedure. The argument of α s in the virtual correction is always set to M 2 Z . For the real correction, all emissions except for the first one (in BEV) are set by the parton shower according to its current q 2 value. For the first emissions, we consider two choices: q 2 and M 2 Z , both in the real correction and in the Sudakov form factor, to keep the shower unitary. The difference between the results corresponding to these two choices will be indicative of the size of  beyond-NLO terms.

Initial results for qq channel only
The first round of numerical tests of KrkNLO is done through comparisons with MCFM and MC@NLO for the qq channel only and stopping the parton shower after first emission. The corresponding results for the total cross section are presented in Table 4 and for the transverse momentum and rapidity distributions in Fig. 9. As we can see, the total cross section and rapidity distributions from KrkNLO agree very well with those of MCFM and MC@NLO. The difference between the KrkNLO results and the pure fixed-order ones is below 1% and comes from the O(α s 2 ) contamination due to the MC PDFs, c.f. Eq.(4.4). Also the differences between the α s (q 2 ) and α s (M 2 Z ) are at the per-mille level. On the other hand, there are remarkable differences in the p T,Z distributions, especially at the lower end. They compensate in the total cross section between the first and the following bins. However, the choice of the running or non-running α s is numerically important for the p T distribution at small and moderate values. The perfect agreement of MC@NLO and MCFM at higher p T is of course enforced by construction. This is not the case in KrkNLO, where some part of higher-order effects, beyond NLO, is included (which is also the case in POWHEG). We see, however, that both MC@NLO and KrkNLO coincide below ∼ 80 GeV, in the case with α s (q 2 ).
In the results of Fig. 9 the parton shower was artificially stopped after the first emission, in order to make the meaningful comparisons with MCFM and limit the effects of the recoil due to subsequent emissions. In Fig 10 we lift this limitation and the parton shower goes to the very end, as in the normal operational mode of PSMC. As we see, switching on to the full PS influences the low p T part of the spectrum quite a lot, in spite of a negligible effect on the rapidity distribution (and hence, on the total cross section).  Figure 10: Comparisons of the transverse momentum and rapidity distributions from MC@NLO and two versions of KrkNLO for the qq channel only as in Fig. 9, but the parton shower backward evolution runs to the end as normally.
Again, KrkNLO with the running α s agrees very well with MC@NLO at low and moderate p T , which results from the domination of the parton shower in that region and from the fact that the influence of the NLO-correcting weight of KrkNLO vanishes there. How-ever, the two KrkNLO results stay above MC@NLO at higher p T , which we again attribute to the admixture of some NNLO terms in the former. The difference between KrkNLO α s (q 2 ) and KrkNLO α s (M 2 Z ) at high p T comes purely from the running of the coupling.

All channels
In the second round of the numerical tests, all the initial-state parton combinations contributing at NLO, that is the qq and qg channels are included. The total cross sections from KrkNLO, MCFM and MC@NLO are compared in Table 5. As we can see, the dif- ference between MCFM/MC@NLO and KrkNLO, which comes from the partial inclusion of the higher-order effects in the latter, is about twice as big as in the the qq channel only. As discussed in Section 4.2, this is related to the large gluon luminosity, which leads to non-negligible differences between the MC and MS PDFs and part of that difference shows up in the total cross section. The results in Table 5 are generally consistent with the differences staying below 5%. In Fig 11 we show the distributions of the transverse momentum and rapidity of the Zboson for the case of both production channels, but with the parton shower stopped after the first emission. The lower panel contains the ratio with respect to MCFM. Similarly to the case of the pure qq channel, we see the expected agreement between MC@NLO and MCFM at large p T,Z . The KrkNLO result with the fixed α s is not far from MCFM in that region as well. On the other hand, all PS-matched results differ from MCFM at smaller p T because of the lack of the Sudakov resummation in MCFM. Each curve behaves slightly differently but the variations are moderate. The rapidity distributions shown in Fig. 11 (right) are close to each other for all calculations. The ∼ 5% difference between KrkNLO and MCFM is the same as for the total cross sections in Table 5. Fig 12 shows the similar distributions but, this time, the parton shower is allowed for arbitrary number of emissions. There, we also show the results from POWHEG, as implemented in Herwig++. Even thought the evolution variable is different (k T rather than q of Eq. (2.7) in all the other PS-matched results), this difference is negligible in practice apart from a few first bins in the p T,Z distribution. Only the first emission is upgraded  Figure 11: Comparisons of the transverse momentum and rapidity distributions from MCFM, MC@NLO and two versions of KrkNLO, now for both the qq and qg channels, with α s set as in the previous figures and the parton shower backward evolution stopped after the first emission as in Fig. 9.
to NLO in the PS-matched results. The distributions of the transverse momentum and rapidity look quite similar to the previous results for the qq channel in Fig 10. The transverse momentum distribution from KrkNLO agrees with MC@NLO at high p T,Z (up to a constant factor in the case with fixed coupling) and most of the differences are seen in the low p T region. The KrkNLO α s (q 2 ) version stays overall very close to MC@NLO. Even better agreement is found when the KrkNLO α s (q 2 ) result is compared with POWHEG. Here, the two p T,Z distributions almost coincide in the range shown in Fig. 12 (left, bottom panel). This is related to a similar, multiplicative way of applying of the NLO correction to PS. In both results, this correction is applied over the entire p T,Z range, hence, at large transverse momenta one becomes sensitive to mixed real-virtual terms. The rapidity distributions shown in Fig. 13 agree very well between all the PS-matched results, up to the normalization differences, the same as in Table 5.  Figure 12: Comparisons of the transverse momentum distribution from MC@NLO, POWHEG and two versions of KrkNLO for both the qq and qg channels as in Fig. 11, but the parton shower backward evolution runs to the end as normally.
in Section 3. The change of the factorization scheme allows one to eliminate troublesome z-dependent terms from the coefficient function and it effectively amounts to creating the MC PDFs. In Section 4, we have discussed how such PDFs can be obtained and how they differ from the standard MS parton distributions. There, we have also validated the MC factorization scheme by studying the Drell-Yan process at the fixed-order NLO level and showing that the MS and MC scheme results are identical up to the order O(α s ).
We have implemented the KrkNLO method on top of the Catani-Seymour type of parton shower in the Sherpa event generator for the case of production of the electroweak boson (hence, the initial-state parton shower). In Sectin 5, we have presented comparisons of the NLO-PS matched results obtained with our technique with the fixed-order NLO results from MCFM and with other matched results, namely those of MC@NLO and POWHEG. In particular, we have demonstrated that the KrkNLO results recover the fixed-order NLO predictions (up to sub-percent differences for the qq channel only and ∼ 5% for both channels, coming from beyond-NLO terms).
As for the comparisons of KrkNLO with MC@NLO and POWHEG at the level of differential distributions, all three methods turn out to give essentially identical results for the y Z spectrum. The p T,Z distributions look somewhat different for each method and the exact features depend on the initial channels and the recoil schemes. In general, the KrkNLO method provides similar predictions to the other two well-established approaches. In par-  Figure 13: Comparison of the rapidity distribution from MCFM, MC@NLO and two versions of KrkNLO for qq and qg channels as in Fig. 11, but parton shower backward evolution runs to the end as normally.
ticular, the results with both channels and the full parton shower stay very close to those from the POWHEG method implemented within Herwig++. Residual differences come from spurious O(α s 2 ) terms, which are different in each of the three matching methods.

APPENDIX A First emission in backward evolution
In case of two LO parton showers in the BEV algorithm the cross section with exactly zero gluons in both showers includes two non-emission ∆-functions which coincides exactly with the n = 0 result in the forward evolution, starting from q 2 F = q 2 B = q 2 s . NB. In the case n = 0, we identifyx F = x F andx B = x B . In the following we shall simplify the notation omitting theŝ argument in D MC (q 2 , x) = D MC (ŝ|q 2 , x) and S MC (q 2 2 , q 2 1 ) = S MC (ŝ|q 2 2 , q 2 1 ), wherever it is unambiguous. The distributions from the BEV algorithm with at least one gluon (n F +n B = 1, 2, 3...) reads as follows: where z 1 integration limits are imposed byŝ(1 − z 1 ) 2 /z 1 > q 2 1 in the K kernels.
The "unitarity" sum rule σ LO 0 +σ LO 1 + = σ LO is of course automatic in the BEV algorithm. Nevertheless, proving it provides an interesting cross-check. The differentiation over the lower boundary in the q 2 integral defining the ∆-function leads to:  (A.7) 29 We denote P(z) = C F αs πP (z) 1−z .
Finally, the longitudinal integrations are streamlined with the help of the substitutions to the x-variables before the emission, x F /z 1 → x F , x 1 B = x B or x B /z 1 → x B , x 1 F = x B , respectively in each shower, and the convolution structure is made manifestly symmetric: The above formula is for the FEV algorithm but for one gluon next to the hard process and any number of trailing gluons down to q 2 s . The same formula in Eq. (3.33) was obtained by means of factorizing such a gluon from the fully exclusive FEV formula of Eq. (3.24).