On the transverse momentum in Z-boson production in a virtuality ordered parton shower

Cross sections for physical processes that involve very different momentum scales in the same process will involve large logarithms of the ratio of the momentum scales when calculated in perturbation theory. One goal of calculations using parton showers is to sum these large logarithms. We ask whether this goal is achieved for the transverse momentum distribution of a Z-boson produced in hadron-hadron collisions when the shower is organized with higher virtuality parton splittings coming first, followed successively by lower virtuality parton splittings. We find that the virtuality ordered shower works well in reproducing the known QCD result.


Introduction
Parton shower algorithms with hadronization models provide a way of generating simulated events according to approximations based on QCD. Since complete final states are generated, one can generate completely exclusive cross sections in this way. By summing over variables that one chooses not to examine, one can also make predictions for inclusive observables. Of special interest are predictions that, in a perturbative expansion, involve large logarithms of ratios of different momentum scales in the physical problem. An important example is the distribution of the transverse momentum p ⊥ of Z-bosons produced in hadron-hadron collisions at some fixed rapidity Y , dσ/(dp ⊥ dY ). When p 2 ⊥ M 2 Z , the perturbative expansion of this cross section contains two powers of the large logarithm log(p 2 ⊥ /M 2 Z ) per power of α s . The large logarithms spoil the usefulness of fixed order perturbation theory for this observable and for observables that contain similar large logarithms. A parton shower algorithm sums contributions to the desired cross section that contain arbitrarily high powers of α s ; thus it sums the accompanying logarithms. For this reason, one can hope that a parton shower calculation will do better than a fixed order perturbative calculation in the region p 2 ⊥ M 2 Z and analogous regions for other processes. Indeed, the basic approximation in a parton shower is that one parton splits into two daughter partons with a probability that matches the singularities of the QCD matrix element when the two daughter partons are collinear or one of them is soft. It is just these soft/collinear configurations that give rise to the large logarithms. Thus one can hope that the cross section generated by a parton shower will be a good approximation to the true QCD result.
In many cases, including the Z-boson transverse momentum distribution, there are predictions based on the full field theory. That is, we know that the cross section "exponentiates" in a sense that one can state precisely and we know some of the leading coefficients that appear in the exponent of the formula that expresses the QCD result.
When, for a particular process, one knows the summation of large logarithms in full QCD, then it is of significant interest to investigate whether a given shower algorithm produces matching results. To do this, one needs to derive the corresponding summation in the shower model, deriving the appropriate evolution equation for the observable in question from the general evolution equation for the shower algorithm.
We have argued above that, because parton shower algorithms are generally based on parton splitting probabilities that have the proper soft/collinear singularities, these algorithms may provide good approximations to the full QCD result in particular cases involving summing large logarithms. However, a parton shower algorithm contains several ingredients beyond the parton splitting probabilities. Among theses are the momentum mapping, the color and spin treatment, and the choice of evolution variable. Depending on the choice of these ingredients, one may obtain agreement with full QCD for a particular observable or one may fail to obtain agreement. Certainly, it is widely understood that a virtuality ordered parton shower without a proper inclusion of the effects of quantum interference can get results for some observables that do not match QCD. In contrast, there is an extensive body of literature [1,2,3,4,5,6,7,8,9] that suggests that parton showers based on ordering in parton emission angles does better for many observables. In this paper, we will indeed see that, for the Z-boson transverse momentum distribution, the choice of ingredients matters.
Because the choice of ingredients matters, we think it important to validate shower evolution schemes against known results for summations of large logarithms. We believe that such a validation program could be useful for understanding the range of validity of current parton shower event generators and could provide important guidelines for the current and future development of such programs.
In ref. [10], in response to ref. [11], we investigated the parton energy distribution in electron-positron annihilation as predicted by the virtuality ordered, color-dipole based parton shower algorithm of refs. [12,13,14] and by k T -ordered variants of this. We concluded that the predictions of these parton shower algorithms are consistent with the field theory result represented by the well know DGLAP evolution equation [15]. This conclusion was confirmed by ref. [16], which included numerical studies. The parton energy distribution is an example of the general program of summing large logarithms, namely the logarithm of the resolution scale for finding the partons (∼ jets) divided by the electron-positron energy. However, this is a rather simple example in that there is only one logarithm per power of α s . Observables for which there are two logarithms per power of α s present a much more subtle problem.
In this paper we investigate the transverse momentum distribution of a Z-boson produced in hadron-hadron collisions, the Drell-Yan process. Here, there are two logarithms per power of α s . We use a slightly modified version 1 of the shower evolution of ref. [12]. This evolution equation describes the evolution of the partonic states in a fully exclusive way. We manipulate the shower equation to produce an evolution equation for the trans- 1 We will find that we need to change the momentum mapping for initial state radiation from that of ref. [12]. In addition, one of the choices given in ref. [14] for a certain function A lk that was left unspecified in ref. [12] gives a satisfactory result, but another of the choices does not work. verse momentum distribution of the Z-boson and we compare the result to the well known field theory prediction that is given in ref. [17]. We find that the virtuality ordered shower works well in reproducing the full QCD result.
The shower evolution equation of ref. [12] includes quantum interference among colors and spins. We note that this shower evolution equation is not directly practical for generating events. A simple approximation to this evolution equation is to average over spins and take the leading color approximation, which yields an evolution equation [13] that can be written as a Markov process and is thus directly practical for generating events. We treat the full evolution equation, but we will see that the same results with some small adjustments hold for the Z-boson transverse momentum distribution produced by the spin-averaged, leading-color shower evolution.
The investigation that we carry out in this paper, and that we believe would be useful for other observables and other shower algorithms, is analytical. That is, we want to test whether the transverse momentum distribution produced by the shower algorithm has the proper exponential form and, if so, whether the coefficients in the exponent are correct. A followup study, not addressed in this paper, would be numerical: how well does an actual implementation of the parton shower algorithm produce results that match numerical results given by the QCD formula. In this case, the summed QCD results, including nonperturbative parameters that are fit to experiment, could be obtained from the Resbos code [18]. We believe that one should start with an analytical study rather than a numerical study for the following reason. The parton shower, for the case of the production a Z-boson with small p ⊥ , represents the physics on "hard" scales from M Z to a few GeV. This is perturbative physics that is adjustable only to a limited extent. In particular, one can modify the parton splitting probabilities in a fashion that does not change them in the soft and collinear limits. 2 A full parton shower event generator also contains elements, such as a hadronization model and a model for the underlying event, that represent soft scale physics. One can adjust the models for the soft physics by tuning various parameters. For very large M 2 Z /Λ 2 QCD and M 2 Z /p 2 ⊥ , only the part of the shower algorithm that cannot be tuned should matter. However, for realistic values of these parameters, the tunable parts of the parton shower can make a numerical difference. One would not want to tune the parameters in order to repair a numerical disagreement that was actually caused by a mismatch between the parton shower and full QCD with respect to the hard scale physics that is fixed. The way to check the match of hard scale physics is to compare analytically. Once this is checked, the numerical comparison is appropriate and needed.
We organize this paper as follows. In sections 2, 3, and 4, we review briefly the needed features of the shower evolution of ref. [12] and set up the notation for our analysis. Then in section 5 we outline the derivation to come and state the nature of the approximations that we will need. The derivation is given in sections 6, 7, 8, 9, and 10. The solution of the evolution equations is given in section 11 and the result is compared to the full QCD result in section 12. We discuss what one would get with other sorts of shower evolution in section 13 and summarize the results in section 14. In appendix A, we summarize results from refs. [12,13,14] that are used in this paper. In appendix B we prove a certain property of integrals of the J 0 Bessel function.

Shower states and shower evolution
We analyze the transverse momentum distribution of a Z-boson as generated by the parton shower evolution equations of ref. [12]. The organizing principle of the shower is that one starts at the hard interaction q +q → Z and moves to softer interactions, always factoring the softer interaction from previous harder interactions.
In the notation of ref. [12], states in the sense of statistical mechanics are represented by ket vectors ρ , while possible measurements are represented by bra vectors F . Thus F ρ is the cross section that one obtains a particular result F from a measurement on an ensemble of systems represented by ρ . We use basis states labelled by lists {p, f, s , c , s, c} m of parton quantum numbers for two initial state partons and m final state partons. As the shower progresses, m increases. The momenta of the final state partons are {p 1 , . . . , p m }. The first final state parton is the Z-boson, with momentum p 1 ≡ p Z . The momenta of the initial state partons are specified by giving their momentum fractions, η a and η b . Then their momenta are where d{p, f, s , c , s, c} m indicates integrating and summing over all of the quantum numbers. A much more detailed specification of our notation is provided in ref. [12]. In ref. [12], the quarks can have non-zero masses, but in this paper we take all of the partons except for the Z-boson to be massless. Let ρ(t) = U(t, 0) ρ(0) be the statistical state at shower time t. Here t = 0 at the beginning of the shower, which starts with the hard process q +q → Z, and t = t f gives the statistical state at the end of the shower, before hadronization. (We do not discuss a hadronization model in this paper.) The total cross section for producing a Z-boson is expressed using the vector 1 , which represents the totally inclusive measurement, The total Z-production cross section is independent of t: the shower evolution maintains 1 U(t, 0) = 1 .
We are interested in the differential cross section dσ/(dp ⊥ dY ) as obtained in the shower at the shower final time t f , Here, the measurement function p ⊥ , Y measures the cross section for the Z-boson to have transverse momentum p ⊥ and rapidity Y . 3 The definition is (2.7) The starting point of shower evolution is the Born cross section, which is proportional to δ(p ⊥ ). (See eq. (4.17) in section 4.5 for details.) The shower evolution is specified in ref. [12] in the form (2.9) 3 In general, we denote vectors in two transverse dimensions by boldface letters like p ⊥ . The transverse part of a four-vector p is p ⊥ . Then p 2 ⊥ = −p 2 ⊥ < 0; however, to avoid confusion, we avoid writing p 2 ⊥ .
Here H I (t) is the splitting operator, which takes a basis state with m final state partons and changes it to a state with m + 1 final state partons. Next, V(t) is a "virtual splitting" operator that leaves number of partons and their momenta, flavors, and spins unchanged. 4 The operator V(t) is determined from H I (t) in such a way that so that the total cross section to produce a Z-boson remains the Born cross section, even though the Z-boson momentum changes as a result of recoils against parton splittings in the shower. In particular, the Z-boson transverse momentum distribution broadens as the shower develops from t = 0 to t = t f . We need to follow the shower evolution to find how the transverse momentum distribution broadens.

Initial state splitting kinematics
We will first need some kinematics for the description of an initial state splitting. For notational convenience, we suppose that it is parton "a" that splits. The initial state parton with momentum p a splits, in backward evolution, to a new initial state parton with momentump a and a final state parton with momentump m+1 , as illustrated in figure 1. The other initial state parton has momentum p b before the splitting and momentump b , equal to p b , after the splitting. In this paper, 5 we describe the splitting 4 In general, the operator V(t) is a non-trivial operator on the partonic color space. In the leading color approximation, valid for Nc → ∞, it is diagonal in the color space, as described in ref. [13]. The derivation in this paper is given for the exact treatment of color, but works also in the leading color approximation. 5 In ref. [12], the splitting operator is expressed in terms of momenta rather than splitting variables y, z, φ and we did not specify a choice of splitting variables. In ref. [13], we did make a choice, but that choice is not the same as the choice used here.
using splitting variables (y, z, φ) or, alternatively, (y, k ⊥ ), defined bŷ Here k ⊥ is the part ofp m+1 that is orthogonal to both p a and p b .
Virtuality. The variable y is We will use a shower time t based on virtuality, 6 where M is the Z-boson mass. We will typically use t as our virtuality variable instead of y, so that y is Momentum fraction. The fraction of the momentump a in the direction of p a that is carried away by the emitted final state parton m + 1 iŝ The variable z must be in the range 0 < z < 1. The momentum fraction η a of parton "a" has a new value after the splitting. From eqs. (2.1) and (3.1), we have, using 2 These are the exact relations. In our applications, we will generally neglect M 2 e −t compared to η a η b s. That is, we consider the virtuality M 2 e −t of a splitting to be small compared to the momentum scale of the hard process, M 2 , which, in turn, is always smaller than η a η b s.
Transverse momentum. The angle φ is the azimuthal angle of k ⊥ . The magnitude of k ⊥ is related to z and y: (3.7) Thus If we use z and φ (along with y or t) as our splitting variables, then k 2 ⊥ is a derived variable. Alternatively, we can use k ⊥ (and thus k 2 ⊥ and φ) as splitting variables. Then z is a derived variable.
Lorentz transformation. If y = 0, momentum differencep a −p m+1 is not exactly equal to p a . In order to maintain momentum conservation at each step in the shower, we must, therefore, take some momentum from the partons in the final state at the time of the splitting. Each parton, with momentum p j , j ∈ {1, 2, . . . , m}, then gets a new momentum p j after the splitting. This includes the momentum p Z ≡ p 1 of the Z-boson. Following the shower algorithm of ref. [12], the momentap j are determined by a Lorentz transformation, (3.9) However, we use a different Lorentz transformation from that chosen in ref. [12]. 7 If thenp = Λp is given byp An equivalent form that is useful if α = 0 iŝ (3.12) 7 The choice of Lorentz transformation in ref. [12] takes the needed total momentum from the final state partons, but it does not properly absorb the transverse momentum recoil onto the Z-boson. The transformation defined here satisfies eq. (3.9), leaves invariant any vector that is orthogonal to pa, p b , and pm+1 and, in addition, transforms p b into a multiple of itself.
The momentap a andp b are not the Lorentz transformed versions of p a and p b . It is, however, of interest to know what the Lorentz transformation does to p a and p b . We have (3.13) The important feature of this is that the emitted parton m + 1 has transverse momentum k ⊥ and the momentum Λp a of parton "a" caries the recoil transverse momentum −k ⊥ .
We can understand what happens to this recoil transverse momentum by thinking of the shower as proceeding forward in time (oppositely to the way it is generated). Initial state parton "a" with momentum Λp a can emit more daughter partons. But a share of the transverse momentum −k ⊥ is finally transmitted to the Z-boson.
Let us now look directly at the transformation of the momentum p Z of the Z-boson. We start with with Thus the Z-boson gets a share x a /η a of the recoil transverse momentum. The new rapidity isŶ This transformation law is complicated. However it simplifies in the limit that we need for this paper. Let us denote by P the Z-boson momentum without its transverse part, We are interested in the development of the Z-boson transverse momentum in the region p 2 Z,⊥ M 2 . Therefore we take Furthermore, the development of the Z-boson transverse momentum distribution is controlled by splittings with y 1. Therefore, we neglect p 2 Z,⊥ /M 2 ,p 2 Z,⊥ /M 2 , and y compared to 1 in Eq. (3.18), givingŶ = Y (3.21) in each splitting. With these approximations, x a and x b are fixed: With the approximations eq. (3.22), we have Although x a and x b are fixed, the momentum fractions η a and η b can change if there is a collinear splitting of an initial state parton. 8 Although we neglect the transverse momentum of the Z-boson in computing its mass, we track changes in the transverse momentum as the Z-boson recoils against emissions from the initial state partons. For an emission from initial parton "a", the new Z-boson transverse momentum isp as stated in eq. (3.17).

Analysis framework
The equations of ref. [12] specify quite precisely the evolution of a certain kind of parton shower. In order to analyze what the parton shower thus defined produces for the transverse momentum distribution of a Z-boson, we develop in this section some theoretical structures beyond those presented in ref. [12].

Measurement operators Q
We are interested in the differential cross section as obtained in the shower at the shower In this paper, we will find it useful to represent the desired measurement with the aid of an operator Q on the space of statistical states, (4.2) Here Q(p ⊥ , Y ) is defined by In the subsequent subsections, we will extend this notation in which an operator Q determines a measurement on the statistical state.

Fourier transformation
As was noted by Parisi and Petronzio [20], it is useful to analyze the evolution of the Fourier transform of the transverse momentum distribution. Thus we study We will refer to 1 Q(b, Y ) ρ(t) as the b-space hadronic cross section.

Tracking the momentum fractions and parton flavors
For our analysis, we will want to keep track of the parton momentum fractions, η a and η b , and the flavors, a and b, of the incoming partons. Thus we consider the function (4.6)

Evolution of the perturbative statistical state
The statistical state vector ρ(t) according to our definition contains a factor This factor gives the parton-parton luminosity. Here f a/A (η a , µ 2 ) and f b/B (η b , µ 2 ) are parton distribution functions and n c (a) and n c (b) are the number of colors carried by partons of flavors a and b, namely 3 for quarks and 8 for gluons. We define an alternative state vector ρ pert (t) in which this non-perturbative factor is removed: (4.8) A convenient notation for this is where F(t) multiplies by the parton distribution factor, The evolution equation for ρ pert (t) can be determined from the evolution equation (2.9) for ρ(t) . We have (4.12) We can write this as  Here the revised real and virtual splitting operators are H pert (4.14) In the last line here, we have noted that F(t) commutes with V(t) since V(t) does not change momenta or flavors. For the perturbative real splitting operator, we have × {p,f ,ŝ ,ĉ ,ŝ,ĉ} m+1 H pert  The factors in the second line appear in the definition of the matrix element of H I in ref. [12]. Thus to obtain the corresponding matrix element of H pert I (t), we simply omit these factors. We will present detailed formulas for H pert I (t) and V pert (t) at the point that we need them.

The function to study
The physics that we want to study is contained in the function defined in section 4.3, in which b and Y are measured and, in addition, we measure η a , η b , a and b. As noted in the previous subsection, this function contains a nonperturbative factor that specifies the parton luminosity. We remove this factor and study We will refer to this function as the b-space partonic cross section. Once we have the b-space partonic cross section, we can obtain the b-space hadronic cross section 1 Q(b, Y ) ρ(t) by convolving it with parton distribution functions according to Our aim is to study how the b-space partonic cross section develops as the shower time t increases. At shower time 0, it is determined from the Born cross section, as in eq. (2.8), and Q ab = 0 a = g or b = g , Q ab = δ a,b [1 − 4 |e a | sin 2 (θ W )] 2 + 1 16 sin 2 (θ W ) cos 2 (θ W ) a = g and b = g . Note that the partonic cross section at t = 0 vanishes unless a is a quark or antiquark flavor and b is the corresponding antiflavor. There is no dependence on b because the corresponding transverse momentum dependent cross section is proportional to a delta function of the transverse momentum. As the shower evolves, we expect 1 We begin the study of the evolution of this function in the next subsection by outlining some key ideas that will go into the derivation.

Outline of the derivation
We are now in a position to outline the derivation that follows, leaving out most of the details and certain subtle points. One of the subtle points is the running of α s . For the purposes of this section, we consider α s to be constant.

Emission kinematics
We will find that the transverse momentum of the Z-boson is the primarily the result of recoils against emission of soft gluons from an initial state quark or antiquark. As we saw in section 3, an initial state splitting can be described by splitting variables {t, z, φ}.
The shower time t gives the virtuality of the splitting. The momentum fraction of the emitted gluon is (1 − z), so that for soft gluon emission we have (1 − z) 1. The transverse momentum k ⊥ of the emitted gluon has azimuthal angle φ and magnitude given by eq. (3.8), k 2 ⊥ = (1 − z)M 2 e −t . It is useful to represent splittings by points in the plane of t and log(k 2 ⊥ /M 2 ), as in figure 2. An allowed splitting has (1 − z) ≤ 1, so log(k 2 ⊥ /M 2 ) ≤ −t. That is, allowed splittings are represented by points below the line labelled "(1 − z) ∼ 1" in figure 2. One should view this line as having a thickness of order 1. As we will discuss in some detail, for a given t, the splitting probability dP from ref. [12] has a term with including a factor 2 due to having two incoming partons that can radiate. This means that splittings are not simply concentrated along the line (1 − z) ∼ 1, but are spread over the region below this line. In fact, a splitting probability proportional to dt dz/(1 − z) would give splittings uniformly distributed in log(k 2 ⊥ /M 2 ) and t. However, there is an effective cutoff at the line (1 − z) ∼ e −t . Thus the splitting probability is approximately constant in the region e −t < (1 − z) < 1 indicated in figure 2.
We are interested in the partonic b-space cross section. For this quantity, a real emission produces a factor exp(ik ⊥ ·b), which simply comes from taking the Fourier transform to get to b-space. For this reason, the line k 2 ⊥ ∼ 1/b 2 in figure 2 is significant. If we integrate the real emission probability over an interval of k ⊥ in the region k 2 ⊥ 1/b 2 , the factor exp(ik ⊥ · b) averages to zero. That is to say, maintaining a non-zero 1 Q(b, Y ; η a , η b , a, b) ρ pert (t) requires not emitting gluons with k 2 ⊥ 1/b 2 . In the region k 2 ⊥ 1/b 2 , gluons can be freely emitted into the final state. We represent gluons that might be emitted in a typical event contributing to 1 Q(b, Y ; η a , η b , a, b) ρ pert (t) by filled circles in figure 2.
Of course, if there is an approximately uniform probability of emitting gluons in any differential unit of area d log(k 2 ⊥ /M 2 ) dt, then the probability that no gluons are emitted for k 2 where A is the area of the triangle in figure 2. This Sudakov factor gives, approximately, the b-dependence of the partonic b-space cross section that we seek. In the following sections, we fill in the details of this argument and make it more precise. We will find that the more precise analysis leads to a Sudakov factor similar to that in eq. (5.2), but with running coupling effects included and an extra "subleading" term.

Strategy
We study the evolution of the b-space partonic cross section defined in section 4.5. Using eq. (4.13), the b-space partonic cross section evolves according to Our aim is to use suitable approximations to turn this equation into a differential equation of the form This differential equation has the solution Here the initial b-space partonic cross section has the simple b-independent form given in eq. (4.17). We will see that the evolution of 1 Here t c is approximately the shower time at which the lines (1 − z) ∼ 1 and k 2 ⊥ ∼ 1/b 2 in figure 2 meet; we will see later the reason for the adjustment factor e 2γ E /4. Then is the Sudakov factor, for which eq. (5.2) is a simple approximation.

Approximations
We will need certain approximations to turn eq. (5.3) into the differential equation (5.4). We describe these approximations in general terms here. First, we note that the behavior of the Z-boson p ⊥ distribution for p 2 ⊥ M 2 is controlled by the b-space partonic cross section for large b 2 . Thus we are interested in the b-space partonic cross section in the region 1/(b 2 M 2 ) 1. Therefore, we simply neglect 1/(b 2 M 2 ) compared to 1 everywhere.
Second, we neglect e −t compared to 1. To justify this, imagine letting the system evolve from time 0 to a time t 1 and then from t 1 to t c . Let t 1 be large enough so that e −t 1 1, but small enough that we can treat t 1 as not being a large logarithm. Then evolution from 0 to t 1 is an approximate version of perturbation theory and gives order α n s corrections to the Born cross section with no large logarithms. We ignore these corrections. For the evolution from t 1 to t c , the approximation e −t ≈ 0 is justified. Furthermore, we can add back the evolution from 0 to t 1 using the approximation e −t ≈ 0, adding more order α n s corrections with no large logarithms. Then we have evolution from 0 to t c with the approximation e −t ≈ 0 at the cost of changing the result by α n s terms with no large logarithms.
For the same reason, we neglect y, eq. (3.4), compared to 1 and k 2 ⊥ /M 2 compared to 1.
Finally, in section 9.2, we will analyze the structure of the splitting function near the line For this analysis, we will make what might be called a low density approximation. For initial state emissions, according to eq. (5.1), the density of emission points per unit dt and d log(k 2 ⊥ /M 2 ) is proportional to α s . We treat α s as small. Consider, then, two emissions, one with parameters {t 1 , log(k 2 ⊥,1 /M 2 )} and the other with parameters {t 2 , log(k 2 ⊥,2 /M 2 )}. Suppose that t 2 > t 1 and that each of these emissions is not far from the line k 2 ⊥ = 1/b 2 . Since the density of points is small, the distance in the {t, log(k 2 ⊥ /M 2 )} plane between any two points is typically large. This suggests that we can neglect e −(t 2 −t 1 ) compared to 1. To see whether this is justified, suppose that we modify the shower algorithm so that it is not allowed to have two splittings be close to the line k 2 ⊥ = 1/b 2 and close to each other. More precisely, choose a distance parameter d 0 and require that no two splittings have | log(k 2 ⊥ b 2 )| < d 0 and simultaneously t 2 − t 1 < d 0 . Choose d 0 such that e −d 0 is small enough to be neglected but d 0 is not so large that it constitutes a large logarithm. Then it is valid to replace e −(t 2 −t 1 ) by zero for any splittings that are not excluded.
The exclusion prescription can be constructed in a different way. Generate 1 → 2 splittings without restriction, but add a new splitting function that describes a 1 → 3 splitting: incoming parton "a" emits partons 1 and 2 in the region | log(k 2 ⊥ b 2 )| < d 0 and t 2 − t 1 < d 0 . The probability for this new splitting should be negative and just big enough to cancel the probability for parton "a" to first emit parton 1 and then emit parton 2 in this region. One then needs to add similar 1 → n + 1 splittings for n > 2 to cancel the probabilities to have more than two splittings that are to close the line k 2 ⊥ = 1/b 2 and to one another, but we limit this discussion to the simple 1 → 3 case. The new 1 → 3 splitting probability is proportional to α 2 s and the excluded area, d 2 0 . With our choice of d 0 , this area is not large.
The new splitting function will modify the Sudakov exponent by adding a term proportional to α 2 s and to the length of the line k 2 ⊥ = 1/b 2 between the two limits in figure 2, We conclude that replacing iterated splittings by zero in the region in which the approximation e −(t 2 −t 1 ) ≈ 0 is not good results in modifying the Sudakov exponent by terms of order α 2 s log(b 2 M 2 ). What we actually do is to replace these splittings by what the inaccurate approximation e −(t 2 −t 1 ) ≈ 0 gives. As long as this approximation results in a fractional change of order 1 in the iterated splitting probability, we also modify the Sudakov exponent by terms of order α 2 s log(b 2 M 2 ). The true QCD Sudakov exponent, as discussed in section 12, has an expansion in powers of α s . In the coefficient of α 2 s , the term with the most powers of log . The next-to-leading term is a constant times α 2 s log 2 (b 2 M 2 ). Thus terms of order α 2 s log(b 2 M 2 ) are third-to-leading. The low density approximation discussed here changes these terms in an uncontrolled way.
We have not analyzed here the effect of 1 → n + 1 splittings for n > 2 that are induced by this approximation. However, it should be clear that these induce α n s log(b 2 M 2 ) changes in the Sudakov exponent.

Evolution of the partonic cross section
We can now begin to simplify eq. (5.3), which gives the evolution of the b-space partonic cross section.
The operator V pert (t) acting on a partonic basis state {p, f, s , c , s, c} m does not add a new parton or change the parton momenta or flavors. For this reason, Here H FS (t) is the part of H pert I (t) that generates the splittings of final state particles. We have dropped the "pert" notation here because, in the definition of ref. [12], a final state splitting does not involve a factor of ratios of parton distributions, so that the part of H pert I (t) that creates a final state splitting is the same as the part of H I (t) that creates the same final state splitting. In the second term in eq. (6.2), we let H pert a (t) be the part of H pert I (t) that generates the splittings of the incoming parton from hadron A. We have decomposed this operator further as an integral and sum of H pert a (t; z, φ, f ), which generates the splittings of the incoming parton from hadron A in which the momentum fraction and azimuthal angle of the splitting are z and φ, respectively and the flavor of the emitted parton is f . (For our analysis, f = g is the most important choice.) Similarly, H pert b (t; z, φ, f ) generates the splittings of the incoming parton from hadron B. In a similar way, we can divide the virtual splitting operator into three parts,

Final state splittings
Consider first the effect of a final state splitting. Using the definitions of ref. [12], we find that a final state splitting replaces one final state parton by two daughter partons, but its effect on the rapidity and transverse momentum of the Z-boson is negligible. When parton l with momentum p l splits to form daughter partons l and m + 1 with momentap l and p m+1 , an amount of momentum ∆p =p l +p m+1 −p l must be taken from the other final state partons. The splitting can be characterized by the virtuality variable y =p l ·p m+1 /p l ·Q, where Q = p a + p b . With the splitting kinematics of ref. [12], the needed momentum is Note that 1 − λ ∝ y for y 1. The momentum ∆p is supplied by applying a Lorentz transformation to each final state parton i with i / ∈ {l, m + 1}:p µ i = Λ µ ν p ν i . In particular, the momentum p Z of the Z-boson is transformed withp µ Z = Λ µ ν p ν Z . As discussed in section 5.3, it suffices to consider only splittings with y 1. For y 1, ∆p is proportional to y and hence Λ µ ν − δ µ ν is also proportional to y. Thus the change in the Z-boson momentum is small. In particular, the rapidity of the Z-boson changes very little, by an amount proportional to y. We can neglect this small change.
Evidently, the change in the Z-boson transverse momentum must also be small, but this statement is not helpful because we are trying to track small changes in the Z-boson transverse momentum. To see what happens, we note that the needed transverse momentum is A fraction of this transverse momentum will come from the Z-boson. Now, we are studying the evolution of the probability that the Z-boson transverse momentum p Z,⊥ is small and remains small. For this to happen, the transverse part of p l must be small, of order p Z,⊥ or smaller. 10 Thus, recalling that (1 − λ) is of order y, we havê where C is of order 1 or smaller. That is, the fractional change in the Z-boson transverse momentum due to a final state splitting is negligible. This discussion can be summarized by saying that, to a sufficient approximation, This gives This is useful because the definition of V FS (t) is based on the requirement, designed to insure that the shower conserves probabilities, that Thus the first term in eq. (6.9) vanishes. We must analyze initial state splittings, but we can ignore final state splittings entirely.

Initial state splittings
We now turn to initial state splittings. We relate Q applied after the splitting to Q applied before the splitting. The relation is To understand this, first look at the b dependence, using the definition (4.6) of Q. When we apply the operator Q after the splitting, it produces a factor exp(−ib·p Z,⊥ ) wherê p Z,⊥ is the Z-boson transverse momentum after the splitting, which is related to the Zboson transverse momentum after the splitting, p Z,⊥ , and the transverse momentum in the splitting, k ⊥ , by eq. (3.24),p Recall that k ⊥ is specified by z and φ: it has azimuthal angle φ and square k 2 The factor exp(−ib·p Z,⊥ ) is generated by the Q operator before the splitting but the second factor in eq. (6.13) must be supplied. The dependence on the momentum fractions is simple. According to eq. (3.6), the momentum fraction η b is unchanged by the splitting, while he momentum fractionη a after the splitting is related to the momentum fraction η a before the splitting by η a ≈ zη a . Thus the η a argument of Q before the splitting is zη a and there is a jacobian factor z because Q is defined with a delta function. Finally, the flavor a after the splitting is related to the flavor a before the splitting by a =â − f , where we use the notation u − g = u, g −ū = u, etc.
With these observations, our equation for the variation b-space partonic cross section with shower time is where K a describes a splitting of an initial state parton "a" from hadron A and is given by (6.15) The operator K b for a splitting of the initial state parton from hadron B is the same with the roles of "a" and "b" are interchanged.

The real splitting function
At this point, we need to know some details about the operation of H pert a (t; z, φ) on a general partonic state {p, f, s , c , s, c} m . Fortunately, we need only the inclusive splitting Formulas from refs. [12,13,14] for this quantity are reviewed in appendix A. The most important case to consider is that of a q → q +g splitting or aq →q +g splitting. However, we include all flavor choices. We treat separately two kinematic regimes: (1 − z) 1 and (1 − z) ∼ 1 since the results in these two regimes have rather different structures. Our subsequent analysis of evolution will make use of the results in these two regions.
The operator H pert a contains a factor α s . In ref. [12], the argument of α s was denoted by µ 2 R and left unspecified. In general, µ 2 R can be a function µ 2 R (z, t) of the kinematic variables that describe the splitting. Our default choice in this paper is where In section 12, we will see why choosing µ 2 R approximately proportional to k 2 ⊥ is useful and we will see why the choice given in eq. (7.3) for the proportionality constant is also useful.
Leaving these points for later, we can immediately understand why a factor (1 − z + y) in eq. (7.2) is preferable to the simpler choice (1 − z). As we will see, having a running scale µ 2 R as in eq. (7.2) with either a factor (1 − z + y) or a factor a factor (1 − z) affects the Sudakov exponent that we obtain in section 12. When α s (µ 2 R ) is expanded in powers of α s (M 2 ) one obtains terms proportional to logarithms of (1 − z + y) or (1 − z) times extra powers of α s . These terms improve the matching between the Sudakov exponent obtained with the shower and the true QCD Sudakov exponent. We will find that, with the accuracy of matching that we can obtain, logarithms of (1 − z + y) or of (1 − z) are equivalent. However, we can still ask which is more desirable in general. The running α s (µ 2 R ) incorporates some features of the singularity structure of higher order graphs into the leading order splitting. We note that the leading order splitting kernel from ref. [12] has a singularity 1/(1 − z + y). This suggests that the higher order contributions might naturally contain a logarithm of this same variable, (1 − z + y). In contrast, the leading order splitting kernel does not have a singularity when (1 − z) → 0 at fixed y, so it would not be natural to introduce a logarithm of (1−z) into the expansion of α s (µ 2 R ). Indeed, soft gluon emissions correspond to y → 0 and (1 − z) → 0 together, while (1 − z) → 0 at fixed y corresponds to an anticollinear emission in which incoming parton "a" emits a gluon in the direction of incoming parton "b". There are such singularities, but they are associated with emissions from incoming parton "b" rather than from incoming parton "a". For this reason, our default choice (7.2) for µ 2 R has a factor (1 − z + y) rather than (1 − z).
We use the splitting operator for an initial state splitting from ref. [12], using the splitting variables t, z, φ defined in section 3 of this paper. When the emitted parton is a gluon, the splitting probability has a "soft gluon emission" singularity that corresponds to a factor 1/(1 − z) when y (1 − z) 1. In this section, we extract the terms that have this soft gluon factor; other terms will be included in the following subsection, where we study the regime (1 − z) ∼ 1. We note, in particular, that contributions in which the emitted parton is not a gluon do not give a 1/(1 − z) contribution.
We use the results of ref. [12], particularly eqs. (12.20), (12.21), and (12.22), as reviewed in appendix A. 11 Part of eq. (12.21) is a function A lk . There is some freedom in choosing this function. We use the definition in eqs. (7.2) and (7.12) of ref. [14], which are equivalent to eqs. (A.6) and (13.4); we discuss the reason for this choice in section 13.2. We neglect y compared to 1 and (1 − z) compared to 1. However, we do not neglect y compared to (1 − z). In the shower formalism of ref. [12], when parton "a" splits by emitting a gluon, we include interference graphs. The soft gluon is emitted from parton "a" in the amplitude and by some other parton k in the complex conjugate amplitude, or by parton "a" in the complex conjugate amplitude and by parton k in the amplitude. Here parton k could be parton "b" or could be a final state parton. For this reason, there is a sum over the indices k of helper partons, with k = a.
After some calculation, we find in appendix A that Here r k is the rapidity of the helper parton k relative to the rest frame of p a + p b and φ k is its azimuthal angle; then f is Because we are calculating an inclusive splitting probability, indicated by the measurement function 1 , there is a quantum inner product {s } m {s} m between the spin states in the amplitude and the conjugate amplitude. 12 There is also a color inner product, which is non-trivial because it contains a color matrix T c a for emitting a gluon from line "a" and a color matrix T c k for emitting the gluon from line k, summed over the eight colors c of the emitted gluon. This same factor appears in the Catani-Seymour dipole subtraction formalism for next-to-leading order perturbative calculations [21].
This splitting function is complicated because of the function f . To understand f , denote the rapidity of the emitted gluon relative to the rest frame of p a + p b by r. For a soft gluon, (1 − z) 1 and y 1, we have For a given splitting time t, y is fixed and we integrate over z. There is a near singularity in this integral for (1 − z) → 0, but this near singularity is cut off when (1 − z) becomes comparable to y. That is, r 1 is favored in the integration. Writing e 2r y for (1 − z) in We see that when r ∼ r k , all three terms in f (z, y, φ; r k ) are comparable, so that we have quite a complicated function. However, We thus see that when the rapidity of the emitted gluon is large compared to the rapidities of previously emitted gluons, the splitting function simplifies to Now the only dependence on the helper parton index k is through the color factor. This enables us to perform the color sum as described in ref. [12], That is, using eq. (2.5), Now consider the collinear limit (1 − z) ∼ 1. We use the general result from ref. [12], reviewed in appendix A. We neglect e −t compared to 1 as discussed in section 5.3. As described in appendix A, the result is a simple color structure that multiplies the standard (unregulated) DGLAP splitting kernels P a,a (z) and the ratio of the number of colors for parton flavor a = a + f to the number of colors for parton flavor a,

The virtual splitting function
We now look at the virtual splitting function. We need Eq. (4.14) gives Following ref. [12], we define the virtual splitting operator V(t) from the requirement which guarantees that shower evolution does not change the cross section when we integrate over all final states that start from the hard scattering. We break H I into pieces as in eq. (6.2), We can similarly write V(t) as a sum and integral in the form With this notation, the definition in ref. [12] is The needed matrix element involving the derivative of F(t) is × 1 {p, f, s , c , s, c} m .
The parton distribution functions obey the DGLAP equation 13 Here P a,â (z) are the standard (unregularized) DGLAP kernels, C a is either C F or C A as in eq. (7.11), and Following our default choice (7.2) for the argument of α s , we have used in the evolution equation for the parton distribution functions. This is not the standard choice, but it can be accommodated without changing the parton distribution functions by modifying the evolution kernel at next-to-leading order and beyond. That is, the terms indicated by O(α 2 s ) are modified from what they would have been had we used α s M 2 e −t . These terms do not affect our analysis.
With these results, we can write the complete V pert (t) as a sum and integral in the form used in eq. (6.3), This simplifies in two limits, which we now discuss. 13 If the parton distributions used in shower generation obey the second or higher order DGLAP equation, then there are more terms with extra powers of αs. We do not display possible higher order terms in order to keep the notation from becoming complicated, but in the end we will find that their inclusion would change the Z-boson transverse momentum distribution at a level that is beyond the accuracy that we aim for.

(1 − z)
1 Consider the limit (1−z) 1. Then the leading terms in V pert a contain a factor that equals 1/(1 − z) as long as 1 − z y; contributions without this factor can be neglected. There is no such factor unless f = g. Thus only f = g is important. In the factors with ratios of parton distributions, a + f = a and η a /z ≈ η a ; thus these factors are well approximated by 1. This gives In the second line on the right hand side of this equation, P a,a (z) contains a term proportional to 1/(1 − z). Other terms, which do not contain this factor, can be ignored. Similarly, we can ignore the term −γ a since it has no 1/(1 − z). In fact, the term in P a,a (z) proportional to 1/(1 − z) is 2C a /(1 − z), which cancels the term −2C a /(1 − z). This leaves the very simple result, Now consider the collinear limit (1 − z) ∼ 1. For the matrix element of H pert a in this limit, we use eq. (7.14). This gives The two terms involving ratios of parton distributions cancel, leaving the very simple result, Here γ E is the Euler constant. We will see later the reason for the factor exp(2γ E )/4. In this section, we make use of the results from the previous section to analyze the evolution before t c , e −t e −tc . Recall from eq. (6.14) that that the evolution is specified by operators that we called K a and K b , The operator K a is defined in eq. (6.15). We will analyze K a ; the analysis of K b is the same. In the end, K a is very simple. However, we will have to analyze it in stages, pruning away complications at each stage.

Distribution of final state partons
We are concerned with the evolution of the b-space partonic cross section, This is an inclusive quantity, but we first note something about the structure of the partonic final states that contribute to it. Recall from its definition (4.6) that Q contains a factor where p j,⊥ is the transverse momentum of the jth final state parton (other than the Zboson). When we form the b-space partonic cross section, we integrate over all of these transverse momenta. The integration region in which one or more partons have p 2 j,⊥ 1/b 2 gives a negligible contribution because the exponential factor exp(ip j,⊥ · b) averages to zero. It is as if Q contained a factor j θ(p 2 j,⊥ < C/b 2 ), where C is a constant that is large compared to 1 but with log C not large. That is, Q effectively projects onto partonic final states in which no partons have been emitted in the region above the line k 2 ⊥ ∼ 1/b 2 in figure 2.
One simple consequence of this concerns any real parton splitting that has the potential to change the flavor a or momentum fraction η a of parton "a". Such a splitting must be close to the line (1 − z) ∼ 1 in figure 2. To see this, note first that in order to significantly change η a it is necessary that (1 − z) not be tiny compared to 1. Second, in order to change the flavor a, the parton emitted into the final state must not be a gluon, but the splitting functions for these splittings do not have (1−z) → 0 singularities and hence have negligible probabilities of occurring with (1 − z) 1. We now note from eqs. (3.8) and (9.1) that Thus any emission with (1 − z) ∼ 1 and e −t e −tc has k 2 ⊥ 1/b 2 . We conclude that all emissions that contribute to the b-space partonic cross section for t < t c leave the flavor a of parton "a" unchanged and leave its momentum fraction η a approximately unchanged. Thus the initial conditions remain (approximately) true and the flavors a and b of the incoming partons do not change.
With eq. (9.5), eq. (3.23) for the splitting variable y reads In the following section, we note another consequence of the fact that Q effectively projects onto partonic final states in which no partons have been emitted with k 2

Angular ordering
In the range e −t e −tc and k 2 ⊥ ∼ 1/b 2 , both real and virtual emissions contribute to the evolution. For e −t e −tc and k 2 ⊥ 1/b 2 , only virtual emissions contribute. In either case, it appears that the analysis of the evolution will be very complicated. The virtual emission operator V pert a (t; z, φ, g) is given in eq. (8.12) in terms of the real emission operator H pert a (t; z, φ, g). The real emission operator is given in eqs. (7.4), (7.13), and (7.14). Here is where the complications are. There is a certain algebraic complexity and there is a non-trivial color structure. Worse, there is a sum over helper partons with labels k, with a different splitting function for each k. Thus, in order follow the evolution of the inclusive quantity b, Y ρ pert (t) , it appears that we need to track the structure of the complete partonic state. 14 Fortunately, there is a simplification available. We consider the real or virtual emission of a gluon with k 2 ⊥ 1/b 2 at shower time t with e −t e −tc . The rapidity r in the Z-boson rest frame of the gluon is approximately given by as long as e −t 1. Using eq. (9.5), this becomes Curves of constant r are shown in Figure 3. The rapidity of a previously emitted gluon (that is, a real gluon in the final state) is 14 This would not be quite so bad in the leading color approximation generally used for parton showers.
Then, we would need only the momentum of the parton that is color connected to the incoming quark or antiquark, which would generally be the gluon previously emitted from that incoming parton line. Here we have omitted the log (1/z k ) because any previously emitted gluon would have had (1 − z k ) 1, so log (1/z k ) ≈ 0. Since the gluon was previously emitted, we have t k < t. Since it was a real emission, the transverse momentum satisfied k 2 ⊥,k 1/b 2 . On the other hand, k 2 ⊥ 1/b 2 . Additionally, log(1/z) ≥ 0. From these inequalities, we conclude that r k r . In the event that k 2 ⊥ and k 2 ⊥,k are both of order 1/b 2 , this strong inequality may not hold. However, the emission probability per unit t and per unit log(k 2 ⊥ /M 2 ) is small, of order α s . Thus we can apply the low density approximation developed in section 5.3 to conclude that in this event we can approximate e 2t e 2t k (9.12) at the cost of affecting terms in the Sudakov exponent that we will obtain at the level of third-to-leading terms and beyond. Eq. (9.12) then implies eq. (9.11).
For this reason, we can make the approximation eq. (9.11) in the expression (7.4) for H pert a (t; z, φ, g). This gives the much simpler formula eq. (7.13), for which the details of the state of the shower at time t are not needed.
We are now prepared to write an evolution equation for 1 Q(b, Y ; η a , η b , a, b) ρ pert (t) for t < log(b 2 M 2 ). We consider first z 1, then z ∼ 1. Then we combine these two cases.

Contribution from (1 − z)
1 In this subsection, we consider K a in the region That is, we consider splittings that are represented by points well below the line (1 − z) ∼ 1 in figure 2 and well to the left of t = t c that passes through the lower right vertex of the triangle. In this region, gluon emission from line "a" is important because the splitting function has a term with a 1/(1 − z) factor. However, splittings with f = g lack this factor and may be omitted. With this approximation, eq. (6.15) becomes (9.14) Since 1 − z 1, we can also replace z by 1 in the exponential factor and in the argument of Q(b, Y ; zη a , η b , a, b). We can also use eq. (9.5) to replace x a /η a in the exponent by 1. Finally, we use eq. (8.14) for V pert a . With these changes, we have We can make one more simplification. The right hand side of eq. (9.15) vanishes when k 2 ⊥ 1/b 2 . In the alternative case that k 2 ⊥ 1/b 2 , we can apply the lesson of section 9.2 so as to use the simple form of 1 zH pert a given in eq. (7.13),  Here we have replaced y by e −t according to eq. (9.6). Thus 1 and e −t e −tc . Notice the key feature of eq. (9.17) that the object for which we seek an evolution equation, appears on the right hand side of this result. We can make one more simplification in this. Recall that we need the average over the emission angle φ of 1 K a ρ pert (t) . The only φ dependence is in the factor exp(i b·k ⊥ ). When we take the average over φ, we get a Bessel function: We will want to integrate this over z. Consider the integral between limits z 1 and can then be approximated by simply (1 − z). Similarly, the argument of α s is, using eq. (9.6) and then we can change the integration variable to k 2 ⊥ . This gives (9.20) Thus we encounter an integral of the form where the function f (k 2 ⊥ ) allows for the dependence of α s on k 2 ⊥ . As we have previously argued, if Q 2 1 and Q 2 2 are both much smaller than 1/b 2 , then we can replace [J 0 (|k ⊥ ||b|)−1] by zero since J 0 (|k ⊥ ||b|) ∼ 1 for |k ⊥ ||b| → 0. Similarly, we have argued that if Q 2 1 and Q 2 2 are both much larger than 1/b 2 , then we can replace [J 0 (|k ⊥ ||b|) − 1] by −1 because J 0 (|k ⊥ ||b|) is a rapidly oscillating function in the integration range. But what happens if Q 2 There is a simple approximation that is accurate provided that f (k 2 ⊥ ) is a slowly varying function. We can replace [J 0 (|k ⊥ ||b|) − 1] by −θ(|k ⊥ ||b| > 2e −γ E ), where γ E is the Euler constant. Thus the integral becomes The integrals in eq. (9.21) and (9.22) differ by terms proportional to a power of 1/(Q 2 2 b 2 ) or a power of Q 2 1 b 2 and by a term proportional to the second derivative of f (k 2 ⊥ ) with respect to log(k 2 ⊥ ) in the region near k 2 ⊥ ∼ 1/b 2 . The exact statement can be found in appendix B. Note that the second derivative of α s λ R k 2 ⊥ with respect to log(k 2 ⊥ ) is of order α s λ R k 2 ⊥ 3 , so that this is a good approximation for our purposes.
With this approximation, where the ∼ in this case indicates that this approximation works inside the integration over φ and z. Specifically, the derivation so far covers a range of z between limits e −t (1 − z 1 ) < (1 − z 2 ) 1. We have in mind that (1 − z 2 ) is chosen small enough that the approximation (1−z) 1 applies within the integration range. The region (1−z) ∼ 1 is treated separately in the following section.
Similarly, we choose (1 − z 1 ) = C 1 e −t , where C 1 1 but log C 1 is not large. Then the integration range 0 < (1 − z) < C 1 e −t needs a separate treatment, which we now provide. For that range, the denominator 1 − z + e −t in eq. (9.18) provides a lower cutoff on (1 − z) at around (1 − z) ≈ e −t . Thus the integration range that needs a separate treatment is really e −t (1 − z) < C 1 e −t . What happens here depends on the factor [J 0 (|k ⊥ ||b|) − 1] in eq. (9.18). Since In the remaining case, e −t (1 − z b ) C 1 e −t , we will use this same replacement. This is not an accurate approximation. However, this inaccuracy occurs only near the point at the intersection of the line (1 − z) = e −t and the line k 2 ⊥ = 1/b 2 in figure 2. Following an argument like that in section 5.3, we recognize that this inaccuracy does not matter because it does not lead to contributions with a large logarithm in the Sudakov exponent. Thus eq. (9.23) is also a sufficient approximation when applied inside an integration over z that includes (1 − z) → 0.

Contribution from
The approximations in the previous section cover the integration region (1 − z) 1 for e −t e −tc . For the integration region (1 − z) ∼ 1 with e −t e −tc , only virtual splittings contribute: For the matrix element of V pert a , we can use eq. (8.16). This gives As in the previous section, we notice that the object for which we seek an evolution equation, 1 Q(b, Y ; η a , η b , a, b) ρ pert (t) , appears on the right hand side of this result.

(9.28)
This gives the evolution equation, for e −t e −tc and assuming that parton "b" caries the opposite flavor from parton "a", b =ā, There is a factor 2 here because emissions from both initial state lines contribute equally when b =ā. This analysis could apply when a = b = g, which would be relevant for computing the transverse momentum distribution in Higgs boson production. However, the case of interest in this paper is that a is a quark flavor and b is the corresponding antiquark flavor, or vice versa. In this case, C a = C F and γ a = (3/2) C F . Then

Evolution for t > t c
After shower time t c , the evolution changes character. To see what happens, it is simplest to This result gives , does evolve for t > t c because initial state emissions can be collinear and thus change the momentum fractions and flavors η a , η b , a, b. What happens is that the partonic cross section evolves and the parton distribution functions evolve with opposite evolution kernels, so that the net change of 1 Q(b, Y ) ρ(t) with t vanishes.

Solution at t = t c
We can solve eq. (9.30) for the b-space partonic cross section with initial condition (4.17), evolving from shower time t = 0 to shower time t = t c . This gives where the Sudakov exponent S 0 is It is useful to change variables from t to k 2 We can now approximate S 0 (M 2 , b 2 ) in a fashion that makes it simpler. We are interested in the behavior of S 0 (M 2 , b 2 ) for M 2 b 2 1. In particular, we are interested in terms with logarithms log(M 2 b 2 ). The region of the k 2 ⊥ integration that dominates for large b 2 is k 2 ⊥ M 2 . In contrast, the region k 2 ⊥ ∼ M 2 cannot contribute logarithms log(M 2 b 2 ) and is thus not of interest for us. For this reason, we can make the approximation k 2 ⊥ M 2 inside S 0 (M 2 , b 2 ). Additionally, it is useful to expand α s (µ 2 R ) in powers of α s (k 2 ⊥ ). Then, after performing the z-integral, we obtain contributions to the integrand of the k 2 ⊥ integral of the form The leading terms in this counting are those with m = n. We need to keep track of those, along with the next-to-leading terms with m = n − 1. However, we will not be interested in the coefficients of terms with m ≤ n − 2. For this reason, we drop all such terms in the approximations below.
To proceed with this program, we note that the renormalization scale in α s is, using the definition (7.2) and then eq. (9.6), Expanding α s (µ 2 R ) in powers of α s (k 2 ⊥ ) and displaying only the order α s and α 2 s terms, we have (11.5) Here The z-integral that multiplies the order α s term is The most important term here is the logarithm log M 2 /k 2 ⊥ . The term −3/2 is next-toleading, so we keep it. The third term in eq. (11.7) contributes to the k 2 ⊥ integral only for k 2 ⊥ ∼ M 2 but is power suppressed in the dominant region k 2 ⊥ M 2 . Consequently, we neglect this term. This same integral multiplies the order α 2 s term proportional to log(λ R ). Here, we keep the large logarithm log M 2 /k 2 ⊥ but neglect the term −3/2 in which the extra power of α s is not multiplied by a large logarithm. For the remaining α 2 s term, we need the integral Here there is no logarithm, log M 2 /k 2 ⊥ , multiplying the extra power of α s , so we neglect this contribution entirely. An analogous analysis shows that we can neglect entirely the order α 3 s and higher order terms in the expansion eq. (11.5).

Result
We have seen that the b-space hadronic cross section 1 Q(b, Y ) ρ(t) evolves from t = 0 to t = t c ≡ log b 2 M 2 e 2γ E /4 and then stops evolving. Therefore the Fourier transform of the physical p ⊥ -space cross section is given by is given as a convolution of parton distributions and the b-space partonic cross section 1 Q(b, Y ; η a , η b , a, b) ρ pert (t c ) by eq. (4.16), We have found that the b-space partonic cross section at t = t c is given by where the Sudakov exponent S given in eq. (11.9). This has been obtained with the approximations discussed in section 5.3 and in section 11. In particular, terms in S of the form α n s log n−1 (M 2 b 2 ) are affected by the approximations. We can compare this to the result of ref. [17] for QCD. In ref. [17], there are two arbitrary parameters, C 1 and C 2 , that do not affect the result summed to all orders of perturbation theory, but do affect the perturbative expansion. Here, we take the simplest choice, C 1 = 2e −γ E and C 2 = 1. Additionally, the result is given for an arbitrary choice of the scale µ 2 = M 2 e −t that appears in the parton distribution functions. Here t should be close to t c . We choose t = t c . With these choices, the result of ref. [17] is that the hadronic b-space cross section is given by eq. (12.1) with ) .

(12.3)
Here A, B, and C have perturbative expansions in powers of α s of the indicated arguments. The first terms in these expansions are [17,22,23,24,25] A(α s ) = 2 C F α s 2π Here we follow the notation of [7] in defining We begin the comparison of the shower result (12.2) with the QCD result (12.3) by noting that they agree at the Born level, order α 0 s . This agreement was built in when we chose the Born cross section as the starting point for shower evolution.
We next note that the QCD cross section exponentiates in b-space. This is a statement about 1 Q(b, Y ; η a , η b , a, b) ρ(t c ) QCD . We consider the logarithm 15 of this function and expand it in powers of α s (M 2 ) and log(M 2 b 2 ). To obtain this expansion, we need to expand α s (k 2 ⊥ ) in powers of α s (M 2 ) and log(M 2 /k 2 ⊥ ) using the evolution equation for α s , then perform the integration over k 2 ⊥ in the Sudakov exponent S. Similarly, in the coefficients C, we need to expand α s 4e −2γ E /b 2 in powers of α s (M 2 ) and log(M 2 /b 2 ). The general form of the perturbative expansion would be We know that the maximum power m of log(M 2 b 2 ) must be no larger than 2n because each loop can give at most two logarithms, one from a collinear singularity and one from a soft singularity. If we look at the perturbative expansion of the partonic b-space cross section instead of its logarithm, then terms with m = 2n actually occur. By the statement that the shower cross section exponentiates in b-space we mean that, for all n, D nm = 0 for m > n + 1 . (12.7) 15 In taking the logarithm, we consider the functions C to be operators on the space of parton distribution functions. Thus the first term, δ a a δ(1 − z), represents the unit operator. The coefficients Dnm(Y ) in eq. (12.6) are then, in general, operators on the space of parton distribution functions. One could diagonalize these operators by taking moments, but there is little reason to do so.
This property of the QCD result (12.3) is shared by the shower result (12.2). For the shower, it is a simple consequence of having a differential equation in which the derivative with respect to shower time of the partonic b-space cross section is a kernel times this same partonic b-space cross section, where the kernel has one α s and one logarithm. One sometimes finds discussions of summations of large perturbative logarithms expressed in terms of the directly observed cross section, which can be written as a convolution of parton distribution functions and a partonic p ⊥ -space cross section defined analogously to the partonic b-space cross section. Thus one may write One can then discuss the coefficients E nm . However, this is not very useful. For instance, no condition on the E nm for m = 2n and m = 2n − 1 will imply eq. (12.7). It is for this reason that we compare the shower result and the QCD result with respect to the logarithm of the partonic b-space cross section.
Let us now compare the n = 1 terms in eq. (12.6) between the QCD result (12.3) and the shower result (12.2). We first note that the terms D 10 do not match between the two results: the order α s contribution to C a a (z, α s ) in the QCD result is lacking in the shower result. These terms arise in part from the one loop virtual graphs, which are not included in the shower splitting functions. For that reason, the shower does not get these terms right. In our analysis we have have, accordingly, completely neglected contributions to the evolution kernel that have a factor α s with no logarithms log(M 2 b 2 ).
We now compare the α 1 s terms with one or two powers of log(M 2 b 2 ). For the shower result (12.2), we have This is the same as D 12 and D 11 in the QCD result. It is not surprising that the leading coefficient, D 12 , matches. The contribution D 12 log 2 (M 2 b 2 ) is 4C F times the area of the triangle in figure 2. The contribution D 11 log(M 2 b 2 ) is more subtle. This contribution is associated with the edges of the triangle. The bottom edge of the triangle is associated with the lower limit on the k 2 ⊥ integral, which is 4e −2γ E /b 2 instead of simply 1/b 2 . In the shower evolution, this value arises as a property in integrals involving Bessel functions J 0 (|k ⊥ ||b|). In the analysis of the contribution to evolution associated with k 2 ⊥ near this limit, it was important that the rapidity of a potential gluon emission was large compared to the rapidities of previous gluon emissions. In the analysis of the contribution associated with the left hand edge of the triangle, it was important that the denominator of the splitting function is 1 − z + e −t and not, say, 1 − z + 2e −t . This denominator is related to the treatment of the soft gluon interference diagrams. Finally, the right hand edge of the triangle is associated with the term −3/2 in eq. (11.9) for the Sudakov exponent. This term comes from the −γ a in the virtual part of the DGLAP evolution equation for the parton distribution functions, eq. (8.9).
Let us now look beyond the order α s contribution to the Sudakov exponent. We note that the QCD result and the shower result share the contribution inside the In the shower result, this form reflects our choice to use µ 2 R defined in eq. (7.2) as the argument of α s ; this is a choice that is not required by the logic of the shower, in which we factor softer interactions from harder interactions. With this choice, we match the contribution in eq. (12.10) between the shower result and the QCD result. Expanding in powers of α s (M 2 ), we see that We note also that the QCD result and the shower result share the contribution inside the k 2 ⊥ integral, This term contributes to the coefficients D nm for n = m. However, the QCD result has a contribution α s k 2 which is to be compared to (12.14) It is not a surprise that these do not match for λ R = 1, since we have used a leading order shower and this is a second order contribution. However, this term in the Sudakov exponent contributes to the coefficients D nm for n = m and n ≥ 2. If we would like to match these terms, we can, by means of a trick introduced in ref. [7]. We choose λ R such that −4β 1 C F log(λ R ) = 2C F K . This scale choice applied to q → q + g splittings produces the second order term in A(α s ) for Z-boson production. The same formalism applies to the transverse momenta of Higgs bosons produced in hadron collisions. In that case, the incoming partons are two gluons instead of a quark and an antiquark. The Sudakov exponent has the same form, but with different functions A, B, and C. For A, the first two terms are now 2C A (α s /(2π)) + 2C A K(α s /(2π)) 2 , but the ratio of the first and second coefficients in A(α s ) is unchanged [4]. Thus the same factor λ R applied to g → g + g splittings produces the second order term in A(α s ) for the Higgs boson transverse momentum distribution. Of course, this is a trick that applies for the particular purpose at hand. It is not equivalent to using a second order splitting kernel.

Other choices
In this section, we briefly investigate the Z-boson transverse momentum distribution that would result if we were to make other choices for the construction of the shower.

Spin-averaged, leading-color shower
We have organized this paper as an analysis of the shower evolution equation in ref. [12], with some minor modifications. This evolution equation contains the effects of spin correlations and applies to arbitrary color states. It contains the physics that we believe should be approximately represented in a parton shower event generator. However, the nature of the evolution is such that finding a way to represent the equations in a practical computer program is challenging. To start with a base approximation, one can average over the spins and take the standard leading color approximation. We analyzed this case in ref. [13]. With a spin average and a leading color approximation for each splitting, the evolution equations have the right form to be implemented as a Markov process, as is commonly used for parton showers. 16 We are thus led to ask whether the results of this paper still hold if one averages over spins and takes the leading color approximation at each parton splitting step. The answer to this question is that the results do still hold.
Spin does not matter because the observable of interest, the Z-boson transverse momentum, is independent of parton spins. The spin dependence is represented in the spin factor in eq. (7.4), {s } m {s} m : the parton spins {s} m in the quantum amplitude must equal those in the quantum conjugate amplitude, {s } m , but the splitting probability is independent of what the spin values are. If we average over spins, this factor goes away.
Color is of some importance. However if we use the exact color dependence then eq. (7.10) allowed us to reduce the color dependence to a simple factor of C F or C A . If we use the leading color approximation, then each gluon is treated as carrying color 3 ⊗3 instead of color 8. The allowed parton pairs forming a dipole are then the color connected partners. Then eq. (7.10) still holds with C F replaced by C A /2. One can get back to C F by simply adjusting the quark-quark-gluon couplings in the splitting functions. With this understanding, the results of the previous section calculated with full color remain unchanged in the leading-color approximation.

Alternative choice of dipole partitioning function
Within the framework of a virtuality ordered shower as presented in ref. [12], there is a choice to make. Consider the emission of a soft gluon from a color dipole composed of the initial state parton "a" and one other parton, k. The total emission probability is proportional to Here we use the eikonal approximation, which is valid for very smallp m+1 . The shower algorithm divides this into pieces, H ak and H ka , defined by Here A ak and A ka are both positive and obey The contribution proportional to A ak is considered to be an emission from parton a and is treated using the momentum mapping associated with emission from parton a. The contribution proportional to A ka is considered to be an emission from parton k and is treated using the momentum mapping associated with emission from parton k. 17 Thus the function A ak tells how to partition the soft emission from a dipole into two parton splitting contributions.
In this paper, we use the A ak defined in eq. (7.12) of ref. [14], , (13.4) with A ka given by the same expression with a ↔ k. In thep a +p b rest frame, this is . (13.5) We have seen that with this choice the shower will produce the proper structure for the Drell-Yan transverse momentum distribution. One can consider using instead the A ak defined in eq. (7.13) of ref. [14], This is the form used by Catani and Seymour and by others for the purpose of partitioning the dipole splitting probability in defining perturbative dipole subtractions. The most important difference between the two choices is that A ak as given by eq. (13.4) is invariant under the scalingp k → λp k . Thus this A ak depends on the angle between the emitted gluon and the momentum of parton k, but not on the energy of parton k.
We have investigated the possibility of using A ak defined by eq. (13.6) to generate the parton shower. We find that if we do so, the argument given in sections 7.1 and 9.2 fails. Now the probability of emitting a gluon m + 1 from the initial state parton "a" with virtuality smaller than that of a previously emitted gluon k can depend on the momentum p k . As a result, we no longer obtain a differential equation for the b-space partonic cross section that has only a kernel times the b-space partonic cross section on the right hand side. The Z-boson will get less transverse momentum recoil from emitted gluons, but we cannot say how much less.
We can draw a lesson from this: a parton shower built on splitting functions that reflect the collinear and soft singularities of QCD does not automatically sum large logarithms that appear in a physical cross section of interest. Seemingly minor details like the choice of A ak matter.

Catani-Seymour dipole shower
The shower discussed in this paper is the virtuality ordered shower of ref. [12] with an improved momentum mapping for initial state radiation and with the renormalization scale for α s given in eq. (7.2). A close relative is the Catani-Seymour (CS) dipole shower, proposed in ref. [26] and implemented in refs. [27,28]. Here one takes the splitting functions and momentum mappings of the Catani-Seymour algorithm [21] for defining subtractions in NLO perturbative calculations and uses them instead to define a shower algorithm. We examine this choice here, with the evolution variable taken to be as in refs. [26,27,28]. The most important difference between the shower of this paper and the CS dipole shower lies in its treatment of momentum conservation. Let us consider the (backward) evolution of an initial state parton by emission of a gluon into the final state. The initial state parton is on shell with zero transverse momentum. It is replaced by a prior onshell parton with zero transverse momentum and the final state gluon that has non-zero transverse momentum. Thus there is a transverse momentum imbalance. In the shower of this paper, the recoil from the final state gluon is taken by all of the other final state partons collectively, by means of a Lorentz transformationk i = Λk i . With this choice, we have seen that the recoil transverse momentum is predominately taken by the Z-boson. In the CS shower, for each splitting there is a partner parton, the other parton in the CS dipole. Let us suppose that we use the leading color approximation, which is commonly applied in parton showers. Then the partner of a radiating initial state quark is the parton that is color connected to the quark. For the first gluon emission, the partner of incoming parton "a" is the other incoming parton "b". In this case, the CS scheme assigns the recoil transverse momentum to the Z-boson. However, once one gluon has been emitted, the partner of "a" is the previously emitted gluon. Then the CS scheme assigns the recoil transverse momentum to the partner parton. That is, the previously emitted gluon recoils against the newly emitted gluon and the Z-boson does not take any recoil.
This leads to a curious situation. With a transverse momentum ordered shower, the shower first evolves by not emitting any gluons. The probability not to emit any gluons up to a certain shower time t ⊥ is a Sudakov form factor. At some point, one real gluon is emitted. Then evolution of the Z-boson transverse momentum distribution stops because the recoil from further gluons emitted is taken by one of the gluons previously emitted. This problem has been analyzed and addressed in ref. [29].
We can formulate the issue in the style of our analysis of the previous sections if we supply one additional variable in the measurement operator Q, a variable recoil that takes values T and F. Including this variable, our operator expressing possible measurements is Q(b, Y ; η a , η b , a, b; recoil). If recoil = T, then when a gluon is emitted from either of the initial state partons, the Z-boson will share in the recoil against the gluon transverse momentum. If recoil = F, then, when a gluon is emitted from an initial state parton, the Z-boson does not share in the recoil.
The evolution equation of the recoil = T cross section is This does not involve the recoil = F cross section. Using the approximations developed in the preceding sections, we can solve this equation. We find (13.9) Note that for very large t ⊥ , the recoil = T cross section tends to zero. The evolution equation of the recoil = F cross section has only the recoil = T cross section on the right-hand side: (13.10) The total partonic b-space cross section is For very large t ⊥ , the recoil = T contribution vanishes and we are left with the recoil = F contribution, which we can obtain by solving eq. (13.10) with the use of eq. (13.9). The result for large t ⊥ is We see that the partonic cross section does not exponentiate in b-space. In fact, it is simplest in p ⊥ -space instead: if we take a Fourier transform back to p ⊥ , the factor J 0 (|k ⊥ ||b|/z) provides a delta function 2(2π) 2 δ(p 2 ⊥ − k 2 ⊥ /z 2 ). This leaves a relatively simple expression consisting of a Sudakov factor and a splitting function for the single allowed splitting.
Again, we see that a parton shower built on splitting functions that reflect the collinear and soft singularities of QCD does not automatically sum large logarithms that appear in a physical cross section of interest. Seemingly minor details like the momentum mapping matter.

Transverse momentum ordered shower
The shower examined in this paper is based on virtuality ordering. What would happen if we kept everything else the same but used transverse momentum ordering? That is, what if we replace the evolution time t of eq. (3.3) by t ⊥ defined in eq. (13.7)? We keep the momentum mapping the same as in the virtuality ordered shower, so that the transfer of recoil transverse momentum to the Z-boson is the same as in the main body of this paper. Transverse momentum ordering is used in Pythia [30,31], but the algorithm that we consider in this section differs in some respects from that of Pythia.
With k ⊥ ordering, the first stage of shower evolution is simple. For t ⊥ log(M 2 b 2 ), only virtual splittings contribute to the evolution of the b-space partonic cross section; the rapidly oscillating J 0 (|k ⊥ ||b|) factor that multiplies H pert I eliminates the real emission contribution.
The later stages of shower evolution are also simple. The argument of section 10 implies that the rate of change of the b-space hadronic cross section vanishes for t ⊥ log(M 2 b 2 ). We are left with the region t ⊥ ∼ log(M 2 b 2 ). This is the bottom of the triangle in figure 2. With a virtuality ordered shower, this region is sampled from left to right, which is in order of increasing rapidity of emitted gluons -that is, decreasing emission angles. Then the splitting function in eq. (7.4) took the simple form in which the complicated function f was replaced by 1 because the rapidity of the real or virtual gluon was much larger than the rapidity of previously emitted gluons. With t ⊥ as the evolution variable, there is no guarantee that a previously emitted gluon had smaller rapidity than a new real or virtual daughter gluon. Thus the form of the evolution equation changes and the evolution of the partonic b-space cross section depends on the previous emission history.
Because the evolution is more complicated than we have dealt with, we do not know what the result is. We do note that the difficulty arises only for a limited region, t ⊥ ∼ log(M 2 b 2 ).

Angle ordered shower
What would happen if we used angular ordering instead of virtuality ordering? That is, what if we replace the evolution time t of eq. (3.3) by t ∠ = − log(tan(θ/2)) where θ is the splitting angle. Then ordering in t ∠ is ordering in splitting angle, starting from large angles and progressing to smaller angles as t ∠ increases. At the same time, we would use a revised version of the dipole partitioning function A lk , as specified below. This gives us two of the main features of Herwig [32,33], although if we start with the shower described in this paper and change only these features, we certainly do not have exactly Herwig.
To be more precise, let P = x a p a + x b p b be the Z-boson momentum approximated as having no transverse part but as obeying P 2 = M 2 Z . Consider the splitting of a parton with label l and momentum p l and let n l be the lightlike vector Let the daughter partons have momentap l andp m+1 . Then define t ∠ = 1 2 log p m+1 ·n l P ·p l p m+1 ·p l P ·n l . (13.14) For an emission from initial state line "a", t ∠ is the rapidity of the splitting that we used in section 9.2. From eq. (9.8) we have, for y 1 and η a ≈ x a , Curves of constant t ∠ are shown in figure 3. Over most of the region of the figure, (1−z) 1 so log(1/z) ≈ 0. This gives straight lines in the t-log(k 2 ⊥ /M 2 ) plane. As these lines approach the line (1 − z) ∼ 1 in the figure, the log(1/z) term becomes important and the t ∠ curves bend.
The curves of constant t ∠ extend down to arbitrarily small virtuality, M 2 e −t . In a virtuality ordered shower, one would simply stop the evolution when M 2 e −t reaches some infrared cutoff at which application of a hadronization model is more appropriate than continued use of perturbative showering. For the angle ordered shower, we should impose this same cutoff as a upper limit on integrations over t at fixed t ∠ . Equivalently, if we parameterize the curves of fixed t ∠ using k 2 ⊥ , then we impose a lower limit on the integrations over k 2 ⊥ .
One should do more than simply switching the evolution variable. One can also modify the dipole partitioning function A ak as used in eq. (13.2). For an angle ordered shower, it is standard [6] to begin with where the angles are defined, for instance, in the Z-boson rest frame. These functions add to one, but are not everywhere positive. One can then replace the functions H a,k and H k,a defined in eq. (13.2) by their averages over the azimuthal angles defined by rotatingp m+1 about, respectively, the directions ofp a andp k . As shown in ref. [6], this gives a simple result, The procedure is equivalent to replacing A  A ak = θ(θ a,m+1 < θ a,k ) 1 − cos θ m+1,k 1 − cos θ a,k , A ka = θ(θ k,m+1 < θ a,k ) 1 − cos θ m+1,a 1 − cos θ a,k . (13.18) These functions are positive. They do not sum to 1, but, as just noted, one gets the right result after suitable averaging over azimuthal angles. It is now easy to see how the partonic b-space cross section Q(b, Y ; η a , η b , a, b) evolves as t ∠ increases. At each step dt ∠ , there are contributions from all values of k 2 ⊥ down to the infrared cutoff. We divide the integration over k 2 ⊥ into three distinct regions: k 2 In the region k 2 ⊥ 1/b 2 , the real emission term with a factor J 0 (|k ⊥ ||b|) cancels the virtual emission term with a factor 1 because J 0 (|k ⊥ ||b|) − 1 ≈ 0, as in our other analyses. Thus this region does not contribute to the evolution of Q(b, Y ; η a , η b , a, b). In the region k 2  ; η a , η b , a, b). The rapidly oscillating J 0 (|k ⊥ ||b|) factor that multiplies H perp I eliminates the real emission contribution. In the region k 2 ⊥ ∼ 1/b 2 , both real and virtual emissions contribute.
We have by no means carried out a detailed investigation of this angle ordered scheme. However, it appears to us that it will give equivalent results for the Z-boson transverse momentum distribution as does the virtuality ordered shower.

Dipole antenna shower
In this paper so far we have discussed dipole parton showers in which the creation of a new gluon is attributed to the splitting of one of the previously existing partons. There is an ambiguity because, for soft gluon emissions, interference diagrams are important: one has to include emission from a parton l squared, emission from a parton k squared, and the l-k interference graphs. Using the eikonal approximation, the sum of these is simple, as in eq. (13.1). We have partitioned the total emission probability into a fraction A lk associated with splitting of parton l and a fraction A kl = 1 − A lk associated with the splitting of parton k, as in eq. (13.2). In the nomenclature of section 16 of ref. [34], this constitutes a partitioned dipole shower. 18 There is a separate momentum mapping P l for each parton l that splits. Here P l is an operator, defined in ref. [12], on the space of parton states. It is not necessary that this momentum mapping depends on the choice of k, and the momentum mapping used in this paper and ref. [12] does not depend on k. Thus, in a highly compressed notation, we can write that the emission from the l-k dipole is partitioned into two terms, In section 13.2, we have seen that the summation of large logarithms can be sensitive to the choice of the partitioning function A lk ; a choice based on the angles between parton l and k and the emitted soft gluon is preferred.
One might be tempted to get rid of this ambiguity, particularly in the leading color approximation where one can use pairs of color connected partons as the dipoles. The idea is to consider each l-k dipole as a unit that can emit a gluon. Thus the basic building blocks are 2 → 3 parton splittings. A shower based on this approach may be called a dipole antenna shower. The pioneering development along these lines is the final state shower of Ariadne [35]. More recent examples include those in ref. [19]. There is a corresponding subtraction scheme for next-to-leading order calculations, antenna subtraction [36].
In a dipole antenna shower, there is no A lk . There is a separate momentum mapping P lk for each dipole l-k. Thus eq. (13.19) is replaced by The freedom to choose A lk now resides in the freedom to choose P lk . It must be symmetric under the interchange of partons l and k. It is usually defined in such a way that momentum is conserved locally in the 2 → 3 splitting, without taking momentum from any of the other partons. For a final state dipole, this means Such a mapping is defined in ref. [37]. We note that it is possible to have a shower that is simultaneously a partitioned dipole shower and an antenna dipole shower. We get that 18 Note that from the point of view adopted in section 13.5, an angle ordered shower also can be considered as a partitioned dipole shower.
case if we define P lk = θ(ϑ l,m+1 < ϑ k,m+1 ) P l + θ(ϑ k,m+1 < ϑ l,m+1 ) P k , (13.22) where the ϑ ij is the angle between momenta i and j and P l represents the momentum mapping that is used in this paper. Then A lk = θ(ϑ l,m+1 < ϑ k,m+1 ). As far as we can see, a dipole antenna shower can reproduce the proper summation of large logarithms for the Z-boson transverse momentum distribution studied in this paper. However, the momentum mapping for initial state splittings cannot be local in the sense of eq. (13.21). For a dipole consisting of the two initial state partons, eq. (13.21) becomes (13.23) We want p a,⊥ =p a,⊥ = 0 and also p b,⊥ =p b,⊥ = 0. This leaves no place for the transverse momentum of parton m + 1 to go. One might hope to use a non-local momentum mapping for this case, but use a local mapping in the case of a dipole consisting of one initial state parton, say "a", and a final state colored parton k. In this case, eq. (13.21) becomes However, that would mean thatp k,⊥ − p k,⊥ = −p m+1,⊥ : the recoil transverse momentum from the emission of parton m + 1 is taken up by parton k instead of by the Z-boson. Thus a non-local momentum mapping is needed for all dipoles that include an initial state parton.

Conclusions
In this paper, we have examined the differential cross section dσ/(dp ⊥ dY ) for producing a Z-boson in hadron-hadron collisions as a function of the transverse momentum p ⊥ and rapidity Y of the Z-boson. For p 2 ⊥ M 2 Z , the perturbative expansion of dσ/(dp ⊥ dY ) contains large logarithms, log(p 2 ⊥ /M 2 Z ). There are known QCD results for the summation of these logarithms, quoted in section 12. We have asked to what extent a virtuality ordered parton shower algorithm of the sort defined in ref. [12] (with some small modifications as discussed in this paper) correctly reproduces the known results.
The parton shower evolves in shower time t, defined to be minus the logarithm of the virtuality of a parton splitting, so that hard splittings come first, soft splittings last. The state of the shower at time t, in the sense of the distribution of parton configurations, is represented by a vector ρ(t) . In this notation, dσ/(dp ⊥ dY ) as it has developed by shower time t is denoted by p ⊥ , Y ρ(t) . Then, at the end of the shower at time t f , the predicted cross section is dσ dp We have examined whether the cross section obtained from the shower algorithm has the structure of the QCD result by examining how p ⊥ , Y ρ(t) evolves with t.
To find d p ⊥ , Y ρ(t) /dt, we used eq. (2.9), which expresses how the completely exclusive partonic state evolves. This expresses d p ⊥ , Y ρ(t) /dt in terms of the complete state of all of the partons at time t. Fortunately, we found that, with appropriate approximations, d p ⊥ , Y ρ(t) /dt is expressed as a convolution of a certain kernel with p ⊥ , Y ρ(t) . That is, the detailed configuration of all of the partons does not matter, only the inclusive distribution p ⊥ , Y ρ(t) matters. Thus we obtain a closed form differential-integral equation for p ⊥ , Y ρ(t) .
The result can be simply stated in terms of the Fourier transform of transverse momentum distribution, defined in eq. (4.4), We found that 1 Q(b, Y ) ρ(t) stops evolving at a certain evolution time t c , with t c = log b 2 M 2 e 2γ E /4 . We found in eq. (12.1) that 1 Q(b, Y ) ρ(t c ) is a convolution of parton distribution functions evaluated at the appropriate scale times a function that does not involve parton distribution functions, 3) The important structure is thus expressed in the function 1 . On general grounds, one would expect that this function has two powers of the large logarithm log(M 2 b 2 ) per power of α s (M 2 ): However, we found that the partonic cross section produced by the shower algorithm exponentiates in b-space in the sense that D nm = 0 for m > n + 1 .
This property is shared with the full QCD result. This exponentiation is not guaranteed by simply having a parton splitting probability that matches QCD in the soft and collinear limits, including the proper interference between emissions of a soft gluon from different partons. We found that two other features of the shower algorithm are important.
First, partons in a parton shower are treated as being on-shell, but it is not possible for an on-shell parton to split into two on-shell partons and still conserve momentum. There has to be a momentum mapping that takes a small amount of momentum from somewhere and supplies it to the daughter partons. In the case of an initial state splitting, the problem is more pronounced because the initial state partons are treated as having zero transverse momentum. This means that when a soft gluon is emitted from an initial state parton, some transverse momentum has to come from somewhere. We found that the momentum mapping has to be such that the Z-boson gets most of the recoil transverse momentum. In fact, the momentum mapping proposed in ref. [12] did not have this property and thus needed to be modified. Similarly, we noted that the momentum mapping used for perturbative subtractions in the Catani-Seymour dipole subtraction scheme does not give the desired exponentiation for 1 Q(b, Y ) ρ(t c ) in a dipole based shower.
Second, in a shower like that of ref. [12] that is based on color dipoles, there is a function that we call A lk that specifies how much of the radiation from dipole l-k is treated as being emitted from parton l and how much is treated as being emitted from parton k. We found that, for initial state emissions, the choice of A lk matters. A choice based on the angles between the momenta of partons l and k and an emitted soft gluon is satisfactory, while the choice used in the Catani-Seymour dipole subtraction scheme fails to give the proper exponentiation in a dipole based shower.
After having checked for the proper exponentiation of the large logarithms, the next question was which non-zero coefficients D nm are correctly produced by the shower algorithm.
We found that the leading coefficient, D 12 , is correctly reproduced. This result has a simple physical interpretation in the parton shower picture. The quantity 1 Q(b, Y ) ρ(t f ) is the cross section to produce the Z-boson with rapidity Y but not radiate initial state gluons with transverse momentum bigger than approximately 1/b 2 . The Sudakov exponent S is the integral of the differential probability to radiate such a gluon, so exp(−S) is the probability not to radiate any such gluons. The leading term in the Sudakov exponent is C F α s /π times twice the area of the triangle in figure 2. This area is total phase space for gluon emission with k 2 ⊥ > 1/b 2 from one of the two initial state partons; doubling the area gives the total phase space for emission from either line. Quantum coherence plays a role here. The initial emission is from a color dipole consisting of the two incoming quarks. Emissions with large positive rapidity count as emission from incoming quark "a" while emissions with large negative rapidity count as emission from incoming quark "b".
We found also that the next-to-leading log coefficient, D 11 , is correctly reproduced. This coefficient corresponds to the behavior of the parton splitting near the three edges of the triangle in figure 2. It thus represents physics that is more subtle than the physics behind the leading coefficient, D 12 .
We note that once the choices of momentum mapping and A lk have been made, no adjustment of parameters is needed to get D 12 and D 11 to match the QCD result.
One could argue that matching the leading order coefficients D 12 and D 11 is all that one should expect from a leading order parton shower. Nevertheless, we found that it is possible to do better (as in ref. [7]). With a suitable choice of the argument of α s in the splitting function, it is possible to match all of the coefficients D nm with m = n + 1 and with m = n. With this level of matching, the Sudakov exponent S produced by the shower algorithm is a good approximation to the full QCD exponent S QCD in the limit log(M 2 b 2 ) → ∞ and α s (M 2 ) → 0 with α s (M 2 ) log(M 2 b 2 ) fixed.
We have investigated how well the virtuality ordered parton shower algorithm defined in refs. [12,13,14] and this paper performs in summing large logarithms for the transverse momentum of a vector boson produced in the Drell-Yan process in hadron-hadron collisions. The same analysis would apply to the transverse momentum distribution of Higgs bosons produced in hadron-hadron collisions. There are other sorts of summations of large logarithms that are of importance for understanding experiments. We believe that the investigation of other cases, perhaps using methods developed in this paper, is an important goal, both for the virtuality ordered parton shower algorithm used in this paper and for other parton shower algorithms.

A. Structure of inclusive splitting
In this appendix, we review formulas from refs. [12,13,14] for inclusive splitting, that is the splitting operator H I (t) times the inclusive measurement function 1 . We start with in eq. (12.20) of ref. [12] The inclusive splitting function is diagonal in spin, but it has a non-trivial color structure. The structure of the splitting is contained in the operator h, which is given in eq. (12.21) of ref. [12]. With a slightly modified notation and with one error corrected, this equation reads Here possible quark masses are included. There is a sum over the index l of the initial state or final state parton that splits. There is a sum over flavor choices ζ f for the splitting and an integration over the momentum choices, dζ p . Once we choose splitting variables t, z and φ, this becomes an integration over the splitting variables. There is a delta function that expresses the definition of the shower time t. Next, there is a ratio of parton luminosities, containing parton distribution functions, as in eq. (4.8). The main momentum and color dependence is contained inside the braces that follow. The second term inside the braces applies when the newly emitted parton with label m + 1 is a gluon. There is a sum over partons k, where k = l. Partons l and k constitute a color dipole that can emit the gluon. There are functions w ll and w lk that describe the direct terms and the interference terms, respectively, for gluon emission. These functions are defined in ref. [12]. Convenient expressions for w ll and w lk as rational functions of dot products of the momenta involved are given in ref. [13]. 19 There is also a function A lk that expresses how the interference term is partitioned between a part considered to be a splitting of parton l and a corresponding part considered to be a splitting of parton k. There is a color operator T k · T l that inserts a gluon color matrix T a on line k and another on line l and sums over the gluon color index a. This operator was denoted g lk ({f } m+1 ) in the original equation.
The first term inside the braces applies when the newly emitted parton with label m + 1 is not gluon. Then there is no interference term. There is a color factor The cases withf m+1 = g do not appear in eq. (A.2), but we will need these cases later. In eq. (12.21) of ref. [12], we failed to note the possibility of an initial state splitting with {f l ,f m+1 } = {g, q} or {g,q}, so we listed this factor simply as T R . We now manipulate eq. (A.2) so as to obtain a more useful form. For thef m+1 = g term, we divide the splitting function into two parts by adding and subtracting the eikonal approximation to w ll . (The eikonal approximation is the limiting form when the gluon momentum tends to zero. The function w lk is already constructed in the eikonal approximation in ref. [12].) Thus we write The first term is important in only in the collinear limit. The momentum dependent factor w ll − w eikonal ll is independent of k and multiplies {c } m T k · T l {c} m . When we sum over k, we can use  19 In ref. [13] and also in ref. [14], we use the notations W ll = w ll and W lk = 2A lk w lk .
For the second term, we take A lk to be of the form specified in eq. (7.2) of ref. [14], in which A lk is expressed using another function A lk . This gives the result of in eq. (7.10) of ref. [14], We can now specialize to initial state splittings with massless partons. In the sum over the label l of the parton that splits, we take l = a. The integration measure dζ p , with the choice of variables {t, z, φ} used in this paper, is dζ p = 2p a · p b (2π) −2 dt dz dφ 2π y 4z =p m+1 ·p a 8π 2 z dt dz dφ 2π . (A.10) The integration over t is performed using the delta function that defines t. We have defined the matrix element of H pert a (t, z, φ, f ) in eqs. (4.15) and (6.2) to be the coefficient of