Markovian MC simulation of QCD evolution at NLO level with minimum k_T
IIFJPAN-IV-2008-6
Markovian MC simulation of QCD evolution at NLO level with minimum k T (cid:63)
P. Stok(cid:32)losa a , W. P(cid:32)laczek b , and M. Skrzypek a a Institute of Nuclear Physics IFJ-PAN,ul. Radzikowskiego 152, 31-342 Cracow, Poland. b Marian Smoluchowski Institute of Physics, Jagiellonian University,ul. Reymonta 4, 30-059 Cracow, Poland.
Abstract
We present two Monte Carlo algorithms of the Markovian type which solve themodified QCD evolution equations at the NLO level. The modifications with respectto the standard DGLAP evolution concern the argument of the strong couplingconstant α S . We analyze the z -dependent argument and then the k T -dependentone. The evolution time variable is identified with the rapidity. The two algorithmsare tested to the 0 .
05% precision level. We find that the NLO corrections in theevolution of parton momentum distributions with k T -dependent coupling constantare of the order of 10 to 20%, and in a small x region even up to 30%, with respectto the LO contributions. IFJPAN-IV-2008-6 (cid:63)
The project is partly supported by the EU grant MTKD-CT-2004-510126, realized in the partnershipwith the CERN Physics Department and by the Polish Ministry of Science and Information SocietyTechnologies grant No 620/E-77/6.PR UE/DIE 188/2005-2008. a r X i v : . [ h e p - ph ] O c t Introduction
With the first LHC results coming soon the precision of the QCD calculations becomesan important issue. As usuall, one finds two approaches to the calculations: fixed-ordercalculations and resumed to all orders Monte Carlo (MC) parton showers. The fixed-order results are impressive. Let us give a few examples of such calculations: the splittingfunctions, calculated up to NNLO level [1, 2] or various differential and semi-inclusivedistributions, calculated in the NNLO and even NNNLO approximation see e.g. [3–5]. Onthe other hand, the MC parton shower approach, indispensable in describing complicatedexperimental signal selection procedures, does not achieve such a high precision. In fact,the complete
NLO level has not been reached yet by any of the available MC partonshower codes. In most of the cases these codes are based on some form of an improved
LOapproximation, see e.g. [6–10]. Another approach is based on combining the NLO matrixelement with the LO parton shower [11, 12].Some time ago some of us started the MC parton shower project with the ambitiousgoal of reaching the NLO precision level. The project is based on the DGLAP evolutionequation [13] at the NLO level [14–17]. We began by creating two different algorithms,and then constructing two MC programs, that solve the evolution equation and simulatethe one-hemisphere collinear parton shower.
First algorithm , called
EvolFMC (MarkovianMonte Carlo), is based on the principle of a Markovian process. It generates the DGLAP-type evolution of Parton Distribution Functions (PDFs) up to the NLO level includingboth gluons and quarks. In Ref. [18] we have shown for the first time that with theaid of modern CPU such a MC code can solve the LO evolution equations with theprecision below 10 − , i.e. comparable or even better than the other numerical methods.This study has been extended to the NLO evolution in Ref. [19]. The second algorithm called CMC is an entirely new type of an algorithm, an example of a wider class that wenamed ”Constrained Markovian Monte Carlo”. CMC is designed to solve the Markovianevolution with superimposed constraint on both x and flavor type of the outgoing parton[20,21]. Such an evolution with predefined end-point, mandatory for the simulation of theresonant processes, is an alternative for the commonly used workaround solution of thisproblem based on the so-called ”backward evolution” [22]. The ultimate goal is to combinetwo initial-state non-collinear NLO constrained evolutions (in two hemispheres) togetherwith the hard process into one NLO parton shower for the Drell-Yan-type processes, eitherin full agreement with the DGLAP evolution or including effects of a modified-DGLAPtype.The EvolFMC code includes also the modified DGLAP evolutions in which the argumentof the coupling constant is replaced by more complicated functions of x , x (cid:48) and Q [23–27].In particular, the use of k T as the argument is known to effectively resum some of thehigher order soft corrections [23, 26, 28]. For this reason it is the preferred choice of mostof the MC parton shower codes. The k T is used as an argument also by the CCFMevolution equation [29]. This equation is designed to effectively interpolate between theDGLAP evolution in the large x region and the BFKL-type behavior in the low x region.The CCFM equation uses, apart from the modified arguments of coupling constant and1ngular ordering, also the ”non-Sudakov form factor”, important in the small x region.In the case of the developed by us EvolFMC code the option of modified argument ofthe coupling constant has been so far implemented only in the LO approximation [30, 31].This is a significant limitation since both the NLO kernels and the k T -dependent couplingconstant are important for the final precision of the parton shower. In the present paperwe will fill in this gap and present the NLO extension of the k T -dependent EvolFMC code.In this scheme the evolution variable is understood to be the rapidity of emitted partons,ensuring the angular ordering. We will discuss also another choice of the argument: Q (1 − z ). It is based on the variant of the CCFM equation, called the ”one-loop” CCFM[32, 33]. Of course, Q (1 − z ) is only one of a few z -dependent functions that have beenused in the literature [24, 25, 34].The primary role of the presented here Markovian NLO code would be to serve as apowerfull testing device capable of reproducing independently any of the above evolutiontypes. In particular, the emulation of the CCFM equation would be possible, if thecollinear evolution is supplied with the transverse momenta and the ”non-Sudakov formfactor” is added. Another area of application of such a Markovian MC would be inperforming fits of the PDFs to the F data. We have demonstrated recently [35] thatusing MC codes for the purpose of fitting is feasible with the modern computer power.Finally, for the simulation of the final-state parton shower, the Markovian algorithmwithout constraints is directly applicable.Let us digress here that the standard evolution kernels of the collinear factorisation,either LO or NLO, are used in the normal integrated form, i.e. as functions of the x and t variables only. This is of course the well-established and extremely successful DGLAPapproach. Nonetheless, for more exclusive observables at the LHC, it will be necessaryto include the effects of transverse momenta beyond the collinear approximation. Someattempts in this direction have already been made, for example by the above mentionedMC@NLO [11, 12], the GR@PPA project [36, 37] or the CCFM-based [29] CASCADEproject [38].Another interesting approach is based on the idea of ”exclusive” DGLAP kernels inwhich the internal degrees of freedom are not integrated out analytically, but insteadsimulated (i.e. integrated) by the MC parton shower itself [39].Finally, let us comment on yet another class of MC algorithms solving the evolutionequations – the veto algorithms [40–42]. They are based on a Markovian evolution withadditional internal pseudo-emission loop. Thanks to this feature it is enough to computean approximate Sudakov form factor, instead of the exact one, for each of the generatedevents and the algorithm becomes faster. On the negative side, the average multiplicity ofthe generated partons is much higher. We did not use the veto scheme in this work. Themain reason for this is that, as discussed earlier, we consider Markovian-type algorithmsprimarily as a testing device for the constrained Markovian algorithms. In these lattercases the veto algorithm does not apply and the Sudakov form factors must be calculated.For that reason we preferred to have the Sudakov form factors implemented (and mutuallytested) in both approaches.This paper is organized as follows. In the next Section we introduce all necessary2otation, briefly present the general formalism of the Markovian MC solutions of theevolution equations and define two particular evolution schemes that we are interestedin. In the following two Sections we describe in detail the MC algorithms that solve theabove two evolutions by means of the collinear Markovian parton shower. We present twoalgorithms, main and auxiliary (mostly for tests), for each of the schemes. As the simplerones we discuss in detail the Q (1 − z )-type algorithms, whereas in the description of the k T -dependent ones we focus mostly on the modifications with respect to the Q (1 − z )-type algorithms. In Section 5 we present numerical results: first we discuss the choiceof the counter term and input parameters, then we briefly describe a variety of technicaltests that we performed in order to establish the technical precision of the EvolFMC code,and finally we compare and discuss various types of evolution. The Summary Sectionconcludes the paper.
In this paper we will follow the general formulation of the evolution equation and itssolutions in terms of the Markovian MC presented in Ref. [42]. In order to establish thenotation let us recall here basic definitions and formulas. For more complete presentationwe refer the reader to the Ref. [42]. We write the evolution equation in the compactmatrix notation ∂ t D ( t ) = K ( t ) D ( t ) (1)and in the explicit representation ∂ t D f ( t, x ) = (cid:88) f (cid:48) (cid:90) x dw K ff (cid:48) ( t, x, w ) D f (cid:48) ( t, w ) , (2)where D f ( t, w ) denotes the parton density function of the parton f and K ff (cid:48) ( t, x, w )denotes the generalized evolution kernel. Note that contrary to the DGLAP case, itdepends on two variables, x and w , instead of their ratio z = x/w only. The kernel isthen decomposed into real and virtual parts K ff (cid:48) ( t, x, w ) = K Vff (cid:48) ( t, x, w ) + K Rff (cid:48) ( t, x, w ) , K Vff (cid:48) ( t, x, w ) = − δ ff (cid:48) δ x = w K vff ( t, x ) . (3)The Markovian solution of eq. (1) is conveniently expressed with the help of the operator¯ E defined as ¯ E D ( t ) ≡ (cid:90) dx (cid:88) f x D f ( t, x ) , i.e. ¯ E f ( x ) ≡ x. (4)This operator in a formal way shows that we use an ”unconstrained” evolution – allvalues of x and f will be generated and that the evolution will be done for momentumdistributions x D ( t ). The solution of eq. (1) and the master formula for the Markovian3C is then¯ ED ( t ) = (cid:90) tt dt (cid:18) (cid:90) tt dt (cid:20) (cid:90) tt dt (cid:26) . . .. . . (cid:90) tt N − dt N (cid:26) ¯ EK R ( t N ) G K V ( t N , t N − ) + ¯ EG K V ( t, t N − ) δ t N = t (cid:27) ... × K R ( t ) G K V ( t , t ) + ¯ EG K V ( t, t ) δ t = t (cid:21) × K R ( t ) G K V ( t , t ) + ¯ EG K V ( t, t ) δ t = t (cid:19) D ( t ) , (5)where the diagonal matrix G K V is the solution of the evolution equation with the virtualkernel K Vff (cid:48) ( t, x, w ) parametrized in terms of the Sudakov form factor Φ f ( t, t (cid:48) | x ): { G K V ( t, t (cid:48) ) } ff (cid:48) ( x, w ) = δ ff (cid:48) δ x = w e − Φ f ( t,t (cid:48) | w ) , Φ f ( t, t (cid:48) | x ) = (cid:90) tt (cid:48) dt (cid:48)(cid:48) K vff ( t (cid:48)(cid:48) , x ) . (6)The above formula (5) follows from the iteration of the momentum conservation principle t (cid:90) t i − dt i (cid:110) ¯ EK R ( t i ) G K V ( t i , t i − ) + ¯ EG K V ( t, t i − ) δ t i = t (cid:111) = ¯ E . (7)Eq. (7) defines also the probabilities of the single step forward in t i , f i and x i variables:1 = 1 x i − t (cid:90) t i − dt i (cid:110) ¯ EK R ( t i ) G K V ( t i , t i − ) + ¯ EG K V ( t, t i − ) δ t i = t (cid:111) f i − ( x i − )= e − Φ fi − ( t,t i − | x i − ) + (cid:90) e − Φ fi − t,ti − | xi − d (cid:16) e − Φ fi − ( t i ,t i − | x i − ) (cid:17) × (cid:34)(cid:88) f i ∂ t i Φ f i f i − ( t i , t i − | x i − ) ∂ t i Φ f i − ( t i , t i − | x i − ) (cid:90) dx i ∂ t i Φ f i f i − ( t i , t i − | x i − ) x i x i − K Rf i f i − ( t i , x i , x i − ) (cid:35) , (8)4here the virtual Sudakov form factor is expressed in terms of the real emission part ofthe evolution kernelΦ f i − ( t i , t i − | x i − ) = t i (cid:90) t i − dt K vf i − f i − ( t, x i − )= (cid:88) f i t i (cid:90) t i − dt x i − (cid:90) dx i x i − x i K Rf i f i − ( t, x i , x i − ) = (cid:88) f i Φ f i f i − ( t i , t i − | x i − ) . (9)In the case of a weighted algorithm with the global correcting weight, one uses the simpli-fied kernel K Rf i f i − ( t i , x i , x i − ) → ¯ K Rf i f i − ( t i , x i , x i − ). This simplification is compensatedby the global weight w ( n ) = e ¯Φ fn ( t,t n | x n ) − Φ fn ( t,t n | x n ) (cid:32) n (cid:89) i =1 K Rf i f i − ( t i , x i , x i − )¯ K Rf i f i − ( t i , x i , x i − ) e ¯Φ fi − ( t i ,t i − | x i − ) − Φ fi − ( t i ,t i − | x i − ) (cid:33) , (10)where n denotes the number of emissions generated within a given MC event, and theform factor ¯Φ f i − ( t i , t i − | x i − ) is constructed from ¯ K Rf i f i − ( t i , x i , x i − ) in complete analogyto eq. (9). The last quantity to be defined here is the exact shape of the kernels, includingthe definition of the argument of the coupling constant in terms of the t and x variables.Following Ref. [42] we will discuss two schemes of modified-DGLAP type, denoted inRef. [42] as (B’) and (C’). The novelty with respect to Ref. [42] is that we will performthe calculation and construct the Markovian MC code at the NLO level, whereas inRef. [42] only the LO case has been discussed. To be specific, the schemes are defined as x K R ( B (cid:48) ) f (cid:48) f ( t, x, w ) == (cid:18) α NLO (ln(1 − z ) + t )2 π zP R (0) f (cid:48) f ( z ) + (cid:16) α NLO (ln(1 − z ) + t )2 π (cid:17) zP R (1) B (cid:48) f (cid:48) f ( z ) (cid:19) θ − z>λe − t ,x K R ( C (cid:48) ) f (cid:48) f ( t, x, w ) == (cid:18) α NLO (ln( w − x ) + t )2 π zP R (0) f (cid:48) f ( z ) + (cid:16) α NLO (ln( w − x ) + t )2 π (cid:17) zP R (1) C (cid:48) f (cid:48) f ( z ) (cid:19) θ w − x>λe − t , (11)where z = x/w and λ is the cut-off on the argument of the coupling constant. Note that λ is not an infinitesimal IR cut-off, as used in the DGLAP case. On the contrary, λ isfinite, typically of the order of a few GeV. Note also that the factor 2 in front of thekernels is due to our definition of the evolution time, which in the DGLAP case is chosenas t = ln Q rather than t = ln Q .The NLO parts of the kernels, P R (1) B (cid:48) f (cid:48) f ( z ) and P R (1) C (cid:48) f (cid:48) f ( z ), consist of the universalpart P R (1) f (cid:48) f ( z ) [16, 17] and the evolution scheme dependent counter terms ∆ P R (1) B (cid:48) f (cid:48) f ( z ) and5 P R (1) C (cid:48) f (cid:48) f ( z ) P R (1) B (cid:48) f (cid:48) f ( z ) = P R (1) f (cid:48) f ( z ) + ∆ P R (1) B (cid:48) f (cid:48) f ( z ) ,P R (1) C (cid:48) f (cid:48) f ( z ) = P R (1) f (cid:48) f ( z ) + ∆ P R (1) C (cid:48) f (cid:48) f ( z ) . (12)These counter terms are necessary to remove the double counting introduced by the shiftof the arguments of the coupling constants. There is some freedom in defining thesecounter terms. The algorithms presented here work for any (reasonable) choice of thecounter terms. For further details on the choice we used in this work we refer the readerto Section 5.1.We will frequently be using also the representation of the universal parts of the kernelsbased on their structure in the z -variable zP R (0) f (cid:48) f ( z ) = 11 − z δ f (cid:48) f A (0) ff + F (0) f (cid:48) f ( z ) ,zP R (1) f (cid:48) f ( z ) = 11 − z δ f (cid:48) f A (1) ff ( z ) + F (1) f (cid:48) f ( z ) . (13)The functions P Rf (cid:48) f ( z ), A ff ( z ) and F f (cid:48) f ( z ), written in the notation used here, are collectedin Ref. [19].The coupling constant at the NLO level has the standard form α LO ( t ) = 2 πβ ( t − ln Λ ) ,α NLO ( t ) = α LO ( t ) (cid:18) − α LO ( t ) β ln(2 t − )4 πβ (cid:19) . (14)We close this general introduction with a brief explanation why we state that theargument ln( x i − − x i ) + t i can be regarded as k Ti of the emitted real parton with thefour momentum k i . It follows immediately from the kinematical mapping of the evolutionvariables into four momenta, provided that the evolution time is identified with the rapid-ity of the emitted parton and the x -variable, as usual, with the light-cone plus componentof the virtual parton t i = ln(2 E h ) −
12 ln k + i k − i , k + i = 2 E h ( x i − − x i ) (15)where E h is an arbitrary reference energy of the incoming hadron. As a consequence ofthe maslessness of k i we obtain k Ti = √ k + k − = e t i ( x i − − x i ) . (16)Let us now proceed with the description of the novel NLO algorithms for the twoschemes. 6 Markovian algorithm for scheme (B’)
In this section we will present two schemes of solving the evolution B’ in a Markovian wayat the NLO level. We begin with the efficient one and then we present the other scheme,devised mostly for the testing purposes.
The main algorithm is based on the simplified LO DGLAP kernel in which the couplingconstant is used in the NLO approximation. Specifically we follow the eq. (3.21) ofRef. [30], see also [43], and we extend it to the NLO case x ¯ K R ( B (cid:48) ) f (cid:48) f ( t, x, w ) ≡ α NLO ( t + ln(1 − z ))2 π z ¯ P R ( B (cid:48) ) f (cid:48) f ( z ) θ − z>λe − t , (17) z ¯ P R ( B (cid:48) ) f (cid:48) f ( z ) = 11 − z ( δ f (cid:48) f A (0) ff + max z F (0) f (cid:48) f ( z ) + M f (cid:48) f ) . (18)We remind the reader that z = x/w . Note that the singular factor 1 / (1 − z ) is artificiallyintroduced into the F -part in order to achieve analytical integrability of the formula, seeRef. [30] for more discussion. The constant M f (cid:48) f is defined as M f (cid:48) f = (cid:40) , if P R (0) f (cid:48) f ( z ) (cid:54) = 0 ,η, if P R (0) f (cid:48) f ( z ) = 0 , (19)where η is a dummy technical parameter. The reason behind introduction of M f (cid:48) f isvery simple – we want to remove all zeroes in the LO transition matrix, because some ofthem might become non-zero at the NLO level and cause infinite weights. The constant η is therefore added to all kernels that are zero at the LO level. Of course this is purelytechnical, dummy, operation, later on corrected by means of a proper rejection weight.The corresponding Sudakov form factor, necessary to generate the time t i , is thendefined as¯Φ B (cid:48) f ( t i , t i − | x i − ) = (cid:90) t i t i − dt (cid:90) d (cid:16) xw (cid:17) (cid:88) f (cid:48) x ¯ K R ( B (cid:48) ) f (cid:48) f ( t, x, w )= 1 π (cid:90) t i t i − dt (cid:90) t − t λ du α NLO ( t − u )( A (0) ff + (cid:88) f (cid:48) (max z F (0) f (cid:48) f ( z ) + M f (cid:48) f ))= (cid:90) t i t i − dt (cid:16) ρ ( t, t − t λ ) − ρ ( t, (cid:17)(cid:16) A (0) ff + (cid:88) f (cid:48) max z ( F (0) f (cid:48) f ( z ) + M f (cid:48) f ) (cid:17) =( ¯ ζ ( t i ) − ¯ ζ ( t i − ))( A (0) ff + (cid:88) f (cid:48) (max z F (0) f (cid:48) f ( z ) + M f (cid:48) f )) , (20)7here ρ ( t, u ) = 1 π (cid:90) du α NLO ( t − u )= 1 π (cid:90) du β G ( t, u ) (cid:18) − β ln(2 G ( t, u )) β G ( t, u ) (cid:19) = − (cid:18) β + β β G ( t, u ) (cid:19) ln(2 G ( t, u )) − β β G ( t, u ) G ( t, u ) = t − u − ln Λ , u = − ln(1 − z ) , t λ = ln λ, (21)and ¯ ζ ( t ) = a (ln ( t − ln Λ )) + ( a t + a ) ln ( t − ln Λ ) + a t, (22) a = 12 β β ,a = 2 β ,a = − β ln Λ − β ln 2 − β β ,a = − β + 2 β ( t λ − ln Λ ) β ( t λ − ln Λ ) ln( t λ − ln Λ ) − β ( t λ − ln Λ ) + β (1 + ln 2) β ( t λ − ln Λ ) . Let us note that the coefficients a ij need to be calculated only once during initializationof the algorithm. The procedure of inverting the function ζ ( t i ), necessary to generate thetime t i , has to be done numerically.In the next step we generate the flavor index f i based on the probability p f i = ∂ t i ¯Φ B (cid:48) f i f i − ( t i , t i − | x i − ) ∂ t i ¯Φ B (cid:48) f i − ( t i , t i − | x i − ) = δ f i f i − A (0) f i − f i − + max z F (0) f i f i − ( z ) + M f i f i − A (0) f i − f i − + (cid:80) f i (max z F (0) f i f i − ( z ) + M f i f i − ) , (cid:88) f i p f i = 1 , (23)where ∂ t i ¯Φ B (cid:48) f i f i − ( t i , t i − | x i − ) = (cid:90) d (cid:16) xw (cid:17) x ¯ K R ( B (cid:48) ) f i f i − ( t i , x, w )= (cid:16) ρ ( t i , t i − t λ ) − ρ ( t i , (cid:17)(cid:16) δ f i f i − A (0) f i − f i − + max z F (0) f i f i − ( z ) + M f i f i − (cid:17) (24)and ∂ t i ¯Φ B (cid:48) f i − ( t i , t i − | x i − ) = (cid:88) f i ∂ t i ¯Φ B (cid:48) f i f i − ( t i , t i − | x i − ) . (25)8s the last variable we generate z i . The normalized density distribution dz i p B (cid:48) ( z i ) isgiven by dz i p B (cid:48) ( z i ) = dz i α NLO ( t i + ln(1 − z i )) π Θ − z i >λe − ti − z i (cid:18)(cid:90) dz α NLO ( t i + ln(1 − z )) π Θ − z>λe − ti − z (cid:19) − = du i α NLO ( t i − u i ) π Θ t i − t λ >u i > (cid:16) ρ ( t i , t i − t λ ) − ρ ( t i , (cid:17) − = du i ∂ u i ρ ( t i , u i )Θ t i − t λ >u i > (cid:16) ρ ( t i , t i − t λ ) − ρ ( t i , (cid:17) − . (26)In the actual generation of z i we use the method of inverse cumulative. The ρ ( t i , u i )function has to be inverted numerically with respect to the u i variable.The last part of the algorithm to be discussed here is the correcting weight, definedin eq. (10), compensating the simplification of the kernel done in eq. (17). The mostcomplicated part of this weight is related to the exact Sudakov form factor. It has theform of the double integral. Numerical evaluation of such a double integral would signif-icantly slow down the algorithm. Therefore it is essential to perform at least one of theintegrations analytically. In the following we will show how this can be done in the NLOcase.Let us define the full Sudakov form factor Φ B (cid:48) f ( t i , t i − | x i − ) of the B’ evolutionΦ B (cid:48) f ( t i , t i − | x i − ) = (cid:90) t i t i − dt (cid:90) d (cid:16) xx i − (cid:17) (cid:88) f (cid:48) x K R ( B (cid:48) ) f (cid:48) f ( t, x, x i − )= (cid:90) t i t i − dt (cid:90) dz (cid:88) f (cid:48) (cid:18) α S ( t + ln(1 − z ))2 π zP R (0) f (cid:48) f ( z )+ (cid:16) α S ( t + ln(1 − z ))2 π (cid:17) zP R (1) B (cid:48) f (cid:48) f ( z ) (cid:19) θ − z>λe − t , (27)and let us write down the desired virtual component of the weight∆ B (cid:48) f ( t i , t i − | x i − ) =Φ B (cid:48) f ( t i , t i − | x i − ) − ¯Φ B (cid:48) f ( t i , t i − | x i − )= (cid:90) t i − t i dt (cid:90) dzθ − z>λe − t (cid:88) f (cid:48) (cid:18) α S ( t + ln(1 − z ))2 π z (cid:0) P R (0) f (cid:48) f ( z ) − ¯ P R (0) f (cid:48) f ( z ) (cid:1) + (cid:16) α S ( t + ln(1 − z ))2 π (cid:17) zP R (1) B (cid:48) f (cid:48) f ( z ) (cid:19) . (28)The integral over dz for general form of the kernel cannot be performed analytically evenat the LO level (cf. Ref. [30]). However, as in the LO case of Ref. [30], the dt integralcan be done analytically also for the NLO case. The calculation looks as follows. Weintroduce the usual variable u = − ln(1 − z ) and then change order of integration. Theresulting integral can be expressed as a sum of two integrals over the two regions of the9 II t i-1 t i t ! t i-1 - t ! t i - tu t ! u = t - Figure 1:
The integration domain in the tu space for the case B’. Regions I and II correspondto two integrals in eq. (29). tu space shown in Fig. 1 : (cid:90) t i t i − dt (cid:90) dzθ − z>λe − t = (cid:90) t i t i − dt (cid:90) t − t λ du e − u = (cid:90) t i − t λ t i − − t λ due − u (cid:90) t i u + t λ dt + (cid:90) t i − − t λ due − u (cid:90) t i t i − dt. (29)The dt integral, which depends on the coupling constant only, can be done analyticallywith the help of the integrals1 π J ( t, u ) = 1 π (cid:90) dtα NLO ( t − u ) = − ρ ( t, u ) (30)and 1 π J ( t, u ) = 1 π (cid:90) dtα NLO ( t − u )= (cid:90) dt (cid:0) − β ( t − u − ln Λ ) + β ln 2 ( t − u − ln Λ ) (cid:1) β ( t − u − ln Λ ) = a G ( t, u ) ln (2 G ( t, u )) + (cid:18) a G ( t, u ) + a G ( t, u ) (cid:19) ln(2 G ( t, u ))+ a G ( t, u ) + a G ( t, u ) + a G ( t, u ) , (31) The similar change of the integration order was also exploited by the authors of HERWIG MC [9,44]. a = − β β , a = − β β , a = 4 β β , a = − β β , a = β β , a = − β . Collecting together the above results we can rewrite eq. (28) as∆ B (cid:48) f ( t i , t i − | x i − ) == (cid:90) t i − t λ t i − − t λ due − u z (cid:88) f (cid:48) (cid:18) ( J ( t i , u ) − J ( u + t λ , u ))2 π (cid:0) P R (0) f (cid:48) f ( z ) − ¯ P R ( B (cid:48) ) f (cid:48) f ( z ) (cid:1) + ( J ( t i , u ) − J ( u + t λ , u ))4 π P R (1) B (cid:48) f (cid:48) f ( z ) (cid:19) + (cid:90) t i − − t λ due − u z (cid:88) f (cid:48) (cid:18) ( J ( t i , u ) − J ( t i − , u ))2 π (cid:0) P R (0) f (cid:48) f ( z ) − ¯ P R ( B (cid:48) ) f (cid:48) f ( z ) (cid:1) + ( J ( t i , u ) − J ( t i − , u ))4 π P R (1) B (cid:48) f (cid:48) f ( z ) (cid:19) (32)where z = 1 − e − u . Note the cancellation of the leading singularity, A (0) ff / (1 − z ), in theabove formula due to (cid:88) f (cid:48) z (cid:16) P R (0) f (cid:48) f ( z ) − ¯ P R ( B (cid:48) ) f (cid:48) f ( z ) (cid:17) = (cid:88) f (cid:48) (cid:32) F (0) f (cid:48) f ( z ) − − z (cid:18) max z F (0) f (cid:48) f ( z ) + M f (cid:48) f (cid:19)(cid:33) . (33)Combining together the real and virtual components we arrive at the final formula for theglobal weight of eq. (10) adopted to the case of the scheme B’ w ( n )( B (cid:48) ) = e − ∆ B (cid:48) fn ( t,t n | x n ) n (cid:89) i =1 (cid:32) K R ( B (cid:48) ) f i f i − ( t i , x i , x i − )¯ K R ( B (cid:48) ) f i f i − ( t i , x i , x i − ) e − ∆ B (cid:48) fi − ( t i ,t i − | x i − ) (cid:33) . (34) The main algorithm described in previous section is a fairly complicated one, especiallydue to the presence of numerical integrations and numerical inversions of various functions.Therefore, it is obligatory to devise a way of testing it down to the sub-per-mille precisionlevel. We have not found any non-Monte-Carlo program that solves the modified-DGLAP-type evolutions at the NLO level and therefore we constructed another, independentMonte Carlo algorithm for the purpose of cross-checks. The algorithm is less efficient butat the same time it is simpler.The algorithm is closely based on the LO algorithm described in Ref. [30]. The entireNLO correction is introduced as a weight. We will only briefly sketch it here and we referthe interested reader to Ref. [30] for details.11he algorithm is based on the simplified kernel of the form x ¯ K R ( B (cid:48) )(2) f (cid:48) f ( t, x, w ) ≡ α LO ( t + ln(1 − z ))2 π z ¯ P R ( B (cid:48) ) f (cid:48) f ( z ) θ − z>λe − t . (35)As compared to simplified kernel (17) of the previous algorithm in this one the couplingconstant is taken in the LO approximation.In complete analogy to eq. (20) the simplified Sudakov form factor reads¯Φ B (cid:48) (2) f ( t i , t i − | x i − ) = (cid:90) t i t i − dt (cid:90) d (cid:16) xx i − (cid:17) (cid:88) f (cid:48) x ¯ K R ( B (cid:48) )(2) f (cid:48) f ( t, x, x i − )= 1 π (cid:90) t i t i − dt (cid:90) t − t λ duα LO ( t − u )( A (0) ff + (cid:88) f (cid:48) (max z F (0) f (cid:48) f ( z ) + M f (cid:48) f ))= (cid:90) t i t i − dt (cid:16) ρ LO ( t, t − t λ ) − ρ LO ( t, (cid:17)(cid:16) A (0) ff + (cid:88) f (cid:48) (max z F (0) f (cid:48) f ( z ) + M f (cid:48) f ) (cid:17) =( ¯ ζ LO ( t i ) − ¯ ζ LO ( t i − ))( A (0) ff + (cid:88) f (cid:48) (max z F (0) f (cid:48) f ( z ) + M f (cid:48) f )) , (36)where ρ LO ( t, u ) = 1 π (cid:90) du α LO ( t − u ) = − β ln(2 G ( t, u )) . (37)The function ¯ ζ LO ( t ) reads¯ ζ LO ( t i ) = 2 β ( t i − ln Λ )(ln( t i − ln Λ ) −
1) (38)In order to generate t i the function ¯ ζ LO ( t i ) must be inverted. This inversion can bedone analytically. However, in the actual Monte Carlo, for the purpose of tests we haveimplemented also the numerical inversion procedure, similar as in the case of the firstalgorithm (there is no visible computing time overhead related to this numerical inversion).The generation of the flavor index f i is based on the probability identical to eq. (23) p (2) f i = ∂ t i ¯Φ B (cid:48) (2) f i f i − ( t i , t i − ) ∂ t i ¯Φ B (cid:48) (2) f i − ( t i , t i − ) = δ f i f i − A (0) f i − f i − + max z F (0) f i f i − ( z ) + M f i f i − A (0) f i − f i − + (cid:80) f i (cid:18) max z F (0) f i f i − ( z ) + M f i f i − (cid:19) ≡ p f i . (39)It is the case because the coupling constants cancel in both eqs. (23) and (39). Thegeneration of the z -variable is based on the LO version of eq. (26) dz i p B (cid:48) (2) ( z i ) = dz i α LO ( t i + ln(1 − z i )) π Θ − z i >λe − ti − z i (cid:18)(cid:90) dz α LO ( t i + ln(1 − z )) π Θ − z>λe − ti − z (cid:19) − = du i ∂ u i ρ LO ( t i , u i )Θ t i − t λ >u i > (cid:16) ρ LO ( t i , t i − t λ ) − ρ LO ( t i , (cid:17) − . (40)12ontrary to the previous algorithm, now the function ρ LO ( t, u ) can be inverted analyti-cally. However, for the purpose of testing various components of the program, we haveimplemented also an option of the numerical inversion.As indicated earlier, the novelty of this algorithm, i.e. the NLO effect, is hidden in theglobal weight. The most complicated part of the weight is the virtual component builtout of the Sudakov form factors. For the calculation of the exact form factor we can usethe results derived for the previous algorithm, and the whole virtual part of the weightfollows from eq. (28)∆ B (cid:48) (2) f ( t i , t i − | x i − ) = Φ B (cid:48) f ( t i , t i − | x i − ) − ¯Φ B (cid:48) (2) f ( t i , t i − | x i − )= (cid:90) t i − t i dt (cid:90) dzθ − z>λe − t (cid:88) f (cid:48) (cid:18) α NLO ( t + ln(1 − z ))2 π zP R (0) f (cid:48) f ( z )+ (cid:16) α NLO ( t + ln(1 − z ))2 π (cid:17) zP R (1) B (cid:48) f (cid:48) f ( z ) − α LO ( t + ln(1 − z ))2 π z ¯ P R (0) f (cid:48) f ( z ) (cid:19) , (41)with the final result∆ B (cid:48) (2) f ( t i , t i − | x i − ) == (cid:90) t i − t λ t i − − t λ due − u z (cid:88) f (cid:48) (cid:18) − ( J LO ( t i , u ) − J LO ( u + t λ , u ))2 π ¯ P R ( B (cid:48) ) f (cid:48) f ( z )+ ( J ( t i , u ) − J ( u + t λ , u ))2 π P R (0) f (cid:48) f ( z ) + ( J ( t i , u ) − J ( u + t λ , u ))4 π P R (1) B (cid:48) f (cid:48) f ( z ) (cid:19) + (cid:90) t i − − t λ due − u z (cid:88) f (cid:48) (cid:18) − ( J LO ( t i , u ) − J LO ( t i − , u ))2 π ¯ P R ( B (cid:48) ) f (cid:48) f ( z )+ ( J ( t i , u ) − J ( t i − , u ))2 π P R (0) f (cid:48) f ( z ) + ( J ( t i , u ) − J ( t i − , u ))4 π P R (1) B (cid:48) f (cid:48) f ( z ) (cid:19) , (42)where z = 1 − e − u and1 π J LO ( t, u ) = 1 π (cid:90) dtα LO ( t − u ) = − ρ LO ( t, u ) . (43)We conclude this section by presenting the complete formula for the global weight inthis algorithm: w ( n )( B (cid:48) )(2) = e − ∆ B (cid:48) (2) fn ( t,t n | x n ) n (cid:89) i =1 (cid:32) K R ( B (cid:48) ) f i f i − ( t i , x i , x i − )¯ K R ( B (cid:48) )(2) f i f i − ( t i , x i , x i − ) e − ∆ B (cid:48) (2) fi − ( t i ,t i − | x i − ) (cid:33) . (44) Having completed the presentation of the algorithms for the scheme B’ we proceed nowto the most important scheme C’. As discussed in Ref. [42], it can be interpreted as13volution in the rapidity variable with k T as the argument of the coupling constant and itis of great physical importance. The NLO algorithms for the scheme C’ are quite similarto the algorithms for the B’ scheme. Therefore in the following sections we will skip someof the details and concentrate on the differences with respect to the scheme B’. As beforewe begin with the efficient one and then we present the other algorithm used for tests. The main algorithm is based on the simplified kernel similar to the one used in the schemeB’ (i.e. the LO kernel with the NLO coupling constant): x ¯ K R ( C (cid:48) ) f (cid:48) f ( t, x, w ) ≡ α NLO ( t + ln w + ln(1 − z ))2 π z ¯ P R ( C (cid:48) ) f (cid:48) f ( z ) θ (1 − z ) w>λe − t , (45) z ¯ P R ( C (cid:48) ) f (cid:48) f ( z ) ≡ z ¯ P R ( B (cid:48) ) f (cid:48) f ( z ) = 11 − z ( δ f (cid:48) f A (0) f (cid:48) f + max z F (0) f (cid:48) f ( z ) + M f (cid:48) f ) . (46)Note that now the kernel depends truly on both the x and w variables (through the θ -function). In the B’ case it depended only on the ratio z = x/w . The simplified Sudakovform factor, necessary to generate time t i , is then defined as¯Φ C (cid:48) f ( t i , t i − | x i − ) = (cid:90) t i t i − dt (cid:90) d (cid:18) xx i − (cid:19) (cid:88) f (cid:48) x ¯ K R ( C (cid:48) ) f (cid:48) f ( t, x, x i − )= 1 π (cid:90) t i t i − dt (cid:90) t +ln x i − − t λ du α NLO ( t + ln x i − − u ) (cid:16) A (0) ff + (cid:88) f (cid:48) (max z F (0) f (cid:48) f ( z ) + M f (cid:48) f ) (cid:17) = (cid:90) t i t i − dt (cid:16) ρ ( t + ln x i − , t + ln x i − − t λ ) − ρ ( t + ln x i − , (cid:17)(cid:16) A (0) ff + (cid:88) f (cid:48) (max z F (0) f (cid:48) f ( z ) + M f (cid:48) f ) (cid:17) = (cid:16) ¯ ζ ( t i + ln x i − ) − ¯ ζ ( t i − + ln x i − ) (cid:17)(cid:16) A (0) ff + (cid:88) f (cid:48) (max z F (0) f (cid:48) f ( z ) + M f (cid:48) f ) (cid:17) . (47)The functions ρ and ¯ ζ are the same as in the algorithm B’. The only difference is in theshifted t -argument: t → t + ln x i − or, equivalently, shifted integration/generation limits t i ( i − → t i ( i − + ln x i − .The generation of the flavor index f i is done, identically as in the previous algorithms,by means of the probabilities p f i defined in eqs. (23) or (39).As the last variable we generate z i using the integrand of the density function dz i p C (cid:48) ( z i )14iven by the formula dz i p C (cid:48) ( z i ) = dz i α NLO ( t i + ln x i − + ln(1 − z i )) π Θ (1 − z i ) x i − >λe − ti − z i × (cid:18)(cid:90) dz α NLO ( t i + ln x i − + ln(1 − z )) π Θ (1 − z ) x i − >λe − ti − z (cid:19) − = du i α NLO ( t i + ln x i − − u i ) π Θ t i +ln x i − − t λ >u i > × (cid:16) ρ ( t i + ln x i − , t i + ln x i − − t λ ) − ρ ( t i + ln x i − , (cid:17) − = du i ∂ u i ρ ( t i + ln x i − , u i )Θ t i +ln x i − − t λ >u i > × (cid:16) ρ ( t i + ln x i − , t i + ln x i − − t λ ) − ρ ( t i + ln x i − , (cid:17) − . (48)As in the case of t i -variable we observe that the whole difference with respect to the caseB’ is in the shift t i → t i + ln x i − .Let us define now the full Sudakov form factor Φ ( C (cid:48) ) f ( t i , t i − | x i − )Φ ( C (cid:48) ) f ( t i , t i − | x i − ) = (cid:90) t i t i − dt (cid:90) d (cid:16) xx i − (cid:17) (cid:88) f (cid:48) x K R ( C (cid:48) ) f (cid:48) f ( t, x, x i − ) (49)As in the case B’ we are able to perform analytically only one of the integrations in Φ ( C (cid:48) ) f .To this end we rearrange order of integrations and decompose the integrals as follows (cid:90) t i t i − dt (cid:90) dzθ (1 − z ) x i − >λe − t = (cid:90) t i t i − dt (cid:90) t +ln x i − − t λ du e − u = θ t λ − ln x i − The integration domain in the tu space for the case C’. Regions I (triangle) andII (rectangle) correspond to two integrals in eq. (51), respectively. Left frame: the case of t λ − ln x i − < t i − . Right frame: the case of t λ − ln x i − (cid:62) t i − . calculate the virtual component∆ C (cid:48) f ( t i , t i − | x i − ) = Φ C (cid:48) f ( t i , t i − | x i − ) − ¯Φ C (cid:48) f ( t i , t i − | x i − )= (cid:90) t i − t i dt (cid:90) dzθ (1 − z ) x i − >λe − t (cid:88) f (cid:48) (cid:18) α S ( t + ln x i − + ln(1 − z ))2 π z (cid:0) P R (0) f (cid:48) f ( z ) − ¯ P R (0) f (cid:48) f ( z ) (cid:1) + (cid:16) α S ( t + ln x i − + ln(1 − z ))2 π (cid:17) zP R (1) C (cid:48) f (cid:48) f ( z ) (cid:19) = (cid:90) t i +ln x i − − t λ t i − +ln x i − − t λ due − u z (cid:88) f (cid:48) (cid:32) J ( t i + ln x i − , u ) − J ( u + t λ , u )2 π (cid:0) P R (0) f (cid:48) f ( z ) − ¯ P R ( C (cid:48) ) f (cid:48) f ( z ) (cid:1) + J ( t i + ln x i − , u ) − J ( u + t λ , u )4 π P R (1) C (cid:48) f (cid:48) f ( z ) (cid:33) + θ t λ − ln x i − EvolFMC version 2. The first version of this program has been presented in the Ref. [18] for thecase of standard DGLAP LO evolution. Subsequently, the NLO evolution has been addedin Ref. [19] and the LO modified-DGLAP evolution in Ref. [30]. The presented here EvolFMC v.2 includes three of the described above four algorithms for NLO modified-DGLAP evolution (the auxiliary algorithm for the B’ evolution is not implemented at themoment). In order to accommodate these new evolution schemes, the overall structure18f the program has been modified, see Ref. [39] for details. As a result, it was importantto perform a number of technical comparisons of the new code in order to establish itstechnical precisions. We very briefly describe these tests in the following subsection.Later on we present numerical results regarding the actual new evolution schemes. Beforeshowing the results we discuss the choice of the counter terms and we list the inputparameters. As indicated in Section 2, modification of the argument of the coupling constant in theevolution kernels is equivalent to adding some higher order (i.e. NLO and higher) terms.Therefore, one has to make sure that there is no double counting with the NLO part ofthe kernel. In the implementation presented here we follow the approach as discussed inRef. [34]. Namely, we take the expansion of the α NLO ( t ) in the form α NLO ( t + ln φ ) = α NLO ( t ) − β π α NLO ( t ) ln φ + O (1 /t ) , (61)where φ represents any arbitrary change in the argument. Then the extra term in thekernels is of the form − β ln φ (cid:16) α NLO ( t )2 π (cid:17) zP R (0) f (cid:48) f ( z ) (cid:39) − β ln φ (cid:16) α NLO ( t + ln φ )2 π (cid:17) zP R (0) f (cid:48) f ( z ) + O (1 /t ) , (62)and consequently the counter terms from eqs. (12) could be defined as follows∆ P R (1) B (cid:48) f (cid:48) f ( z ) = β ln(1 − z ) P R (0) f (cid:48) f ( z ) , for scheme B’ , (63)¯∆ P R (1) C (cid:48) f (cid:48) f ( z ) = β (cid:0) ln w + ln(1 − z ) (cid:1) P R (0) f (cid:48) f ( z ) , z = x/w, for scheme C’. (64)However, a more detailed inspection of eq. (64) reveals that this formula over-subtracts the double counting. Namely, in the DGLAP evolution the kernels are functions of “local” z -variables only. In eq. (64) there is a corresponding ln(1 − z ) term. The other term, theln w , interconnects two emissions, and as such it is of genuinely beyond-DGLAP origin.It means that it is absent in the DGLAP kernel and there is no reason to remove it. Asa result, the counter term in the C’ case becomes identical to the counter term in the B’case. In the following we will present results for both variants in the C’ case, with theunderstanding that the preferred choice is the∆ P R (1) B (cid:48) f (cid:48) f ( z ) = ∆ P R (1) C (cid:48) f (cid:48) f ( z ) = β ln(1 − z ) P R (0) f (cid:48) f ( z ) , (65)for both schemes B’ and C’. 19 .2 Input parameters The set-up of the EvolFMC code is the same for all the presented results. We use the gluon( G ) and quark singlet ( Q ) PDFs with three massless quarks D Q = (cid:88) i (cid:16) D q i + D ¯ q i (cid:17) . (66)As the initial distributions of the evolution we take D G ( x ) = 1 . · x − . (1 − x ) . ,D sea ( x ) = 0 . · x − . (1 − x ) . ,D u val ( x ) = 2 . · x − . (1 − x ) . ,D d val ( x ) = 1 . · x − . (1 − x ) . (67)and D u ( x ) = D u val ( x ) + 16 D sea ( x ) ,D d ( x ) = D d val ( x ) + 16 D sea ( x ) ,D s ( x ) = D u ( x ) = D d ( x ) = D s ( x ) = 16 D sea ( x ) ,D Q ( x ) = D sea ( x ) + D u val ( x ) + D d val ( x ) . (68)The QCD constant Λ = 0 . λ = 1 GeV, N f = 3 and the dummy parameter η = 0 . We have performed three different sets of the technical comparisons of the code EvolFMCv.2 :1. With the semianalytical code APCheb40 [42,45] based on the expansion in Chebyshevpolynomials.2. With the previous version of the EvolFMC code: the EvolFMC v.1 , which was ex-tensively tested in the past.3. Between different algorithms within the EvolFMC v.2 . It is the most important testof the new NLO modified-DGLAP evolutions, which are not available in any othercode.The overall conclusion of all the tests is that the technical precision of the program EvolFMC v.2 is at least 5 × − (half of a per mille). For the details of the tests we referthe reader to Ref. [39]. 20 .4 Comparison of evolutions Having established the technical precision of the EvolFMC code let us now proceed withthe comparison of the various evolutions discussed earlier. Before getting into details letus remind the reader that in this paper for all of the results (i.e. for all of the evolutions)we use the same parameter setup: the same initial distributions at Q = 1 GeV and thesame Λ . In a more realistic study one should perform fits to the experimental datafor each of the evolutions separately and then use different initial distributions for eachevolution.We organize the numerical comparisons in the following way: We begin by showingthe ”reference” result, i.e. the standard DGLAP in the LO and NLO approximations(Fig. 3). Next, we show the new results for the B’-type and C’-type evolutions. Thesetwo plots (Figs. 4 and 5) are the main numerical results of this paper, showing the NLOcorrections in the modified DGLAP evolutions. The next two plots (Figs. 6 and 7) areof technical character and show in a more detail certain aspects of the C’-type evolution.We conclude the section by comparing in a common plot all three evolutions in the NLOapproximation (Fig. 8). . In the Fig. 3 we show the case of the standardDGLAP. We show these well known results because DGLAP will serve us as a referencepoint in discussing the modified evolutions of the B’- and C’-type. We present the gluonand quark momentum distributions in LO and NLO approximations as well as their ratios.Three evolution time limits are shown: 10, 100 and 1000 GeV. The characteristic featureof the plots is that DGLAP NLO corrections are systematically bigger in the large x region and they show the tendency of diverging in the x → x regionthe NLO corrections are small. Results for the other evolutions will be presented in asimilar way. The B’-type evolution is shown in the Fig. 4. It is one of the two main new numericalresults of this paper. As compared to the DGLAP case we can notice that the NLOcorrections are a bit bigger in the small x region, of the size up to 20%. In the large x region, on the contrary, the corrections are much smaller and showing less divergentbehavior. This is in agreement with the general principle of discussed here modificationsof the DGLAP equation. These modifications are supposed to improve the description ofthe emission of soft partons [23, 26, 28] and it is the x → The C’-type evolution. In the Fig. 6 we compare the LO and NLO evolutions in thecase of modified-DGLAP C’-type. We consider this figure as the most important new numerical result of this paper. We show the gluon and quark momentum distributions fortwo evolution time limits: 10, 100 and 1000 GeV. We see that at 100 and 1000 GeV theNLO corrections in the small- x region are somewhat bigger than in the case B’, reachingeven 30%. In the large x region the NLO corrections seem to be even milder and slightlyless divergent than in the B’ case, showing the same improvement over DGLAP. For the21 x) log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O Figure 3: Left frames: the DGLAP-type evolutions in the LO and NLO approximations. Uppercurves (LO: magenta and NLO: blue): the gluon xD G ( x ) distr.; lower curves (LO: black andNLO: red): the quark xD Q ( x ) distr. Right frames: the ratio of the NLO to LO distributionsfor the gluon (magenta) and quark (black) distributions. Top frames: the evolution up to 10GeV. Bottom frames: the evolution up to 100 GeV. shorter time of 10 GeV the effects are much smaller. In fact, in this evolution type, dueto the cut-off k T > λ , much less of the evolution happens before 10 GeV, and both theLO and NLO curves are close to the initial condition as well as to each other. Note thatin the case of both DGLAP and B’, the evolution at 10 GeV is already well developed. Inorder to reduce the evolution and get closer to the initial condition, in the DGLAP andB’ cases the evolution time must be much shorter, below 2 GeV at least. Technical details related to C’-type evolution: The C’-type evolution with the disfavoured counter term ¯∆ P from eq. (64). For the sake ofcomparison, in the Fig. 6, we present also the other, disfavored, choice of the counter termin the C’ evolution, given in eq. (64). The plots clearly confirm that the NLO correctionsare much bigger, and, in addition, strongly divergent in the small x limit. The C’-type evolution with no counter term. As the last exercise, in Fig. 7, we completely22 x) log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O Figure 4: Left frames: the modified-DGLAP B’-type evolutions in the LO and NLO approxi-mations. Upper curves (LO: magenta and NLO: blue): the gluon xD G ( x ) distr.; lower curves(LO: black and NLO: red): the quark xD Q ( x ) distr. Right frames: the ratio of the NLO toLO distributions for the gluon (magenta) and quark (black) distributions. Top frames: theevolution up to 10 GeV. Bottom frames: the evolution up to 100 GeV. switched off the counter term in the C’-type evolution. This way we try to understandbetter what actually causes the changes in the shape of the NLO corrections in themodified schemes: is it the genuine change of the argument of the coupling constant orrather the cut-off λ ? In a crude approximation the counter term can be regarded as anestimate of the size of the pure effect of shifting the argument in the coupling constant.As one can see, the shapes in the right hand side plots in Fig. 7 have changed significantlyin the large x region, and have remained similar in the small x region, as compared tothe complete C’ plots from Fig. 5. This demonstrates that indeed it is the change of theargument that drives the effect in the region of soft emission. On the other hand, fromthe comparison of Figs. 5 and 7 to the DGLAP evolution, Fig. 3, we inferr that the higherorder terms, resumed in the coupling constant, contribute as much as the counter term. Remark on the small x limit of the C’-evolution. As we have mentioned, C’ evolution23 x) log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O Figure 5: Left frames: the modified-DGLAP C’-type evolutions in the LO and NLO approxi-mations. Upper curves (LO: magenta and NLO: blue): the gluon xD G ( x ) distr.; lower curves(LO: black and NLO: red): the quark xD Q ( x ) distr. Right frames: the ratio of the NLO toLO distributions for the gluon (magenta) and quark (black) distributions. Top frames: theevolution up to 10 GeV. Central frames: the evolution up to 100 GeV. Bottom frames: theevolution up to 1000 GeV. is motivated by the CCFM equation. This equation treats the x → x evolution. It depends on the transverse momenta of all the emitted partons.It is in principle not difficult to include such an effect into the Monte Carlo evolution ofthe C’-type by means of the rejection mechanism, provided the transverse momenta areproperly generated within the cascade, and we plan to do it in the future [39].Let us summarize the comparisons: 1. The modifications of the argument of thecoupling constant decrease the size of the NLO corrections and reduce the divergences24 x) log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O Figure 6: Left frames: the modified-DGLAP C’-type evolutions in the LO and NLO approx-imations with the modified disfavoured (C’-type) counter term . Upper curves (LO: magentaand NLO: blue): the gluon xD G ( x ) distr.; lower curves (LO: black and NLO: red): the quark xD Q ( x ) distr. Right frames: the ratio of the NLO to LO distributions for the gluon (magenta)and quark (black) distributions. Top frames: the evolution up to 10 GeV. Central frames: theevolution up to 100 GeV. Bottom frames: the evolution up to 1000 GeV. in the large x region. 2. The NLO corrections in schemes B’ and C’ are significant, upto 30% in small x region and must be included in the MC parton showers for the LHC.3. One must remember, however, that these comparisons are to some degree artificialbecause, as discussed earlier, we use the same input PDF distributions for all evolutiontypes. In a more complete study each of the evolutions should be separately fitted tothe experimental data and obtained this way different input PDF distributions should beused. 25 x) log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V f x D log-3 -2.5 -2 -1.5 -1 -0.5 0 E v o l F M C . v2 N L O / L O Figure 7: Left frames: the modified-DGLAP C’-type evolutions in the LO and NLO approx-imations without the counter term. Upper curves (LO: magenta and NLO: blue): the gluon xD G ( x ) distr.; lower curves (LO: black and NLO: red): the quark xD Q ( x ) distr. Right frames:the ratio of the NLO to LO distributions for the gluon (magenta) and quark (black) distribu-tions. Top frames: the evolution up to 10 GeV. Central frames: the evolution up to 100 GeV.Bottom frames: the evolution up to 1000 GeV. As the last exercise we compare the three analyzed types of the evolution. In the Fig.8 we show simultaneously all three evolutions in the NLO approximation. It is clearlyvisible that the difference between DGLAP and B’-type evolution (with α ( Q (1 − z ))) israther small and of the quantitative form. On the contrary, the C’-type evolution (with α ( k T )) looks very different, both in the shape and in the magnitude. There is a visibleflattening of the C’-type distributions around the value of x ∼ λ/Q caused by the cut-off λ . This follows from the condition x i − (1 − z i ) ≥ λ/Q which stops any cascade as soonas the x i variable falls below the λ/Q threshold. In particular no evolution at all candevelop for any x ≤ λ/Q . 26 x) log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V g x D log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V q x D log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V g x D log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V q x D log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V g x D log-3 -2.5 -2 -1.5 -1 -0.5 0 ( x ) , Q = G e V q x D Figure 8: Comparison of three NLO evolutions: DGLAP (magenta), B’-type (blue) and C’-type (black). Left frames: the gluon xD G ( x ) distr. Right frames: the quark xD Q ( x ) distr.Top frames: the evolution up to 10 GeV. Central frames: the evolution up to 100 GeV.Bottom frames: the evolution up to 1000 GeV. In this paper we have presented a series of Markovian MC algorithms that solve theQCD evolutions of the modified-DGLAP type in the NLO approximation. One of thetwo discussed modifications of the DGLAP evolution is of high practical importance. Inthis evolution, as the argument of the coupling constant the transverse momentum of theemitted parton is used, and the evolution time is identified with the rapidity variable.Such a modification is known to describe better the emission of soft partons. It can serveas a first step towards incorporating the complete CCFM effects into the evolution as well.In this paper we have called this scenario the C’-type evolution. The other scenario, calledthroughout the paper the B’-type evolution, uses Q (1 − z ) as the argument of the couplingconstant, modelling the ”one-loop” CCFM evolution equation. It is one of many possible z -dependent modifications of the argument discussed in the literature. Proper counter27erms have been added to the evolution kernels in order to remove double counting at theNLO level.The algorithms have been implemented into the new version of the MC program EvolFMC and extensively tested by comparisons with the semianalytical code APCheb40 ,with the previous version of EvolFMC and by comparisons of independent algorithmswithin the new version of EvolFMC itself. As the overall conclusion of the tests we claimthe technical precision of the EvolFMC to be at least 5 × − .The comparison of the modified-DGLAP evolutions at both LO and NLO level showsthat the NLO corrections are in general smaller in the modified evolutions and moreimportantly, the divergent behavior of the NLO corrections in the large x region is limited,as expected. The only exception is the very low x region where the NLO corrections inthe modified schemes are larger. This is however the region where the DGLAP equationbecomes less accurate anyway.Quantitatively, the NLO corrections are relatively modest: in the B’ case they aresmall, of the order of up to 20% of the LO terms. For the k T -dependent evolution (C’)they are, in most of the parameter space, of the order of 10%, but in some limitedregions of small x they grow up to 30%. These results, on the one hand, show thatthe NLO contributions are numerically significant and should be taken into account inthe construction of the parton shower MC for the LHC experiments. This is especiallyimportant in the case of the physically well motivated k T -dependent evolution. On theother hand, the convergence of the QCD perturbative expansion looks reasonably well,better than in the DGLAP case.The main limitations of the presented in this paper MC algorithms are: missing effectsdue to the non-zero masses of the quarks in all of the discussed types of evolution andlack of the dedicated fits to the data for the modified-DGLAP type evolutions. Anotherinteresting development line of the modified-DGLAP type evolutions would be to includea non-perturbative parametrization of the behavior of the coupling constant below theLandau pole. We hope to address some of these issues in the future. The authors would like to thank S. Jadach, K. Golec-Biernat and Z. Was for numerousdiscussions and comments. References [1] S. Moch, J. A. M. Vermaseren, and A. Vogt, Nucl. Phys. B688 , 101 (2004), hep-ph/0403192.[2] A. Vogt, S. Moch, and J. A. M. Vermaseren, Nucl. Phys. B691 , 129 (2004), hep-ph/0404111. 283] C. Anastasiou, K. Melnikov, and F. Petriello, Nucl. Phys. B724 , 197 (2005), hep-ph/0501130.[4] M. Grazzini, (2008), arXiv:0801.3232.[5] A. Gehrmann-De Ridder, T. Gehrmann, E. W. N. Glover, and G. Heinrich, JHEP , 094 (2007), arXiv:0711.4711.[6] T. Sjostrand, S. Mrenna, and P. Skands, (2007), arXiv:0710.3820.[7] T. Sjostrand et al. , Comput. Phys. Commun. , 238 (2001), hep-ph/0010017.[8] S. Gieseke et al. , (2006), hep-ph/0609306.[9] G. Corcella et al. , JHEP , 010 (2001), hep-ph/0011363.[10] L. Lonnblad, Comput. Phys. Commun. , 15 (1992).[11] S. Frixione and B. R. Webber, JHEP , 029 (2002), hep-ph/0204244.[12] P. Nason, JHEP , 040 (2004), hep-ph/0409146.[13] L.N. Lipatov, Sov. J. Nucl. Phys. (1975) 95;V.N. Gribov and L.N. Lipatov, Sov. J. Nucl. Phys. (1972) 438;G. Altarelli and G. Parisi, Nucl. Phys. (1977) 298;Yu. L. Dokshitzer, Sov. Phys. JETP (1977) 64.[14] E. G. Floratos, D. A. Ross, and C. T. Sachrajda, Nucl. Phys. B129 , 66 (1977).[15] E. G. Floratos, D. A. Ross, and C. T. Sachrajda, Nucl. Phys. B152 , 493 (1979).[16] G. Curci, W. Furmanski, and R. Petronzio, Nucl. Phys. B175 , 27 (1980).[17] W. Furmanski and R. Petronzio, Phys. Lett. B97 , 437 (1980).[18] S. Jadach and M. Skrzypek, Acta Phys. Polon. B35 , 745 (2004), hep-ph/0312355.[19] K. Golec-Biernat, S. Jadach, W. Placzek, and M. Skrzypek, Acta Phys. Polon. B37 ,1785 (2006), hep-ph/0603031.[20] S. Jadach and M. Skrzypek, Comput. Phys. Commun. , 511 (2006), hep-ph/0504263.[21] S. Jadach, W. Placzek, M. Skrzypek, P. Stephens, and Z. Was, (2007), hep-ph/0703281.[22] T. Sjostrand, Phys. Lett. B157 , 321 (1985).[23] D. Amati, A. Bassetto, M. Ciafaloni, G. Marchesini, and G. Veneziano, Nucl. Phys. B173 , 429 (1980). 2924] S. J. Brodsky, G. P. Lepage, and P. B. Mackenzie, Phys. Rev. D28 , 228 (1983).[25] W. K. Wong, Phys. Rev. D54 , 1094 (1996), hep-ph/9601215.[26] G. Sterman, Nucl. Phys. B281 , 310 (1987).[27] B. I. Ermolaev, M. Greco, and S. I. Troyan, Phys. Lett. B522 , 57 (2001), hep-ph/0104082.[28] B. I. Ermolaev and S. I. Troyan, Phys. Lett. B666 , 256 (2008), 0805.2278.[29] M. Ciafaloni, Nucl. Phys. B296 (1988) 49;S. Catani, F. Fiorani and G. Marchesini, Phys. Lett. B234 Nucl. Phys. B336 (1990) 18;G. Marchesini, Nucl. Phys. B445 (1995) 49.[30] K. Golec-Biernat, S. Jadach, W. Placzek, P. Stephens, and M. Skrzypek, Acta Phys.Polon. B38 , 3149 (2007), hep-ph/0703317.[31] M. Skrzypek, S. Jadach, K. Golec-Biernat, and W. Placzek, Acta Phys. Polon. B38 ,2369 (2007).[32] G. Marchesini and B. R. Webber, Nucl. Phys. B349 , 617 (1991).[33] J. Kwiecinski, Acta Phys. Polon. B33 (2002) 1809;A. Gawron, J. Kwiecinski and W. Broniowski, Phys. Rev. D68 (2003) 054001;A. Gawron, J. Kwiecinski, Phys. Rev. D70 (2004) 014003.[34] R. G. Roberts, Eur. Phys. J. C10 , 697 (1999), hep-ph/9904317.[35] P. Stoklosa and M. Skrzypek, Acta Phys. Polon. B38 , 2577 (2007).[36] S. Tsuno, T. Kaneko, Y. Kurihara, S. Odaka, and K. Kato, Comput. Phys. Commun. , 665 (2006), hep-ph/0602213.[37] Y. Kurihara et al. , Nucl. Phys. B654 , 301 (2003), hep-ph/0212216.[38] H. Jung and G. P. Salam, Eur. Phys. J. C19 , 351 (2001), hep-ph/0012143.[39] S. Jadach et al. , in preparation.[40] M. H. Seymour, Comp. Phys. Commun. , 95 (1995), hep-ph/9410414.[41] T. Sjostrand, S. Mrenna, and P. Skands, JHEP , 026 (2006), hep-ph/0603175.[42] K. Golec-Biernat, S. Jadach, W. Placzek, and M. Skrzypek, Acta Phys. Polon. B39 ,115 (2008), 0708.1906.[43] W. Placzek, K. Golec-Biernat, S. Jadach, and M. Skrzypek, Acta Phys. Polon. B38 ,2357 (2007), arXiv:0704.3344. 3044] G. Marchesini and B. R. Webber, Nucl. Phys. B310 , 461 (1988).[45] K. Golec-Biernat,