Factorization for groomed jet substructure beyond the next-to-leading logarithm

Jet grooming algorithms are widely used in experimental analyses at hadron colliders to remove contaminating radiation from within jets. While the algorithms perform a great service to the experiments, their intricate algorithmic structure and multiple parameters has frustrated precision theoretic understanding. In this paper, we demonstrate that one particular groomer called soft drop actually makes precision jet substructure easier. In particular, we derive a factorization formula for a large class of soft drop jet substructure observables, including jet mass. The essential observation that allows for this factorization is that, without the soft wide-angle radiation groomed by soft drop, all singular contributions are collinear. The simplicity and universality of the collinear limit in QCD allows us to show that to all orders, the normalized differential cross section has no contributions from non-global logarithms. It is also independent of process, up to the relative fraction of quark and gluon jets. In fact, soft drop allows us to define this fraction precisely. The factorization theorem also explains why soft drop observables are less sensitive to hadronization than their ungroomed counterparts. Using the factorization theorem, we resum the soft drop jet mass to next-to-next-to-leading logarithmic accuracy. This requires calculating some clustering effects that are closely related to corresponding effects found in jet veto calculations. We match our resummed calculation to fixed order results for both $e^+e^-\to$ dijets and $pp\to Z+j$ events, producing the first jet substructure predictions (groomed or ungroomed) to this accuracy for the LHC.


Introduction
The high luminosity proton collisions at the Large Hadron Collider (LHC) enable an unprecedented sensitivity to rare and high scale physics. The cost of such high luminosities is the presence of significant amounts of pile-up radiation present in every event, arising from numerous secondary proton collisions per bunch crossing. Pile-up is truly uncorrelated with the hard scattering and can contaminate any potential measurement. This is particularly important for measurements made on jets, for which pile-up can effect a large systematic bias in observables like the jet mass. In searches for resonances that decay to boosted electroweak objects which have definite masses, pile-up can significantly degrade the ability to separate signal from background. Over the past several years, numerous methods [1][2][3][4][5][6][7][8][9][10][11] have been developed for grooming jets and events for pile-up mitigation and removal, and are now standard experimental tools at both ATLAS and CMS experiments. Especially in analyses of jets, measurements made at the LHC often involve some form of grooming.
With this motivation, it is imperative to understand these jet grooming techniques from first principles QCD. There have been a few studies of the theoretical aspects of jet groomers [6,8,12,13], with predictions for jet-observable distributions calculated to next-to-leading logarithmic (NLL) accuracy for two widely used jet groomers: the modified mass drop tagger (mMDT) and soft drop. These explicit analytic studies showed that these jet groomers not only have desired experimental properties, but can also dramatically simplify theoretical calculations as compared to their ungroomed counterparts. Non-global logarithms (NGLs) that arise from correlations between in-and out-of-jet scales have proven to be a significant obstruction to resummation of ungroomed jet observables to NLL accuracy and beyond. In particular, it was demonstrated by explicit calculation in Refs. [6,8,12] that mMDT and soft drop groomers eliminate the leading non-global logarithms in jet mass distributions [14]. mMDT and soft drop pave the way for systematically improvable resummed predictions of jet observables.
In this paper, we open the door to systematically improvable jet substructure calculations by presenting an all-orders factorization theorem for the soft-drop [8] groomed observables using soft-collinear effective theory (SCET) [15][16][17][18]. An overview of the method we discuss here and some of our results were presented recently in Ref. [19]. This paper provides a more detailed presentation of those results as well as a derivation of the factorization formula and its remarkable properties.
The soft drop groomer walks through the branching history of a jet, discarding soft branches until a sufficiently hard branching is found. This is enforced by effectively requiring where E i and E j are the energies of the particles in that step of the branching, θ ij is their relative angle, and R is the radius of the jet. z cut is a parameter that sets the scale of soft, wide angle emissions in the jet; the typical value is z cut = 0.1. β is a parameter that controls the aggressiveness of the groomer: β = ∞ removes the groomer, β = 0 coincides with mMDT and is simply an energy cut, and β < 0 removes all soft and collinear singularities. We will consider β ≥ 0. If Eq. (1.1) is not satisfied, the softer of the two branches is removed from the jet, and the grooming procedure continues on the harder branch. When Eq. (1.1) is satisfied, the procedure terminates and the groomed jet is returned. For concreteness, on this groomed jet, we measure the two-point energy correlation functions e (α) 2 with angular exponent α > 0 [20][21][22].
In e + e − → dijets events, the factorization formula we derive in this paper for soft-drop groomed left and right hemisphere jets is: This factorization theorem applies when z cut 1 and the left-and right-hemisphere energy correlation functions are asymptotically small: e (α) 2,L , e (α) 2,R z cut 1. We illustrate the physical configuration corresponding to this factorization theorem in Fig. 1. In Eq. (1.2), H(Q 2 ) is the hard function for e + e − → qq. S G (z cut ) is the global soft function, which is only sensitive to the scale set by z cut since all of its emissions fail soft drop. S C (z cut e (α) 2,L ) is a soft function that is boosted along the direction of the jet in the left hemisphere; its corresponding modes are referred to as collinear-soft [23][24][25][26][27][28]. Emissions in S C (z cut e 2 ) denotes the jet modes.
As we will explain in detail, there are several important consequences of this factorization formula. Because the formula depends on the observables e (α) 2,L , e (α) 2,R only through collinear objects each of which has a single scale, there are no non-global logarithms. The elimination of the purely soft contribution also makes the shape of soft-drop groomed jet shapes largely independent of what else is going on in the event. For example, the shape of the left hemisphere jet mass is independent of what is present in the right hemisphere. Additionally, the scale associated with the collinear-soft mode is parametrically larger than the soft scale associated with ungroomed masses, so non-perturbative corrections such as hadronization are correspondingly smaller.
This factorization theorem allows us to go beyond NLL accuracy to arbitrary accuracy. In this paper, we show that next-to-next-to-leading logarithmic (NNLL) accuracy is readily achievable. We focus on α = 2 where the two-point energy correlation function is equal to the squared jet mass (up to a trivial normalization). This lets us extract most of the necessary two-loop anomalous dimensions from the existing literature. For β = 0, the global soft function S G (z cut ) is closely related to the soft function with an energy veto [28,29] which is known to two-loop order. There are additional clustering effects from the soft drop algorithm, but these are straightforward to calculate. Interestingly, we find that the clustering effects in the soft drop groomer are intimately related to similar effects observed in jet veto calculations [30][31][32][33][34]. For β = 1, we compute the two-loop anomalous dimension of S G (z cut ) numerically using the fixed-order code EVENT2 [35]. 1 We thereby achieve full NNLL resummation for the soft-drop groomed jet mass. 2 It is straightforward to generalize from e + e − to pp collisions, since the distribution is determined by collinear physics within the jet, independent of the initial state. The main new ingredient in pp collisions is that jets may be initiated by quarks or gluons. As we will show, soft-drop grooming the jet enables an infrared and collinear safe definition of the jet flavor at leading power in e (α) 2 and z cut by simply summing the flavors of partons in the groomed jet. Using this procedure, we are able to match our NNLL resummed distribution of soft-drop groomed jet mass to fixed order results for pp → Z + j events (including relative O(α 2 s ) corrections to the Born process). The outline of this paper is as follows. In Sec. 2, we review the definition of the soft drop grooming algorithm and the energy correlation functions. In Sec. 3, we present the factorization theorem for soft-drop groomed energy correlation functions in e + e − → dijets events. In this section, we will also present a detailed power-counting analysis of soft-dropped observables to determine the range of validity of the factorization theorem. Our factorization theorem has many non-trivial consequences, which we review in Sec. 4. These include absence of non-global logarithms, process independence, and small hadronization corrections. In Sec. 5, we describe and present the ingredients necessary for NNLL resummation. Here, we also describe our method for extracting anomalous dimensions from EVENT2. We then match our NNLL results with fixed-order calculations for e + e − collisions in Sec. 6 and for pp → Z + j events in Sec. 7, comparing with Monte Carlo simulations in each case. In Sec. 8, we summarize and conclude. The calculational details for NNLL resummation are collected in appendices.

Observables
In this section, we review the soft drop grooming algorithm and the energy correlation functions. Although previous work has focused on jets produced in pp collisions, we will provide definitions for both lepton and hadron collider environments.

Soft Drop Grooming Algorithm
Given a set of constituents of a jet with radius R, the soft drop grooming algorithm [8] proceeds in the following way: for jets in e + e − collisions is where E i , E j are the energies of particles i and j and θ ij is their relative angle. p is a real number that defines the particular jet algorithm. For jets produced in pp collisions, the k T clustering metric is where p T i , p T j are the transverse momenta of particles i and j with respect to the beam and R 2 ij is their relative angle in the pseudorapidity-azimuth angle plane. While the original implementation of soft drop was restricted to reclustering with the Cambridge/Aachen algorithm (p = 0) [48][49][50], we will also briefly consider reclustering with the anti-k T algorithm (p = −1) [51] in Sec. 5.2.
2. Sequentially step through the branching history of the reclustered jet. At each branching, check the soft drop criterion. For e + e − collisions, we require This is known as the soft drop criterion. If the branching fails this requirement, then the softer of the two daughter branches is removed from the jet. The soft drop groomer then continues to the next branching in the remaining clustering history. For pp collisions, the soft drop criterion is 3. The procedure continues until the soft drop criterion is satisfied. At that point, soft drop terminates, and returns the jet groomed of the branches that failed the soft drop criterion.
Once the jet has been groomed, any observable can be measured on its remaining constituents.

Energy Correlation Functions
On jets that have been groomed by soft drop, we measure the two-point energy correlation functions [20][21][22]. We do this mainly for concreteness; the general properties of the factorized formula we will present apply for a much broader class of observables. For jets in e + e − collisions, the two-point energy correlation function e where E J is the sum of the energies of particles in the jet, the sum runs over distinct pairs i, j of particles in the jet, p i is the four-vector momentum of particle i, and the angular exponent α is required to be greater than 0 for IRC safety. For α = 2 and a jet that has massless constituents, the two-point energy correlation function reduces to the normalized, squared jet mass: The energy correlation functions have the nice property that they are insensitive to recoil effects [20,40] and do not include explicit axes in their definition. For jets produced in pp collisions, the energy correlation functions are appropriately modified by replacing spherical coordinates with cylindrical coordinates: where p T J is the transverse momentum of the jet and R ij is the separation of particles i and j in the pseudorapidity-azimuthal angle plane. For jets at central rapidities and in the limit that all particles in the jet are collinear, Eq. (2.7) reduces to Eq. (2.5). This property in particular will enable us to recycle results calculated in e + e − collisions to the case of pp collisions.

Factorization Theorem
In this section, we derive the factorization formula for energy correlation functions measured on soft-drop groomed jets in the region of phase space where e (α) 2 z cut 1 using SCET [15][16][17][18]. We begin with a power counting analysis based on the scales e (α) 2 and z cut relevant to soft-drop groomed jets. This enables us to identify all modes and their momentum scalings that contribute at leading power. Using these scales and the associated modes, we derive the factorization formula. We then show that the jet function in the factorization formula can be re-factorized due to a collinear-soft mode which decouples from the collinear-but-not-soft modes as a result of soft drop.

Power Counting and Modes
For jets on which the soft drop groomer is applied and the energy correlation functions are measured, there are three relevant dimensonless scales: the jet radius R, the soft drop parameter z cut , and e (α) 2 . Typically, jet radii are R ∼ 1 . We are interested in the singular region e (α) 2 → 0 for a fixed value of z cut . Thus we can assume e (α) 2 z cut . We will also assume z cut 1 to refactorize the jet function. The limits R 1 or z cut ∼ 1 could be considered as well, but are beyond the scope of our analysis.
We will use scaling arguments to identify the regions of phase space that are present at leading power and then take the limit where each region becomes a separate sector, that no longer interacts with the other regions. For a jet to have e 1, all particles must be either soft or collinear to the jet axis. In particular, a particle with energy E = zE J at an angle θ from the jet axis must satisfy This is a line in the log(1/z)-log(1/θ) plane, as shown in Fig. 2. Anything below the dashed line in this figure is too hard to be consistent with a given value of e (α) 2 . The soft drop criterion is that This is the region below the solid line in Fig. 2.
To find the relevant modes for the factorized expression, we need to identify the distinct characteristic momentum scalings that approach the singular regions of phase space in the limit e (α) 2 z cut 1. For a particular scaling, the constraints in Eqs. (3.1) and (3.2) will either remain relevant or decouple. We can characterize the relevant regions by their scalings in light-cone coordinates. Defining n µ as the jet direction andn µ as the direction backwards to the jet, then light-cone coordinates are triplets p = (p − , p + , p ⊥ ) where p − =n · p, p + = n · p and p ⊥ are the components transverse to n. On-shell massless particles have p + p − = p 2 ⊥ . The energy fraction is z = p 0 /Q = 1 2 (p + + p − )/Q and the angle to the jet axis in the collinear limit is θ = p ⊥ /p 0 .
We start with the soft modes, emitted at large angles θ ∼ 1, but still within the jet. If such radiation were to pass soft drop, with energy fraction greater than z cut , it would set e (α) 2 z cut ; this contradicts our assumed hierarchy e (α) 2 z cut . Therefore, soft wide-angle radiation is removed by soft drop and is not constrained by e (α) 2 . These modes thus have momenta that scale like p s ∼ z cut Q(1, 1, 1) .
They contribute only to the normalization of the distribution, not to its shape. Next consider the collinear radiation, emitted at small angles θ 1. All collinear radiation has p − p + . Then, from Eq. (3.1), we find Collinear modes can either have z ∼ 1 or be parameterically soft z 1. For modes with z ∼ 1, we have z z cut . Thus p − ∼ Q and p + ∼ Q(e (α) 2 ) 2/α independent of z cut . Their scaling is We call these modes collinear modes, although strictly they are not-soft collinear modes.
Collinear radiation that can have z ∼ z cut 1 we call collinear-soft. In this case, p − ∼ zQ and p + ∼ θ 2 zQ. These modes are simultaneously compatible with Eqs. (3.1) and (3.2). Their scaling is determined by saturating these parametric relationships, which leads to This is the point in phase space labeled S C (z cut e (α) 2 ) in Fig. 2.

Factorization and refactorization
With the relevant scalings identified, we proceed to derive the factorization formula. For simplicity, we focus on the case of e + e − → hemisphere jets, with e (α) 2 measured on each hemisphere. Jets at hadron colliders can be treated similarly, as we discuss in Sec. 4. Fig. 3 illustrates the relevant modes and their scales.
We begin with the usual SCET factorization formula, in the absence of soft drop grooming. The hard, collinear and soft modes are separated in the limit of small observables. This leads to [39,52,53] for the ungroomed hemispheres in e + e − → dijets events, provided e S G (z cut ) Figure 3: Illustration of the multi-stage matching procedure to derive the soft drop factorization theorem. As discussed in the text, we first match QCD to SCET, then factorize the jet function into collinear and collinear-soft modes. Canonical scales of all modes in the factorization theorem are shown on the right, ordered in virtuality where we assume that α > 1 and β ≥ 0.
full QCD to get the hard function, then decouple the soft and collinear degrees of freedom to pull the jet and soft functions apart [15][16][17][18]. Alternatively, one can use the method of regions approach [54,55], or the on-shell phase space approach [56][57][58]. Importantly, e (α) 2 is insensitive to recoil effects from soft emissions that displace the jet axis from the direction of hard, collinear particles [20,40], and so the jet and soft functions are completely decoupled.
Next we write down the hard-soft-jet factorization formula in the presence of soft drop grooming, assuming the hierarchy e (α) 2 z cut 1. With this assumption, soft radiation emitted at large angles must necessarily fail the soft drop criterion. Thus, all wide angle soft radiation in the jets (in this case, the hemisphere jets) is groomed and cannot contribute to the observable. All that remains of the global soft function is a z cut -dependent normalization factor S G (z cut ). This leads to (3.8) S G (z cut ) gives the cross section for the radiation from a set of Wilson lines that fails the soft drop criterion. An explicit calculation of S G for hemisphere jets at one-loop is given in Appendix C. With the collinear and soft modes decoupled, we can lower the virtuality of the collinear modes without further matching. The jet function J ze still depends on multiple scales, so to resum all the large logarithms it must be re-factorized. To see that it refactorizes, note first that in addition to being collinear, radiation in the jet function that is sensitive to the scale set by z cut must also be soft, by the assumption that z cut 1. Equivalently, emissions with order-1 energy fractions are not constrained by the scale z cut . We can thus factorize the jet function into two pieces depending on their energy fraction: Here, J(e 2 ) is the jet function that only depends on e (α) 2 and only receives contributions from emissions with order-1 energy fraction. S C (z cut e (α) 2 ) is the soft limit of the unfactorized jet function J ze (z cut , e (α) 2 ). The scaling of the collinear and collinear-soft modes are given in Eqs. (3.5) and (3.6) as discussed above. Note that, importantly, because the collinear-soft mode arises from refactorization of a jet function, it is a color singlet and only depends on two back-to-back directions. Because the jet function only depends on e (α) 2 , it is sensitive to a single infrared scale.
The step in Eq. (3.9) is the most unusual and important in the derivation. That the collinear-soft function depends on only a single combination of z cut and e (α) 2 is absolutely critical to being able to resum all the logs of e The pieces of the factorization theorem are: • H(Q 2 ) is the hard function for production of an e + e − → dijets event.
• S G (z cut ) is the global soft function. It integrates the radiation coming from Wilson lines in the jet directions that fails the soft drop criterion. Its modes fail soft drop and have momenta that scale as determined in Eq. (3.3).
2 ) is the jet function describing the emission of collinear radiation from a jet. Its modes parametrically pass soft drop and have momenta that scale as determined in Eq. (3.5).
is the collinear-soft function describing the emission of soft radiation boosted along the direction of a jet. Its modes may or may not pass soft drop and have momenta that scale as determined in Eq. (3.6). We denote the single scale that the collinear-soft function depends on as z cut e We present the operator definitions and explicit one-loop results for all of these functions in the appendices.
The appearance of collinear-soft modes in this factorization theorem has some similarities and differences with respect to the identification of other collinear-soft modes in the literature [23][24][25][26][27][28]. The original construction of a collinear-soft mode in Ref. [23] followed from boosting two jets in an event far from their center-of-mass frame in an effective theory of collinear dijets called SCET + . The collinear-soft mode in SCET + is sensitive to three Wilson line directions: the two from the collinear jets and the backward direction from boosting all other jets in the event. This collinear-soft mode was also exploited in Ref. [26] in the resummation of jet observables that are sensitive to multi-prong substructure.
The collinear-soft mode in the factorization theorem presented here, however, is more similar to modes identified from the measurement of multiple observables on a jet, each of which is only sensitive to radiation about a single hard core [24][25][26][27][28]. For example, Ref. [24] presented a factorization theorem for jets on which two angularities [37][38][39][40] are measured. At leading power, angularities are only sensitive to the hard jet core, and so the collinear-soft modes only know about two Wilson line directions: the jet axis and the backward direction. More recently, collinear-soft modes of this type have been used to resum NGLs [25,27] and logarithms of the jet radius [28].

The Single Scale of the Collinear-Soft Function
To demonstrate explicity that the collinear-soft function only depends on a single scale, we can make the following scaling argument. The collinear-soft function has the following form: Here, n is the number of final state collinear-soft particles, dΠ n is on-shell Lorentz-invariant phase space in d = 4 − 2 dimensions: µ is the renormalization scale, and M n is the amplitude for the production of the final state. Θ SD represents the soft drop grooming algorithm, which applies constraints on the final state and δ e (α) 2 represents the measurement of e (α) 2 on the final state: where the sum runs over the set of final state particles {i} that remain in the jet after grooming. To write this expression, we have used the definition of e (α) 2 from Sec. 2.2 and expanded in the collinear-soft limit, as in Eq. (3.4). Now, we rescale the momenta in light-cone coordinates that appear in the phase space integral in the following way: 14) At leading power in exactly d = 4, the phase space measure dΠ n and the squared matrix element |M n | 2 scale exactly inversely. Therefore, in d dimensions, under this rescaling, we have Next, look at how the measurement functions Θ SD and δ e (α) 2 change under the rescaling of Eq. (3.14). First, consider the soft drop groomer Θ SD . This consists of two parts: one, the reclustering with the Cambridge/Aachen algorithm and the second, the energy requirement on the clustered particles. The clustering metric of the Cambridge/Aachen algorithm is just the pairwise angle d This is invariant under the rescalings of Eq. (3.14). If in the clustering history we compare the angle between a collinear-soft particle i to the jet axis and the angle between two collinear-soft particles j and k, we have the constraint which is also invariant under the rescalings of Eq. (3.14). Finally, we can compare the angle between a pair of collinear-soft particles i and j to the angle between another pair of collinearsoft functions k and l, this leads to This too is invariant under Eq. (3.14). Therefore, for all possible clustering structures, the Cambridge/Aachen algorithm is invariant under the rescalings of Eq. (3.14).
The soft drop energy requirement on any number of particles that have been reclustered takes the form: where z i is the energy fraction of particle i and θ is the angle that the cluster of particles {i} makes with the jet axis. In terms of light-cone coordinates, this can be written as: Applying the rescalings of Eq. (3.14), this constraint becomes Note that the low scale z cut has been removed from this constraint. Under the rescaling, the measurement constraint δ e (α) Therefore, the low scale e (α) 2 has been removed from this constraint. Putting this all together, the collinear-soft function can be rewritten as (3.26) We have used the notation that Θ zcut=1 SD is the soft drop grooming algorithm with z cut = 1 and δ e (α) 2 =1 is the measurement with e (α) 2 = 1. All low scales have been explicitly removed from the phase space integral. This function is now seen to be a function only of the single scale "z cut e (α) The quotes just mean that the left-hand side is our abbreviation for the unwieldly quantity on the right-hand side. This proves that, to all orders, the collinear-soft function has dependence on only a single infrared scale, defined by this combination of z cut and e

Consequences of Factorization Theorem
Before we use the factorization theorem of Eq. (3.10) to make predictions for the cross section, we discuss consequences of this formula in some detail. Because the factorization theorem was derived without respect to any fixed order, these results hold to all orders.
Many of these consequences follow from the fact that soft wide angle radiation does not contribute to the shape of the soft-drop groomed e (α) 2 distribution for e (α) 2 z cut 1, a property that persists even for jets at hadron colliders. For example, it follows immediately from this fact that the shape of such a distribution is insensitive to contamination from pile-up and underlying event.
In this section, we will furthermore prove that at leading power, there are no NGLs that affect the shape of the soft-drop groomed e (α) 2 distribution in this regime. This was explicitly shown at O(α 2 s ) in Refs. [6,8,12] and plausibility arguments were presented for all orders, but this is the first proof. The factorization theorem also exhibits sample independence to a large degree, because the shape of the distribution is only sensitive to collinear physics. We will also demonstrate that soft-drop groomed energy correlation functions are less sensitive to hadronization than their ungroomed counterparts.

Absence of Non-Global Logarithms to All Orders
NGLs in cross sections of observables measured on individual hemispheres in e + e − collisions arise from a parametric separation of the scales in the hemispheres. Their leading effects are exclusively non-Abelian and quantify the correlation between the two hemispheres. Clearly, for a correlation to be present, there must be correlated radiation emitted into both hemispheres. If we measure the energy correlation functions e 1, then the radiation in the event must be soft wide-angle or collinear. At leading power, it is not possible to have correlations between different collinear directions (beyond total momentum conservation) as this would violate the collinear factorization of gauge theory amplitudes. Therefore, correlations and NGLs can only arise from soft, wide-angle radiation in the event with these assumptions.
The factorization theorem for ungroomed hemisphere energy correlation functions is [39,52,53] where "ug" denotes ungroomed. The cross section explicitly depends on soft wide-angle radiation through the soft function S(e 2,L , NGLs will be present in this factorization theorem. Because the soft function depends on two scales, all of the singular dependence cannot be determined by renormalization group invariance. More generally, non-global structure present in the soft function has been studied at O(α 2 s ) [59][60][61] and beyond [62,63] and recently, methods have been developed to control all-orders behavior [25,27,64,65]. However, NGLs represent an obstruction to resummation of the cross section to NLL and beyond. For groomed hemisphere energy correlation functions, our factorization theorem instead takes the form of Eq. (3.10): All soft, wide angle radiation throughout the event is described by S G (z cut ), which is sensitive only to the single scale z cut . Therefore, there are no NGLs present in this factorization theorem. Even with a hierarchy between e (α) 2,L and e (α) 2,R , these observables are completely decoupled at leading power in e (α) 2 and z cut . Additionally, the shape of the distribution is also independent of jet radius effects and the precise way in which the hemispheres are defined.
When we discuss soft-drop groomed jets in pp collisions in Sec. 7, we will place no constraint on global soft radiation throughout the event, unlike the case of e + e − → hemisphere jets. Nevertheless, the shape of the soft-drop groomed e (α) 2 distribution will still have no NGLs, jet radius effects, etc., due to universality of the collinear limit of QCD amplitudes. The normalization, however, will in general be sensitive to scales both in the jet (set by z cut and the jet radius R) and scales outside of the jet (set by the partonic collision energy). To eliminate these effects in pp collisions, we can normalize the cross section, say, to integrate to unity.

Process Independence
Strictly speaking, the factorization theorem of Eq. (3.10) depends on the process. It includes the hard function, which is process dependent, and a soft function, that knows about all hard jet directions. Nevertheless, there is a sense in which the factorization theorem is process independent. Normalizing the cross section completely removes the hard and soft function dependence. Then, by the universal collinear factorization of QCD amplitudes, if we are completely inclusive over the right hemisphere, then the differential cross section of the softdrop groomed energy correlation function in the left hemisphere is given by where we assume that e (α) 2,L z cut 1 and N is some normalization factor. That is, in the deep infrared where e (α) 2,L z cut 1, all radiation in the groomed jet is constrained to be collinear. Therefore, in this limit and for a fixed jet energy, the shape of the distribution for quark jets is independent of the process that created the quark jets, due to the universality of QCD matrix elements in the collinear limit.
This collinear factorization property of soft-drop groomed observables can be exploited for jets in pp collisions. Unlike the dominant case in e + e − collisions, jets at a pp collider can be either quark or gluon. Of course, on a jet-by-jet level, we cannot determine whether a jet was initiated by a quark or gluon. However, for a given process, we can determine the relative fraction of quark and gluon jets in the sample. For jets produced at a pp collider, the process independence manifests itself in the cross section as where D q (D g ) is proportional to the fraction of quark (gluon) jets in the sample. The relative fraction of quark and gluon jets can be determined from fixed-order calculations, using a simple algorithm for determining the flavor of a groomed jet. We will describe this in detail in Sec. 7 when we match our resummed distribution to fixed order in pp → Z + j events.

Hadronization Corrections
With a factorization formula, one can estimate the size and importance of non-perturbative corrections to the cross section. We will only consider non-perturbative corrections to the shape, as the normalization can be set by hand. Therefore, non-perturbative corrections can only enter into our factorization theorem via the jet or collinear-soft functions. From Eqs. (3.5) and (3.6), the scales appearing in the jet and collinear-soft functions are If either of the scales approaches Λ QCD , then we expect there to be large corrections to the perturbative cross section due to non-perturbative physics. We can estimate when nonperturbative corrections become large by setting these scales to be Λ QCD . For α > 1 and β ≥ 0, the collinear-soft mode has a lower virtuality than the collinear mode, so it will probe the non-perturbative region of phase space first. The value of e (α) 2 at which the collinear-soft mode becomes non-perturbative is This estimate can be compared with the Monte Carlo analysis of hadronization corrections to the soft-drop groomed energy correlation functions from Ref. [8]. In particular, the estimate of Eq. (4.7) of when non-perturbative corrections become important for α = 2 as a function of β agrees exceptionally well with Fig. 10(a) of Ref. [8].
For β < ∞, the soft drop groomer reduces the effect of non-perturbative corrections with respect to the ungroomed observable. This can be simply seen from Eq. (4.7), in which the prefactor Table 1: α s -order of ingredients needed for resummation to the accuracy given. Γ cusp is the cusp anomalous dimension, γ is the non-cusp anomalous dimension, and β is the QCD β-function.F (∂ ω ) are the logarithms in the low-scale matrix elements that have been Laplace transformed and c F are constants in the low-scale matrix elements. The final column shows the relative order to which the resummed cross section can be matched to fixed-order.
approaches unity as the grooming is removed (β → ∞). This factor is less than 1 for β < ∞, provided α > 1 and Λ QCD < z cut Q. For high energy jets, this suppression can be substantial. For example, for α = 2 (corresponding to jet mass) and β = 0 (corresponding to mMDT groomer) non-perturbative effects become important at This agrees with the estimate of the size of nonperturbative corrections for the mMDT groomer from Ref. [6].

Achieving NNLL Accuracy
In this section, we determine the anomalous dimensions necessary to resum the large logarithms of soft-dropped energy correlation functions through NNLL accuracy. The practical details of how one assembles these ingredients, in the framework of SCET, to construct a resummed cross section are given in App. F. We will discuss matching to fixed-order and demonstrate our ability to make phenomenological predictions in subsequent sections. Resummation in SCET is accomplished with renormalization group evolution. Solving the renormalization group equations to a given logarithmic accuracy requires anomalous dimensions to a particular fixed order. The anomalous dimensions of the functions in the factorization theorem must sum to zero, because the cross section is independent of the renormalization scale.
Recall that, for e + e − → hemisphere jets, the factorization theorem for soft-drop groomed energy correlation functions is (5.1) Table 1 presents the order to which anomalous dimensions and constants of the functions in this factorization theorem must be computed for particular logarithmic accuracy (see, e.g., Ref. [66]). The cusp anomalous dimension Γ cusp and the QCD β-function are known through three-loop order [67][68][69][70][71][72] and we present them in App. A. The hard function H(Q 2 ) for e + e − → qq is known to high orders and its non-cusp anomalous dimension γ H is known at three-loop order [73,74]; we present the relevant pieces in App. B. For arbitrary angular exponents α and β, little else in the factorization theorem is known at sufficiently high accuracy to resum to NNLL. The goal of this section is to fill in the rest of the table, to achieve full NNLL accuracy. We start in Sec. 5.1 restricting to α = 2 (jet mass) and β = 0 (mMDT groomer). For this case, all of the missing ingredients can be determined by recycling results from the literature, up to calculable clustering effects from the soft drop algorithm. In Sec. 5.3, we consider α = 2 and β ≥ 0 and demonstrate that one can extract unknown two-loop non-cusp anomalous dimensions with EVENT2. It is possible to extend our analysis to angular exponents for the energy correlation functions beyond α = 2, but we do not do it in this paper. 3

NNLL for
We first consider angular exponents α = 2 and β = 0. In this case, the soft drop requirement enforced at every branching reduces to an energy cut On the soft-drop groomed jets we then measure where the subscript g denotes that the mass and energy are measured on the groomed jet. The jet functions in the factorization theorem are independent of the soft drop groomer, so we are able to use results from the literature for these. The inclusive jet function has been calculated to two loops [75][76][77][78] and the non-cusp anomalous dimension of the inclusive jet function is known to three loops [79,80]. We present the relevant expressions in App. D. This leaves the soft function S G (z cut ) and the collinear-soft function S C (z cut e 2 ) to be determined. Their one-loop expressions are easily calculable, and we present the results in App. C and App. E. To determine their two-loop non-cusp anomalous dimensions, we exploit the renormalization group consistency of the factorization theorem. The sum of the anomalous dimensions must vanish at each order: where γ F denotes the anomalous dimension of function F in the factorization theorem, and we have used the symmetry of the left and right hemispheres of the event. Therefore, only one unknown anomalous dimension remains, which we take to be γ S .

Two-Loop Soft Function
To calculate the two-loop non-cusp anomalous dimension γ S we need to calculate the soft function S G (z cut ) with two real emissions. The two-loop expression for the soft function is and |M(k 1 , k 2 )| 2 is the squared matrix element for two soft emissions from a qq dipole. The explicit expression for |M(k 1 , k 2 )| 2 can be found in Ref. [81]. Θ SD is the phase space constraint imposed by the soft drop groomer. Recall that, for consistency with the assumed hierarchy e (α) 2 z cut , soft modes must fail soft drop. If the particles in the hemispheres are reclustered using the Cambridge/Aachen algorithm, Θ SD can be written as The first line of Eq. (5.7) corresponds to particles 1 and 2 lying in different hemispheres (opposite rapidity with respect to the qq dipole), and so each particle individually must fail soft drop. Q is the center of mass energy and so Q/2 is the energy in one hemisphere. The second and third lines correspond to the configuration where both particles lie in the same hemisphere. θ 12 is the angle between the particles and θ iJ (for i = 1, 2) is the angle particle i makes with that hemisphere's axis. If θ 12 is less than both θ 1J and θ 2J then, according to the Cambridge/Aachen algorithm, the soft particles are clustered first. Therefore, the sum of the energies of particles 1 and 2 must fail soft drop. If instead one of the particles is closer to the jet axis, then they are clustered separately and must individually fail soft drop.
To proceed, we separate the squared matrix element into Abelian and non-Abelian pieces, according to their color coefficient. At this order, the squared matrix element takes the form Here, "n-A" denotes the non-Abelian component of the squared matrix element, which includes the C F C A and C F n f T R color channels. The Abelian contribution is just the symmetrized product of the one-loop result, with a color factor of C 2 F . We will consider these two pieces separately, starting with the non-Abelian term.

Non-Abelian Clustering Effects
Note that except for the effects from Cambridge/Aachen clustering, soft drop is just imposing a soft energy veto on each hemisphere. The two-loop soft function with a soft energy veto was calculated in Ref. [29]. That calculation showed that the two-loop Abelian piece (proportional to C 2 F ) to the energy vetoed soft function satisfies non-Abelian exponentiation. The two-loop non-cusp anomalous dimension for a hemisphere energy vetoed soft function is then purely non-Abelian and was extracted in Ref. [28]. The non-Abelian part of the soft function with an energy veto at two-loops is The phase space cut Θ veto is where Λ is the veto scale. We can then write the two-loop soft function for soft drop as where the veto scale is set to Λ = z cut Q/2. The difference between the soft drop and energy veto phase space constraints is purely a clustering effect, given by Eq. (5.10) enables us to calculate much more simply the two-loop non-cusp anomalous dimension of the soft function. The anomalous dimension can then be written as γ veto is the two-loop non-cusp anomalous dimension of S veto extracted in Ref. [28]: where β 0 is the one-loop β-function coefficient: Then, we only need to determine the contribution to the anomalous dimension from residual Cambridge/Aachen clustering effects, γ C/A . The non-Abelian clustering effects are contained in The squared non-Abelian matrix element does not have collinear singularities when the angle of the particles from the jet axis is strongly ordered. Therefore, in this integral there is only a collinear divergence when the two emissions become collinear to the jet axis in a non-strongly ordered way. The coefficient of this divergence is proportional to the correction to the twoloop anomalous dimension due to clustering effects in the non-Abelian color channel. The divergence can be extracted with the standard plus-function prescription and the correction to the anomalous dimension can be found. While we were unable to find an analytic expression, its approximate numerical value is 4 The contribution to the anomalous dimension is then

Abelian Clustering Effects
The Abelian contribution can be calculated similarly. However, unlike the non-Abelian contribution, the exponentiation of the one-loop result will describe at least some of the two-loop Abelian piece. If the square of the one-loop result does not account for all of the two-loop result, then non-Abelian exponentiation breaks down. This does not mean that exponentiation breaks down or that the cross section cannot be resummed, just that the anomalous dimension of the purely Abelian piece will need to be corrected at every logarithmic order. So, for the two-loop non-cusp anomalous dimension, we need to determine the part of the soft function that is not accounted for by non-Abelian exponentiation.
To do this, we start from the full expression for the Abelian term at two-loops: We then add and subtract the one-loop phase space constraints: The difference between the phase space constraints is a clustering effect, given by As with the non-Abelian term, this phase space constraint completely removes all soft divergences and the strongly-ordered collinear limit. The remaining divergence can be isolated by standard plus-function techniques. For the two-loop Abelian Cambridge/Aachen clustering term, we find the numerical result for the second integral in Eq. (5.19). The contribution to the anomalous dimension is then

Two-Loop Anomalous Dimension and Comparison with EVENT2
Combining Eqs. (5.13), (5.17) and (5.22), the total two-loop non-cusp anomalous dimension for the soft function is (5.23) The two-loop non-cusp anomalous dimension for the collinear-soft function is found by consistency using Eq. (5.4). Note that this anomalous dimension has no log z cut terms. Therefore, the anomalous dimensions of no functions in the factorization theorem have log z cut dependence. This is a consequence of the fact that each function of our factorization theorem in Eq. (3.10) depends on a single infrared scale, allowing NNLL resummation of all logarithms of z cut and e We can verify this result by comparing the resummed distribution, truncated at O(α 2 s ), with the singular region of the full QCD result, computed to the same fixed order. For the full QCD result, we have implemented soft drop into EVENT2 [35], a Monte Carlo code that generates fixed-order results up to O(α 2 s ) in e + e − collisions. Our specific implementation is as follows. We generate e + e − collisions at 1 TeV center of mass energy and identify event hemispheres with the exclusive k T algorithm [45]. We then recluster each hemisphere using the Cambridge/Aachen algorithm and apply soft drop with β = 0. On each of the soft-drop groomed hemispheres, we then measure the energy correlation function e (2) 2 and record the larger of the two values, which we denote by e (2) 2,H and refer to as the heavy groomed mass. This is simply related to the cross section of our factorization theorem: In Fig. 4, we compare EVENT2 results to the prediction of the factorized expression at NNLL expanded to O(α 2 s ). For soft drop with β = 0, soft logarithms are removed, which means that at O(α 2 s ), the cross section has the schematic form where C 0 and C 1 are constants. We plot the cross section separated into the three color channels (C 2 F , C F C A , and C F n f T R ). We set z cut = 0.001 to suppress power corrections of z cut . Excellent agreement between our factorization theorem and EVENT2 is observed at small e (2) 2 , demonstrating that we have captured all singular terms of the full QCD result in our factorization theorem to O(α 2 s ).

Reclustering with anti-k T
It is illuminating to study the clustering effects in the soft function in more detail. In this section, we re-calculate the clustering effects with the anti-k T algorithm, instead of the standard Cambridge/Aachen algorithm. We find that the clustering effects with the anti-k T algorithm are intimately related to the corresponding effects calculated in jet veto calculations. This can be understood relatively simply by re-expressing the clustering conditions in a form analogous to the clustering metric of the longitudinally-invariant k T algorithm.
To calculate the two-loop soft function for soft drop defined with anti-k T reclustering, we only need to calculate the clustering effects unique to this algorithm. We will denote the phase space constraints for the anti-k T reclustering as Θ ak T SD , but we will not explicitly present them here. The two-loop soft function is The relevant phase space constraints can be written as With this, we can calculate the divergent part of the two-loop soft function from clustering effects and extract the anomalous dimension. As with Cambridge/Aachen, we can write the two-loop non-cusp anomalous dimension as where γ ak T is the part of the anomalous dimension purely from clustering effects. We find This anomalous dimension is fascinating. First, note that there is no C 2 F term, implying that non-Abelian exponentiation holds for anti-k T reclustering, in contrast to what we found for the Cambridge/Aachen algorithm. That is, all logarithms at O(α 2 s ) with color factor C 2 F are accounted for by exponentiating the one-loop result. This is to be expected: Since the anti-k T algorithm clusters soft gluons (with energy fractions of order z cut ) one-by-one with the hard jet core unless two soft gluons have angular separation ∆R z cut , clustering effects are merely a power correction for Abelian gluons. Also, unlike the case for Cambridge/Aachen reclustering, there is explicit log z cut dependence in the anomalous dimension of Eq. (5.29). This shows that we do not resum logarithms of z cut to full NNLL accuracy when anti-k T clustering is used in soft drop. The coefficient of the log z cut term is identical to the coefficient of the logarithm of the jet radius R found from clustering effects in jet veto calculations [30][31][32][33][34]. This connection between soft drop and jet veto calculations can be made clearer by a simple rewriting of the clustering metric.
The k T class of clustering metrics for e + e − collisions can be written as for particles i and j, with p an integer that defines the jet algorithm. In the soft function, soft particles are either clustered with each other or with the jet axis. For β = 0, these soft particles have characteristic energy fraction z cut 1. In terms of energy fractions, the clustering metric of two soft particles is and for a soft particle i with the jet axis it is where θ i is the angle between particle i and the jet axis. Consider p < 0. In this case, the two soft gluons are (parametrically) clustered together when 33) or equivalently, when The effective clustering metric in this case is then With p = −1, this is the clustering metric for the inclusive anti-k T algorithm with effective jet radius R = z cut 1. There will now be logarithms of the jet radius that arise. The log z cut term in the anomalous dimension has the identical coefficient as the log R term in jet veto calculations because z cut and R act as the angular scale for collinear splittings in the respective soft functions.
In summary, while we could use anti-k T to recluster the jet for soft drop grooming, we could not resum all large logarithms to the same precision without a different factorization theorem. Therefore, reclustering in soft drop with the Cambridge/Aachen algorithm is preferred from a theory perspective.

NNLL for
For soft drop with angular exponent β > 0, we cannot recycle results from the literature to reach NNLL precision. Instead, a completely new two-loop calculation of either the soft or collinear-soft function is needed. But without such a calculation, we can perform NNLL resummation for particular values of β > 0, using numerical simulations to estimate the ingredients we lack. We will demonstrate this explicitly in the case of β = 1, and the result will allow us to study features of NNLL distributions for energy correlation functions with less aggressive grooming.
The same method we used to validate anomalous dimensions for β = 0 can be used to extract the anomalous dimension for β > 0. This method relies on the fact that all ingredients necessary for NNLL resummation with α = 2, β > 0 are known except the two-loop noncusp anomalous dimensions of the soft and collinear-soft functions. As mentioned above, renormalization group invariance determines one of these, say γ (1) S C , in terms of the other anomalous dimensions. So only one unkown, γ (1) S , remains and we can extract it at fixed order.
To do this for a given β > 0, we can use EVENT2 to obtain numerical results at O(α 2 s ) for the groomed e (2) 2,H distribution with several moderately small values of z cut . From each of these distributions, we can subtract the known terms, which we get by expanding the NNLL distribution to fixed order. This leaves a term proportional to the unknown γ 2,H ). At O(α 2 s ) it is appropriate to use 0 ≤ n + m ≤ 3, though in practice we found terms with m ≥ 2 to be difficult to fit. With the non-negligible power corrections thus removed, we can then extract the remaining anomalous dimension.
While the procedure outlined above is straightforward, an explicit calculation of γ S or γ (1) S C for β > 0 is of course desirable. On practical time scales, numerical extractions are limited to rough approximations, due to inadequate numerical precision in the deep infrared. Nevertheless, an estimate is sufficient for our purposes here, which are to demonstrate the advantages of resumming jet substructure observables to NNLL, and to examine various levels of grooming. Thus, we will test the above procedure on β = 0, and learn about the associated uncertainties by comparing with our direct calculation, Eq. (5.23). Then we will move to β = 1, and extract γ (1) S in that case. In Fig. 5a we show numerical results at O(α 2 s ) from EVENT2 with β = 0 and z cut = 0.1. Also shown is the NNLL distribution, expanded to fixed order, but without the γ (1) S term. The discrepancy between the curves is thus due to the missing γ (1) S term and z cut power corrections. Using several distributions like this one, with values of z cut between 10 −4 and 10 −1 , we fit the z cut power corrections. Fig. 5b shows the remaining offsets between our analytical curves and the results of EVENT2, after z cut log n (z cut ) log(e (2) 2,H ) power corrections have been subtracted. On each point in this plot, the error bar represents the standard deviation in EVENT2 output, across e 2,H bins. The offset that remains as z cut → 0 is the γ (1) S we would extract using this method. One can see from Fig. 5b that agreement with our analytical calculation, Eq. (5.23), is quite good, with some discrepancy in the C F channel. Table 2 lists the numerical results of γ S using this method. 5 The uncertainties quoted for the β = 0 extraction in the table come from the standard deviation in EVENT2 output across e 2,H bins, which introduces an error in the identification of constant offsets. These should 5 In carrying out the procedure just described, we tuned EVENT2 parameters to favor the infrared. In particular, we use of order 1 trillion events, with CUTOFF = 10 −15 and phase-space sampling exponents NPOW1 = NPOW2 = 5. This procedure corresponded to centuries of CPU time.
be compared with our direct calculation in the second line of the table. The discrepancy in the C F channel gives us a sense of additional numerical uncertainties, which are significant. Similar disagreements have been encountered before, e.g. in Ref. [82], in the context of C F channel extractions from EVENT2, and to resolve it might require significantly longer run times.
As stated above, a rough estimate of γ S for β > 0 is sufficient for our purposes, so we have applied the method described above to the case of β = 1. See the third line of Table 2 for the results of the extraction. Uncertainties quoted in this line of the table have two sources: (i) variance in EVENT2 output, and (ii) additional numerical precision issues, which we took to be the difference (both absolute and relative) between extraction and direct calculation in the β = 0 test. In each color channel, we took the maximum of these uncertainties and inflated it by a factor of 2.
The estimate in Table 2 allows us to study NNLL distributions of e 2,H groomed with β = 1. In the resulting distributions, the uncertainties associated with the imperfect extraction are relatively small; e.g. see Fig. 6b below. Still, a direct calculation of either γ  S C for β > 0 would of course be preferred, but we leave this to future work.
6 Matching NNLL to Fixed Order in e + e − → dijets Using the results calculated in the previous sections, here we match our resummed differential cross section for soft-drop groomed energy correlation functions to fixed-order for hemisphere jets produced in e + e − collisions. We first match resummed results at NLL and NNLL to O(α s ) and O(α 2 s ), respectively, using EVENT2 and demonstrate that theoretical uncertainties are greatly reduced at NNLL. We then compare several Monte Carlo parton shower simulations to our matched NNLL results. We compare both parton and hadron level Monte Carlo to our perturbative analytic results, and include a simple model of hadronization in our calculation. We leave a detailed understanding and justification of incorporating hadronization into the resummed and matched cross section to future work.

Matching Resummation to Fixed-Order
With the explicitly calculated and extracted two-loop non-cusp anomalous dimensions of the soft function in the soft drop factorization theorem Eq. (3.10), we are able to resum the differential cross section through NNLL accuracy in the region where e z cut 1, and will not provide an accurate description of the cross section outside this region. To accurately describe the cross section throughout the full phase space requires matching the resummed result to fixed-order.
While there are many ways to do this at various levels of sophistication, we choose to use simple additive matching. That is, we construct matched distributions according to Here, dσ resum is the resummed cross section, calculated to the appropriate logarithmic accuracy. dσ FO is the fixed-order differential cross section calculated to a particular order in α s . dσ resum,FO is the resummed cross section truncated at the same accuracy as the fixed-order cross section. In the infrared phase space region, this term will exactly cancel the singularities in the fixed-order cross section, only leaving the resummed cross section plus power corrections. In the hard phase space region, this term cancels the resummed cross section, up to higher orders in α s . The logarithmic accuracy of the resummed cross section was defined in Sec. 5, and here we specify the fixed orders that we use in the matching procedure. We additively match the analytic NLL distributions to O(α s ) fixed order results, which include one real emission from the qq dipole. We match NNLL distributions to O(α 2 s ) results, which include up to two real emissions. EVENT2 is able to generate e + e − collisions through O(α 2 s ), except for the two-loop virtual contribution. The two-loop virtual term only contributes at e (2) 2 = 0, so our differential distributions are unaffected by this omission.
In Fig. 6, we plot the resummed and matched differential cross sections for the larger e (2) 2 of the two hemispheres at NLL and NNLL with various levels of soft drop grooming. Here, we consider dijet production in e + e − collisions at 1 TeV center-of-mass energy and identify hemispheres with the exclusive k T algorithm [45]. The parameters of soft drop are z cut = 0.1 and we show both β = 0 and β = 1. We also show the ungroomed heavy hemisphere e (2) 2 distribution. In these plots, we include estimates of theoretical uncertainties represented by the lighter bands about the central curve. While more sophisticated methods for estimating uncertainties exist, we simply vary the natural scales that appear in the functions of the factorization theorem up and down by a factor of two. We then take the envelope of these scale variations as an estimate of theoretical uncertainties. This simple prescription is sufficient for our main purpose in showing uncertainty bands: to demonstrate the reduction in theoretical uncertainty in moving from NLL to NNLL. Included in these uncertainty estimates is a variation in our treatment of the Landau pole of the strong coupling α s . For scales µ > 1 GeV, α s is evaluated according to its perturbative running. For µ < 1 GeV, we freeze α s to its value at the scale µ = 1 GeV. This is not intended to be a model for hadronization or non-perturbative physics, but is just intended to maintain finite cross section predictions at small e 2 values. To estimate the sensitivity of our results to the scale at which we freeze the coupling, we vary this 1 GeV scale by a factor of two, and include the effect in the uncertainty bands of Fig. 6 as well.
Finally, we have shown the uncertainty bands around the β = 1 curves at NNLL with and without the uncertainty in our estimate of the two-loop non-cusp anomalous dimension of the soft function. One can see from the figure that this imperfect extraction has only a relatively small effect on the overall uncertainty at this order.
Importantly, we allow the normalization of the cross section to change under these scale variations. That is, the curves in Fig. 6 are constructed according to Eq. (6.1). The normalization of each distribution displayed is meaningful, since we resum all large logs in both the shape and the normalization. While the central value curves don't change much in going from NLL to NNLL, the uncertainties are dramatically reduced, and this is partly due to the increased accuracy in the normalization.

Comparison to Monte Carlo
In this section, we compare our NNLL resummed and matched soft-drop groomed e (2) 2 distributions to the output of several standard Monte Carlo simulations. We generate e + e − → dijets events at 1 TeV center-of-mass collision energy with Herwig++ 2.7.1 [83,84], Pythia 8.210 [85,86], and Vincia 1.2.02 [87][88][89][90]. While the Herwig++ and Pythia events are showered from the leading order process e + e − → qq, we consider Vincia with and without fixed-order matching included. The matched Vincia results are accurate effectively through O(α 2 s ). For the most direct comparison of the simulations to our NNLL matched results, we run α s at two loops in all Monte Carlos (in the CMW scheme [91,92]) and we fix α s (m Z ) = 0.118, which is the same value used in our analytic calculations. We include Monte Carlo events both at partonic level and after hadronization. These events are then clustered into hemispheres using the exclusive k T algorithm [45] using FastJet 3.1.3 [93]. The soft drop grooming and subsequent measurement of the two-point energy correlation functions of these e + e − events is implemented in FastJet with custom code.
We compare the Monte Carlo distributions to our NNLL resummed and matched calculations in Figs. 7 and 8. In these plots, all distributions integrate to the same value over the range e   The uncertainty bands for the analytic curve includes the uncertainty in the two-loop noncusp anomalous dimension.
As a more direct comparison of the Monte Carlos, Fig. 9 displays the relative difference between each of the hadron-level Monte Carlos and our matched NNLL predictions. Again, soft drop is performed with z cut = 0.1, and both β = 0 and β = 1 are shown. All the Monte Carlo curves lie within our shaded band of theoretical uncertainty, but discrepancies between the different simulations are visible.
One striking feature in these plots, especially for β = 0, is the presence of additional structure in the hadron-level Monte Carlo distributions at small e (2) 2 . It is clear that this feature is due to non-perturbative physics, and so is therefore not included in our NNLL calculation. Nevertheless, we can include a simple model of hadronization into our calculation to see if this structure is easily explained.
For additive IRC safe observables, like thrust or jet mass, it can be shown from general principles that hadronization corrections can be incorporated in perturbative distributions by convolution with a model shape function [94,95]. In general, the energy correlation functions are additive observables, so we should be able to use shape functions to model hadronization corrections. However, once soft drop is applied on the jet, emissions in the jet may or may not contribute to the energy correlation functions, so the observable is no longer strictly additive. We leave a more careful study of whether shape functions can be used to model hadronization effects in groomed observables to future work. Here we convolve our matched results with a simple shape function to see if qualitative agreement with the Monte Carlos can be achieved.  Figure 10: Perturbative NNLL results for soft-drop groomed e 2 with z cut = 0.1 and β = 0 (left) and β = 1 (right), compared to analytic results that include the shape function of Eq. (6.2) for modeling hadronization, and compared to hadron-level Monte Carlo. The parameter Ω = 1 GeV. Note that, qualitatively, the shape function produces a hadronizationbump similar to those seen in the Monte Carlos.
Because shape functions describe non-perturbative physics, they only have support for energies comparable to Λ QCD . The shape function we choose is the parametrization suggested by Ref. [96]: and has first moment equal to Ω. As discussed in Sec. 4.3, of all modes present in our factorization theorem, the collinear-soft mode has the lowest virtuality, so it will have the largest sensitivity to non-perturbative physics. We thus convolve our perturbative distribution with the shape function, assuming non-perturbative effects are primarily associated with the collinear-soft mode. That is, we include hadronization corrections in the soft drop groomed e (α) 2 distribution according to where the argument of the perturbative distribution is shifted by the virtuality of the collinearsoft mode, Eq. (4.7). In Fig. 10, we compare the matched NNLL distribution of e 2 with and without convolution with the shape function of Eq. (6.2), in which we set Ω = 1 GeV to be comparable to the scale of hadron masses. We show this comparison for soft drop grooming with β = 0 and β = 1. The peak at small e (2) 2 for β = 0 agrees qualitatively with the structure of the hadronized Monte Carlo distributions. Similarly, the shape at small e (2) 2 for β = 1 agrees with the simulations as well. This suggests that there might exist a shape function for describing hadronization effects in groomed jet observables, though we leave a detailed discussion and justification for such a model to future work.
7 Matching NNLL to Fixed Order in pp → Z + j In this section, we present predictions for soft-drop groomed e (α) 2 distributions as measured on the jet in pp → Z + j events at the LHC. The definitions of soft drop and energy correlation functions appropriate for jets in pp collisions are given in Sec. 2. As with jets from e + e − collisions, we match our NNLL resummed distribution to fixed-order results that include relative O(α 2 s ) corrections to the Born process. There are two complications we must deal with. First, at pp collisions, the jets will be both quark and gluon initiated. Second, because we only measure the observable within the jet and do not constrain radiation throughout the rest of the event, the simple hard-soft-jet factorization that we employed for e + e − → hemisphere jets will not apply here. Nevertheless, while the normalization of the jet-observable distribution will thus be complicated and sensitive to multiple scales, the shape of the distribution will still be controlled exclusively by collinear physics. To address these complications, we first show how soft-drop groomed quark and gluon jets can be unambiguously defined order-by-order in perturbation theory. Then we discuss how the normalization of the distribution can be obtained by matching to full QCD at fixed order. The discussion will focus on the Z + j sample for concreteness, but these ideas apply equally well to any process with hard jets at a hadron collider.

Resummed Cross Section in pp → Z + j
We define our observable on soft-drop groomed jets in pp → Z + j events in the following way. First, we cluster the final state according to a jet algorithm with some jet radius R ∼ 1. Of the jets with pseudorapidity |η J | < η max , we then identify the jet with the largest transverse momentum p T J and require that p T J > p min T . We groom this jet with soft drop and measure e (α) 2 according to the definitions given in Sec. 2 for jets in pp collisions. In this procedure, we remain inclusive over all other hadronic activity in the final state: we only care about the hardest jet.
For this process, the relevant factorization formula is 2 ) . z cut 1. There will be logarithms of z cut (and other scales in the events) in the D k prefactor that we do not resum. When referring to calculations of this cross section, we will specify the accuracy to which logarithms of e (α) 2 are resummed (i.e., NLL or NNLL). We now explain the components of this formula in detail.
In Eq. (7.1), S C,k (z cut e D k is a matching coefficient that can be extracted from fixed-order calculations, and it sets the normalization and relative contributions from the different jet flavors. In addition to the dependence explicitly shown, D k also depends implicitly on parton distributions, as different initial states produce different flavors of final state jets. Unlike the case in e + e − collisions, where the jet energy was (almost exactly) half the center-of-mass energy, due to the non-trivial parton distributions, the distribution of the jet p T has a finite width and depends on the cut, p min T . For a true precision prediction, we would compute the matching coefficient D k as a function of p T J and include an integral in Eq. (7.1) convolving the jet and collinear-soft functions with D k (p T J ). An approach to doing this in a semi-automatic manner was discussed recently in Refs. [97,98]. But, for simplicity we instead employ the following approximation: we evaluate the jet and collinear-soft functions at p T J , the average p T J .
The average jet transverse momentum p T J can be estimated by using the fact that the cross section for a jet with transverse momentum p T J takes the power-law form: This distribution is normalized and the mean value of p T J is The typical exponent is n ∼ 5, and we take n = 5 in our numerical computations. The full cross section for soft-dropped e (α) 2 (including power corrections) can be expressed as dσ Here, the right-most term includes all power corrections suppressed by e (α) 2 or z cut . The functions S C,k and J k should be evaluated at p T J but we have suppressed their arguments for brevity. We will use this form of the cross section to define the matching coefficient D k at fixed-order. For NNLL resummation, the relative O(α s ) corrections to D k are required.
First, at leading order in α s , Eq. (7.4) becomes where the superscript (0) denotes the leading order in α s . Here, we have used J 2 ). Also, since a jet has only one constituent at this order, the distribution has no support away from e At the next-to-leading order in α s , the extraction of D k requires separating the jets by flavor. Since D k is defined in each flavor channel, we need to determine the flavor of the hardest jet in each pp → Z + j event included in our sample. Ordinarily, any definition of jet flavor based on the constituents of the jet is infrared-unsafe and ill-defined at leading power, because soft wide-angle emissions into a jet can change its flavor. 6 Soft drop eliminates this problem at leading power in e (α) 2 and z cut by removing soft wide-angle radiation from the jet. This allows for an infrared and collinear safe definition of jet flavor at leading power in e (α) 2 and z cut . We define the jet flavor f J as the flavor sum of the constituents of the groomed jet: where f q = 1, fq = −1 and f g = 0. The subscript on J g means that one only sums over the jet constituents that remain after grooming with soft drop. If f J = ±1, then the jet is quark-type, while if f J = 0, it is gluon-type. With this jet flavor identification, we are able to determine the total fixed-order cross section for each jet flavor channel in pp → Z + j. We will denote the next-to-leading order term in the cross section for a jet of flavor k as σ (1) k , defined according to the phase space cuts described at the beginning of this section.
Then, at next-to-leading order in the k flavor channel, Eq. (7.4) becomes Here, S C,k and J We computed σ (1) k using MCFM [100,101] with settings detailed in the next section. We computed the power corrections according to (7.10) For the first term in the integrand, we use a numerical distribution obtained with MCFM.
Since we do not have access to this distribution at arbitrarily small values of e This completes our extraction of the matching coefficient D k through relative O(α s ). With it, the resummed cross section of Eq. (7.1) is complete and ready to be matched to relative O(α 2 s ) fixed-order results.

Matching Resummation to Fixed-Order
With the resummed differential cross section for soft-drop groomed e (α) 2 defined in Eq. (7.1), we next match to fixed order for pp → Z + j. Our matching procedure will be identical to the procedure we used for e + e − collisions; we add the difference between the exact fixed order and the expansion of the resummed distribution to fixed order: We match the analytic NLL resummed distributions to fixed-order results that include the relative O(α s ) corrections to the Born process for pp → Z + j. We match NNLL distributions to fixed-order results including relative O(α 2 s ) corrections and up to 3 partons in the jet. We use MCFM v. 6.8 [100,101] to generate the fixed-order cross sections for soft-drop groomed e (α) 2 in pp → Z + j events. Currently, MCFM can only generate fixed-order corrections at O(α s ) relative to a Born-level process, and so we will have to use some properties of the observable to be able to calculate to relative O(α 2 s ) accuracy. For e (α) 2 > 0, as we did in e + e − collisions, we can ignore the purely two-loop virtual contribution to pp → Z + j, as it has no effect on the differential distribution away from e (α) 2 = 0. MCFM can generate both inclusive pp → Z + j and pp → Z + 2j processes through relative O(α s ) accuracy. Therefore, we can use pp → Z + 2j at relative O(α s ) in MCFM to calculate the relative O(α 2 s ) distribution for pp → Z + j, in the region where e (α) 2 > 0. In practice, this procedure requires some care. To define the cross section for pp → Z +2j in MCFM, we must set a minimum p T for the two jets as identified by MCFM. This is set by the parameter ptjet min within MCFM. To compute the fixed-order cross section correctly for e (α) 2 as measured on the soft-drop groomed jet in pp → Z + j events, ptjet min should be set to 0; this would of course produce infinity because pp → Z +2j lacks the virtual corrections of pp → Z + j. To regulate this divergence, we set ptjet min = 1 GeV and have verified that for jets with p T J > 500 GeV, this choice has a negligible effect on the differential cross section of e (α) 2 until deep in the infrared region, well beyond the point where resummation dominates. Additionally, we have verified that the distribution of e (α) 2 as measured in pp → Z+j at relative O(α 2 s ) is identical to that measured in pp → Z + 2j at Born level with ptjet min = 1 GeV, up to differences deep in the infrared. Using this procedure, we are therefore able to match to relative O(α 2 s ) with MCFM. We generate pp → Z +j events through relative O(α 2 s ) accuracy at the 13 TeV LHC using MSTW 2008 NLO parton distribution functions [102]. We require that the p T of the Z boson is greater than 300 GeV and the absolute value of its pseudorapidity is less than 2.5. Jets are clustered with the anti-k T algorithm with radius R = 0.8. We study the hardest jet in these events that satisfies p T J > 500 GeV and |η J | < 2.5. On these identified jets, we then softdrop groom and measure e (α) 2 using custom code. This is an exceptionally computationally demanding procedure at relative O(α 2 s ), due to the complicated phase space of real emissions and the small width of the bins required to calculate the e (α) 2 distribution. This precision jet substructure study is only possible because of the development of highly efficient methods for generating fixed-order corrections.
In Fig. 11 we plot matched distributions for soft-drop e 2 with z cut = 0.1 and both β = 0 and β = 1 at NLL and NNLL. Here, we show both the distributions normalized to the total cross section and normalized over the range e 2 ∈ [0.001, 0.1]. The shaded bands represent estimates of theoretical uncertainties due to residual infrared scale sensitivity. 7 We show these bands mainly to allow comparison of the uncertainty remaining at different levels of formal precision. For the collinear-soft and jet functions in the resummed cross section, we vary the low scales by a factor of two. To estimate the scale dependence of the matching coefficient D k in the resummed cross section is more complicated, and we discuss this in detail in App. G. To estimate scale uncertainties in the fixed-order cross section, we vary the factorization and 7 The relatively large size of the uncertainty bands for e (2) 2 0.1 is an artifact of our simplistic additive matching. Additionally, due to the large K factor, the absolute scale of the matched NNLL distribution in Fig. 11b is roughly twice as large as the matched NLL distribution in Fig. 11a. ∈ [0.001, 0.1]. Note the reduction in uncertainties as one moves from NLL to NNLL, and also as one considers the normalized distribution. renormalization scales in MCFM by a factor of 2 about 500 GeV p T J . We then take the envelope of all of these scale variations to produce the shaded bands in Fig. 11. For β = 1 at NNLL, we have also explicitly shown the additional uncertainty due to the two-loop non-cusp anomalous dimension of the collinear-soft function. In going from NLL to NNLL accuracy, the relative size of the scale uncertainty bands decreases by about a factor of 2 or 3 for both choices of normalization of the distributions. However, normalizing the distributions over the range e (2) 2 ∈ [0.001, 0.1] dramatically reduces residual scale uncertainties; at NNLL, these normalized distributions have residual scale uncertainties at the 10% level and smaller.

Comparison to Monte Carlo
We now compare our NNLL resummed and matched calculation of soft-drop groomed e (2) 2 distributions to Monte Carlo simulations. We generate pp → Z + j events at the 13 TeV LHC with Herwig++ 2.7.1 and Pythia 8.210. To improve statistics somewhat, we have turned off Z/γ interference in the Monte Carlos. The Z boson is forced to decay to electrons, and we require that the invariant mass of the electrons is within 10 GeV of the mass of the Z boson. We then require that the identified Z boson has p T Z > 300 GeV and |η Z | < 2.5. Jets are clustered with FastJet 3.1.3 using the anti-k T algorithm with radius R = 0.8 and we identify the hardest jet in the event with p T J > 500 GeV and |η J | < 2.5. We then soft-drop groom this jet and measure e (2) 2 . Both soft drop and the energy correlation functions are implemented using FastJet contrib v. 1.019 [93,103].
We have generated two samples from both Herwig++ and Pythia to study the effect of hadronization and underlying event. One sample is purely parton level: As observed with jets in e + e − collisions, there is good agreement between our precision calculation and the Monte Carlos over a wide dynamic range. Importantly, this measurement of the soft-drop groomed e (2) 2 is very different from the case in e + e − . In e + e − collisions we calculated the heavy groomed and ungroomed jet masses. By measuring the heavier of the two jet masses, both masses have to be small, and the observable is global. For pp → Z + j events, we want to make no restrictions on the out-of-jet radiation. Thus although the soft drop jet mass is still free of non-global contributions, the ungroomed mass will not be. That is, we do not have control over all the large logarithms of ungroomed jet mass in pp → Z + j events, and thus cannot predict them using our factorized expression, although other approaches are possible. 8 For this reason, we only show distributions of soft-drop groomed e  10 −3 . That hadronization effects are small is expected from our e + e − analysis, but this also demonstrates that underlying event effects are negligible. A similar observation was made in Ref. [8], though at a much higher jet p T (p T > 3 TeV). As in e + e − collisions, we expect that the hadronization effects that are observed in the Monte Carlo can be explained by a shape function, though we leave this to future work. That the shape of the resummed distribution is both completely determined by collinear dynamics and is insensitive to underlying event suggests that by grooming jets with soft drop, we are able to completely isolate factorization-violating effects into an overall normalization. Therefore, we conjecture that the shape of the leading-power distribution of soft-drop groomed observables as measured in hadron collision events completely factorizes, just like the p T spectrum in Drell-Yan events [104]. We leave a proof of this conjecture to future work. 9

Conclusions
In this paper, we presented the first calculation for an observable measured exclusively on the constituents of a jet to NNLL accuracy and matched to fixed-order results at O(α 2 s ) relative to the Born process. The ability to do this calculation required grooming the jet with the soft drop algorithm, which eliminates the complications due to non-global logarithms that afflict ungroomed jet measurements. The soft drop groomer also significantly reduces nonperturbative effects from hadronization and underlying event, rendering the perturbative calculation of energy correlation functions accurate over several decades. The insensitivity of soft-drop groomed jet observables to underlying event suggests that the normalized cross section fully factorizes in hadronic scattering events. 9 Due to the presence of the complicated object D k (p min T , ηmax, zcut, R), Eq. (7.1) is not strictly a factorization theorem. It may not be possible to factorize D k to all orders due to the presence of so-called Glauber modes [104] in the cross section. While it is beyond the scope of this paper, recent work suggests that Glaubers can be included into the cross section directly [105], and our numerical work indicates that the effect may be absorbable into the normalization.
To complete the resummed calculation to NNLL accuracy required determining the twoloop non-cusp anomalous dimension for the soft function for which all emissions are removed. For β = 0, we were able to use results from the literature to extract the non-cusp anomalous dimension, up to calculable clustering effects. While not used for results in this paper, the clustering effects when using the anti-k T algorithm with soft drop are closely related to similar effects found in jet veto calculations. For soft drop angular exponent β > 0, we demonstrated a numerical procedure for determining the anomalous dimension using EVENT2. This was sufficient to approximate the non-cusp anomalous dimension, but a full calculation of the two-loop soft function for soft drop with β ≥ 0 is desired.
With a complete calculation of the two-loop soft function, including constants, we would be one step closer to resumming to next-to-next-to-next-to-leading logarithmic accuracy (N 3 LL). Up to the unknown four-loop cusp anomalous dimension (whose effects have been shown to be small [106][107][108]), the only other piece to get to N 3 LL would be the three-loop non-cusp anomalous dimension of the soft-dropped soft function. Without an explicit threeloop calculation, this anomalous dimension could in principle be estimated using a technique similar to what we used at two loops, using a fixed-order code like EERAD3 [109]. If this is possible, then resummation to this accuracy would potentially reduce residual scale uncertainties to the percent-level, assuming a scaling of uncertainties like observed in going from NLL to NNLL.
For our complete predictions, it was vital to match our resummed calculations to high precision fixed-order distributions. Fixed-order calculations have been traditionally used for observables that are inclusive over soft and collinear radiation, like total cross sections or p T spectra. The generation of fixed-order differential distributions for the plots in this paper required CPU-centuries, which we attained only by running on thousands of cores. For calculations of more complicated jet observables, precise fixed order computations are likely infeasible with presently available tools. As jet substructure pushes to higher precision, it will be necessary to have fixed-order calculations that more efficiently sample the infrared regions of phase space.
The calculations in this paper represent a new frontier of precision QCD. While jet substructure techniques have been used for some time in experimental analyses at the LHC, they are just now approaching the level of theoretical precision that can be meaningfully compared to data. By soft-drop grooming jets, we greatly reduce the theoretical challenges, enabling the calculation of a wide range of jet substructure observables to full NNLL accuracy.
Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.

A Three-Loop β-function and Cusp Anomalous Dimension
The β-function is defined to be For NNLL resummation, we need the β-function to three-loop order [71,72]. The first three coefficients are For NNLL resummation, we need the cusp anomalous dimension to three-loop order. The first three coefficients of the cusp anomalous dimension are [67][68][69][70]: T R n f , (A.7)

B Hard Function
The hard function for dijet production in e + e − collisions is defined by the Wilson coefficient for matching the full QCD current onto the SCET dijet operator. For e + e − → qq events, the Wilson coefficient C Q 2 , µ is Here,ψΓψ is the QCD current for the production of a qq pair from the vacuum. χn is a quark jet operator collinear quark operator defined in the light-like directionn in SCET. For calculations at leading power, χ n = W † t ψ, with W t a Wilson line pointing in some direction t not collinear to n and ψ is an ordinary quark field. The soft Wilson lines Y n and Yn point in the n andn directions respectively. Γ represents a generic Dirac matrix. We have ignored contraction with the leptonic tensor for simplicity. The Wilson lines Y n is defined as where P denotes path-ordering. Yn and W t are defined similarly withn µ and t µ replacing n µ . In SCET, the gluon fields in the Wilson line are soft gluons for the Y 's and collinear gluons for the W 's, but once the sectors are decoupled one can treat any of these gluons simply as a gluon field of full QCD. The hard function is the square of the Wilson coefficient: While we do not present its expression here, the hard function for e + e − → gg events is defined analogously, by matching the Higgs current F µν F µν onto SCET.
The one-loop hard function for the process e + e − → qq is [23,39,110,111] The cusp anomalous dimension of the hard function to all orders is where Γ cusp is the cusp anomalous dimension defined in Eq. (A.5). Similar to the cusp anomalous dimension, we define the coefficients of the non-cusp anomalous dimension γ via Through two-loops, the non-cusp anomalous dimension coefficients of the hard function are [73,74] γ (0) In the infinite top quark mass limit or with a finite Yukawa coupling, e + e − scattering can produce final state gluon jets. The hard function for such a process can be extracted from gg → H calculations. To all orders, the cusp anomalous dimension of the e + e − → gg hard function is where Γ cusp is the cusp anomalous dimension in Eq. (A.5). Through two-loops, the coefficients of the non-cusp anomalous dimension are [112][113][114] γ (0)

C The Global Soft Function
For arbitrary exponent β in the soft-drop groomer, the soft function can be calculated by requiring that soft gluons in measured jets fail the soft drop criterion. For hemisphere jets in e + e − → qq events, for example, the soft function is defined by the forward matrix element of soft Wilson lines: Here, n andn are the light-like directions of the qq dipole, T denotes time ordering, and Θ SD denotes the soft drop groomer operator which requires the final state to fail soft drop. The action ofΘ SD on soft final states cannot be written in a closed form for an arbitrary final state due to clustering effects, though it can be defined order-by-order. For example, the matrix element ofΘ SD for β = 0 on a final state with two soft particles was presented in Sec. 5.1.1. At one-loop for hemisphere jets in e + e − collisions, the soft function S G can be calculated from where n,n are back-to-back light-like vectors with n ·n = 2. The requirementn · k > n · k restricts the radiation to lie in one hemisphere, while the requirement restricts the soft gluon to fail soft drop. We find where C i is the appropriate color factor (C F for e + e − → qq; C A for e + e − → gg) and To all orders, the cusp anomalous dimension of the hemisphere wide-angle soft function is where Γ cusp is the cusp anomalous dimension from Eq. (A.5). To one-loop order, the non-cusp anomalous dimension is 0: γ (0) S = 0 . For NNLL resummation, we need the non-cusp anomalous dimension to two-loop order. As discussed in Sec. 5.1.1, for soft drop with angular exponent β = 0, this can be extracted from energy veto calculations, up to clustering effects that we calculated. For soft drop with β = 0 and Cambridge/Aachen reclustering, we find the two-loop non-cusp anomalous dimension to be (C.7)

D Jet Functions
Here, we present the quark and gluon jet functions on which the energy correlation function e (α) 2 is measured. The quark jet function, for example, is defined by the forward matrix element: Here, the jet is collinear to the light-like direction n, the operator δ(Q −n · P) restricts the large light-cone component of momentum to be equal to the center-of-mass collision energy Q, and δ (2) ( P ⊥ ) restricts the jet function to have zero net momentum transverse to the n direction. The measurement operator is defined by its action on an n-particle collinear final state |X n as:ê To write this expression, we have expanded the definition of the energy correlation function from Sec. 2.2 to leading power with collinear momenta. The gluon jet function is defined similarly: where B µ ⊥ is the collinear-gauge invariant operator in SCET that creates physical collinear gluons.
The following expressions will be presented in Laplace space, where renormalization is multiplicative and the Laplace space conjugate is ν. That is, The one-loop quark and gluon jet functions were first calculated in Ref. [25] for jets on which the two-point energy correlation functions with arbitrary angular exponent are measured.

D.1 Quark Jets
To one loop, the Laplace-space quark jet function is To all orders, the cusp anomalous dimension of the quark jet function is where Γ cusp is the cusp anomalous dimension from Eq. (A.5). For all α, the one-loop non-cusp anomalous dimension is γ For NNLL resummation, we also need the two-loop non-cusp anomalous dimension. For α = 2, corresponding to jet mass or thrust, this is known exactly. In that case, the non-cusp anomalous dimension is [79] γ q,(1) C α=2 (D.9)

D.2 Gluon Jets
To one-loop, the Laplace-space gluon jet function is To all orders, the cusp anomalous dimension of the gluon jet function is where Γ cusp is the cusp anomalous dimension from Eq. (A.5). For all α, the one-loop non-cusp anomalous dimension is γ For NNLL resummation, we also need the two-loop non-cusp anomalous dimension. For α = 2, corresponding to jet mass or thrust, this is known exactly. In that case, the non-cusp anomalous dimension is [115] γ g,(1) C 14)

E Collinear-Soft Function
The final piece in the factorization theorem is the collinear soft function, defined from soft radiation that is collinear to the jet. As it describes soft radiation, the collinear-soft function is defined as a forward matrix element of Wilson lines: The Y and W Wilson lines are the same as the ones in the soft and jet functions respectively, but depend on collinear-soft fields (which, like any of the others, can be treated as full QCD fields at leading power). Now, collinear-soft modes only contribute to e (α) 2 if emissions pass the soft drop groomer: this is denoted by 1 −Θ SD in the measurement function. (Recall thatΘ SD removes emissions from the jet according to soft drop.) Again, this operator cannot be written in closed form for an arbitrary final state due to clustering effects, but below, we will calculate it explicitly at one-loop. Theê (α) 2 measurement operator is defined by its action on an n-particle collinearsoft final state |X S,n : This follows from expanding the definition of the energy correlation function from Sec. 2.2 to leading power with collinear-soft momenta. This can be calculated at one-loop accuracy from where C i is the color factor of the jet. In the second equality, we have rearranged the phase space constraints and explicitly removed scaleless integrals. For this collinear-soft function, at one-loop in Laplace space we find where To all orders, the cusp anomalous dimension of the collinear-soft function is where Γ cusp is the cusp anomalous dimension from Eq. (A.5). To one-loop order, the non-cusp anomalous dimension is 0: γ S C = 0 . For NNLL resummation, we need the non-cusp anomalous dimension to two-loop order. For α = 2 and β = 0, this can be determined by renormalization group consistency of the cross section directly, using either the e + e − → qq or the e + e − → gg process. For soft drop with Cambridge/Aachen reclustering, the two-loop non-cusp anomalous dimension is γ (1) S C α=2,β=0 = C i −17.00 C F + −55.20 + 22π 2 9 + 56ζ 3 C A + 23.61 − 8π 2 9 n f T R . (E.7)

F Resummation
Because we work in Laplace space, defined according to where the anomalous dimension of F is γ. The anomalous dimension can be written as where Γ F (α s ) is the cusp part of the anomalous dimension, µ 1 is the infrared scale in the logarithm and γ F (α s ) is the non-cusp part of the anomalous dimension. The solution 10 to the renormalization group equation can be written more conveniently as an integral with respect to α s , by using the definition of the β-function as where µ 0 is a reference scale. The exponentiated kernels can be explicitly evaluated to any logarithmic accuracy given the anomalous dimensions. The cusp-part of the anomalous dimension, Γ F (α s ), is proportional to the cusp anomalous dimension, Γ F (α s ) = d F Γ cusp , where d F includes an appropriate color factor. The cusp anomalous dimension has an expansion in α s given by Eq. (A.5). The non-cusp anomalous dimension has a similar expansion defined in Eq. (B.7). For resummation to NNLL accuracy, we need the γ 0 and γ 1 coefficients, corresponding to computing the anomalous dimensions of the functions in the factorization theorem to two-loops.
With these expansions, we are able to explicitly evaluate the exponentiated kernel to NNLL accuracy. We have: The other exponentiated factor is Then, we can write the solution in Laplace space to the renormalization group equation in Eq. (F.5) as .

(F.8)
Because the hard function and the wide-angle soft function are independent of the observable e .

(F.9)
Note that the logarithms that appear in the low-scale jet function J(ν, µ 0 ) have the same argument as the factor that is raised to the ω J power. Therefore, using the relationship (noted by Ref. [80]) ∂ n ∂q n ν q = ν q log n ν , (F.10) we can re-write the jet function as .

(F.11)
Here J(L → ∂ ω J ) means that the logarithms in the low-scale jet function J(ν, µ 0 ) are replaced by derivatives with respect to the exponentiated factor ω J (µ, µ 0 ). The exact same replacement can be made for the collinear-soft function. In that case, we have S C (ν, z cut , µ) = e K S C (µ,µ 0 ) S C (L → ∂ ω S C ) µ 2 0 (ν e γ E ) .

(F.12)
This re-writing of the jet and collinear-soft functions allows for very straightforward inverse Laplace transformation. In Laplace space, the total differential cross section for left and right hemisphere jets in e + e − collisions is σ(ν) (F.13) = exp K H (µ, µ H ) + K S (µ, µ S ) + K S C (µ, µ (L) Note that the inverse Laplace transform commutes with the derivatives, and we have 2 ) −q−1 Γ(−q) .
(F.14) Therefore, the differential cross section in real space can be written as:

G Renormalization Group Evolution of D k
In this appendix we discuss in detail the renormalization group evolution of the jet flavor coefficient D k and explain the procedure we used to estimate the scale uncertainty introduced by neglecting higher-order terms.
The cross section for soft-drop groomed jets in pp → Z + j events factorizes in the limit e The fact that D k depends on multiple scales prohibits its resummation to all orders. Nevertheless, its renormalization scale dependence is completely determined by renormalization group invariance of the cross section. We can improve our prediction by solving the following renormalization group equation, which holds at leading power in z cut : where Q = 2p T J . The anomalous dimensions Γ D k and γ D k are Here, C k is the color Casimir for the jet of flavor k. The anomalous dimension has log z cut dependence, which means Q is not a natural scale of D k where all logarithms are minimized.
Nevertheless, we can still formally evolve D k from a scale µ 0 ∼ Q to a renormalization scale µ common to the jet and collinear-soft function. Solving the renormalization group evolution Eq. (G.2), the improved D k takes the form Here, µ f represents the factorization scale; i.e., the scale at which the parton distribution functions in D k are defined. The ω D k and K D k functions are defined in App. F. To estimate uncertainties from higher-order corrections due to residual scale dependence in D k , we will vary both µ 0 and µ f over the values For evaluating D k at fixed-order, we keep the full leading and next-to-leading terms as well as singular terms at the next-to-next-to-leading order in the following expansion of the solution to the renormalization group equation, Eq. (G.5). Expanding D k in powers of α s as D k (α s , µ, µ f ) = The non-singular terms c D k term cannot be determined without the next-to-next-to-leading pp → Z + j cross section. Note that the size of c (2) D k is no greater than O(α 3 s log 4 z cut ). Thus we can estimate that the size of uncertainty introduced by the unknown higher-loop non-singular term c (2) D k is roughly a factor of 1 ± α 2 s log 4 z cut e α n s L n+1 +··· , which is beyond NNLL accuracy.