New insights on an old problem: resummation of the D-parameter

The D-parameter is one of the oldest and most experimentally well-studied hadronic observables for e+e− collisions. Nevertheless, unlike other classic observables like the C-parameter or thrust, the D-parameter has never been resummed throughout its entire singular phase space. Using insights and techniques motivated by modern multi-differential jet substructure calculations, we are able to predict the D-parameter distribution with no additional phase space cuts. Our approach is to measure both the C- and D-parameters on hadronic final states in e+e− collisions. We can tune the value of the C-parameter with respect to the D-parameter to specify simple, physical configurations of final state particles in which to perform calculations. There are three parametric regions that exist: D ≪C2 ∼1, D ≪C2 ≪1, and D ∼C2 ≪1, and we calculate the D-parameter in each region separately. In the first two of these three regions, we present all-orders factorization theorems and explicitly demonstrate resummation to next-to-leading logarithmic accuracy. The region in which D ∼ C2 ≪ 1 corresponds to the dijet limit and where the D-parameter loses the property of additivity. In this region we introduce a systematically-improvable procedure exploiting properties of conditional probabilities and resum to approximate next-to-leading logarithmic accuracy. The contributions from these regions can be consistently combined, and the value of the C-parameter integrated over to produce the cross section for the D-parameter. With these results, we match to leading fixed order as proof of principle and compare our resummed and matched prediction to data from LEP.


Introduction
The legacy of the event shape measurements from experiments at the Large Electron-Positron Collider (LEP) on precision QCD calculations is profound. Long after the experiments at LEP collected their last data, new frontiers of fixed-order and resummed perturbation theory have been established and used to study that data. Programs for numerical evaluation of e + e − → three-jet scattering to next-to-next-to-leading fixed order (NNLO) now exist [1][2][3], automated next-to-next-to-leading logarithmic order (NNLL) resummation for general observables is available [4,5], and resummation for particular observables has even been extended to next-to-next-to-next-to-leading logarithmic order (N 3 LL) [6][7][8][9][10]. These advances have enabled extractions of the strong coupling at the Z-pole, α s (m Z ), with resummation at NNLL or N 3 LL matched to NNLO [7,11,12] that are now included in world averages provided by the Particle Data Group [13]. All the observables that have been calculated at this high precision, such as thrust [14], C-parameter [15][16][17], and broadening [18][19][20], have the property that they are first nonzero at tree-level for e + e − → qqg (at O(α s )). This is vital, especially for resummation. In the limit where resummation for thrust τ becomes important, τ → 0, all radiation in the event is constrained to be collinear with the hard final state quarks or are soft (low energy). By momentum conservation, the hard quarks are back-to-back and the soft radiation is emitted off of the dipole formed from these quarks. This physical configuration is simple and it is not necessary to include hard gluons sourcing jet structure for precision resummation. Also, for these types of observables, including diagrams up through O(α 3 s ) corresponds to NNLO. Indeed, while it can calculate in principle arbitrary observables, the program EERAD3 by default only calculates observables which are first nonzero for e + e − → qqg [21].
In addition to these O(α s ) observables, the experiments at LEP measured a number of observables that are sensitive to more exotic particle configurations. In particular, because e + e − scattering occurs in the center-of-mass frame, several observables that are sensitive to the aplanarity of the final hadronic system have been defined and measured. These include observables like the thrust-minor [22] and the D-parameter [15][16][17]. Because these observables are only non-zero if the momentum of all final state particles is not coplanar, they are first non-zero at O(α 2 s ). As such, fixed-order calculations of these observables cannot yet reach NNLO precision, and their resummation requires analysis of many different phase space regions that can contribute. For example, a measurement of the D-parameter in JHEP02(2019)104 the resummation regime D 1 only constrains out-of-plane emissions. The configuration of the particles in the plane are completely unconstrained, as long as they remain in the plane. Because of this, theoretical studies of these observables had been restricted to fixedorder [23,24] or resummation in the region in which there are three planar, well-separated jets, with soft or collinear emissions out of the plane [25][26][27][28]. This can be accomplished by, for example, requiring that the three-jet resolution variable y 3 [29] is relatively large. However, such a procedure does not produce an inclusive prediction for these observables because the measured value of y 3 restricts radiation. Therefore, these predictions cannot be directly compared to the LEP measurements that had no such restriction.
In this paper, we present the first inclusive resummation calculations for these observables sensitive to aplanarity. Our analysis will focus exclusively on the calculation of the D-parameter, though the procedure should be applicable to a wide range of similar observables. This result is possible due to recent advances in resummation of multi-differential cross sections, which has seen extensive development in the field of jet substructure (see ref. [30] for a review). A fundamental problem in the field of jet substructure is the construction of optimal observables for discrimination between jets of different origins; for example, discrimination of QCD jets initiated by light partons from jets from hadronic decays of highly-boosted W , Z, or H bosons. One difference between these jets is their prong structure: jets from decays of electroweak bosons will dominantly have a two-prong structure, while QCD jets will dominantly have only one prong. There are many ways that a jet can manifest having two prongs, and so one must carefully identify the different phase space regions and then sum them up consistently. This is the approach we take here to the resummation of the D-parameter: we will isolate the distinct phase space regions that can contribute to the D-parameter and then sum them up.
As mentioned earlier, to identify the distinct regions that contribute to the D-parameter, it is not sufficient to just measure the D-parameter. This is in contrast to thrust, for example, where the value of thrust regulates all singular configurations. So, with insight from jet substructure calculations, we will measure the D-parameter and another observable on the final state which will control the planar configuration of particles. As this auxiliary observable, we choose to measure the C-parameter. This choice is due to the relation of the C-and D-parameters as defined by the spherocity tensor, but one could choose another observable to accomplish the same task. For the resummation regime when D 1, we will demonstrate that there are three regions defined by the relationship between the values of C and D: 1. Region 1 (Trijets): D C 2 ∼ 1. This is the region where the C-parameter is well described at tree-level fixed-order because the final state of e + e − → qqg is planar, in the center-of-mass frame. This region was studied and resummed in ref. [28], though using the three-jet resolution variable y 3 to define the planar configuration.
2. Region 2 (Hierarchical Trijets): D C 2 1. The radiation that sets the Dparameter in this region is still parametrically separated from the radiation that sets the C-parameter. However, now, the value of the C-parameter must be resummed as well, as there will be large logarithms of C in addition to D in this region.

JHEP02(2019)104
3. Region 3 (Dijets): D ∼ C 2 1. In this region, logarithms of C must be resummed, but there will be no additional logarithms of the D-parameter. The physical configuration of this region corresponds to the dijet limit, for which the D-parameter is no longer an additive observable.
There are no regions with more than three resolved jets that contribute at leading power in the limit that D 1. With more than three jets in the final state, they are no longer constrained to lie in a plane in the center-of-mass frame. The existence of more jets in the final state does not produce any large logarithms, and enforcing D 1 requires that all the jets lie within an angle of order the value of D of the plane.
To accomplish resummation of the C-and D-parameters (as appropriate) in the first two of these three regions, we construct effective theories in each region with the formalism of soft-collinear effective theory (SCET) [31][32][33][34]. The general framework for resummation with multi-jet hierarchies was presented in ref. [35] with an effective theory called SCET + . We factorize the cross section differential in C and D into hard, collinear, and soft contributions that are each defined at natural scales. Resummation is then accomplished by renormalization group evolution between the hard, collinear, and soft scales. While our factorization will be valid to any logarithmic accuracy, for concreteness we will present results at next-to-leading logarithmic accuracy (NLL). By "NLL" we mean that we are able to successfully resum logarithms of the form (α s log D) n in the exponent of the cumulative distribution of the D-parameter. In the language of the factorization theorem, we calculate and resum with one-loop anomalous dimensions.
While we calculate all necessary pieces in regions 1 and 2 for complete resummation to NLL accuracy, for the plots we present in this paper, we will combine them with a simple prescription that resums all but one source of logarithms relevant for NLL. The factorized cross section valid in regions 1 and 2 where D C 2 can be expressed in the compact form of 1 σ Here, H(Q 2 ) is the hard function for production of dijets in e + e − collisions, J i (x i , D) is the jet function describing collinear radiation off of jet i in the final state, and S 123 (x 1 , x 2 , D) is the soft function that describes low-energy gluon emission off of the dipoles formed from pairs of final state jets. H (2→3) (x 1 , x 2 ) is the hard function describing the formation of the third jet; its presence and resummation in this factorization theorem is what enables a consistent combination of regions 1 and 2. The x i are the three-body phase space variables which in the center-of-mass frame are just energy fractions. The C-parameter in these phase space regions is just a function of the x i . The jet and soft functions are implicitly convolved with one another. This expression is valid for all values of x 1 and x 2 such that they enforce D C 2 and only misses contributions from logarithms of the D-parameter from collinear-soft radiation [36].
Because the D-parameter is non-additive in the third region, resummation via factorization is confounded and so we take a different approach. As there is no hierarchy between C and D, only large logarithms of C need to be resummed which suggests a way JHEP02(2019)104 to split the double differential cross section in this region. Recognizing the double differential cross section as a joint probability distribution, we can express it with the definition of a conditional probability: Here, dσ(C)/dD is the differential cross section of D conditioned on the value of C. As written this is an identity: the equality is exact if everything is calculated to the same accuracy. However, in this form, this enables the cross section of C and the conditional cross section of D to be calculated at different accuracy for an approximate result. We use results from the literature for the NLL resummed cross section of the C-parameter, and calculate the conditional cross section to lowest non-trivial order, O(α 2 s ). This produces an approximate NLL accurate double differential cross section in this region which, at any rate, is formally accurate to fixed-order and systematically improvable.
With the resummed cross section calculated in each of these regions, we then just need to match them to construct a cross section that is differential in D 1, for any value of C. Region 3 can then be included by a simple additive matching, appropriately subtracting the overlap with regions 1 and 2. Once this super-cross section has been constructed, then the cross section for the D-parameter exclusively is found by integrating over the value of C: The upper bound of the integral over C is 3/4, as that is the maximum value for the Cparameter when D 1. The lower bound of this integral is 4 3 D 1/2 as this is the minimum value that the C-parameter can take given a value of D. We will show both of these limits in the next section. This procedure then produces a resulting resummed cross section for the D-parameter that is completely inclusive of any other restrictions on the event.
The outline of this paper is as follows. In section 2, we define the C-and D-parameter observables, and present useful, but perhaps non-standard, expressions for them. In sections 3 and 4, we present the factorized expressions for the cross section in regions 1 and 2. By appropriately defining scales in the functions of these factorization theorems, we can combine these regions, which is done in section 5. The region 3 calculation is presented in section 6, which has a unique form because the D-parameter is not additive in this region. Our procedure for consistently combining the cross sections in the three regions is described in section 7. Section 8 describes the scaling of non-perturbative corrections to the Dparameter and the regions for which they are most important. The factorization theorems and completeness of the three regions are validated in section 9 by numerical comparison of our predictions to the fixed-order code EVENT2 [37] at leading-order, O(α 2 s ). The matching of fixed-order to our resummed cross section for D is presented in section 10, along with a comparison to data from the experiments at LEP. We conclude in section 11. Appendices contain explicit calculations and other technical details of the factorization theorems.

JHEP02(2019)104 2 Observable definitions
The C-and D-parameters are defined from the spherocity tensor (assuming all massless particles) [15][16][17]: The sum runs over all particles in the final state of e + e − collisions with Q the center-ofmass collision energy, E i the energy of particle i, and p iα is the α component of particle i's three-momentum. The C-and D-parameters are defined by the eigenvalues of the spherocity tensor λ i , with the ordering λ 1 ≥ λ 2 ≥ λ 3 and tr Θ = λ 1 + λ 2 + λ 3 = 1. The C-parameter is The normalization factors of 3 and 27, respectively, are such that the C-and D-parameters range from 0 to 1. The D-parameter can also be written in terms of energies and relative angles as The expression for D on the second line of this equation will be especially useful, where θ ij is the angle between particles i and j. Note that these identities hold when considering exclusively massless particles; we will assume this throughout this paper. The possible hierarchical relationships between the eigenvalues of the spherocity tensor can be exploited to identify the different regions that we need to consider for resummation of the D-parameter. First, if the eigenvalues are all of the same order: then both C and D are of order-1, and are well-described at fixed-order. Therefore, there is no resummation necessary with this scaling. If there is one eigenvalue that is parametrically small than the others: this corresponds to the physical configuration when the final state is nearly planar. In this limit, the C-parameter reduces to

JHEP02(2019)104
and has an order-1 value and is described at fixed-order. The D-parameter, on the other hand, is parametrically smaller: Therefore D C 2 ∼ 1 which corresponds to region 1 as defined above. The resummation of the D-parameter in this phase space region was done to NLL accuracy in ref. [28]. Instead of using the value of the C-parameter, ref. [28] used the three-jet resolution variable y 3 [29] to constrain the system to have three jets. Note that the maximum value of the C-parameter in this region of phase space occurs when λ 1 = λ 2 = 1/2, and so C ≤ 3/4.
If the three eigenvalues are all parametrically separated: this corresponds to the physical configuration when the final state is nearly planar and approaching the dijet limit. The C-parameter in this limit reduces to as λ 1 → 1 in this limit. The D-parameter is related to it as: Therefore, D C 2 1, which corresponds to region 2 as defined above. Restricting D C 2 1 means that there are three jets in the final state, but one of them is soft, or collinear to another jet. This configuration can be resummed using the factorization theorems of ref. [36] (for the collinear jets) or ref. [38] (for the soft jet). Ref. [39] described how to combine them consistently for the particular application of the D 2 observable [40,41]. Note that because there are two hierarchies, one must resum both C and D/C 2 , which is accomplished in the factorization theorems referenced. Also, in this limit the D-parameter is additive, so its resummation is "simple." There is one final hierarchical relationship between the eigenvalues: (2.12) In this limit, the C-parameter becomes while the D-parameter is D = 27λ 2 λ 3 . (2.14) That is, D ∼ C 2 1 which corresponds to the dijet limit, region 3. Because of this relationship, the D-parameter is no longer additive, but there are no large logarithms of D/C 2 ∼ 1. Therefore, in this phase space region, we only need to resum logarithms of C 1, which is still an additive observable. Resummation of this form is like the "soft haze" region described in ref. [39].

JHEP02(2019)104
In this final region, it is useful to determine the precise relationship between the C-and D-parameters. Note that the D-parameter can be arbitrarily smaller than the C-parameter as λ 3 → 0. However, the maximum value that the D-parameter can take occurs when λ 2 = λ 3 . In this limit, we have the relationship This upper limit will be important to remember when constructing the contributions to the D-parameter in this region. To summarize: in region 1, we only resum D and calculate C at fixed-order; in region 2, we need to resum both C and D; and in region 3, we only need to resum C and calculate D/C 2 at fixed-order. We will now discuss in detail the form of the factorized cross section in each region.

Region 1 (trijets):
We begin by constructing the cross section in region 1, when D C 2 ∼ 1. Because C ∼ 1 in this region, the C-parameter is well-described at fixed-order by the production of three well-separated, energetic jets in the final state. Additionally, because we assume that D 1, the radiation in the final state that is out of the plane defined by those three jets is constrained to be at small angle with respect to the plane, or have low energy. Those emissions that are at small angle to the plane could in principle be at any angle with respect to the energetic jets that set the value of the C-parameter. However, the emission probability for this small-angle radiation is only enhanced by a large logarithm if it is additionally collinear to one of the energetic jets. Therefore, to leading power, the emissions that set the value of the D-parameter are soft (low energy) or collinear to one of the three energetic jets that set the value of the C-parameter.
The cross section in this region of phase space factorizes just assuming hard-collinearsoft factorization of QCD. The form of the factorized cross section in this region is 1 σ To write this cross section, we have used the three-body phase space variables x 1 , x 2 , and x 3 where E i is the energy of the ith jet in the event, which can be operationally defined through an exclusive jet algorithm that finds three final-state jets. However, this multi-differential cross section is not directly measurable in experiment because it requires identification of the quark and gluon jets that compose the final state. Because D 1, the masses of these jets are very small compared to their energies and so, to leading power, their phase space is completely defined by these energy fractions. Q is the total center-of-mass energy and x 1 + x 2 + x 3 = 2. Because of conservation of energy, there are only two independent phase JHEP02(2019)104 space variables, which we take to be x 1 and x 2 , the energy fractions of the quark and antiquark in the final state. The tree-level expression for the 3-jet hard function H 3−jet (x 1 , x 2 ) is simply the leading-order matrix element for e + e − → qqg divided by the cross section for e + e − → qq: 3) The C-parameter as measured on this final state has the value We refer the reader to refs. [35,42] for matrix element definitions of hard functions for multi-jet production and soft functions with multiple Wilson lines. The other functions in the factorization theorem describe emission of collinear radiation or soft radiation that contributes to the measured value of D. In this region of phase space, the D-parameter can be expressed in terms of the out-of-plane momentum of a final state particle as [28] where the sum over j runs over all final state particles. Note that this is directly dependent on the value of C. Importantly, this demonstrates that the D-parameter is an additive observable in this region; this will enable resummation to be accomplished in a relatively simple way. At the level of the factorized cross section, additivity is imposed by convolution (denoted by ⊗) of the various functions that depend on the value of D.
The jet function J i (x i , D) describes the production of radiation that is collinear to energetic jet i and contributes to the value of the D-parameter. Using the expression for the D-parameter as in eq. (3.5), the contribution to D from collinear radiation off of this jet is Here, the sum runs over those collinear particles in jet i, z j is the energy fraction with respect to the jet energy, θ j is the angle of the particle from the jet axis, and φ j is the azimuthal angle about the in-plane jet axis. This angle is 0 or π if the particle is in the plane. For a jet with two particles (sufficient for resummation to NLL accuracy), this simplifies to where now z and 1 − z are the energy fractions of the particles in the jet, and θ is their relative angle. The soft function S 123 (x 1 , x 2 , D) describes low-energy radiation emitted off of the dipoles formed from any pair of the three final state jets. It depends only on the relative angles between the jets, and therefore by momentum conservation on the energy fractions x i . The contribution to the D-parameter from soft emissions can be expressed as: In the first equality, we have written it in terms of the energy E j of the soft emission j and its angle θ j from one of the hard jets and its angle φ above the plane of the three hard jets. In the second equality, we have changed variables to express D in terms of the momentum k ⊥,j of soft emission j out of the plane, angle φ j out of the plane, and rapidity η j with respect to one of the hard jet directions.
To resum this factorized cross section to NLL accuracy in the D-parameter, we must calculate the one-loop anomalous dimensions of the functions in the factorization theorem. The calculation of the anomalous dimensions in this phase space region are presented in appendix A. The relative scales of the functions appearing in the factorization theorems are displayed in figure 1. Because this factorization theorem only required hard-soft-collinear factorization, this scales plot doesn't encode much information. However, when considering the other relevant phase space regions it will be important, and so is a benchmark for additional scales that appear in those cases. With this resummed and factorized form of the cross section, we can then calculate the cross section differential in the C-and D-parameters by marginalizing over x 1 and x 2 : (3.9) While the C-parameter δ-function integral can be done in principle, in practice, this integral is easiest to do by Monte Carlo. The hierarchal jets region, region 2, corresponds to the relationship D C 2 1. Because now both C and D are parametrically less than 1, we must resum both observables. This is accomplished through a re-factorization of the cross section presented in the previous section. There are two possible configurations that must be considered in this region. At leading power, this corresponds to either the gluon jet becoming collinear with one of the quark jets, or the gluon jet becoming soft as compared the center-of-mass energy. Because we still require D C 2 , this means that the final state is still parametrically more planar

JHEP02(2019)104
than it is dijet-like. This will enable a direct connection with the factorization theorem of region 1, which we will discuss in the next section. We will address these two configurations separately. The form of the factorization theorems presented in this section and the method for their consistent combination was first discussed in ref. [39] and the general jet hierarchy configuration was analyzed in ref. [35].

Collinear gluon jet
For concreteness, consider the configuration where the gluon becomes collinear with the anti-quark. In terms of the three-body phase space variables, this corresponds to the limit x 1 → 1. In this limit, the expression for the C-parameter simplifies greatly: The expression for the D-parameter is modified appropriately in this limit, and still has contributions from emissions collinear to and at low energy with respect to the hard(er) jets in the final state. The factorized cross section in this region of phase space was first studied in ref. [36]. Here, we present the derivation of the factorization theorem in this phase space region as a re-factorization of the factorization theorem of eq. (3.1). Starting from that factorized cross section, we then take the limit that x 1 → 1. In this limit, the jet functions J i are not re-factorized; they just take their leading-power value with x 1 → 1. The hard function H 3−jet , however, must be re-factorized because there is a new hierarchy introduced by x 1 → 1. In this limit, the hard function factorizes as Here, H 2-jet is the hard function for dijet production in e + e − collisions with center-of-mass energy Q which is just 1 at tree-level. H 1→2 is the quark collinear splitting function, which is the x 1 → 1 limit of the tree-level e + e − → qqg cross section: We have denoted the energy fraction of the gluon in the splitting as z where The soft function also re-factorizes because the angle between two of the jets that source soft radiation is becoming small. Honest wide-angle soft radiation can only resolve the total qq dipole, while soft radiation boosted along the collinear anti-quark and gluon direction can resolve their splitting. Therefore, the soft function factorizes into these two components:

JHEP02(2019)104
Here, S 12 (C, D) is the soft function for radiation off of the qq dipole and C s (C, D, x 2 ) is the collinear-soft function that describes soft radiation that can resolve the anti-quark-gluon collinear splitting. This cross section is more naturally expressed in terms of the energy fraction z and the splitting angle θ of the gluon off of the anti-quark. In terms of the three-body phase space variables, note that the invariant mass of the anti-quark-gluon system is In the limit where x 1 → 1, the splitting angle θ is therefore This enables an equivalent representation of the C-parameter in this limit: The fully factorized cross section is therefore (4.10) In this expression, we have used eqs. (4.1) and (4.9) to express the cross section exclusively in terms of the gluon energy fraction z, C, and D. In figure 2 we show the re-factorization and hierarchy of scales in this factorization theorem. The anomalous dimensions and treelevel expressions of the hard and soft functions are given in appendix B.1. The anomalous dimensions of the jet functions are found from taking limits of the anomalous dimensions presented in appendix A.5. There is an essentially identical region when x 2 → 1 which is just found with the replacement x 1 ↔ x 2 in this factorization theorem.

Soft gluon jet
The phase space region where the gluon becomes soft corresponds to the limit in which the three-body phase space coordinates become: Note that this corresponds to the limit in which the gluon energy fraction x 3 → 0. In this limit, the C-parameter becomes Just like in the collinear jet region, the expressions for the D-parameter are appropriately modified from its value in region 1, but we won't show the explicit expressions here. The factorization of this soft jet region of phase space was first studied in ref. [38]. As with the collinear jets, we will present the refactorization necessary to derive the cross section in this region. In this soft gluon region, the jet functions J i and soft function S 123 are just modified according to the leading expression in the x 1 , x 2 → 1 limit. This is true because of the assumed hierarchy of D C 2 1, and so the final state is more co-planar than it is dijet-like. The hard function H 3−jet factorizes into the e + e − → dijets function H 2−jet and the function that creates a soft gluon at an arbitrary angle from the dijets, H s :

JHEP02(2019)104
Here, z is the energy fraction of the soft gluon with respect to the total collision energy, (4.14) The tree-level expression for H s is found by taking the appropriate limit of the e + e − → qḡ matrix element In the second line, we have expressed the cross section in terms of the energy fraction z of the gluon and the angle from the quark-anti-quark pair, θ.
The cross section in this soft gluon region of phase space then factorizes as Here, we have expressed the cross section in terms of the value of the C-and D-parameters, and the energy fraction of the gluon, z. Figure 3 displays the hierarchal scales of this factorization theorem. The anomalous dimension of the new hard function H s appearing in this factorization theorem is presented in appendix B.2. The forms of the factorization theorems for either the collinear gluon jet or the soft gluon jet in region 2 both followed from a refactorization of the factorization theorem of region 1. That is, a consistent combination of the cross sections in regions 1 and 2 will result in a cross section for the region where D C 2 , with no constraint on the absolute scale of the C-parameter. One way to do this combination is to match the cross sections of the different regions and then subtract their overlap. This was effectively the approach taken in the calculations of the D 2 observable [40,41] in the two-prong region of phase space in refs. [39,43,44] or the approach advocated in the general hierarchical jet analysis in ref. [35]. With one-loop anomalous dimensions, such a combination of regions would correctly resum all large logarithms for D C 2 to NLL accuracy. However, for our purposes in this paper, there is a significantly simpler way in which to combine phase space regions 1 and 2. We present this procedure which accomplishes approximate NLL accuracy, only missing logarithms from the collinear-soft function C s (C, D, z). We can get away with approximate NLL accuracy for making plots because a complete calculation of the D-parameter requires including region 3 as well, but there currently exists no factorization theorem for this region, as will be discussed in the next section.

JHEP02(2019)104
Recall the factorization theorem for region 1, when D C 2 ∼ 1: In region 2, because we still assume a hierarchy between the values of D and C, the functions set by the measured value of the D-parameter (the jet functions and soft function) are identical to the calculation of the corresponding function in region 1. In the evaluation of region 2, we have just taken the appropriate collinear or soft limit of the emission that

JHEP02(2019)104
sets the value of the C-parameter. This requires refactorization of the hard function and refactorization of the soft function for collinear subjets. In region 1, because C ∼ 1, there is no hierarchy between the center-of-mass energy and the scale of the C-parameter. However, in region 2, there is a hierarchy, and this is resummed by either the function H 1→2 or H s , as appropriate to the collinear or soft limits, respectively. In the following, we present a single formula that resums all logarithms present in the hard and jet functions of regions 1 and 2, and only misses those logarithms that arise from refactorization of the soft function.
To completely describe the region where D C 2 , we need to account for these different regions. To our approximate NLL accuracy, we can do this in the following way. We construct a new factorized form of the cross section: is the hard function for e + e − → qq production, which to lowest order is just 1 and has an anomalous dimension of Using the results presented in appendix A, we then imply the anomalous dimension of the "antenna function" H (2→3) (x 1 , x 2 ) by renormalization group invariance of the cross section to be By the renormalization group invariance of the cross section, we can fix the renormalization scale µ = Q, the center-of-mass energy. This enables us to consistently resum over the entire range when D C 2 . Note that the natural scale of the hard function H(Q 2 ) is set to Q also: µ H = Q. That is, there is no running of the hard function from its natural scale. Therefore, with this prescription H(Q 2 ) = 1. Then, one needs to appropriately set the natural scale of the antenna function H (2→3) (x 1 , x 2 ) in order to resum all logarithms of the C-parameter. When C ∼ 1, we want the natural scale of H (2→3) (x 1 , x 2 ) to be µ H (2→3) = Q so that it takes its tree-level value: That is, when C ∼ 1, the hard function H (2→3) has no large logarithms because we choose µ = Q. When C 2 1, the natural scale of this function should correspond to resumming the hierarchy between the center-of-mass scale Q and the value of the C-parameter. From the results presented in appendix B, this scale is

JHEP02(2019)104
So, we just need to design a scale that interpolates between these two regimes, as a function of the value of the C-parameter.
Our prescription for doing this will be very simple. We just multiply the natural scale of H (2→3) (x 1 , x 2 ) by a factor so that it matches onto the hard scale µ H = Q at the kinematic endpoint. The maximum value of the C-parameter at O(α s ) is 3/4 at which point all phase space variables x i = 2/3. At this kinematic endpoint, the scale of eq. (5.6) is therefore So, we rescale µ H (2→3) by a factor of √ 12: and rearrange the anomalous dimension of this function appropriately: We will use this form of the anomalous dimension to produce our results that we plot later. This scale setting approach exactly reproduces the factorization theorems of region 1 and the soft jet of region 2, so it will successfully resum all of those logarithms. It misses, however, the re-factorization of the soft function of the collinear jet of region 2, and therefore does not completely capture all logarithms of this region, which is why we refer to it as approximate NLL.
6 Region 3 (dijets): The third region defined by the simultaneous measurement of the C-and D-parameters is the most subtle. In regions 1 and 2, because D C 2 , the emissions that set the Cparameter also set the event plane. The emissions that set the D-parameter, however, had too low of a relative k T to affect the event plane, and so the D-parameter is additive. This property is vital for the factorization of the cross section in eq. (5.2), for example. The convolutions over the value of the D-parameter are a consequence of this additivity, and so resummation in regions 1 and 2 is straightforward.
In region 3, where D ∼ C 2 1, the D-parameter is no longer additive. Emissions that set C also affect D, and vice-versa, so the event plane is not well-defined. As such, this seems to suggest that additional emissions need to know about all previous emissions to set the event plane, and therefore the value of D. This is somewhat of a similar situation to that of recoil-sensitive observables such as broadening [18][19][20], and may suggest that resummation requires cross-talk between soft and jet functions in addition to just their measured values of D [45]. At the very least, the resummation in this region will be nonstandard. So, we will approach its calculation from a different perspective, rather than attempting factorization of the cross section.

JHEP02(2019)104
Let p(C, D) be the joint probability distribution of the C-and D-parameters. This is just the (normalized) double differential cross section. Then, from the definition of conditional probability, this can be re-written as where p(C) is the probability distribution of the C-parameter and p(D|C) is the conditional probability of D given C. Expressed in terms of differential cross sections, we have where σ(C) denotes that C is a parameter of the cross section, and not a random variable. So far, this is an identity: as long as everything on both sides of this equation are calculated to the same accuracy, the equality is exact. However, written in this form enables another approach that will allow us to make progress. Our goal will be to calculate the double differential cross section of C and D in region 3 to NLL accuracy. The cross section of the C-parameter dσ/dC has been calculated to high perturbative accuracy both at fixed-and resummed-order [1][2][3]9], and so getting it to NLL accuracy is no problem. The conditional cross section, dσ(C)/dD, on the other hand, contains all of the subtleties of region 3, and (currently) cannot be resummed. However, it can be calculated to fixed-order in region 3, which is what we will do here. As the D-parameter is first non-zero at O(α 2 s ), we will calculate the conditional cross section to that order. So, the expression for our resummed double differential cross section for the C-and D-parameters is The superscripts denote the accuracy to which each factor is calculated. The fixed-order conditional cross section is defined as Note that this conditional cross section is technically O(α s ), because it is formed from the ratio of a cross section at α 2 s to one at α s ; i.e., the first non-trivial order for the numerator and denominator, respectively. The lowest-order cross section of the C-parameter in the where σ 0 is the Born-level cross section for e + e − → qq scattering. In the evaluation of eq. (6.3) for comparison to LEP data in later sections, we will include the effects of a running coupling in the cross section for the C-parameter. The expression including the running coupling is provided in appendix C. This approach to constructing the cross section has been used to study observables that are formally infrared and collinear unsafe, and yet have finite cross sections when all-orders JHEP02(2019)104 effects are accounted for. Such observables are referred to as Sudakov safe [47,48]. In the analogy with the C-and D-parameters, D would represent the unsafe observable, while C would be a safe companion. By measuring the safe companion, the conditional probability of the unsafe observable is calculable in fixed-order perturbation theory. Then, one can marginalize over the safe companion and a finite result is obtained when a Sudakov factor is included to exponentially suppress the singular region of phase space. Here, however, the Dparameter is itself infrared and collinear safe, so the analogy only goes so far. Nevertheless, measuring the C-parameter was necessary to isolate this dijet region of phase space.
As written, eq. (6.3) is not restricted to region 3, where D ∼ C 2 1. There is a non-trivial contribution from the region where D C 2 1, which is already covered by regions 1 and 2. Therefore, we need to eliminate this contribution, which we will do by simple subtraction. That is, to restrict to region 3, we just subtract the limit when D C 2 from the conditional cross section: We calculate the subtracted double differential cross section of region 3 to α 2 s at singlelogarithmic accuracy in appendix C. At single-logarithmic accuracy, we only need to consider the soft and collinear limit to α 2 s , which simplifies the calculation. We were unable to find exact expressions for the results in region 3, and we use numerical interpolation to perform the integrals over C to determine the inclusive D-parameter cross section.
The running coupling here needs to be evaluated at the relative transverse momentum of the emissions. This accomplishes a resummation of some of the collinear logarithms that arise at higher perturbative orders. The exact scales and integrals are worked out in appendix C, but to NLL accuracy, the running couplings can just be evaluated as: Here, Q is the center-of-mass collision energy, and the one-loop running coupling is As mentioned above, this procedure does not strictly account for all logarithms that may be present at NLL accuracy. However, we believe that it only misses potential contributions from soft emissions at angles comparable to the harder emissions that set the value of C (or D). The effect of collinear emissions off of the emissions that set the Cand D-parameters are included at NLL accuracy by evaluating the coupling α s in eq. (6.6) at the relevant k T scale. The resummed cross section of C accounts for logarithmicallyenhanced emissions off of the original quark-anti-quark dipole. Soft emissions from dipoles formed from the emissions that set the value of D, on the other hand, resolve the structure JHEP02(2019)104 of the event plane. We do not expect them to, in general, be correctly accounted for in eq. (6.3). Nevertheless, eq. (6.3) enables a concrete prediction, formally accurate at least O(α 2 s ), and is systematically improvable by simply calculating the conditional probability to higher orders. However, when logarithms of D are comparable in size to the inverse of the coupling, the accuracy of the expression breaks down, and honest resummation of all terms is required.
The approximation used here is not equivalent to modified leading logarithmic accuracy (MLLA) [49], in which one just exponentiates those logarithms that arise at leading-order. For additive observables like C-parameter or thrust, MLLA is a sensible approximation, because it is clear how to improve the approximation. For the D-parameter in this region, it is not obvious how to improve MLLA to higher accuracy, because it is not a priori obvious how logarithms exponentiate at higher orders. This is why we advocate for the conditional probability approach. Further, the effect of multiple emissions setting the value of the D-parameter will only first contribute at NNLL order. Emissions that set both the values of C and D at leading power must be non-strongly ordered in both energy and angle. The leading-order contribution therefore scales like α 2 s log 2 D, because there need to be two emissions to be out of the event plane and those two emissions must be soft and collinear to one of the jets. If there was a third emission that contributed at leading power to D, it cannot be strongly-ordered with respect to the first two emissions, and therefore it contributes at order α 3 s log 2 D, or NNLL.
7 Combining factorization theorems of regions 1, 2, and 3 With the double differential cross sections in regions 1, 2, and 3 established, we now want to combine them and integrate over the C-parameter to produce the cross section exclusive in D only. As established above, the double differential cross section of the C-and Dparameters in the limit where D 1 is just the sum of the contributions from regions 1 and 2 and region 3: In region 3, where D ∼ C 2 , we have already removed the overlapping contribution with regions 1 and 2 which is why the cross section is just a sum. Now, we just integrate this over C. As discussed in section 2, the range of the C-parameter when D 1 is This is thus the bounds of integration for the integral over C: With this integral, we are in some sense done, but it is useful to analyze it a bit further.
Let's first focus on the contributions from regions 1 and 2, where D C 2 :

JHEP02(2019)104
Formally, the cross section is only accurate in the region where D C 2 ; however, the lower bound of the integral extends to the region in which D ∼ C 2 . This is perfectly fine: this cross section nevertheless resums some of the large logarithms in this region, but does not capture all of them. To get all logarithms of D where D ∼ C 2 , we of course need to include the contribution from region 3.
The region 3 expression for the cross section is As we have already explicitly subtracted the contribution from the region D C 2 , the upper bound of C ≤ 3/4 strictly only adds power corrections in D. So, we can safely take the upper bound in this region to ∞: Further, let's write the double differential cross section in this region using conditional cross sections, as in the previous section: In the second equality, we have written the differential cross section of the C-parameter as the derivative of the cumulative distribution of the C-parameter, Σ(C). We can then integrate by parts: In writing this expression, we use the fact that the cross sections vanish at C = ∞. Because our approach to this region calculates the cross section for C to NLL accuracy, while the conditional cross section is only calculated to O(α 2 s ), it will be a bit easier to numerically evaluate the final equality, after integration by parts.

Non-perturbative effects
Our analysis thus far has been restricted to perturbation theory, but we also need to determine the regime in which this analysis is valid. In this section, we estimate the size and effect of corrections due to non-perturbative physics on the D-parameter distribution. The JHEP02(2019)104 separation of the cross section into the three regions enables a straightforward estimation of non-perturbative effects in each region. We can then identify the region with the largest non-perturbative correction to specify the range of validity of our perturbative calculation. The regions with the largest non-perturbative effects are regions 1 and 2, in which D C 2 , which we will justify.
In regions 1 and 2, non-perturbative corrections become important when the natural scales of functions in the factorized cross section approach the non-perturbative scale of QCD, Λ QCD . From the results of appendix A, the function with the lowest scale is the soft function; either S 123 in region 1 and the soft jet of region 2 or S 12 for the collinear jet of region 2. The soft scale µ S Setting this equal to the QCD scale Λ QCD , we find that the D-parameter cross section in this region is dominated by non-perturbative effects when The inequality on the right was found by setting the C-parameter to its largest value, C < 3/4. This linear scaling in Λ QCD of non-perturbative corrections to the D-parameter agrees with what was found in ref. [28].
In region 3, the scale of the C-parameter itself is comparable to that of the D-parameter. Therefore, non-perturbative corrections to the D-parameter are potentially further suppressed by the value of the C-parameter. Non-perturbative effects dominate when the lower bound of the integral over the C-parameter enters the nonperturbative regime. This lower bound is D < 3 4 C 2 , and so the non-perturbative value of the D-parameter is where C NP is the value of the C-parameter at which it becomes non-perturbative. This non-perturbative scale is linearly proportional to the QCD scale, and so the non-perturbative value of D in region 3 is suppressed by two powers of the QCD scale: Because the these corrections scale with energy Q parametrically smaller than for regions 1 and 2, they are formally suppressed. Thus, we only will consider non-perturbative corrections from regions 1 and 2. At the Z-pole, Q = m Z , the D-parameter remains perturbative as long as As we will show in section 10, the best data on the D-parameter from LEP has several bins at values of the D-parameter that are below this value.

JHEP02(2019)104 9 Comparison to EVENT2
Before comparison of our resummed predictions for the D-parameter with LEP data, we will compare the fixed-order expansion to the output of the program EVENT2 [37]. EVENT2 calculates the fixed-order cross section for e + e − → qqg to next-to-leading order. As such, it can predict the leading-order expression for the D-parameter, which is first non-zero at O(α 2 s ). With a program like EERAD3 [21] or NLOJET++ [23], which can calculate the fixedorder cross section for e + e − → qqg to next-to-next-to-leading order, we could, in principle, predict the D-parameter at next-to-leading order (O(α 3 s )). However, it will be non-trivial enough to agree at leading-order, as our resummed result consists of multiple moving parts.
The first thing we need to do is to extract the O(α 2 s ) cross section from our resummed predictions in regions 1, 2, and 3. In region 3, numerical estimates for the single-logarithmic contribution are presented in appendix C.2. The fixed-order expression in regions 1 and 2 can be found from the expressions for the anomalous dimensions in appendix A. The anomalous dimensions for the jet and soft functions is presented in appendix A.5 in Laplace space for the D-parameter. With the initial condition that the tree-level functions are just 1 in Laplace space, we can solve the anomalous dimension equations and inverse Laplace transform to find the expressions in real space. Then, we sum the jet and soft functions and multiply by the tree-level matrix element for e + e − → qqg scattering to find the cross section differential in D and the three-body phase space variables x 1 and x 2 . To find the cross section exclusively in D, we then just integrate over x 1 and x 2 .
Solving the anomalous dimension equations and inverse Laplace transforming, the sum of the jet and soft functions at O(α s ) is The superscript (1) denotes that this is the result at O(α s ). Then, the full cross section triply differential in D, x 1 , and x 2 is Recall that the value of the C-parameter in this region is The cross section inclusive over x 1 and x 2 is found by integrating over them. However, the integration only extends to the point when C = 4 3 D, so the leading logarithmic term

JHEP02(2019)104
in eq. (9.2) stays positive. So, we need to integrate over x 1 and x 2 with this constraint. Then in this region, we find This integral can be done numerically and combined with the result from region 3.

Leading-logarithmic cross section
When comparing to the limit in region 3, it will be useful to explicitly evaluate the leadinglogarithmic cross section for the D-parameter in each color channel. To do this, we take the soft and collinear limit of the qqg matrix element, and keep only the leading terms in the x 1 , x 2 → 1 limit. In doing this, we can rewrite the cross section in terms of the gluon's energy fraction z and splitting angle θ 2 : In this expression, we note that It's also useful to change variables to the ratio y = D/C 2 , so the leading-logarithmic limit can be identified with y 1. The cross section then becomes To continue, we can integrate over θ 2 ∈ [0, 1] and find To isolate the leading-logarithmic cross section in each color channel, we can set all numerical factors in logarithms to 1 and remove subleading logarithmic terms. The cross section then further reduces to  To determine the cross section exclusive in only D, we just integrate over y. This produces:

O(α 2 s ) cross section and EVENT2
With the numerical integral of eq. (9.4) representing the contributions from regions 1 and 2, and the results from appendix C.2 for region 3, we can validate our expression for the Dparameter cross section with EVENT2. The results of this comparison are shown in figure 4. Shown here are the differential cross sections plotted logarithmically in the value of D. The individual contributions are separated into the three color channels at O(α 2 s ), C 2 F , C F C A , and C F n f T R . The Born-level cross section σ 0 has been divided out, as well as the appropriate color factor and an overall factor of (α s /2π) 2 . On the left in figure 4, we compare the three different color channels as computed with EVENT2 and from our factorization theorems for the different regions (labeled as NLL). Excellent agreement is observed, with only a relative constant offset between the EVENT2 results and our calculation. This difference is beyond NLL accuracy, and so demonstrates consistency of this calculation.
The plot on the right of figure 4 focuses in on the contribution from region 3 specifically. In this plot, we have subtracted the region 1-2 contribution from the EVENT2 distributions. The result of this subtraction is a residual cross section that is linear in each color channel on this logarithmic plot. This demonstrates that eq. (9.4) does correctly describe the leading-logarithmic contributions that scale like α 2 s log 3 D/D and α 2 s log 2 D/D. The residual α 2 s log D/D terms that remain are described by contributions from region 3. As calculated in appendix C.2, the approximate logarithmic slopes of the different color To evaluate these coefficients, we had to numerically integrate over the non-strongly-ordered soft and collinear region of phase space in each color channel. The C F channel is significantly simple that all but one numerical integral could be done analytically, so its slope is known to high precision. In the n f T R channel, two numerical integrals must be done, but there is no strongly-ordered soft divergence which makes the matrix element regular over a large region of phase space. So, the slope is also known well for this color channel. This is in contrast to the C A color channel in which the strongly-ordered soft and collinear divergences in the matrix element make numerical integration significantly challenging. Nevertheless, these slopes are in good agreement with the result of EVENT2, given these caveats. However, as figure 4 demonstrates, the C F channel dominates the cross section, so any uncertainties on numerical integration of the result from the C A channel will only have a small effect.

Comparison to LEP data
With our resummed predictions for the D-parameter in hand, in this section we compare to data. We start in this section by first comparing our resummed predictions to fixed-order distributions to gain an understanding of how resummation is affecting the cross section. We will compare our resummed results to both leading fixed order from EVENT2 as well as next-to-leading fixed order of the D-parameter from ref. [23] which used the NLOJET++ program. Our resummed prediction is a combination of the region 1, 2, and 3 predictions as discussed in earlier sections, and leading-order matching. This is accomplished by simple additive matching: The leading-order cross section ("LO") is the result from EVENT2, the resummed cross section ("resum") is the combination of regions 1, 2, and 3, and the fixed-order expansion of the resummed result ("resum, α 2 s ") was discussed in section 9. Figure 5 compares these three cross sections plotted logarithmically (left) and linearly (right). Our resummed predictions terminate when scales hit the QCD Landau pole to ensure that they remain in their region of validity. As observed in refs. [23,24], the NLO K-factor is about a factor of 2, demonstrating that higher-order effects are very important for the D-parameter. These plots demonstrate that this large NLO K-factor seems to be largely reproduced by resummation, and so is dominantly described by higher-order logarithms. Resummation appears to become important below a value of about D 0. the cross section by a factor of 2. There are many scales that exist in the resummed cross section, but the most sensitive scale is that of the soft function in the region 1-2 factorization theorem. The scale uncertainty band in figure 5 is exclusively from varying the soft function scale up and down by a factor of 2. We leave a detailed uncertainty analysis to future work. It is also useful to discuss which phase space regions are relevant in these plots. Recall that the upper limit on the value of D in the resummation region is To be able to compare our prediction to data, we also need to include the effects of non-perturbative physics. In section 8, we demonstrated that regions 1 and 2, where D C 2 , has the largest non-perturbative corrections, and the D-parameter is additive in this region. The average size of non-perturbative corrections to the D-parameter in regions 1 and 2 are Because D is additive in these regions, the dominant non-perturbative effects can be incorporated by simply translating the argument of the perturbative cross section in region 1 [50][51][52][53]: Non-perturbative corrections can be included in a more detailed and quantitative way with shape functions [54,55], though such an analysis is beyond the scope of this paper. All we will do to include the non-perturbative effects here is this translation of the regions 1-2 cross section argument. In the following we set Λ QCD = 1 GeV so that D NP 0.038. The most precise measurements of the D-parameter were performed at LEP at a centerof-mass energy of the Z boson mass [56][57][58]. While the D-parameter was measured at several other energies, we restrict our comparisons to the Z pole to be able to make quantitative comparative statements. This LEP data is compared to our resummed, leading-order matched, and including the non-perturbative shift in figure 6. Good qualitative agreement is observed between our prediction and data, demonstrating that resummation is important for describing the data, especially in the turn-over observed in the log-plot in figure 6a in data around D 0.1. Our resummed prediction significantly undershoots the data at large values of D, which we ascribe to a lack of higher-fixed order corrections. We only match our resummed prediction to leading order, and so we miss the large K-factor between leading and next-to-leading order. These plots also show that next-to-leading order predictions do indeed improve agreement with data over leading-order, at large values of D. A complete matching to next-to-leading order that properly resums all logarithms requires next-to-next-to leading order resummation, which we leave to future work.
Fixed-order studies of the D-parameter [23,24] noted that to match next-to-leading order predictions to LEP data requires an excessively small choice of renormalization scale. These studies, however, did not include any estimates of non-perturbative corrections, which was first done in ref. [28]. The very small renormalization scale in the fixed-order studies is a proxy for these non-perturbative corrections. However, as our analysis in this paper and figure 6 demonstrates, these non-perturbative corrections are necessary for a quantitative description of data, in addition to resummation. The D-parameter is sensitive to infrared scales in a very different manner then familiar observables like thrust, and so a complete quantitative description requires the inclusion of fixed-order, resummation, and non-perturbative physics. Figure 6. Comparison of our resummed prediction with non-perturbative corrections and the next-to-leading fixed order prediction to LEP data on a logarithmic scale (left) and a linear scale (right). DELPHI represents data from ref. [56], L3 represents data from ref. [57], and OPAL represents data from ref. [58]. Uncertainties in the data are comparable in size to the plot markers. The termination of the resummed prediction just below D 0.1 is where the non-perturbative corrections become large.

Conclusions
Resummation of the D-parameter throughout its phase space requires techniques that were only recently developed in the context of jet substructure for LHC applications. Unlike familiar observables like thrust, the value of the D-parameter doesn't isolate individual phase space configurations as it is only sensitive to the aplanarity of the final state hadrons in e + e − collision events. To separate phase space configurations, we additionally measure the C-parameter on the final state, and different hierarchical relationships between D and C correspond to different phase space boundaries. One of these regions, the near-toplanar three jet configuration, has been studied before, but the two other regions are novel. We presented all-orders factorization theorems or systematically-improvable calculations in each region and compared resummed and matched predictions including estimates of nonperturbative effects to precision data measured at LEP. This work is just the beginning for analytically understanding complicated event shapes measured at LEP. A precision program at the level of thrust or C-parameter for observables like the D-parameter is not likely any time soon, but because D probes significantly exotic phase space configurations, its uncertainties, both experimental and theoretical, are unique. Therefore, predictions of the D-parameter at any accuracy are a detailed cross check on the consistency of strong coupling values, for example. With improved predictions for the D-parameter, they might feed into Monte Carlo tunings. Four-jet observables at e + e − colliders are the first that are directly sensitive to the self-coupling of the gluon and nontrivial color flow, and therefore are used to set parameters like color reconnections within Monte Carlo simulations.
In the regions where the D-parameter is additive, we are able to derive all-orders factorization theorems for which each function is just convolved with one another. When JHEP02(2019)104 D ∼ C 2 , emissions that set the D-parameter therefore also set the event plane, and this backreaction is responsible for the D-parameter to lose additivity and prohibit any standard factorization. Related issues with factorization breakdown in endpoint regions of phase space have also been observed in studies of observables sensitive to multi-prong jet substructure [39,44,59]. Nevertheless, this does not imply that factorization is impossible. With motivation from several examples, factorization of multi-differential cross sections away from strongly-ordered phase space boundaries would be a major advance and enable precision studies of a large class of observables relevant at both LEP and LHC.
The analysis of this paper also demonstrates that we continue to learn new ways of calculating in perturbative QCD that can be applied to old problems. LEP left a legacy of numerous precision event shapes, of which only a very few are understood at sufficient accuracy to be used for extracting α s , for example. We hope that this rich data inspires more studies and enables more advances in analytically studying QCD in extreme phase space regions.

Acknowledgments
A.L. thanks Andrea Banfi for originally suggesting this problem, for collaboration in early stages of this work, and comments on the manuscript. A.L. also thanks Ian Moult and Duff Neill for comments on the manuscript, and Wouter Waalewijn and Frank Tackmann for addressing an error in the original description of the region 2 factorization theorems.

A.1 Operator definitions of functions in factorization theorem
In this appendix, we provide the SCET operator definitions for the jet and soft functions that appear in the factorization theorem of region 1, eq. (3.1). For a quark jet function J q (x q , D), it is defined as the following matrix element of SCET operators: Here, we assume that the jet direction is defined by the light-like vector n, the opposite direction is defined byn for which n ·n = 2. In this expression, P is the momentum measurement operator which, for a collinear jet, has zero perpendicular component P ⊥ , and whose large light-cone momentum componentn · P is set equal to the center-of-mass energy Q times the energy fraction x q . χ(x) is the gauge-invariant collinear quark field defined as The δ-function in this expression simply sets the large light-cone component of momentum to be ω, ξ n (x) is the n-collinear quark field in SCET that transforms in the fundamental representation of SU(3) color, and W n (x) is a collinear Wilson line that ensures gauge invariance:

JHEP02(2019)104
The sum runs over all permutations of exponential factors, g is the QCD coupling, and A n (x) is the Lie algebra-valued n-collinear gluon field. Finally, the operatorD measures the D-parameter on the collinear quark field, as defined by the leading-power expression for the D-parameter in the collinear limit, eq. (3.6). The gluon jet function, J g (x g , D), is defined similarly, as Now, the gauge-invariant collinear gluon field is The covariant derivative D µ ⊥ is composed of momentum operators and gluon fields as where the ⊥ subscript denotes that the momentum and fields only have non-zero components transverse to the n direction; that is, describing physical, transverse polarized gluons. The soft function that appears in the factorization theorem is defined as a matrix element of soft Wilson lines: The soft Wilson line in the i direction is where n i is the light-like direction that defines i. P is the path-ordering operator and the color representation in which the soft gluon field A(x) is evaluated is implicit, as defined by i's particle type. T andT are the time ordering and anti-ordering operators and the operatorD measures the value of the D-parameter in the soft limit, eq. (3.8). Note that the Wilson line is outgoing, as all colored particles are produced in the final state from the hard collision. Finally, the dependence of the soft function on the energy fractions x 1 and x 2 is really a shorthand for dependence of the soft function on the angles between the Wilson lines, as soft emissions cannot resolve the energy of hard particles.

A.2 Hard function
The leading-order hard function for the process e + e − → qqg is just the tree-level matrix element: Its one-loop anomalous dimension is [42]

JHEP02(2019)104
In terms of the C-parameter, this can be expressed in the suggestive form as

A.3.1 Quark jet
The phase space and matrix element for collinear emissions from a quark is The D-parameter jet function is then Using the +-function expansion we find We only expand to O( −1 ) to extract anomalous dimensions. Laplace transforming the D-parameter, we find Its anomalous dimension is therefore

A.3.2 Gluon jet
The jet function for a gluon jet on which the D-parameter is measured is Using the +-function expansion, we find Laplace transforming the D-parameter, we find Its anomalous dimension is therefore (A.20)

A.4 Soft function
The one-loop calculation of the soft function for an emission off of jets i and j for the D-parameter in this region can be written as We can choose coordinates such that Then, n j · k can be expressed as The soft function is then The remaining integrals can be done with appropriate changes of variables. With x = e −η , the integrals become The integral over y can be changed to range over y ∈ [0, 1]: Here, 2 F 1 (a, b, c; z) is the hypergeometric function. Then, the remaining integral over x is

JHEP02(2019)104
Here, the symbol means that the expressions agree up through 0 order. To do this, we expanded the hypergeometric function as 2 F 1 2 , This can be expressed as an integral on x ∈ [0, 1]: The first integral can be broken into two regions, corresponding to the two divergences at x = 0 and x = 1: Then, putting it together, we find Then, we find in total (in MS) Laplace-transforming in the D-parameter, we have

JHEP02(2019)104
The anomalous dimension is therefore The color matrix products for e + e − → qqg are: We can then sum over the colors to find: This can be expressed in terms of the C-parameter in the suggestive form: .

A.5 Summary of anomalous dimensions
Here, we summarize the anomalous dimensions for this region. We have: .
As required, the sum of these anomalous dimensions is 0.

B Region 2 calculations
There are two configurations that contribute to region 2, when D C 2 1. The emissions that set the value of the C-parameter can be hard and collinear or they can be soft and wide-angle. We will consider each of these in turn.

B.1 Collinear region
The factorization theorem for the collinear contributions is 1 σ The one-loop hard function for e + e − → qq production is [36,42,60,61] H( Its anomalous dimension is thus The leading-order hard splitting function for collinear gluon emission off of the anti-quark is where x 2 is the energy fraction of the anti-quark in the splitting. To write this, we have used the expression for the C-parameter in this region as Its one-loop anomalous dimension is [17,36,62] We can use the results from the end of section A.4 to evaluate the anomalous dimensions of the wide-angle soft function, S 12 (C, D). The product of color factors is and the angle between the two jets in the event is 1 − cos θ 12 = 2. Therefore, its anomalous dimension in Laplace space is The anomalous dimension of the collinear-soft function C s (C, D, x 2 ) can also be found from the results at the end of section A.4. We take the x 1 → 1 limit of the soft function anomalous dimension in region 1: Then, subtracting the wide-angle soft function's anomalous dimension of eq. (B.8), we find

B.1.5 Summary of anomalous dimensions
Here, we summarize the anomalous dimensions for this region. We have: As required, the sum of these anomalous dimensions is 0.

B.2 Soft region
The factorization theorem for the soft contribution is 1 σ The hard function, H(Q 2 ), is defined above in eq. (B.2). The leading-order hard splitting function for soft gluon emission off of a quark and anti-quark dipole H s (x 3 , C) can be found from changing variables in the soft limit of the e + e − → qqg matrix element. We have That is, the hard splitting function is where x 3 = 2 − x 1 − x 2 is the energy fraction of the gluon. The anomalous dimension is [38,63] γ All other functions in this region are just the x 3 → 0 limits of the corresponding functions from region 1.

C Region 3 calculations
C.1 C-parameter at NLL accuracy The C-parameter was first resummed to NLL accuracy in ref. [46], by relating its value in the limit as C → 0 to thrust at NLL from ref. [64]. The cumulative distribution at NLL accuracy is Here, the function R(C, µ) is the radiator, which to NLL is Γ(x) is the Euler Gamma function, and γ E = 0.5772 . . . is the Euler-Mascheroni constant. The function λ(C) is and the function S(C) is β 0 and β 1 are the one-and two-loop β-function coefficients which are T R n f , (C.5)

JHEP02(2019)104
The coupling α s ≡ α s (µ) is evaluated at the scale µ, which is of the order of the center-ofmass collision energy, Q.
The C-parameter has been resummed to N 3 LL accuracy by relating it to thrust by factorizing the cross section with soft-collinear effective theory [12]. While such a factorization is useful for resumming the D-parameter in regions 1 and 2, for compactness here, we prefer to use the form from the earlier literature. The differential cross section of the C-parameter is just the derivative of this expression: For calculating the conditional cross section of the D-parameter given C, we also need the cross section expanded to lowest order in α s . Here, we choose to include resummation of collinear logarithms through one-loop running of the coupling as well, as this will result in successfully resumming more logarithms of D. This can be calculated by simply differentiating the radiator with respect to C and keeping those terms that correspond to one-loop running. That is, the lowest-order differential cross section of the C-parameter including the effects of one-loop running of α s is 1 σ Expanding this to lowest order in α s produces the result in eq. (6.5).

C.2 Double differential cross section for C-and D-parameters
To calculate the C and D parameters in this region of phase space we have It is most convenient to express the phase space integrals in terms of energy E, polar angle θ, and azimuthal angle φ for both emissions. In the soft and collinear limit, the phase space integral is The maximum value of the energy in the soft and collinear limit is arbitrary, but we will find that the choice of E max = Q/2 is convenient. Only the relative angle between the particles 1 and 2 is important, so we can integrate over one the particles' φ angles. Doing this, we find Then, the integrals we need to do are

JHEP02(2019)104
Now, we would like to pull out all of the scales from the δ-functions, so the integrands are purely O(1). First, energy Q dependence is trivial: where z 1 and z 2 are now energy fractions with respect to Q/2. In doing this, we have explicitly written the α s dependence, multiplying by (4πα s ) 2 . Also, to resum collinear logarithms, we evaluate the coupling at the transverse momenta of particles 1 and 2. To NLL accuracy, we only need one-loop running of the coupling, which is (C.14) Here, Q is a reference scale, which we will take to be the center-of-mass collision energy. β 0 is the one-loop coefficient of the QCD β-function Because region 3 corresponds to D ∼ C 2 or that k ⊥1 ∼ k ⊥2 , to NLL accuracy, we can just evaluate the couplings at the same scale. So, we have where k ⊥1 = z 1 θ 1 Q 2 . With this rescaling, the C and D parameters are So, the cross section with the explicit δ-functions is Now, we will make the change of variables

JHEP02(2019)104
The cross section then becomes To write this expression, we note that in the soft and collinear limit, the matrix element |M(k 1 , k 2 )| 2 ∝ z −4 1 θ −4 1 , which we have pulled out and made explicit. Now, using the C-parameter δ-function, we can do the integral over z 1 . We find The integral over θ 1 can then be done as the only non-trivial dependence is in the running coupling. We can simplify the expression for the running coupling as follows. By selection of this region, uv 2 ∼ 1 and so to NLL accuracy, we can just ignore the multiplicative factors in the running coupling scale. That is, we can simplify this expression to To continue, the phase space constraints imposed by the Θ-functions can be manipulated into

JHEP02(2019)104
Then, to NLL accuracy, by integrating over θ 1 , the cross section becomes To evaluate this expression, we have used the behavior of the integral over the running coupling In what follows, we use the shorthand notation To evaluate the remaining δ-function integral, we will integrate over φ. To do this, it is convenient to change variables to x = cos φ , (C.27) so that the cross section becomes

JHEP02(2019)104
Doing the integral over x, we then find Here, A is the solution of the δ-function: It is useful to change variables to D and y, where Then, C = D 1/2 y −1/2 and the Jacobian from this change of variables is With this change, the cross section becomes This is all the further we can go without an explicit form for the matrix element. In the C F color channel, the matrix element for two soft gluon emission is [65] |M C F (k 1 , k 2 )| 2 = 16g 4 C 2 where we have expanded in the collinear limit and used the changes of variables developed above. We have also multiplied by a factor of 2 in the soft and collinear limit to account for emission off of the quark or anti-quark in the final state. In this color channel, the matrix element is then It is then useful to make the change of variables t = uv 2 (C. 36) in place of the u integral. The cross section with this change of variables becomes The integral over v can be done and the integral can be massaged into a simple form, with one integral over t ∈ [0, 1] remaining: To isolate the regime in which y ∼ 1 (D ∼ C 2 ), we need to remove the y 1 limit. In this limit, y and t scale similarly, so we can expand the integral above to first order in y, t 1: We just need to subtract this from eq. (C.38) to determine the cross section in the region where y ∼ 1. It is also useful to verify that this expression agrees with the leadinglogarithmic calculation of section 9. To leading-logarithmic accuracy, we can set all numerical factors in logarithms to 1, yielding 1 σ 0 d 2 σ y 1,LL dD dy = α 2 s C 2 Then, to calculate the leading-logarithmic cross section for the D-parameter, we just integrate over y: This agrees with the results of section 9.1.
We are unable to do this integral exactly, and to implement the resummation of region 3, we choose to just numerically integrate over y. To compare to the O(α 2 s ) results from EVENT2, we need to integrate over y. To good approximation, we find the following result:

C A channel
In the C A color channel, the matrix element for two soft gluon emission in the soft and collinear limit is [65] |M C A (k 1 , k 2 )| 2 = 128 (4πα s ) 2 Q 4 z 4 1 θ 4 where we have expanded in the collinear limit and used the changes of variables developed above. We have also multiplied by 2 to account for collinear emission off of the quark or anti-quark. As everything else in the integral is symmetric in cos φ → − cos φ, we can explicitly do this symmetrization: .
The cross section in the C A color channel is then We need to make the replacement then that sin 2 φ = y (1 + uv 2 ) 2 3uv 2 , (C. 46) as well as the change of variables t = uv 2 . (C.47)

JHEP02(2019)104
The cross section is then To restrict to the region where y ∼ 1, we need to subtract the limit in which y 1. Unlike the C F color channel, this expansion in the C A channel is challenging, so we just extract the limiting behavior numerically. We find The leading logarithmic term agrees exactly with the expression from the region 1/2 factorization theorem derived in section 9.1. Subtracting this from the full expression then restricts to y ∼ 1. We are unable to do this integral exactly, and to implement the resummation of region 3, we choose to just numerically integrate over y. We find the following approximate result: (C.50)
In this color channel, the cross section is then We need to make the replacement then that sin 2 φ = y (1 + uv 2 ) 2 3uv 2 , (C. 54) as well as the change of variables t = uv 2 . (C.55) The cross section is then