Multivariable evolution in parton showers with initial state partons

One can use more than one scale variable to define the family of surfaces in the space of parton splitting parameters that define the evolution of a parton shower. In an earlier paper, we developed this idea for electron-positron annihilation. Here, we use multiple scale variables for a parton shower with initial state partons. Then we need a more sophisticated analysis because the evolution of parton distribution functions must be coordinated with the parton shower evolution. We make the needed connections more precise than in our earlier work, even for the case of just one scale variable. Then we develop an example with three scale variables, which leads to advantages compared to the usual shower formulation with only one scale variable. We provide results for Drell-Yan muon pair production.

In a parton shower event generator, one can view the parton state as evolving according to an operator based renormalization group equation. Starting with a state with just a few partons, the shower evolves as a scale µ s changes from a large value µ h characteristic of the starting state to a low value µ f on the order of 1 GeV. As the shower evolves, more and more partons are emitted.
The function of the shower scale µ s is to divide possible parton splittings into resolvable splittings, with scales µ > µ s , and unresolvable splittings, with scales µ < µ s . There is substantial freedom to choose exactly what this means. The space of possible splittings is divided into the resolvable and unresolvable regions by a surface labelled by µ s . Many different choices are possible for defining this surface. For instance, one can use a measure of the transverse momentum in the splitting to define the surface or one can use a measure of the virtuality in the splitting.
In this paper, we explore the possibility of using more than one variable to define a family of surfaces within the framework of our parton shower program Deductor. Instead of one µ s , we use µ s = (µ s,1 , µ s,2 , . . . ). Then evolution means moving from large values of the component scales µ s,n to small values along a path µ s (t). Defining this path is then part of defining the shower algorithm.
There is an additional freedom available when multiple scales are involved. It may be possible to divide the shower splitting functions into separate terms such that one of the terms is not sensitive to one of the scales in the sense that no singularity is encountered when this scale approaches zero. When this happens, we can modify the definition of the unresolved region for this term in a way that makes this term exactly independent of this scale. This redefinition can simplify the shower evolution.
In this paper, we explore the additional freedom obtained by using three scales instead of one.
This general concept works for proton-proton, e ±proton, and e + e − collisions. In Ref.
[1], we considered the simplest case, e + e − collisions. In this paper, we consider collisions involving incoming hadrons, with an emphasis on hadron-hadron collisions. With incoming hadrons, the needed theoretical construction is more involved than in e + e − annihilation. First, the representation of the cross section needs both a shower operator U(t f , t h ) and an operator U V (t f , t h ) that sums threshold logarithms. Second, the description of initial state splittings using "backward evolution" requires the use of parton distribution functions (PDFs). Since the evolution of these functions must be coordinated with the evolution of the shower, they are not, in general, the MS PDFs used in fixed order perturbation theory. We refer to the PDFs used in the shower as shower oriented PDFs. We review this connection in some detail in this paper, expanding on our previous treatments.
Much of the analysis of this paper concerns what we call the infrared sensitive operator [2], D(µ r , µ s ), where µ r is the renormalization scale and µ s represents one or more shower scales that define the unresolved region. In a first order shower with just one shower scale and with µ r = µ s , we use the first order contribution D [1] (µ s , µ s ). This can be thought of as being the more familiar shower splitting operator, S [1] (µ ) integrated over µ from µ = µ s down to µ = 0. This integral would be infrared divergent, but we perform the integration in 4 − 2 dimensions to regulate the divergence. This makes D [1] seem like a derived quantity, but we consider D [1] to be the starting point of the definition of a shower [2]. Then it is the shower splitting operator S [1] that is the derived quantity. By using D as the starting point, we connect to the definitions of renormalization, factorization, and parton distribution functions used for the calculation of cross sections at any perturbative order in QCD.
We introduce the analysis in this paper with a review in Sec. II of the vector space of parton states used in the Deductor framework. In Sec. III, we sketch the role of the infrared sensitive operator D and, with some important notation established, we provide in Sec. IV a preview of what is technically new in this paper.
We turn in Sec. V to a detailed analysis of the structure of D. There is another operator that we need, the infrared finite operator V. We analyze its structure in Sec. VI. With the operators D and V available, we turn in Sec. VII to the operator U(t , t) that generates the probability preserving parton shower and the operator U V (t , t) that corrects for the mismatch between shower evolution and PDF evolution and thereby generates a summation of threshold logarithms. This leads us to the definition of the shower oriented PDFs in Sec. VIII.
In Sec. IX, we start with one scale that represents k T ordering, Deductor's standard Λ ordering, or angular ordering. Then we add two additional scale parameters, µ e that controls soft splittings and µ iπ that controls the imaginary part of virtual exchange graphs. Next, we see how one can modify the definition of the unresolved region for parts of the generators of U(t , t) and U V (t , t) that are sensitive only to some of the scale parameters but not others. In Sec. X we suggest a special three component choice of path in the space of scale parameters. In Sec. XI we show how the operators used to express an infrared safe cross section behave on the suggested path.
We provide a numerical example in Sec. XII: the Drell-Yan cross section for muon pair production via vector boson production. We offer a summary and some comments in Sec. XIII.
We also include a number of appendices, A, B, C, D, E and F, that document derivations and formulas that are treated only briefly in the main text.

II. PARTONS AND THE STATISTICAL SPACE
This paper generally concerns the definition of the unresolved region in the space of parton momenta in a parton shower and then how the definition of the unresolved region affects the parton shower. We will concentrate on emissions in a first order shower, but we include a discussion of the general case of a shower algorithm with splitting functions defined at an arbitrary order of perturbation theory. We use the general framework presented in Ref. [2]. This general framework allows for substantial freedom in choosing the functions that define a particular parton shower algorithm. We have developed a partic-ular realization of these choices for a first order shower [3], a realization that is in many ways similar to other first order parton shower algorithms. It remains an open problem to realize these choices for a parton shower with splitting functions beyond order α s .
The description that we use is based on linear operators on a vector space that describes the state of the partons in the shower, which we call the statistical space [3]. We begin by recalling the nature of this space.
We consider the description of a cross section in hadron-hadron collisions.
At a particular stage in the shower, there are m final state partons plus two initial state partons with labels "a" and "b." The partons have momenta and flavors {p, f } m = {p a , f a ; p b , f b ; p 1 , f 1 ; p 2 , f 2 , . . . , p m , f m }. We define Q by (1) Then also Q = p a + p b . The theory is expressed using linear operators that act on a vector space that we call the statistical space. Basis vectors for this space have the form |{p, f, c, c , s, s } m ). Here (c, c ) and (s, s ) represent the quantum colors and spins of the initial and final state partons. We use the apparatus of quantum statistical mechanics, with the color and spin part of |{p, f, c, c , s, s } m ) representing the density matrix |{c, s} m {c , s } m |.
In order to properly incorporate quantum mechanics into a parton shower, one must include both quantum color and spin. It is quite straightforward to include both color and spin in the evolution equations for a first order shower. For color, one then needs an approximation to obtain a description that can be implemented in a practical computer program. Deductor uses the "LC+" approximation [4]. Thus our description in this paper uses full color, with the LC+ approximation playing a special role in Secs. IX C, X, and XI. A proper implementation of spin [5,6] is replaced by simple spin averaging in Deductor. Because of this, we omit spin quantum numbers in everything that follows in this paper. This results in a simpler notation without seriously obscuring conceptual issues. Using quantum color but not spin, the basis vectors for the statistical space have the form |{p, f, c, c } m ), where the color quantum numbers represent the density matrix |{c} m {c } m |. We use the trace basis for the vectors |{c} m [3].
In order to use the statistical space to express a cross section, we need the probability associated with a vector |ρ) in the statistical space. We define a vector (1| so that the probability associated with |ρ) is (1|ρ) [3]. We define (1| by using the probability associated with a basis state |{p, f, c, c } m ): 1 {p, f, c, c } m = 1 pf {p, f } m 1 color {c, c } m (2) with 1 pf {p, f } m = 1 ,

III. THE INFRARED SENSITIVE OPERATOR
The shower description in Ref. [2] begins with an operator on the statistical space, D(µ 2 r , µ 2 s ), that we call the infrared sensitive operator. This operator includes a specification of what one means by an unresolved splitting in the shower. In this section, we introduce this operator and the operators that are derived from it.
The operator D(µ 2 r , µ 2 s ) describes the soft and collinear singularities of QCD. Here µ r is the standard renormalization scale and µ s is called the shower scale. In this paper, we contemplate the possibility of having more than one independent shower scale, µ s = (µ s,1 , µ s,2 , . . . ). The infrared sensitive operator is expanded in operators D [n] (µ r , µ s ) that are proportional to α n s : The operator D [n] is further expanded as Here D [nr,nv] creates n r real emissions and n v virtual exchanges. After one of these operators acts on a state |{p, f, c, c } m ), we have partons with momenta and flavors {p,f }m withm ≥ m. The new total momentum ism i=1p i =Q .
The first order contribution to D forms the basis for a first order shower. It consists of In D [1,0] l with the replacement a ↔ b, so it suffices to describe the case l = a.
The infrared sensitive operator D provides the basis for the definition of the probability preserving parton shower evolution operator ( [7], Eqs. (27), (51), and (54), or [2]) Here the scale parameters evolve along a path µ r (t), µ s (t) as a function of the shower time t and T indicates ordering in t. The generator operator S at lowest perturbative order is related to parton distribution functions supplied by an operator F(µ 2 r ) = F a (µ 2 r )F b (µ 2 r ) defined in Eq. (A11) and to the infrared sensitive operator at lowest perturbative order, D [1] .
The first order generator of the probability preserving shower consists of three terms, The operator in the second term, which we denote by S [1,0] (t) P , leaves the number of partons and their momenta and flavors unchanged. As described in Appendix A, it is determined from the real emission operator S [1,0] (t) by integrating over the splitting variables for a splitting that might have happened, but did not. With this definition, U(t 2 , t 1 ) conserves the total probability of the statistical state: 1 U(t 2 , t 1 ) = 1 . The operator in the third term also leaves the number of partons and their momenta and flavors unchanged. It is obtained by differentiating the imaginary part of the virtual exchange operator D [0,1] : This operator does not change probabilities: 1 S [0,1] iπ (t) = 0. We return to these relations in more detail in Sec. VII.
The infrared sensitive operator D also provides the basis for the definition of the inclusive infrared finite operator V(µ r , µ s ) [2,10], which we examine starting in Sec. VI. The structure of this operator is crucial to the definition of shower oriented parton distribution functions, examined in Sec. VIII. It is also needed for the definition of the operator U V (t, t ) that produces the summation of threshold logarithms and is examined in Sec. VII.

IV. NEW FEATURES IN THIS PAPER
The structure of the shower cross section with initial state splittings is more complicated than with just final state splittings because this structure has to include parton distribution functions and the interplay between PDF evolution and shower evolution. For this reason, we provide a rather detailed analysis in what follows. We extend the analysis of Refs. [2,[8][9][10] in several ways. First, we allow for the presence of more than one independent shower scale: µ s → µ s . Second, we present a more precise derivation of the needed relation between standard MS PDFs and the PDFs needed internally in the shower. Third, we provide an improved determination of the shower oriented PDFs from the MS PDFs. Fourth, we write the real emission operator D [1,0] a (µ r , µ s ) in what we think is a clearer notation. Fifth, we provide a simpler definition of the real part of the virtual exchange operator Re D [0,1] a (µ r , µ s ) than appears in Refs. [2,9,10]. Sixth, we provide a simpler definition of the imaginary part of the virtual exchange operator Im D [0,1] (µ r , µ s ) that appears in Refs. [4,11].

V. STRUCTURE OF THE INFRARED SENSITIVE OPERATOR
The infrared sensitive operator includes a specification of what one means by an unresolved splitting in the shower. We cover what this means in several steps in this section.

A. Kinematics of parton splitting at first order
In the real emission operator D [1,0] a (µ r , µ s ), initial state parton "a" splits, in the sense of backward evolution [12]. Before the splitting, the momentum of this parton is where p A is the momentum of the incoming hadron (with the approximation p 2 A = 0) and η a is the momentum fraction carried by the incoming parton. The new initial state parton carries momentum A new final state parton carrying momentump m+1 witĥ p 2 m+1 = 0, is produced. We define the momentum fraction in the splitting by We define a dimensionless virtuality variable for the splitting by We also denote by φ the azimuthal angle ofp m+1 around the direction of p a in the rest frame of Q. The three splitting variables y, z, φ suffice to define the splitting. The momenta {p} m+1 of the partons after the splitting are determined by y, z, φ and the momenta {p} m , as described in Appendix A. It is sometimes useful to define a transverse momentum variable for an initial state splitting, as in Eq. (A.8) of Ref. [9]. This is approximately the absolute value of the square of the part ofp m+1 transverse to p a and Q. 1 We will specify D [1,0] a (µ r , µ s ) in some detail in Appendix A, but for now these details do not matter. What is important is that D [1,0] a exhibits collinear and soft singularities. To describe these, it is useful to define an angular variable This is (1 − cos θ)/2 where θ is the angle betweenp a and p m+1 as measured in the rest frame ofQ. The operator D [1,0] a is singular in the collinear limit ϑ → 0 with fixed z, in the soft limit (1 − z) → 0 with fixed ϑ, and in the soft×collinear limit (1 − z) → 0 and ϑ → 0. In the integrations in D [1,0] a , these singularities are regulated with dimensional regularization.
We can also write ϑ as a function of z and k 2 T , where We will use the variables z and ϑ to describe a splitting. Then y and k 2 T can be determined from (z, ϑ) using B. The unresolved region The idea of the singular operator D(µ r , µ s ) is that, when applied to a parton basis state |{p, f, c, c } m ), it produces approximated cut Feynman diagrams that have the same infrared singularities as actual cut Feynman diagrams. Of course, it does not match the behavior of Feynman diagrams when the momenta are not close to one of these singularities. For this reason, we define an unresolved region, a bounded region in the parton momentum space that surrounds the singularities of interest. For a parton splitting included in D [1,0] a (µ r , µ s ), the unresolved region is a region U ( µ s ) in the space of the splitting variables (z, ϑ). Since the angle variable ϑ defined in Eq. (17) has the range 0 < ϑ < 1, we always take the unresolved region U ( µ s ) to be a subset of the region The shape of the region depends on parameters µ s . We define a family of unresolved regions by letting µ s depend on a parameter t with 0 < t < t f such that the resolved region shrinks as t increases. Often t f = ∞ is a convenient choice. The singular lines 0 < z < 1 at ϑ = 0 and 0 < ϑ < 1 at z = 1 must be included in U ( µ s (t)) for all t. This concept is discussed in some detail, with examples, in Ref. [1].
In this paper, we will impose a fixed cutoff on k 2 T using a parameter m 2 ⊥ that is on the order of 1 GeV 2 . A splitting in which will always be counted as unresolved. This cutoff will serve to keep resolved parton splittings in a range for which perturbation theory is applicable. For an initial state splitting involving a heavy quark or antiquark, we impose an additional restriction on the unresolved region. Here we count c and b quarks as heavy with m 2 ⊥ < m 2 c < m 2 b , while we consider other quarks to be massless. When a massive parton with flavor a turns into a gluon in the sense of backward evolution (so, going forward in time, g → a with an emittedā) we count the splitting as unresolved if k 2 T < m 2 a . This is consistent with switching from a 5 flavor scheme for m 2 b < µ 2 to a 4 flavor scheme for m 2 c < µ 2 < m 2 b and then a 3 flavor scheme µ 2 < m 2 c in MS parton distribution functions. Although the definition of the unresolved region involves the masses of heavy quarks, the shower splitting functions approximate the quark masses by zero. 2 We can summarize this by saying that an initial state splitting with a →â in backward evolution is always counted as unresolved when 2 It would be better to use on-shell charm and bottom quarks with masses in the shower evolution and then use corresponding quark mass dependence in the PDF evolution equations [13]. However, the current version of Deductor treats charm and bottom quarks as massless except in the functions that define the unresolved region for initial state splittings. With this approach, the parton shower misses important effects when the scales µ are not much greater than the quark masses. where In order to further specify the family of unresolved regions, define a function a cut (z, µ s (t)) such that (z, ϑ) ∈ U ( µ s (t)) if 0 < z < 1, 0 < ϑ < 1, and As t increases, a cut (z, µ s (t)) decreases. In this paper, we impose some conditions on this function. We require that a cut (z, µ s (t)) → 0 for t → t f . We further require that for t → t f , a cut takes a simpler limiting form, where µ lim is a function µ lim ( µ s ) of the scales µ s , with for t → t f . We require that the limiting function a lim (z, µ 2 lim (t)) takes the factored form Thus the scale parameter µ lim controls the collinear limit of a cut when ϑ → 0 with fixed z. It is useful to choose for large t.
The function a lim (z, µ 2 lim (t)), together with the choice of µ 2 r (t), specifies the unresolved region and its relation to the renormalization scale in the collinear limit, in which a cut 1 for fixed z. For the definitions of the unresolved region considered in this paper, a lim (z, µ 2 lim (t)) has the very simple form shown in Eq. (28). This simple form enables us to define shower oriented PDFs as a function of a single scale µ 2 r . For not-so-large values of t, or not-so-small values of a cut , we can allow a more general structure. First, we can maintain the relation (29) in the large t limit while modifying µ 2 r (t) away from this limit. Additionally, for small (1 − z) with t not large, we can have a lim (z, µ 2 c (t)) > 1. We therefore allow a more general form for a cut (z, µ s (t)): All that we need is that h(z, µ s (t)) → 0 at least as fast as µ 2 r (t)/Q 2 in the large t, small µ 2 r (t)/Q 2 limit.

C. Simple cases for the unresolved region
Before we examine the relation between D [1,0] a (µ r , µ s ) and the shower, it may be helpful to provide some examples of the unresolved region in cases with only one shower scale, µ 2 c . With only one scale, we let µ lim (µ c ) = µ c and where we let a c (z, µ 2 c ) be given by one of three choices, a ⊥ (z, µ 2 ⊥ ), a Λ (z, µ 2 Λ ), or a ∠ (z, µ 2 ∠ ). That is, in a c (z, µ 2 c ), C =⊥, Λ, or ∠. In each case, we take the corresponding renormalization scale to be µ 2 r = µ 2 c . k T definition: We can base the unresolved region on the transverse momentum variable defined in Eq. (16). A splitting is unresolved if k 2 T < µ 2 ⊥ , where µ 2 ⊥ is the shower scale. We choose µ 2 r = µ 2 ⊥ . Accounting also for the fixed m 2 ⊥ (a,â) cutoff, we adopt the definition that where a ⊥ (z, µ 2 ) is defined in Eq. (19). Λ definition: The default ordering variable in Deductor [14] is where Q 0 is the total momentum of the final state partons at the start of the shower. Thus, letting η (0) a and η (0) b be the momentum fractions of the initial state partons at the start of the shower, where A splitting is considered unresolved if Λ 2 < µ 2 Λ , where µ 2 Λ is the shower scale. We choose µ 2 r = µ 2 Λ . Accounting also for the fixed m 2 ⊥ (a,â) cutoff, we say that (z, ϑ) ∈ U (µ 2 Λ ) if 0 < z < 1, 0 < ϑ < 1 and where a ⊥ (z, µ 2 ) is defined in Eq. (19) and ϑ definition: We can use the angle variable ϑ as the ordering variable. 3 Then we consider a splitting to be unresolved if ϑ Q 2 < µ 2 ∠ , where µ 2 ∠ is the shower scale. We choose µ 2 r = µ 2 ∠ . Accounting also for the fixed m 2 ⊥ (a,â) 3 We could use the angle measured in the fixed reference frame defined by Q 0 , but then the formulas are more complicated.
cutoff, we say that (z, ϑ) ∈ U (µ 2 ∠ ) if 0 < z < 1, 0 < ϑ < 1 and where a ⊥ (z, µ 2 ) is defined in Eq. (19) and In the context of a shower cross section, the state is always accompanied by PDFs that are part of the probability associated with that state. The PDFs are supplied by an operator F a (µ 2 r )F b (µ 2 r ), Eq. (A11). For initial state splittings of parton "a" realized by backward evolution, we need to remove the prior PDF factor by applying the operator F −1 a (µ 2 r ). Then after the splitting we supply the new PDF factor by applying the operator F a (µ 2 r ). Thus it is useful to consider the operator F a D [1,0] a F −1 a . We apply D [1,0] a (µ r , µ 2 s ) with its accompanying PDF operators to an m-parton state and write the result in the form Here D [1,0] a adds one new parton and we integrate over the momenta and flavors {p,f } m+1 of the partons after the splitting. There are dimensionally regulated singularities, so this integration is in 4−2 dimensions for each momentum. ThenD a aâ is a function of the momenta and flavors before and after the splitting and is an operator that maps the m-parton color space to the m+1 parton color space. The index a inD a aâ is the flavor of the incoming parton "a." The operatorD a aâ also depends on the scales µ 2 r and µ s and it depends on , but here we do not make that dependence explicit in the notation.
The operatorD â aa takes the form The pieces in this formula are described in Appendix A.
Here we mention only the most important features. The operatorD a aâ depends on the splitting variable z and on a, which is the flavor of the incoming parton after the splitting (in the sense of backward evolution). The righthand side of Eq. (41) begins with integrations over the other two splitting variables ϑ and φ, with appropriate dependence on the dimensional regulation parameter . There is a theta function that requires (z, ϑ) to be in the unresolved region U ( µ s ) according to the shower scales µ s . After the integrations, there is a delta function that sets {p,f } m+1 to the momenta and flavors obtained from a splitting with variables (z, ϑ, φ,â) applied to partons with momenta and flavors {p, f } m according to Deductor conventions.
Deductor is a dipole shower. There is a sum over the index k ∈ {a, b, 1, . . . , m} of a dipole partner parton for the splitting of parton "a." For the case that k = a, there is a splitting probabilityP aâ (z, ϑ, ), Eqs. (A18), (A19), (A20), and (A21). For = ϑ = 0, this is the standard Dokshitzer-Gribov-Lipatov-Alterelli-Parisi (DGLAP) splitting function. The function N (a,â) provides a color factor, Eq. (A15). For k = a, we have a gluon emitted from parton "a" in the ket state interfering with a gluon emitted from parton k in the bra state, or the same configuration with bra and ket interchanged. The function W 0 depends on ξ ak = p a · p k /p k · Q. This function is described in detail in Appendix A. It has the important property that W 0 /[(1 − z)ϑ] is singular when the gluon is soft, but not when it is collinear to p a .
The final factor in Eq. (41) contains color operators that act on the ket color state |{c, c } m ) to give the linear combination of new color states |{ĉ,ĉ } m+1 ) that one gets after emitting the new parton m+1. In t † a ⊗ t k , t † a acts on the ket color state |{c} m to give a new ket color state t † a |{c} m and t k acts on the bra color state {c } m | to give a new bra color state {c } m |t k . When parton m+1 is a gluon, the color operators obey the identity We have used this identity to rewrite the operator D [1,0] a used in Deductor [9] in what we think is a more transparent form.

E. Inclusive probability for real emission
We have presented the matrix element to obtain a particular state |{p,f ,ĉ,ĉ } m+1 ) produced by D . We also need the inclusive probability for a splitting starting from the state {p, f, c, c } m . Using Eqs. (2) and (3), the probability corresponding to F a D We write this using another operatorP a aâ as where (1 color | times the operatorP a aâ is Here we have used the momentum conserving delta function inD a aâ to eliminate the integration over {p,f } m+1 . In the color factor, we have used the instruction in Eq. (3) to take the trace of the color density matrix after the splitting. The flavors in the color factor are determined by the flavor indices a andâ. The argument {p} m ofP a aâ refers to the dependence of W 0 on the variables ξ ak .
We can simplify the color here. In the case that k = a, where N (a,â) is the Casimir eigenvalue (A15) appropriate to the flavor content of the splitting. When k = a, the emitted parton is always a gluon. The gluon line attaches to line "a" with a color generator matrix T c a in the 8, 3 or3 representation according to the flavor of parton a. The gluon line attaches to line k with the appropriate generator matrix T c k . Then we sum over the gluon color index c. The result can be denoted by T k · T a . Thus for k = a, These simplifications give us This specifies 1 color P a aâ (z; {p} m ) but not the operatorP a aâ (z; {p} m ). We note that Thus we can define the color content ofP a aâ (z; {p} m ) bŷ This particular choice of the color operator in the second term inP a aâ gives this operator the color structure of a virtual gluon exchange between partons "a" and k. . We now would like to define the real part of D [0,1] a (µ r , µ s ). We consider the imaginary part of D [0,1] (µ r , µ s ), for exchanges involving both initial state and final state partons, in the following subsection.
The operator D [0,1] a comes from virtual graphs, in which we integrate over a momentum q that flows around a loop. This operator describes the infrared singularity structure when q → 0 or q becomes collinear with p a [2]. Since D does, however, change colors. First, it contains terms from self-energy insertions on one of the parton legs. These terms are proportional to the unit operator on the color space times C F or C A . Second, there are terms coming from gluon exchanges between two parton legs. These terms are proportional to either [T k · T a ⊗ 1] for a virtual graph on the ket amplitude or [1 ⊗ T k · T a ] for a virtual graph on the bra amplitude. The virtual graphs have 1/ 2 and 1/ poles. By using the identity (42), we can arrange that the terms proportional to the unit operator on the color space have 1/ 2 and 1/ poles, while the terms with [T k · T a ⊗ 1] and [1 ⊗ T k · T a ] color operators have only 1/ poles that arise from the exchange of a soft gluon.
Since Re D where a is the flavor of the incoming parton "a." The operator Γ a ({p} m ) has a soft divergence that should be regulated. As we will see, this divergence is cancelled. Now, we need to define Γ a . The probability associated with Re D The probability associated with Im D [0,1] a vanishes [9]. Thus We can combine the probabilities (44) associated with D where We now propose a definition for the operator Γ a . Both the real and virtual terms in P a aa (z; {p} m ) contain 1/ 2 and 1/ poles. However, as we will see in Appendix B, collinear factorization and real-virtual cancellations lead to This fixes the pole part of Γ a ({p} m , ) but leaves its finite part undefined. It is not straightforward to impose an ultraviolet cutoff on the unresolved region for virtual graphs that matches the cutoff that we used for real emission graphs. In Ref. [9] we proposed a method for this.
Here, we propose a simpler method that gives essentially the same result. We impose a momentum sum rule on P a aa (z; {p} m ): This defines Γ a ({p} m ): We can write the value of P with the virtual contributions determined by the momentum sum rule (MSR) as where [P a ] msr is defined using a subtraction at z = 1 on P a considered as a matrix in aa . For a generic operator A aa (z; {p} m ), the definition is and Here [· · · ] + denotes the usual + prescription, G. Imaginary part of the virtual exchange operator The operator D [0,1] , in which a virtual parton is exchanged between two partons, has an imaginary part when the two partons are either two final state partons or the two initial state partons. Our calculation here builds on that in Ref. [4]. We write the imaginary part of For a final state interaction, illustrated in Fig. 1, let us denote the parton momenta before the exchange by k l and k k . For an initial state scattering, k l and k k denote the parton momenta before the exchange in the sense of backward evolution. The partons with momenta k l and k k connect to a graph representing much harder interactions in which k l and k k are approximated as p l and p k .
To obtain the imaginary part of the graph, we replace the propagators for the partons with momenta k l and k k by delta functions, i/(k 2 l + i ) → 2πδ(k 2 l ) and i/(k 2 k + i ) → 2πδ(k 2 k ). Thus the imaginary part of the diagram corresponds to elastic parton-parton scattering. We integrate over the scattering angle θ in the p l + p k rest frame. The integral is logarithmically divergent at θ → 0. We regulate the divergence using dimensional regularization with a scale µ 2 IR and impose a cut so that we integrate over an "unresolved" interval with boundary x lk that surrounds the singularity. Then, neglecting contributions suppressed by a power of x lk , the integration gives There are two terms, proportional to T l · T k ⊗ 1 and 1 ⊗ T l · T k , which correspond to a gluon exchange in the ket state and in the bra state, respectively. If we had color factors 1 ⊗ 1, these terms would cancel. With the color factors that we have, Im D [0,1] can potentially have an effect. Now we can define a scale parameter µ 2 iπ ( µ s ) that depends on the shower scale or scales µ s (t), so that then µ 2 iπ depends on t. If there is only one scale, µ c , as described in Sec. V C for k T ordering, Λ ordering, or ϑ ordering, then we can make the simplest choice, µ 2 iπ = µ 2 c . We let the scattering angle cutoff x lk be related to µ 2 iπ by where f lk iπ is any function of the external parton momenta that we choose. Then We can use the relation (11) between Im D [0,1] and the corresponding contribution S Notice that S

VI. THE INCLUSIVE INFRARED FINITE OPERATOR
In order to define a parton shower operator U(t f , t h ) starting with the operator D(µ r , µ s ) that describes the infrared singularities of QCD, the construction of Ref. [2] defines, in its Eqs. (79) and (92), an inclusive infrared finite operator V(µ r , µ s ). 4 The operator V leaves the number of partons and their momenta and flavors unchanged 4 In papers prior to Ref. [2], we used the notation V to denote a different operator.
and is related to D by The operator V at first order contains contributions from emissions from each of the two initial state partons: b . There is no contribution from emissions from final state partons because of cancellations between ] is an operator that multiplies by the bare PDFs. This factor, omitting K, would be present in any perturbative calculation of a cross section in which PDFs are factored from the hard scattering cross section [2]. The circles, a • b, represent convolutions in the PDF momentum fraction variables. The operator that performs MS subtractions of 1/ n poles is Z F (µ 2 r ). The factor K(µ 2 r ) transforms from the MS scheme for the PDFs to the shower scheme used in the PDFs represented by F(µ 2 r ), Eq. (A11). We will see in what follows why we need a shower scheme instead of the MS scheme for the PDFs.
The right-hand-side of Eq. (70) is infrared finite, even though D(µ r , µ s ) is singular, for two reasons. First, we form an inclusive probability by multiplying by 1 , which allows real-virtual cancellations. Second, the initial state singularities that do not cancel in this fashion are removed by Z F (µ 2 r ) in the PDF factor. As required in Ref. [2], we need to ensure that V(µ r , µ s ) is not only finite after dimensional regularization is removed, but also that it approaches the unit operator in the limit of small scale parameters. Specifically, we need V(µ r , µ s ) → 1 when µ r → 0, µ s → 0 and the infrared cutoff parameter is removed, m 2 ⊥ (a,â) → 0, while µ r and µ s are related as in Eq. (29) and α s (µ 2 r ) is held constant. Without these conditions, the operator U V that is constructed from V and sums threshold logarithms [2] would be contaminated by effects with a 1 GeV scale. We will see that these conditions can be achieved by making a suitable choice of K(µ r ).
The operator V can operate non-trivially on the parton color space [2]. We will define the color content of the first order contribution to V in what follows.
We will need the detailed formula for the first order contribution to V from emissions from parton "a" in Deductor. We will write this formula in the form Here Here the term proportional to P aa (z)/ is the first order contribution from Z F (µ r ). The term K aa (z, µ 2 r ) allows the shower to use a shower-oriented PDF instead of an MS PDF: Here g ms a/A (η, µ 2 r ) obeys the MS evolution equation. The remaining term in Eq. (72), P a aa , gives (1|FD After some manipulation of P a aa in Eq. (72), we reach a result given in Eq. (C18). Using this result, Eq. (72) becomes The kernelP ( ) aa (z) is given in Eq. (A20). The contribution proportional to P aa (z)/ in Eq. (72) has cancelled an identical singular term in P a aa (z; {p} m ; ) that appears in Eq. (C18). After this cancellation, the second and third terms in Eq. (74) remain and an infrared "non-sensitive" operator P a,NS aa remains. The operator P a,NS aa is an operator on the color space, while the first three terms in Eq. (C18) are proportional to the unit operator [1 ⊗ 1] in color. The operator P a,NS aa is given in Eqs. (C17) and (C19). It has the crucial property that it has no 1/ poles and vanishes in the limit of small scales µ 2 r , µ s and small infrared cutoff m 2 ⊥ (a, a ) when integrated against a fixed test function h(z). This insensitivity to effects at small momentum scales is exactly what we want in V a aa (z; {p} m ). We are left with the second and third terms in Eq. (74). These terms are not singular, but they are not small when the scales µ 2 r , µ s and the infrared cutoff m 2 ⊥ (a, a ) become small. In order to cancel these infrared sensitive terms, we set is insensitive to effects from small momentum scales, as desired. We will see in Sec. VIII how Eqs. (73) and (75) can be used to produce PDFs that can be used in the shower evolution.

VII. THE SHOWER AND THRESHOLD OPERATORS
We now outline, very briefly, how the operators introduced above are used to define a parton shower.
First define ( [7], Eq. (18), or [2],) Then the shower evolution operator is ( [7], Eqs. (27), (51), and (54), or [2]) Here (µ r (t), µ s (t)) defines a path in the space of the scale variables as in Sec. V B and T indicates ordering in the path parameter t, the shower time. The operator S(t) is the generator of shower evolution, defined by To use this in a first order shower, one evaluates S(t) at order α s and then exponentiates it. We also define the threshold operator ( [7], Eqs. (56), (58), and (60), or [2]), where Here t h is the shower time at the start of the shower, corresponding to scale parameters comparable to the scale of the hard interaction considered. Then t f is the shower time at the end of the shower, corresponding to scales on the order of 1 GeV. It is important that V(µ r (t f ), µ s (t f )) ≈ 1. Then the integration in the exponent of Eq. (80) is dominated by scales near the hard scale, while contributions from scales near 1 GeV are power suppressed. To use U V (t f , t h ) in a first order shower, one evaluates S V (t) at order α s and then ex- With the use of these operators, the cross section for an infrared safe observable is, as described in some detail in Ref. [2], Here ρ h is the parton statistical state at the hard interaction and F(µ 2 h ) supplies the parton distribution function factor, using PDFs that match the organization of the parton shower. Then U V (t f , t h ) provides a factor that sums threshold logarithms. It does not change the number of partons or their momenta or flavors. Next, U(t f , t h ) provides a parton shower, creating many partons, while preserving the total probability obtained by summing inclusively over the parton states. Finally, O J makes the desired measurement on the resulting parton state, and 1 specifies an inclusive sum over the parton state variables.

VIII. SHOWER ORIENTED PDFS
We now examine the shower oriented parton distribution functions f a/A (η, µ 2 r ) defined in Eq. (73) in the simple cases of k T ordering, Λ ordering, and ϑ ordering as specified in Sec. V C. The shower oriented PDFs are different for different choices of shower ordering. This difference plays a role in how the cross section sums threshold logarithms, as outlined in Appendices D and E. Ref. [15] presents a different view of the connection between PDFs and parton showers.

A. kT ordering
We define shower oriented parton distribution functions f a/A (η, µ 2 r ; λ), where λ = 0 refers to k T ordering. Rearranging Eq. (73), we have, neglecting terms of order where Consider first the case that µ 2 Thus f a/A (η, µ 2 r ; 0) is close to g ms a /A (η, µ 2 r ), but there is a small change because we do not renormalize the k Tordered parton distribution functions by simply integrating over all k T and subtracting poles. Now consider the case that m 2 ⊥ < m 2 c < µ 2 r < m 2 b . Then there is a small complication. We take the operator Z F (µ 2 r ) in Eq. (70) to be the operator that renormalizes parton distribution functions according to the MS prescription in five flavor QCD. Then g ms a/A (η, µ 2 r ) is the five flavor MS parton distribution function for flavor a. With the conventional definition in which there are no intrinsic b-quarks, g ms b/A (η, m 2 b ) = 0, but (applying five flavor Since P aa (z) is the kernel for PDF evolution in µ 2 , this contribution to K aa (z, η, µ 2 r ; 0) cancels the evolution of the five flavor MS PDF in the region µ 2 r < m 2 b and gives us That is, f a/A (η, µ 2 r ; 0) is the four flavor MS PDF as conventionally defined for m 2 We now generalize the definitions in the previous subsection by defining shower oriented parton distribution functions f a/A (η, µ 2 r ; λ), where λ = 0 refers to k T ordering and λ = 1 refers to Λ ordering. Rearranging Eq. (73) gives Using Eqs. (19) and (37) in Eq. (75), we now have, where Eq. (90) follows from Eqs. (19) and (37) for λ = 0 and λ = 1. For 0 < λ < 1, it is a simple interpolation.
We recall from Eq. (35) that r a = η/η a is the momentum fraction of parton "a" at the start of the shower. The appearance of r a has two effects. First, it makes the kernel K depend on both η and z instead of just z. Second, it introduces an external parameter η (0) a into K.
For Λ ordering, we need f a/A (η, µ 2 r ; λ) at λ = 1. To find this starting with f a/A (η, µ 2 r ; λ) at λ = 0, we consider λ in the range 0 ≤ λ ≤ 1. To find f a/A (η, µ 2 r ; λ) while working around the complications produced by the presence of r a in K, we first define an auxiliary scale variable, Then we define an auxiliary parton distribution function that is, in the end, simpler than f a/A (η, µ 2 r ; λ), From its definition, we havẽ For λ > 0, we use Eqs. (89), (90), and (93) to givẽ Since P aa (z) is the generator of scale changes for g ms , the sum of the first two terms is g ms a/A (η, µ 2 ) up to order α 2 s corrections. In the remaining terms, we can replace the scale r −λ a µ 2 by just µ 2 at leading order in α s . This givesf The sum of the first two terms and the λ = 0 contribution from K aa givesf a/A (η, µ 2 ; λ) at λ = 0, while in the last term, we can replace g ms byf a/A (η, µ 2 ; 0) at leading order in α s . This gives us Now, differentiating with respect to λ gives the differential equation where We can solve Eq. (98) forf a/A (η, µ 2 ; 1) withf a/A (η, µ 2 ; 0) as the boundary value. Then the shower oriented PDF at λ = 1 is Eq. (98) results from the interpolation that we chose in Eq. (90) between K [1] (µ 2 , 0) and K [1] (µ 2 , 1). It is of interest to know how much the solution of Eq. (98) for f a/A (η, µ 2 ; 1) depends on this choice of interpolation. To this end, we note thatf a/A (η, µ 2 ; 1) can be expressed, using a more compressed operator notation, as where we exponentiate using the convolution product • and T represents ordering in λ. We have noted that the evolution kernel is the derivative of the first order term of the operator K. We have included the ordering instruction T because d K(µ 2 , λ)/dλ is a matrix in the parton flavor space and the matrices for different values of λ do not commute. However, at λ = 0 and λ = 1, the most important terms in K are the terms with 1/(1 − z) singularities. These terms are diagonal in flavor and do commute. Assuming the terms with 1/(1 − z) singularities in K for 0 < λ < 1 are diagonal in flavor, we can conclude that to a reasonable approximation, the ordering instruction T is not needed and, for µ 2 > m 2 ⊥ (a, a ), Thus, to a reasonable approximation, only the endpoints matter, not the interpolation. How doesf (µ 2 , 1) evolve under changes of µ 2 ? At leading order in α s , we can omit the ordering instruction T in Eq. (101), so that we can use Eq. (102) forf (µ 2 , 1). Forf (µ 2 , 0), we have the ordinary first order DGLAP equation, Differentiating Eq. (102) then gives where That is, the first order kernel for evolution in µ 2 of f (µ 2 , 1) is the familiar DGLAP kernel but with a cut We write Eq. (104) in detail as The kernel P aa (z, µ 2 ) in Eq. (106) was introduced in Ref. [9], but without enforcing the momentum sum rule. This gives us two ways to determinef a/A (η, µ 2 ; λ) at λ = 1. We can start withf a/A (η, µ 2 ; 0), which is g ms a/A (η, µ 2 ) with a small correction given by Eq. (86). For g ms a/A (η, µ 2 ) one can use a PDF set that is fit to data and uses at least NLO evolution in µ 2 . Then we can solve the λ-evolution equation (98) to findf a/A (η, µ 2 ; 1). We know the kernel for this differential equation only at first order, Eq. (99), but the change in λ is not large. Alternatively, we can note thatf a/A (η, is known and we can use the µ 2 -evolution equation (107) to determinẽ f a/A (η, µ 2 ; 1) at larger values of µ 2 . We know the kernel for this differential equation only at first order, Eq. (106), but the lack of higher order corrections should not be a problem if we are interested in values of log(µ 2 /m 2 ⊥ ) that are not too large.

C. ϑ ordering
We define shower oriented parton distribution functions f a/A (η, µ 2 r ; λ), where λ = 0 refers to k T ordering and λ = 1 refers now to ϑ ordering. Rearranging Eq. (73), we have, up to order α 2 s , Here K is a different kernel than the one used for Λ ordering. Using Eqs. (19) and (39) in Eq. (75), we have, Now, differentiating f with respect to λ gives the differential equation where We can solve this for f a/A (η, µ 2 r ; 1) with f a/A (η, µ 2 r ; 0). As with Λ ordering, we can also write a first order equation for the evolution of f a/A (η, µ 2 r ; 1) as µ 2 r varies with ϑ ordering: Here for ϑ ordering P is

D. Modified argument of αs
In the differential equation (98) for the λ dependencẽ f a/A (η, µ 2 ; λ), we have taken the argument of α s to be the scale µ 2 . We made the same choice in Eq. (107) for the µ 2 evolution off a/A (η, µ 2 ; λ) at λ = 1. Similarly, in Eqs. (110) and (112), we evaluated α s at the scale µ 2 r of the shower oriented parton distribution functions. These are the natural choices for a formalism that is designed to make sense at any perturbative order. However, in this paper, we have only a first order shower. We can modify the argument of α s so as to include some desired higher order corrections in the evolution kernels.
In Eq. (98), we include a factor Here K g is the standard factor [16], with n f set to the number of active quark flavors at scale k 2 T . We note that ∂ K aa (z, µ 2 , λ)/∂λ contains a factor θ(k 2 T > m 2 ⊥ (a, a )), which protects us from a singularity of α s (k 2 T ). This substitution, based on Ref. [16], can help to sum large logarithms [7]. An example, for threshold logarithms, is analyze in Appendix E.

E. Numerical values of the PDFs
In this subsection, we look at the shower oriented PDFs in a proton. We obtain these PDFs from the corresponding MS PDF by first applying the P ( ) transformation (86) to obtain the k T ordered PDF. Then we solve either Eq. (98) to obtain the Λ ordered PDF or Eq. (110) to obtain the ϑ ordered PDF. In the case of Λ ordering, we choose r a = 1 or, alternatively, examine the PDF f a/p (η, µ 2 ) with a modified scale parameter, Eq. (100).
In Fig. 2, we plot the up-quark momentum distribution η f u/p (η, µ 2 ) as a function of η at a large value of the scale, µ ≈ 1 TeV. We show first the MS distribution as a dashed curve. There are quite a lot of up quarks around η ∼ 0.1 because the up quark is a valence quark. The shower oriented distribution for k T ordering is almost the same as the MS distribution. The shower oriented distributions for Λ ordering and for ϑ ordering are somewhat smaller for η < 0.1. It is difficult to see the differences Up-quark distributions η f u/p (η, µ 2 ) at µ = 1007 GeV: the MS distribution and the shower oriented distributions for kT ordering, Λ ordering, and ϑ ordering.
among ordering choices for η > 0.1 because the distributions themselves are small. To make these differences visible, we plot in Fig. 3 the ratios of the three shower oriented PDFs to the MS PDF. We see that for Λ and ϑ ordering, these ratios rise quite steeply as η increases. As described in Appendix E, part of the summation of threshold logarithms is contained in the shower oriented parton distribution functions. This threshold logarithm summation produces the rise at large η that we see in Fig. 3.
In Fig. 4, we plot the up-quark momentum distribution at a value of the scale, µ ≈ 50 GeV, that is much larger than 1 GeV but not nearly as large as in Fig. 2. We plot the distributions down to a smaller minimum value of η, in keeping with the smaller value of µ. We show the MS distribution and the shower oriented distributions for k T ordering, Λ ordering, and ϑ ordering. In Fig. 5, we plot the ratios of the three shower oriented PDFs to the MS PDF. The pattern is similar to what we saw at a 1 TeV scale, but somewhat more pronounced because α s is larger.
In Fig. 6, we plot the bottom-quark momentum distribution η f b/p (η, µ 2 ) at the scale µ ≈ 50 GeV. As before, we show the MS distribution and the shower oriented distributions for k T ordering, Λ ordering, and ϑ ordering. The shape of the MS bottom quark distribution is different from that of the up quark distribution because bottom quarks all come from g → bb splittings. However, the relationships among the curves for η < 0.1 is quite similar to the relationships for up quarks: the distribution k T ordering is slightly larger than the MS distribution, while the distributions for Λ ordering and ϑ Ratio of shower oriented up-quark distributions to the MS up quark distribution at µ = 1007 GeV: shower oriented distributions for kT ordering, Λ ordering, and ϑ ordering. The rise at large η illustrates the PDF part of the summation of threshold logarithms.
ordering are smaller.
In Fig. 7, we plot the ratios of the three shower oriented PDFs to the MS PDF. Here, we see something quite different from what we saw for up quarks. In the cases of Λ ordering and ϑ ordering, the ratios of shower oriented distributions to the MS distribution do not rise with increasing η. Instead, they fall. Here, we must face the fact that we are using a first order evolution kernel to go from k T ordering to Λ ordering or ϑ ordering. This is plausibly justified if the change in the parton distribution is reasonably small, say When the condition (117) is violated, we judge that higher order contributions to the evolution kernel are needed. The regions for which the lowest order evaluation appears to be untrustworthy by this criterion are indicated by dashed curves in Fig. 7.

F. Which PDFs to use
Assume that we base the shower cross section on one of the options in Sec. V B: k T ordering, Λ ordering, or ϑ ordering. (Later in this paper, we propose more complicated choices.) Then we use the corresponding definition of shower oriented parton distribution functions as described above in this section. The PDFs appear as the operator F(µ 2 h ) in the cross section formula (82). We still have a choice to make. One possibility is to define the shower oriented PDFs using evolution in the parameter λ. For this, we begin with the MS PDFs and apply the P ( ) transformation (86) to obtain the k T ordered PDF. Then we use evolution in λ: we solve either Eq. (98) to obtain the Λ ordered PDF or Eq. (110) to obtain the ϑ ordered PDF. This produces the PDFs exhibited in the previous subsection.
Alternatively, we can use µ 2 evolution. We begin with the MS PDFs at the starting scale m 2 ⊥ , adjust these with the P ( ) transformation, and then let the PDFs evolve from µ 2 = m 2 ⊥ to higher values of µ 2 using the evolution equation (106) and (107) for Λ ordering or (112) and (113) for ϑ ordering or just the ordinary DGLAP equation for k T ordering.
The MS PDFs are produced using next-to-leading order evolution in µ 2 . However, both the λ evolution and the µ 2 evolution used to define shower oriented PDFs use just leading order evolution kernels. With this limited accuracy, the results from µ 2 evolution can be noticeably different from the results from λ evolution.
For µ 2 m 2 b , we expect the results using λ evolution to be the most reliable, except for charm and bottom quarks at very large values of η. 5 In the operators F(µ 2 h ) and U V (t f , t h ) in Eq. (82) for the cross section, the PDFs are evaluated at or not far from the scale µ h of the hard scattering. For this reason, we use the λevolution version of the shower oriented PDFs in these operators. In the shower evolution operator U(t f , t h ) in Eq. (82), the PDFs are evaluated at all scales µ 2 down to m 2 ⊥ . At small scales, we find that the PDFs for charm and bottom quarks obtained using λ evolution are badly behaved, while the PDFs obtained from µ 2 evolution are well behaved and nicely matched to the evolution of the shower. For this reason, we use the µ 2 -evolution version of the shower oriented PDFs in U(t f , t h ).

IX. MULTIPLE SCALES
In Eqs. (40) and (41), we have described the action of FD [1,0] a F −1 {p, f, c, c } m resulting from an initial state splitting using an operatorD a aa that is a function of z and the momenta and flavors of the partonic states and is an operator on the partonic color state. In general, this operator has both soft and collinear singularities. Our first task in this section is to divideD a aa into terms with different singularity structures. imation with respect to color, even though the Deductor code then treats color approximately. The default color approximation in Deductor is the LC+ approximation [4], 6 in whichD a aa is replaced by an approximate operatorD a,LC+ aa . The terms inD a aa in which the dipole partner parton is the same as the emitting parton, k = a, have a simple color structure. In the LC+ approximation, these k = a terms are unchanged. In the terms proportional to the function W 0 in Eq. (41), which have k = a, some contributions are dropped. The terms proportional to W 0 inside the integrations are singular in the soft limit, (1 − z) → 0 at fixed ϑ, but not in the collinear limit, ϑ → 0 at fixed (1 − z), as we see explicitly in Eq. (A29). Thus the difference operator is singular in the soft limit but not in the collinear limit [4]. These terms are proportional to δ aa , so a = a . Thus we can use Eq. (118) to decomposeD a aa in the form   In this section, we propose an alternative in which we define different unresolved regions for the "LC+" operators and the "soft" operators based on the differences in their singularity structures. We then use Eq. (119) and (120) as the definitions ofD a aa and P a,NS aa .

B. Scales for LC+ operators
In our proposed choice for the unresolved region for D a,LC+ aa andP a,LC+ aa , the evolution is described by three scales, The third scale, µ iπ , is the scale for the imaginary part of virtual graphs, regarded now as an independent scale. This scale controls the color phase operator, iπS iπ (t), as specified in Eq. (69), but does not affect real emissions. The scale µ c will describe part of the evolution of real emissions. The subscript in µ c denotes one of the scale choices defined previously using a c (z, µ 2 c ) with C = ⊥ for k T as the ordering variable (Eq. (19)), C = Λ for Λ as the ordering variable (Eq. (37)), or C = ∠ for ϑQ 2 as the ordering variable (Eq. (39)). We add an additional scale, µ e , that will provide a cut on the energy of an emitted parton.
With scales µ e and µ c , we define the unresolved region U (µ e , µ c ) to be used for the operatorsD a,LC+ aa andP LC+ aa as follows. We say that (z, ϑ) ∈ U (µ e , µ c ) if 0 < z < 1, 0 < ϑ < 1 and as in Eq. (25), where Here we use the definition (19), (37), or (39) for a c (z, µ 2 c ) and define The inclusion of a e (z, µ 2 e ) in Eq. (122) means that a splitting (z, ϑ) is unresolved when (1 − z) < µ 2 e /Q 2 . Thus µ e controls the approach to the soft limit by putting an upper bound on the momentum fraction (and thus the energy, E) of the emitted parton. The detailed specification of the functions that we need inD a aa andP a aa are given in Appendix F.
The unresolved region in the plane of 1 − z and ϑ is illustrated in Fig. 8 for Λ ordering, C = Λ, with r a µ 2 c /Q 2 = 0.1 and µ 2 e /Q 2 = 0.2. The unresolved region is shaded blue, while the resolved region is shaded yellow. We illustrate a potential m 2 ⊥ (a, a ) cut with m 2 ⊥ (a, a )/Q 2 = 0.002. In this example, the m 2 ⊥ (a, a ) cut has no effect.

C. A special treatment for soft splittings
We adopt a different strategy for the soft operatorŝ D a,soft aa andP soft aa , which have (1 − z) → 0 singularities but no collinear, ϑ → 0, singularities. We define the unresolved region U (soft; µ e ) forD a,soft aa andP soft aa by letting (z, ϑ) ∈ U (soft, µ e ) if 0 < z < 1, 0 < ϑ < 1 and

X. A PATH FOR THE THREE SCALE PARAMETERS
We have defined an unresolved region (z, ϑ) ∈ U (µ e , µ c ) for scale parameters, (µ e , µ c ), for real emission graphs and a separate unresolved region (1 − cos θ) < µ 2 iπ f lk iπ ({p} m ) for the integral that gives the imaginary part of virtual graphs.
We now need to specify a path µ(t) in the space of µ = (µ e , µ c , µ iπ ), with t in the arbitrarily chosen range 0 < t < 3. First, we should specify the endpoints. at t = 0, for k T or Λ ordering, we first choose µ 2 c (0) to be a scale µ 2 h appropriate to the hard scattering that initiates the shower. For ϑ ordering, we choose µ 2 For µ e (0), we need an appropriate matching scale. For this purpose, define z h as the solution of This gives 1 − z h that is not close to 0. For ϑ ordering, Eq. (126) is satisfied for any z h in the range 0 < z h < 1. Thus for ϑ ordering, we can choose 1 − z h = 1.
There are many paths that one might choose that connect (µ e (0), µ c (0), µ iπ (0)) and (µ e (3), µ c (3), µ iπ (3)). We make a particular choice in which we change each of the three scales one at a time as we change t from 0 to 1, then from 1 to 2, and finally from 2 to 3. The path µ(t) is illustrated in Fig. 10. We let µ e (t) decrease from µ e (0) to 0 as t increases from 0 to 1. For 1 < t < 3, µ e (t) remains at 0. For example, we can take µ e (t) = (1 − t)µ e (0) 0 < t < 1 0 We let µ c (t) retain its initial value, µ h , for 0 < t < 1. Then we let µ c (t) decrease from µ h to 0 as t increases from 1 to 2. For 2 < t < 3, we let µ c (t) remain at 0. For example, we can take Finally, we let µ iπ (t) remain at µ h for 0 < t < 2. Then we let µ iπ (t) decrease from µ h to m ⊥ as t increases from 2 to 3. For example, we can take Let us consider the three operators U(1, 0), U(2, 1), and U(3, 2) corresponding to the three segments of this path. In each case, U is determined by dD [1,0] (t)/dt or d Im D [0,1] (t)/dt on that path segment, according to Eqs. (9), (10), and (11).
We find in Eq. (F15) that, for 0 < t < 1, (132) This is independent of t. Thus when we construct the generator S [1,0] (t) of shower evolution by differentiating D [1,0] (t) with respect to t,D a,soft aa contributes butD a,LC+ aa does not.
On the second segment, 1 < t < 2, µ iπ is fixed, so there is no contribution from Im D [0,1] (t). Also, µ e (t) is fixed at µ e (t) = 0. The unresolved region forD a,soft aa depends only on µ e according to Eq. (125), so there is no contribution fromD a,soft aa . On the second segment, µ 2 c decreases from µ 2 h to 0. InD a,LC+ aa , the unresolved region is determined by µ 2 c (t) with µ 2 e set to zero. Thuŝ D a,LC+ aa does contribute to U(2, 1). With one notable qualification, this gives us a shower based on k T ordering (C = ⊥), Λ ordering (C = Λ), or ϑ ordering (C = ∠) from Sec. V C. The qualification is thatD a,LC+ aa uses the LC+ approximation for color.
On the third segment, 2 < t < 3, µ 2 e and µ 2 c are fixed, so there is no contribution from D [1,0] (t). Since µ 2 iπ decreases from µ 2 h to m 2 ⊥ , there is a contribution from Im D [0,1] (t). Thus there is a color phase factor iπS We thus arrive at what seems to us quite a surprising structure. The operator U(2, 1) is the most important of these evolution operators. It represents a full first order shower from the hard scale µ 2 h to zero using k T , Λ, or ϑ ordering and using the LC+ approximation for color. The operator U(1, 0) provides an exponentiated correction to the LC+ color approximation. It is generated by an operator that produces soft gluon emissions and exchanges using the difference between color matrices for full color and the corresponding color matrices for LC+ color. The remaining operator, U(3, 2) produces an exponentiation of the color phase from the imaginary part of virtual graphs, generated by S iπ (t).

XI. THE CROSS SECTION
We now describe the components of an infrared safe cross section when we use three scale parameters µ = (µ e , µ c , µ iπ ) with the special path described in the previous section in which µ e changes first, then µ c , then µ iπ . The scale µ c can correspond to any of k T , Λ, or ϑ ordering. We consider any infrared safe cross section, but concentrate on cross sections for which the beginning hard scattering process is either Drell-Yan muon pair production or the hard scattering of two partons. The measured cross section could involve more than just the beginning hard process. For instance, one could measure the transverse momentum of the muon pair in the Drell-Yan process.
The parton shower representation for the cross section is described in some detail in Ref. [2] and is stated in Eq. (82).
If we decompose U(3, 0) as U(3, 2) U(2, 1) U(1, 0) and decompose U V (3, 0) as At the end of the shower, O J makes the desired measurement on the many-parton state. Then 1 represents an inclusive sum over the parton state variables. At the start of the shower, ρ h is the parton statistical state at the hard interaction. It can be computed beyond the leading perturbative order if we use infrared subtractions that are matched to the shower [2]. Since this matching is not yet implemented in Deductor, we concentrate in this section on the case that ρ h is calculated at lowest order. The operator F(µ 2 h ) supplies the parton luminosity factor needed to form a cross section. It uses the PDFs that match the organization of the parton shower according to the ordering represented by µ c , calculated using evolution in the parameter λ as specified in Sec. VIII. If we use k T ordering, the shower PDFs are very close the MS PDFs. If we use Λ or ϑ ordering, the shower PDFs can differ substantially from the MS PDFs. This change is part of the summation of threshold logarithms [10] that is part of the shower algorithm represented in Eq. (136). (We present a simple example in the following section and an analysis of the structure of the summation of threshold logarithms according to the shower algorithm in Appendix E.) The next operator in Eq. (136) is U V (3, 0), written as three separate factors. This operator is also part of the summation of threshold logarithms. It uses shower oriented PDFs calculated using evolution in the parameter λ as specified in Sec. VIII. In general, U V (3, 0) leaves the number of partons and their momenta and flavors unchanged. The operator U V (t, t ) is the ordered exponential of S V (t), as given in Eq. (80). The generator S V (t) at lowest order is given by the first order version of Eq. (81)   (120): Consider the three operators U V (1, 0), U V (2, 1), and U V (3, 2) corresponding to the three segments of this path. In each case, U V is determined by d V [1] (t)/dt on that path segment, according to Eq. (137).
On the second segment, 1 < t < 2, µ e is fixed at µ e = 0 while µ 2 c varies from µ 2 h to 0. Thus P a,soft aa does not contribute. Only P a,LC+ aa contributes. That is, U V (2, 1), is U V (2, 1) as it would be calculated using the LC+ approximation. The LC+ approximation has the property that the color basis vectors used in the shower algorithm are eigenvectors of U V (2, 1). Thus U V (2, 1) applied to a parton basis state {p, f, c, c } m gives just a numerical factor.
On the first segment, 0 < t < 1, µ c is fixed at µ 2 c = µ 2 h while µ 2 e varies from µ 2 e (0) to 0. As we saw in Eq. (132), the unresolved region that applies in P a,LC+ aa is independent of t on the first path segment, so that P a,LC+ aa does not contribute. Thus we are left with P a,soft aa , which provides a correction in U V (1, 0) to the LC+ approximation used in U V (2, 1). The operator V leaves the number of partons and their momenta and flavors unchanged, but it can change the parton color state. The operator P a,soft aa , in fact, contains nontrivial color operators. Thus U V (1, 0) is, in general, a nontrivial operator on the color space for the partons involved in the hard process, as described in |ρ h ). Fortunately from the point of view of computation, this color space is finite dimensional. Thus it should be possible to calculate U V (1, 0)|ρ h ) numerically. Even more fortunately, in the color space for the qq final state in the lowest order Drell-Yan process, the LC+ approximation is exact, so that P a,soft aa acting in this space vanishes. That is, for the lowest order Drell-Yan process, Let us now consider the three operators U(1, 0), U(2, 1), and U(3, 2) corresponding to the probability conserving shower evolution on the three segments of the path.
On the first segment, 0 < t < 1, because of Eq. (132), D a,LC+ aa does not contribute to U(1, 0). Only D a,soft aa contributes. This operator contains non-trivial color operator and it creates more partons each time it acts. That is, the dimensionality of the color space grows each time S [1,0] acts. This means that, at least as far as we know, one cannot, in general, calculate the action of U(1, 0) exactly. However, we expect that one can expand U(1, 0) in powers of D a,soft aa up to order [D a,soft aa ] n and calculate these contributions numerically for any n that we choose, limited by the available computer power available. It remains a future research project to implement this scheme. 7 There is one case in which U(1, 0) can be computed exactly. In the Drell-Yan process at lowest order, one starts with qq color states. These color states form a one-dimensional space and on this space the LC+ approximation is exact. Thus D a,soft aa acting on the single state in this space vanishes. That is, U(1, 0) = 1 when acting on this space. We can simply ignore U(1, 0).
On the second segment, 1 < t < 2, U(2, 1) creates a full shower using first order splitting functions derived from D a,LC+ aa . The soft splitting functions from D a,soft aa do not contribute because µ 2 e is fixed. This gives us a first order shower based on the LC+ approximation from µ 2 h to zero (with infrared cutoffs based on m 2 ⊥ (a,â)). On the third segment, 2 < t < 3, neither µ are interleaved with ordinary LC+ splittings. We found that the contributions from D a,soft aa are practical to compute and can be substantial in the case of the rapidity gap fraction in hadron-hadron collisions but are small in other cases that we have examined. commutes with S iπ (t), so that However, 1 S iπ (t) = 0 so that 1 U(3, 2) = 1 . This leaves That is, the effect of U(3, 2) cancels and there is no need to calculate U(3, 2). There is a possible exception to this. One may want to add a nonperturbative model of hadronization to the perturbative parton shower. For instance, one can start with a perturbative Deductor event, add a nonperturbative underlying event, and then pass the event to Pythia for hadronization according to the Lund string model, as in Ref. [10]. In this case, the result depends on the color configuration in the Deductor event, so that the effect of U(3, 2) does not automatically cancel.
We can summarize these results for the three segment path. First, for the third segment we can use U V (3, 2) = 1 and, as long as the measurement operator O J is blind to parton color, U(3, 2) → 1. For the second segment, we need both U V (2, 1) and U(2, 1), but these operators are straightforward to calculate with the current Deductor code because they use the LC+ approximation. For the first segment, if we consider the lowest order Drell-Yan process then U V (1, 0) → 1 and U(1, 0) → 1, leaving Alternatively, if we consider parton-parton scattering as the hard process, then U V (1, 0) should be fairly simple to calculate by exponentiating a finite dimensional matrix. We expect that one can calculate shower evolution operator U(1, 0) perturbatively in powers of D a,soft aa . We leave implementation of the calculation of U V (1, 0) and U(1, 0) to future research.

XII. THE DRELL-YAN CROSS SECTION
We can test how the methods of this paper work by calculating the cross section for p + p → µ +μ + X. We use the three segment evolution path presented in this paper. We start with the expression (141) for the cross section. Here the starting statistical state |ρ h ) is obtained from the Born matrix elements. It would be preferable to correct this to NLO, but the required matching is not available in the present version of Deductor and, with matching, the treatment of U V (1, 0) and U(1, 0) would be nontrivial.
We first look at the cross section dσ/dM µμ , in which only the mass M µμ of the muon pair is measured. In this case the measurement operator O J commutes with the shower operator U(2, 1) and we can use 1 U(2, 1) = 1 to eliminate the shower operator, leaving There are three versions of each of the operators U V (2, 1) and F(µ 2 h ), one version for each of k T ordering, Λ ordering, and ϑ ordering.
In Fig. 11, we show the separate effect of the operators U V (2, 1) and F(µ 2 h ). First, we replace U V (2, 1) by 1 and use the k T -ordered, Λ-ordered, and ϑ-ordered versions of the shower oriented PDFs in the operator F(µ 2 h ). These PDFs are determined from the standard CT14 [19] MS PDFs using the λ-evolution equation (98). We find the cross sections shown as solid curves in Fig. 11. Evidently, the choice of ordering variable makes a substantial difference in the results. Second, we fix the PDFs in F(µ 2 h ) to be the MS PDFs and use the full operators U V (2, 1) corresponding to the three choices of ordering variable. We find the cross sections shown as dashed curves in Fig. 11. Again, the choice of ordering variable makes a substantial difference in the results.
In Fig. 12, we calculate the cross section using the product U V (2, 1) F(µ 2 h ) of operators for each of k T ordering, Λ ordering, and ϑ ordering. We find the cross sections shown as solid curves in Fig. 12. Now we see that the choice of ordering variable makes almost no difference in the results. We also see that the effect of sum-  Fig. 3 of Ref. [9]. The curve labeled analytic (BNX) is the analytic summation of threshold logarithms of Becher, Neubert, and Xu, comparable to the NNLO curve of Fig. 8 of Ref. [20]. The curve labeled NLO is obtained from a perturbative calculation using MCFM [21]. ming threshold logarithms using U V (2, 1) F(µ 2 h ) is quite substantial. This is not a surprise because threshold logarithms are known to be quite important for processes that involve partons with large momentum fractions. In Fig. 12, we show as a dashed red curve the result of the analytic summation of threshold logarithms by Becher, Neubert, and Xu [20]. We see that the shower version of the summation of threshold logarithms matches the analytic version quite nicely except that the shower result needs an approximately constant correction factor of about 1.15 to equal the analytic result. This factor would plausibly be supplied by matching the shower result to a perturbative NLO calculation. In Fig. 12, we also show the perturbative NLO result obtained from MCFM [21]. This fixed order result lacks the summation of threshold logarithms, but it does show the approximately 15% correction to the shower result at M µμ ≈ 1 TeV.
In Appendix E, we use a simple leading logarithm approximation to analyze how the operators U V (2, 1) and F(µ 2 h ) combine to sum threshold logarithms. This analysis shows why the curves in Fig. 11 have the qualitative behavior that we see in the figure. This analysis also shows why, at least in the leading approximation, the three curves in Fig. 12 match.
We now look at the cross section dσ/(dM µμ dp T ), in which both the mass and the transverse momentum p T of the muon pair are measured. Now the shower opera- (1/σ) dσ/dpT Drell-Yan pT distribution kT ordered Λ ordered ϑ ordered ResBos

FIG. 13.
The normalized Drell-Yan transverse momentum distribution, (1/σ)dσ/dpT for the LHC at 13 TeV. Here pT is the transverse momentum of the µμ pair, dσ/dpT is dσ/(dMµμdpT) integrated over 2 TeV < Mµμ < 2.1 TeV and σ is this cross section integrated over 0 < pT < 100 GeV, so that the area under the curve is 1. We show curves for kT ordering, Λ ordering, and ϑ ordering of the shower. We also show the corresponding result obtained using ResBos [23,24] as a dashed black curve.
tor U(2, 1) in Eq. (142) does not commute with the measurement operator O J , so that the parton shower has an effect. In Fig. 13, we show dσ/(dM µμ dp T ) integrated in M µμ over the range 2 TeV < M µμ < 2.1 TeV and normalized by this cross section integrated in p T over 0 < p T < 100 GeV. The shower approximately sums logarithms of p T /M µμ . We compare to an analytic summation of these logarithms [22] using the program ResBos [23,24]. The ResBos result includes some nonperturbative broadening that is not present in the shower result. Fig. 13, is comparable to Fig. 4 of Ref. [9], but now the definitions for the shower algorithm are a little different and now we show results for k T ordering, Λ ordering, and ϑ ordering of the shower.
We see that the distributions for k T ordering and Λ ordering agree with each other pretty closely. The Res-Bos result is, as expected, a bit broader than these. The result with ϑ ordering is qualitatively similar to the k T and Λ ordering results, but noticeably narrower.

XIII. SUMMARY AND COMMENTS
In a previous paper [1], we showed how, in the case of electron-positron annihilation, one can formulate a first order parton shower using more than one shower scale, µ = (µ 1 , µ 2 , . . . ). Then one has the freedom to specify an evolution path µ(t). This can lead to a simpler end result for cross sections than one would get with, say, k T evolution. In this paper, we extend the use of multiscale evolution to showers with initial state partons. This paper covers many issues, so we briefly summarize the main points here.
A parton shower with initial state partons in inherently more complicated than a shower with just final state partons. One needs to factor parton distribution functions out of the cross section and one needs to connect the evolution of the shower with the evolution of the parton distribution functions. This has two consequences.
First, the parton distribution functions used internally in the shower should not be the same as the MS parton distribution functions used for calculations in fixed order perturbation theory. We have analyzed this in earlier papers [2,[8][9][10], but in Sec. VIII of this paper, we have presented a more precise derivation. We also present a direct relation between the shower oriented PDFs and the MS PDFs that is better suited than the relationship used in Refs. [2,[8][9][10] when the scale arguments of the PDFs are large.
A second consequence of the connection between shower evolution and PDF evolution is that, in addition to the shower evolution operator U, a second operator, U V appears in the shower cross section. The operator U V creates a summation of threshold logarithms. This operator was included in Refs. [2,[8][9][10], but in this paper in Sec. VII we have presented what we believe is a clearer derivation. We also present an improved expression for U V and the infrared sensitive operator D on which the derivation is based. In particular, we have simplified the definition of the part of D at first order that comes from virtual parton exchanges. This is achieved by introducing what we call the momentum sum rule in Sec. V F.
With initial state partons included, the first order virtual exchange operator D [0,1] acquires an imaginary part. Sec. V G of this paper provides a simpler definition of Im D [0,1] than appears in Refs. [4,11].
Using this structure for initial state showers, we consider a shower with multiple shower scales. We use three shower scales, defined in Sec. IX. The first is µ c , which stands for one of µ ⊥ , µ Λ , or µ ∠ , corresponding to k T ordering, Λ ordering (based on virtuality), or angular ordering. The scale µ c determines the boundary of a region in the space of splitting variables (z, ϑ) that is defined to represent splittings that are unresolved at scale µ c . With k T ordering or Λ ordering, a splitting is unresolved if it is either too soft ((1 − z) too small) or too collinear (ϑ too small.) The second scale is µ e , for which a splitting is unresolved if it is too soft, (1 − z) < µ 2 e /Q 2 . There are some terms in the shower splitting functions that have a soft singularity but not a collinear singularity. We take advantage of having the scale µ e by eliminating the cuts based on µ c for these terms. We introduce a third scale µ iπ that defines an unresolved region for Im D [0,1] but not for other terms in D [1] .
With three scales available, we define a simple, threecomponent path in the space of scales that takes us from large scales to small scales. This choice gives a simpler result for a calculated cross section than one would have with a single scale µ c . First, as long as the measurement applied at the end of the shower is not sensitive to the color state of the partons, the imaginary part of the splitting function has no effect. Second, one can choose to use the LC+ approximation for color on the second segment of the path. Then there is a term in the splitting functions that is proportional to the difference between full color and the LC+ approximation. This term has a soft singularity by no collinear singularity. This difference term appears by itself on the first segment of the path. In the case of the lowest-order Drell-Yan process, the shower operator, U(1, 0), is just the unit operator and the threshold operator U V (1, 0) is also just the unit operator.
This gives us a very simple result for the cross section. One can simply use the result for one scale µ c in which we ignore the imaginary part of virtual exchanges and in which we use the LC+ approximation for color. In Sec. XII we exhibit results using this approach for the Drell-Yan cross section dσ/dM µμ and for the p T distribution of the produced muon pair. These results are obtained for k T ordering, Λ ordering, and ϑ ordering. The results for dσ/dM µμ are remarkably independent of the ordering choice and the results for the p T distribution show a noticeable dependence on the ordering choice in the case of ϑ ordering.
Having summarized the main points of the paper, we can offer some comments.
We view a parton shower as something that is, in principle, defined at all perturbative orders [2]. In particular, the infrared sensitive operator D should be thought of as having an expansion in powers of α s : D = 1 + D [1] + · · · . Then, in the formulation of Ref. [2], the generator S(t) of parton splittings is determined and has an expansion in powers of α s beginning at order α 1 s . In a first order shower, we use the first order term. In fact, it is an important open problem to define a prescription in which the second order contributions to D and thus S(t) can be calculated.
With a first order shower, there is considerable freedom to define D [1] and thus S [1] (t). For instance, one can use a single scale that defines either k T ordering or angular ordering. The resulting first order splitting functions are not the same. However, as analyzed in Ref. [1], the contribution to a cross section resulting from this difference starts at order α 2 s . If we were able to work with order α 2 s splitting functions, S [2] , then the differences in calculated cross sections between prescriptions would be order α 3 s and therefore would be less significant. In addition, if we were able to use splitting functions S [2] , we would have a sensible way to judge the merits of two prescriptions for a first order shower: we would prefer a prescription at first order that makes the second order splitting functions S [2] numerically small in applications.
In this paper, we have used the option to use multiple shower scales µ s in a way that makes the first order shower simpler. The difference between the resulting shower and a more standard one scale shower would then be generated by higher order contributions to the shower splitting functions. We can hope, but cannot now prove, that the first order shower thus defined will not have large contributions to its splitting functions when extended to second order.
We expect the D [2] operator in a second order shower to be quite complicated. In order to deal with the complications, it may be helpful to use multiple splitting scales along the lines used in this paper.

ACKNOWLEDGMENTS
This work was supported in part by the United States Department of Energy under grant DE-SC0011640. The work benefited from access to the University of Oregon high performance computer cluster, Talapas. We thank Thomas Becher for advice about threshold summation for the Drell-Yan cross section and Frank Petriello for advice about the perturbative Drell-Yan cross section. the splitting. There are dimensionally regulated singularities, so this integration is in 4 − 2 dimensions for each momentum. ThenD a aâ is a function of the momenta and flavors before and after the splitting and is an operator that maps the m-parton color space to the m+1 parton color space. The index a inD a aâ is fixed. It is the flavor of the incoming parton "a." The operatorD a aâ also depends on the two scales µ 2 r and µ 2 c , but we have not made that dependence explicit in the notation.
The operator F = F a F b provides the parton distribution functions needed to make a cross section, Eq. (6.3) of Ref. [9], Here n c (a) is the number of colors for flavor a and n s (a) is the number of spins, 2 for a quark and 2(1 − ) for a gluon. The operatorD â aa iŝ The operatorD a aâ depends on the splitting variable z and onâ, which is the flavor of the incoming parton after the splitting (in the sense of backward evolution). We useâ as a splitting variable that specifies the flavor content of the splitting. For instance,â = a corresponds to gluon emission from the incoming line, so that parton m+1 is a gluon. If a is a quark flavor and a = g, then parton m+1 is the corresponding flavor of antiquark. The right-hand side of Eq. (A12) begins with integrations over the other two splitting variables ϑ and φ, with appropriate dependence on the dimensional regulation parameter . The variable φ is a unit vector in the 2 − 2 dimensional transverse momentum space and represents the azimuthal angle ofp m+1 . There is a sum over the index k ∈ {a, b, 1, . . . , m} of a dipole partner parton for a splitting. Then φ k is the azimuthal angle of p k if k / ∈ {a, b}. The integration over φ is an integration over a unit sphere that is a 1 − 2 dimensional surface. The function S(2 − 2 ) is the surface area of this sphere, so that Integration over φ of a function of φ − φ k is a shorthand notation for integration over a unit vector φ in a coordinate system with the unit vector φ k equal to (1, 0, . . . ). After the integrations, there is a delta function that sets {p,f } m+1 to the momenta and flavors obtained from a splitting with variables (z, ϑ, φ,â) applied to partons with momenta and flavors {p, f } m according to Deductor conventions.
The idea of the singular operator D [1,0] is that it integrates over splittings that are arbitrarily close to the soft and collinear limits, but with an ultraviolet cutoff that depends on a scale parameter µ 2 c . The region of (z, ϑ) allowed by the cutoff is called the unresolved region and is denoted by U (µ 2 c ). The subscript c denotes the definition of the unresolved region, as described in Sec. V C. We therefore insert a theta function that specifies that (z, ϑ) lies in the unresolved region.
There is a special feature that applies when the hard scattering process is parton-parton scattering, a + b → 1 + 2, with transverse momentum P ⊥ defined by p 1 = α p a + β p b + P ⊥ . Here the partons are massless and P ⊥ · p a = P ⊥ · p b = 0. In a parton shower, the initial state partons split in backwards evolution into other final state partons besides partons 1 and 2 from the hard scattering, leaving at the end of the shower new initial state partons a and b with larger momentum fractions than the partons that created the hard scattering. In this situation, the scattering identified as the hard scattering should be the hardest of all the scatterings. That is, the k 2 T of each initial state emission should be no greater than |P 2 ⊥ |. 8 For this reason, the unresolved region U (µ 2 c ) in Eq. (A12) should be a subregion of 0 < z < 1, 0 < ϑ < 1, and k 2 T < |P 2 ⊥ |. That is, we define Θ((z, ϑ) ∈ U (µ 2 c )) in Eq. (A12) by Θ (z, ϑ) ∈ U (µ 2 c ) = θ(ϑ < a c (z, µ 2 c )) θ(0 < ϑ < 1) × θ(ϑ < a ⊥ (z, |P 2 ⊥ |)) .

(A14)
Then the resolved region consists of all points (z, ϑ) with 0 < z < 1, 0 < ϑ < 1, and k 2 T < |P 2 ⊥ | that are not in the unresolved region. Only splittings with k 2 T < |P 2 ⊥ | can then be generated in the shower. When the hard process is Drell-Yan muon pair production, this restriction does not apply. When the hard scattering is parton-parton scattering, if the shower is k T ordered, this restriction does not matter, but in the case of Λ or ϑ ordering, this restriction does matter. In order to avoid extra clutter, we omit effects of the factor θ(ϑ < a ⊥ (z, |P 2 ⊥ |)) in our formulas that follow.
In the next factor in Eq. (A12), there is a sum over dipole partner partons k. In the first term, the partner parton is the same as the emitting parton, k = a. This term contains a color factor N (a,â) defined by where q is any quark or antiquark flavor. Then there is a splitting functionP aâ (z, ϑ, ). For the case thatâ = a, this function is related to the function w aa ({p,f } m+1 ) that appears in Eq. (5.7) of Ref. [4] by Here n s (a) is the number of spin states for a parton of flavor a as in Eq. (A11). We calculate w aa ({p,f } m+1 ) and similar functions in the cases below from the definitions in Ref. [3] so that the sums over parton spins are performed in 4 − 2 dimensions. For the case thatâ = a, P aa (z, ϑ, ) is related to the functions w aa ({p,f } m+1 ) and w eikonal aa ({p,f } m+1 ) that appear in Eq. (5.7) of Ref. [4] by The functionsP aâ (z, ϑ, ) are rather simple. We simplify the notation by defininĝ P aâ (z, ϑ, 0) =P aâ (z, ϑ) , P aâ (z, 0, 0) =P aâ (z) , Then we findP where y = (1 − z)ϑ/z and q is any quark or antiquark flavor.
The contributions proportional to (at y = 0) arê Finally, at y = 0 and = 0 we havê That is, these are the familiar unregulated DGLAP kernels. Next, in Eq. (A12) is a term proportional to a function W 0 (ξ ak , z, ϑ, φ − φ k ). This term comes from interference between emission of a gluon from parton "a" and emission from dipole partner parton k with k = a. The variable ξ ak is (1−cos θ a,k )/2 where θ a,k is the angle between p k and p a as measured in the rest frame of Q, The function W 0 is related to the functions w dipole ak ({p} m+1 ) and A ak ({p} m+1 ) that appear in Eq. (5.7) of Ref. [4] by The functions w dipole ak ({p} m+1 ) and A ak ({p} m+1 ) are (from [4], Eq. (5.3) and [6], Eq. (7.12)) The function (1 − z) A ak w dipole ak is singular both in the collinear limit, ϑ → 0 at fixed (1 − z), and in the soft limit, (1 − z) → 0 at fixed ϑ. In the collinear limit, In Eq. (A23), we have subtracted this singular behavior, leaving a function W 0 that vanishes in the collinear limit, although it remains finite in the soft limit.
The explicit expression for W 0 is In our study of the inclusive splitting probability, we will need the azimuthal average of W 0 at = 0, It has a reasonably simple form: where We can rewrite this in a form that shows that W (ξ ak , z, ϑ) is proportional to ϑ 2 as ϑ → 0 with fixed z: We will also need W at z = 1: Finally in Eq. (A12) there is a factor with color operators. The operator t † a (f a →f a +f m+1 ), acting on the ket color state {c} m , gives the new color state {ĉ} m+1 that one gets after emitting the new parton m+1 from parton "a" with flavor f a . This operator is described in some detail in Ref. [3]. Similarly, t k (f k →f k +f m+1 ), acting on the bra color state {c } m , gives the new color state {ĉ } m+1 that one gets after emitting the new parton m+1 from parton k with flavor f k .
When parton m+1 is a gluon, the color operators obey the identity This identity arises from the fact that the parton color state is an overall color singlet, so that attaching a color generator matrix T c k to all of the parton lines k in the state, including k = a, gives zero. We have used this identity to add the same term, proportional to z/[(1 − z)ϑ], to both the k = a term and the k = a terms in Eq. (A12). We have added this term in both places in order to move the soft×collinear singularity from the k = a terms to the k = a term. After this change, the k = a terms, proportional to W 0 , have a soft singularity but not a collinear singularity.
We now turn to the splitting operator. The shower operator U(t 2 , t 1 ) is the ordered exponential of the integral over shower time t of a splitting operator S(t) according to Eq. (8). Consider the first order contribution to the splitting operator for splitting of initial state parton "a," S with respect to the shower time using Eq. (79) specialized to first order, which gives Eq. (10). For the contribution from initial state parton "a," this is We differentiate D [1,0] a as given in Eq. (A10). In Eq. (A10), the ratio of PDFs corresponds to the PDF operators F a (µ 2 r ) · · · F −1 a (µ 2 r ) rather than D [1,0] a , so we do not differentiate the PDF factor. There is a factor of α s (µ 2 r ) that is part of D [1,0] a , but the derivative of α s with respect to its scale argument is of order α 2 s and we want only the first order contribution to S a (t), so we do not differentiate the α s factor. This leaves the derivative of D a aâ , which is given in Eq. (A12). The dependence on µ s (t) is contained in the factor Θ((z, ϑ) ∈ U (µ 2 c )), so it is only the theta functions contained in this factor that we should differentiate, producing delta functions that fix the integration variable ϑ inD a aâ . Then the dimensional regularization is not needed, so we can set = 0. These considerations give one a straightforward calculation of S [1,0] (t). We omit further details.
The virtual exchange part, S (t) can be determined rather simply because the definition of S(t) is arranged so that the shower operator U(t 2 , t 1 ) preserves probability. This implies that 1 S Then the dimensional regularization is no longer needed, so we set = 0. Finally, the color operator t † a ⊗ t k for a real emission is mapped into the operator 1 ⊗ t k t † a for the corresponding virtual exchange and t † k ⊗ t a is mapped to t a t † k ⊗ 1. The result is This is the operator that appears in Eq. (9).

Appendix B: Pole structure
Inclusive splitting at first order is given in Eq. (54) by an operator P a aa (z; {p} m ) that operates on the parton color space. This operator has two parts. The first, P a aa (z; {p} m ) comes from real emission graphs, while the second, Γ a ({p} m ) comes from virtual exchange graphs: In P a aa (z; {p} m ), there are infrared singularities in the integration over splitting variables y, z, φ. The singularities are regulated by working in 4 − 2 dimensions, so that integrals over z ofP a aa (z; {p} m ) have poles 1/ and 1/ 2 . The integrals over y, z, φ are confined to the region of unresolved splittings set by a scale or scales µ s .
For the virtual graphs that contribute to Γ a ({p} m ), there is some freedom available in defining a region of unresolved loop momenta that corresponds to the unresolved region for real emission graphs. However, the pole structure of Γ a ({p} m ) is independent of the definition of an unresolved region and is given by We then show that this leads to Eq. (B3). We call Eq. (B4) the momentum-sum-rule (MSR) ansatz. Of course, Eq. (B4) could be supplemented by adding to Γ a ({p} m ) any operator that has no poles. In Eq. (B1), the inclusive probabilityP a aa for real emissions is related to D [1,0] in Eq. (44). Its structure is given in Eq. (50). We can rewrite this equation by changing the integration variable from y to ϑ = zy/(1 − z). We write the result in terms of three functions, A k,a (z), and C (R) k,a (z; ξ ak ), where the superscript R denotes real emissions, Here whereP aa (z) is the familiar DGLAP kernel, equal toP aa (z, ϑ, ) at ϑ = = 0. Next, Finally, For the virtual exchange operator Γ a ({p} m ), we can define functions A a , and C (V) k (ξ ak ), with superscripts V for virtual exchanges: Here a virtual gluon exchange between the initial state parton "a" and another parton k is represented by C k,a (ξ ak ). It has the same color structure as the C (R) k,a (z; ξ ak ) term inP a aa . We start our investigation with the third term in Eq. (B5) forP a aa , proportional to the function C (R) k,a (z; ξ ak ). This term has a distinctive color structure, {[T k · T a ⊗ 1] + [1 ⊗ T k · T a ]}. This term inP a aa does not have a collinear, ϑ → 0, singularity because W 0 vanishes as ϑ → 0, as one can see directly from the expression (A25) for W 0 . However it has a soft gluon singularity, as seen in the factor 1/(1 − z). This singularity produces a pole 1/ after we integrate over z. For any smooth function h(z) we have The corresponding contribution to virtual exchange integrated against a smooth function h(z) is In order for the two 1/ poles in Eqs. (B10) and (B11) to cancel, we need In order to obtain this relation, we use the definition given by the momentum sum rule: k,a (z; ξ ak ) .
Of course, there are other ways to satisfy Eq. (B12).
k,a is expressed as a dimensionally regulated integral over a 4 − 2 dimensional loop momentum q. We need to cut out the ultraviolet contributions to this integral using the shower scales µ s used for the real emission integrations. However we do this, the infrared pole 1/ will match between the real and virtual graphs. We have seen that if we simply define the virtual integral in terms of the corresponding real integral by using the momentum sum rule, the 1/ contributions will cancel as required. Then the momentum sum rule ansatz fixes the finite part of the virtual contribution as a function of µ s .
Thus the integral of B aa (z) according to the momentum sum rule, it will contain no poles also. We conclude that B aa (z) is free of 1/ poles. Finally, we consider the first term in Eq. (B5), proportional to A (R) aa (z). We can perform the ϑ integration to give where Q (R) Caution is needed becauseP aa (z) has a 1/(1 − z) singularity:P whereP reg aa (z) is finite for z → 1. Treating Q (R) aa (z) as a distribution by integrating against a test function h(z) that vanishes for z → 0, we obtain Here we have used the dimensional regularization factor [(1 − z) 2 /z] − in the singular term. Then instead of an unregulated 1/(1−z) singularity, we obtain a distribution that is well defined as long as = 0, with δ(1 − z) and [z/(1−z)] + terms. Overall, counting the 1/ in Eq. (B15), we have 1/ 2 and 1/ infrared singularities in A We can now add the contributions from virtual graphs. We can use the momentum sum rule as an ansatz to define these contributions. Then we can check whether this ansatz works to remove all 1/ n poles in a physical cross section. The momentum-sum-rule ansatz is to replace Q with This gives Note that the contribution proportional to δ(1 − z) with a 1/ singularity has cancelled. Since P aa (z) obeys the momentum sum rule, we have We can now define A aa (z) = A (R) aa (z) with the result that Of course, we could also add a term proportional to δ(1 − z) with no poles to A (V) aa (z). The construction presented here determines only the pole structure of the contributions from virtual graphs to P a aa (z; {p} m ). We can use Eq. (B23) in Eq. (72). For V in Eq. (72), counting the factorization subtraction, we obtain The 1/ singularity cancels from V.
We conclude that as long as the contribution to P a aa (z; {p} m ) from virtual graphs is defined using the momentum sum rule, we obtain as required for the infrared finite operator V in Eq. (70) to be free of infrared poles in the dimensionally regulated theory. This is the result reported in Eq. (56).
Appendix C: Structure of V at first order In this appendix, we examine the structure of the operator V at first order in α s . We start with Eq. (70), The operators K, Z F , and D are simply unit operators at order α 0 s . Thus at first order, the part of V that applies to hadron A becomes The • symbols in the first two terms stand for convolutions, as in Eq. (63) of Ref. [2], so that We are free to choose what K aa (z, µ 2 r ) should be. In the second term, the factor in the MS scheme that removes the 1/ singularity from an initial state splitting is 9 The operators in Eqs. (C3) and (C4) are proportional to unit operators on the parton color space. In Eq. (54), is expressed as a convolution of the same PDF factor with an operator P a that acts nontrivially on the parton color space.
We write V a as a convolution with PDFs according to Eq. (71), This gives We now need to examine the structure of P a aa (z; {p, f } m ; ) for → 0. Since we define P a aa from its real emission partP a aa using the momentum sum rule, Eq. (59), we can begin withP a aa . For this, we can use Eq. (50): (C7) 9 This is in 4 − 2 dimensions. Recall that αs is dependent and has dimension (mass) −2 in the dimensionally regulated theory [26].

(C10)
Here, in general, µ lim is a function µ lim ( µ s ) of the scales µ s . Then a lim (z, µ 2 lim ( µ s (t)) is the limiting form of a cut (z, µ s (t)) at large shower times t, as defined in Eq. (26).
In the term in Eq. (C7) proportional toP aa , we use the notation of Eq. (A18) to writeP aa (z, ϑ, ) in the form P aa (z, ϑ, ) =P aa (z) − P ( ) aa (z) + [P aa (z, ϑ, ) −P aa (z, 0, )] For the termP aa (z, ϑ, ) −P aa (z, 0, ), it suffices to take the → 0 limit inside the integrations. For the first two terms, we need the integral Performing the integrations gives After some rearrangement, this gives us an infrared sensitive contribution and a non-sensitive (NS) contribution: The infrared nonsensitive contribution iŝ The first term here isP aa (z) times [− log(a limit ) + log(a max )θ(a max < 1)]. For our later purposes, we have written this as an integral. 10 Eq. (C16) givesP a aa , the part of P a aa associated with real parton splittings. To obtain the full P a aa , we use the momentum sum rule, Eq. (59): where Here we have noted that the DGLAP kernel obeys the momentum sum rule: [P aâ (z)] msr = P aâ (z). We note from Eq. (A29) that W (ξ ak , z, ϑ) is proportional to ϑ 2 for ϑ → 0, so that the integration of W over ϑ is not singular. The infrared nonsensitive contribution has no 1/ poles and vanishes in the limit of small scales µ s 10 The form used in Eq. (C17) assumes that a limit ≤ amax, as is the case when we use Eq. (C12). A slightly different form would apply if a limit > amax.
when P a,NS aa (z; {p} m ) msr is integrated against a fixed test function f (z).
We can now insert these results into Eq. (C6), giving The contribution proportional to P aa (z)/ in Eq. (C6) has cancelled an identical singular term in P a aa (z; {p} m ; ) in Eq. (C18). This is the result reported in Eq. (74) in the main text. We now define K aa (z, µ 2 r ) so that it cancels the second and third terms in Eq. (C20), leaving This is the result reported in Eq. (76) in the main text.
Appendix D: Structure of SV at first order We use the operator V(µ r (t), µ s (t)) to define the threshold operator U V (t 2 , t 1 ) according to Eqs. (80) and (81): We let t = 0 be the shower time at the start of the shower, corresponding to a renormalization scale parameter µ r = µ h , where µ h is comparable to the scale of the hard interaction. Then t f is the shower time at the end of the shower. Thus we consider At first order, There are two contributions to S If we use the three segment path from t = 0 to t = t f = 3, then the T instruction tells us to break U V (3, 0) into U V (3, 2) U V (2, 1) U V (1, 0). We saw in Sec. XI that U V (3, 2) = 1. Next, U V (2, 1) is calculated using the LC+ approximation. No T instruction is needed within U V (2, 1) because all of the operators commute. Finally U V (1, 0) is determined by the difference between the full color operators and their LC+ approximate version. This contribution is subleading in threshold logarithms and subleading in color. It is not simple to treat U V (1, 0) in computer code. We leave that to future work. For the leading order Drell-Yan process investigated in this paper, the parton state just after the hard interaction is very simple and this color difference operator applied to this state gives zero, so that U V (1, 0) = 1.
In this appendix, we examine V is not sensitive to the infrared cutoff used. For that reason, we simplify the notation by setting m 2 ⊥ (a,â) = 0 for now. Then a max = a limit = a c (z, µ 2 c ). We set µ 2 r = µ 2 c . This gives for V On the second segment of our chosen path, we use the LC+ approximation for the color operators in Eq. (D6), but we do not specify that choice here. We need the derivative of this with respect to t: When we differentiate with respect to µ 2 c , we differentiate with respect to the upper endpoint of the ϑ integral, but we do not differentiate α s (µ 2 c ), f a /A (η a /z, µ 2 c ), or f a/A (η a , µ 2 c ) because the derivatives of these objects are of order α s and we want only the first order contribution to S V,a (t) Appendix E: Summation of threshold logs In the shower cross section repeated in Eq. (82) from Ref. [2], the first two operator factors applied to the hard scattering statistical state ρ h are F(µ 2 h ) and U V (t f , t h ). (For application to the second segment of our three segment path, the U V factor would be U V (2, 1).) The operator U V (t f , t h ) leaves the number of partons and their momenta and flavors unchanged but multiplies by certain factors that relate to how the shower splitting functions match to the evolution of parton distribution functions [8]. The operator F(µ 2 h ) multiplies by some standard normalization factors and by parton distribution functions f a/A (η a , µ 2 h ) and f b/B (η b , µ 2 h ). These are shower oriented parton distribution functions, which can be different from the normal MS PDFs.
When the cross section to be computed involves a very hard scattering, so that η a and η b are close to 1, the operator U V (t f , t h ) applied to ρ h can be large. Additionally, the ratios f a/A (η a , µ 2 h )/f ms a/A (η a , µ 2 h ) and can be large. In this case, there are big effects because Eq. (82) is summing what are called threshold logarithms. The summation of threshold logarithms was first analyzed by Sterman in 1987 [27]. There is a large literature on this summation using analytical calculations. 11 Ref. [9] cites some of this literature. In an analytical approach, one normally analyzes a suitable Mellin transform of the cross section. Then the transformed cross section is a function of a Mellin moment variable N and we are interested in the limit 11 There has been substantial recent interest in the summation of large logarithms by parton shower algorithms [7,[28][29][30][31][32]. However, to our knowledge the summation of threshold logarithms is not included in parton shower algorithms except for Deductor [2,9,10].
N → ∞. In this appendix, we apply the Mellin transform approach to the operators that appear in one of the simple single scale versions of the shower cross section introduced in Sec. V C. Our aim is to compare the parton shower version of threshold summation to the analytical results, extending the less general analysis in Ref. [9]. We find that U V (t f , t h ) becomes an exponential where E V,a (N ) and E V,b (N ) can be expanded in a series of powers of α s (µ 2 h ) and powers of log(N ). Similarly, the PDF factors become exponentials f a/A (η a , µ 2 h ) f ms a/A (η a , µ 2 h ) → exp(E pdf,a (N )) , → exp(E pdf,b (N )) .

(E2)
The leading approximation obtained with analytical approaches takes the form threshold factor = exp(E a (N ) + E b (N )) .
With the shower approach in this appendix, we find We begin with the operator U V (t f , 0). Using Eqs. (D3) and (D5), we write with contributions from the two initial state parton legs, "a" and "b." For initial parton leg "a," we have where S V,a (t) is given by Eq. (D14). In Eq. (D14), the unresolved region is defined using scale parameters µ 2 c , where C =⊥ for k T ordering (Eq. (19)), C = Λ for Λ ordering (Eq. (37)), and C = ∠ for ϑ ordering (Eq. (39)).
For the purpose of summing the leading threshold logarithms, the only term that matters in S V,a , in Eq. (D14), is the first, proportional to θ(ϑ(z, t) > 1)P aa (z). Furthermore, it suffices to approximateP aa (z) by its most singular term,P aa (z) → δ aa 2C a 1 − z . (E7) We choose the argument of α s to be k 2 T = (1 − z) 2 ϑQ 2 /z as in Eq. (D10), but we can omit the α 2 s term in Eq. (D12). Additionally, we can set the infrared cutoffs m 2 ⊥ (a,â) and m 2 ⊥ to zero with the understanding that we work term by term in the expansion of α s (k 2 T ) in powers of α s (µ 2 h ), where µ 2 h is the scale of the hard scattering.
We use k 2 T as the integration variable instead of t. The theta function in the first term in Eq. (D14) limits the integration to 1 < ϑ < a c (z, µ 2 h ) or We identify the square of the momentum of the partons just after the hard scattering, Q 2 , with µ 2 h . For the leading behavior of the integral, we are interested only in the (1 − z) → 0 behavior of the integrand. Thus we can replace the factor 1/z in Eq. (E8) and the factor z in Eq. (37) for a Λ (z, µ 2 ) by just 1. We can also replace the factor r a = η a /η In order to compare to standard analytical results that work with a Mellin transform of the cross section, we need to relate this to the Mellin transformed parton distributions. We can use the "single power approximation" [] in which we replace Now we can look at the PDF factor, Eq. (E2). We use the Mellin transform of the parton distribution, at the scale µ 2 h of the hard scattering that starts the shower, as a function of the parameter λ in Sec. VIII, g a/A (N, µ 2 h ; λ) = 1 0 dη η η N f a/A (η, µ 2 h ; λ) .
We need the Mellin transform of the shower oriented PDF, f a/A (η, µ 2 ; 1) with λ = 1. The shower oriented PDF with λ = 0 corresponds to k T ordering and differs from the MS PDF only by a small adjustment specified in Eq. (84). We obtain f a/A (η, µ 2 h ; 1) from f a/A (η, µ 2 h ; 0) by solving the differential equation (98) with the first order kernel d K/dλ given by Eq. (99) for Λ ordering or Eq. (111) for ϑ ordering.
We are interested only in the leading threshold logarithms in g a/A (N, µ 2 ; 1), so we make some small adjustments in the differential equation, as we did in the previous analysis of U V . We set the fixed infrared cutoff to m 2 ⊥ (a, a ) = 0. We replace the DGLAP kernel in the differential equation by its leading singularity as z → 1, Since only the (1 − z) → 0 limit affects the threshold logarithms, we replace z A (1 − z) B by just (1 − z) B in the differential equation. For Λ ordering, there is a factor r a in the scale choice that relates f a/A (η, µ 2 h ; λ) to an adjusted functionf . Since r a is just 1 in the hard scattering state for which the PDF is measured, we set r a = 1. The argument of α s in the first order kernel is not determined by strictly lowest order considerations. In order to incorporate leading logarithms beyond leading order in α s (µ 2 h ), we follow the choice generally made in the Deductor code and take the argument of α s to be Here we understand that α s (k 2 T ) is to be expanded in powers of α s (µ 2 h ). With these adjustments, the differential equation for g a/A (N, The solution of this is We can change integration variables from λ to k 2 T specified in Eq. (E17): .

(E22)
The sum of the U V and PDF contributions is .
Notice that the net result is the same for k T ordering (n = 0), Λ ordering (n = 1), and ϑ ordering (n = 2). For k T ordering, the entire result comes from E V,a (N ). For ϑ ordering, the entire result comes from E pdf,a (N ). For Λ ordering, both E V,a (N ) and E pdf,a (N ) contribute. The comparable standard formula given in Ref. [33] is .
Here we have z N −1 instead of z N −2 and 2C a instead of 2zC a , but these differences are not important in the N → ∞ limit.
Appendix F: Structure of operators with more scales In this appendix, we provide details about the opera-torsD a aa andP a aa needed in Sec. IX C when we use three scales, (µ 2 e , µ 2 c , µ 2 iπ ). Here µ 2 e specifies a cut on (1 − z) according to Eq. (124) and µ 2 c specifies the unresolved region for k T ordering for C =⊥ (Eq. (19)) or Λ ordering for C = Λ (Eq. (37)) or ϑQ 2 ordering for C = ∠ (Eq. (39)). We do not need to consider the scale µ 2 iπ sinceD a aa and P a aa do not depend on µ 2 iπ , which controls Im D [0,1] . We first need to specify in detail the two versions of the unresolved region defined in Sec. IX C: the region U (µ e , µ c ) that we use forD a,LC+ aa andP LC+ aa and the region U (soft; µ e ) that we use forD a,soft aa andP soft aa . With scales µ e and µ c , the unresolved region is defined in Eq. (122) and Eq. (123): (z, ϑ) ∈ U (µ e , µ c ) if 0 < z < 1, 0 < ϑ < 1 and ϑ < a max (z, µ 2 e , µ 2 c ) , where we define a max (z,µ 2 e , µ 2 c ) = max a cut (z, µ 2 e , µ 2 c ), a ⊥ (z, m 2 ⊥ (a,â) (F2) with a cut (z, µ 2 e , µ 2 c ) = max a e (z, µ 2 e ), a c (z, µ 2 c ) .
This gives us just the evolution for a k T -, Λ-, or ϑ-ordered shower according to our choice of C ∈ {⊥, Λ, ∠}.

Decomposition ofP a,NS
aa andD â aa We now record the details of the operatorsP a,NS aa andD â aa when we use two scales, µ e and µ c , that affect the definition of the unresolved region. We define a max (z, µ 2 e , µ 2 c ) and a limit (z, µ 2 e , µ 2 c ) using Eq. (C12). We use Eq. (C16) to obtainP a aa . In the third term on the right-hand-side of Eq. (C16), we use Eq. (C12). In the last term, we write 12 The operatorP a,LC+ aa (z; {p} m ) multiplies α s (µ 2 r ). In order to provide an improved summation of large logarithms, when