EFT for Soft Drop Double Differential Cross Section

We develop a factorization framework to compute the double differential cross section in soft drop groomed jet mass and groomed jet radius. We describe the effective theories in the large, intermediate, and small groomed jet radius regions defined by the interplay of the jet mass and the groomed jet radius measurement. As an application we present the NLL$'$ results for the perturbative moments that are related to the coefficients $C_1$ and $C_2$ that specify the leading hadronization corrections up to three universal parameters. We compare our results with Monte Carlo simulations and a calculation using the coherent branching method.


Introduction
The physics of jets and their substructure has seen rapid development during the past decade. Originally designed for boosted object studies and pile-up mitigation at the LHC, the tools of jet grooming now play a key role in defining problems at the forefront of collider phenomenology [1]. Recently, soft drop grooming [2] (generalizing the modified mass drop algorithm [3,4]) has been widely studied, as the special nature of the algorithm enables a high precision description of measurement of infrared and collinear safe (IRC) observables in perturbative QCD. Furthermore, due to reduced sensitivity to nonperturbative effects of hadronization and contamination from the underlying event (UE), groomed jet observables are often more suited for LHC than their ungroomed counterparts. This has fueled a widespread interest in perturbative calculations of a variety of interesting jet based observables (see [5] for a review) such as the groomed jet mass [6][7][8][9][10], groomed angularities [11], soft drop thrust [12,13], groomed multi-prong jet shapes [14], groomed jet radius [15], energy drop [16], as well as calculations related to jets initiated by heavy quarks [17,18].
Groomed jet observables have also found applications in the field of nonperturbative (NP) nuclear physics, such as in improving our understanding of hadron structure by accessing the NP initial state in the form of structure functions such as the PDFs and the polarized/unpolarized TMDPDFs via scattering experiments such as DIS. This program entails that the observables chosen have small final-state nonperturbative effects, while at the same time be sensitive to the initial state hadronic physics. Groomed jets with an identified light/heavy hadron in the jet were proposed as probes of TMD evolution and distribution in [19,20]. Other observables defined with jet grooming such as transverse momentum spectrum of groomed jets [21] also meet this criteria. Soft drop jet observables have also recently been studied in the context of precision measurements at the LHC, such as the strong coupling constant α s [22] and the top mass [18]. To obtain precision collider physics analyses using groomed observables, theoretical control and understanding of final state NP effects is as essential as an accurate description in the perturbative regime, since the NP effects can be as significant as higher order perturbative corrections.
Analyses of NP corrections for jet substructure have been carried out in Refs. [4,6,23,24]. In Ref. [24] a field theory based formalism was developed to analyze nonperturbative (NP) corrections to soft drop jet mass due to hadronization. There the leading nonperturbative modes relevant for the largest hadronization correction to the groomed jet mass were identified, and were used to demarcate the regions modified by grooming in the jet mass spectrum into the soft drop nonperturbative region (SDNP), where the NP corrections are O(1), and the soft drop operator expansion region (SDOE), where perturbative contributions dominate, while NP corrections are small but still relevant for precision physics. It was shown that the leading nonperturbative where Q cut = 2 β Qz cut for hemisphere jets in e + e − collisions and is defined in Eq. (2.2) below for generic jet radius R and for pp collisions. Q stands for the large momentum such that Q = 2E J for the e + e − case and Q = 2p T cosh(η J ) for pp collisions. E J is the energy of the groomed jet, m J be in a reasonable agreement with the parton shower Monte Carlo (MC) simulations. From the point of view of phenomenological relevance, the terms appearing in the LL treatment are likely to suffice to capture the bulk of the NP corrections in the SDOE region, with higher order NP effects being below the level of perturbative uncertainty in parton-level predictions for groomed cross sections. However, since the coefficients C κ 1 and C κ 2 are defined such that they do not contain leading double logarithmic terms, ∼ (α s ln 2 ) k , they are so far not known with even a leading treatment of resummation effects. This requires an analysis beyond LL order, which is a main goal of our work here. At the same time we intend to provide a more realistic assessment of perturbative uncertainties from higher order effects.
In Ref. [15] results were presented for the cumulative soft drop R g distribution at NLL accuracy in SCET, including resummation of non-global logarithms due to CA clustering effects. The factorization for groomed jet mass was developed and studied at NNLL accuracy in Refs. [6][7][8] and N 3 LL accuracy for β = 0 in Ref. [36]. In this paper we build upon these results and describe the joint resummation of the soft drop cross section differential in the jet mass and cumulative in R g at NLL + NLO accuracy in the SDOE region.
Our analysis will involve so-called non-global large logarithms [37,38] that occur due to correlations between emissions across boundaries in phase space, and start at O(α 2 s ln 2 X), with X involving a ratio of scales associated to the measurements or specifying the phase space boundaries. In our case a boundary is introduced by measurement of θ g (as well as the jet boundary R). Another type of logarithm can be introduced by the jet clustering algorithm [39] such as the C/A used in the definition of the soft drop reclustering. These are referred to as abelian clustering logarithms, and are induced by uncorrelated soft emissions close to the jet boundary. Various techniques have been introduced to compute the non-global effects [35,37,38,40,41] and abelian clustering effects [42][43][44][45][46]. For our analysis we will follow Ref. [15], where both effects were studied in detail and accounted for in the cross section single differential in groomed jet radius. We will find the same pattern on non-global effects and clustering effects for the double differential distribution considered here. We also find that these contributions have only a small impact on the calculation of C κ 1 and C κ 2 . As a first application of our double differential cross section framework, we compute the following partonic moments that are related to C κ 1,2 (m 2 J ): modification for the Soft-Drop condition reads Θ sd = Θ(z − z cut θ β g ) → Θ sd (ε) = Θ(z − z cut θ β g + ε) , (1.6) with the obvious generalization to the case of more emissions. This modification preserves the IRC properties of the groomer for ε → 0. The superscript ' ' in M κ −1 (m 2 J ) signifies that the double differential distribution is evaluated at the "boundary" of the soft drop constraint, which is implemented via this ε-derivative.
From Eq. (1.4) we have C 1 = M κ 1 and C 2 = M κ −1 at LL order. We have used the notation for the moments displayed in Eq. (1.5) to stress that these functions can be defined independently of their interpretation at LL accuracy as the Wilson coefficients C 1,2 . Furthermore, higher order resummed results for the corresponding double-differential cross sections immediately translate into perturbative predictions for M κ 1 and M κ −1 at the same accuracy. At NLL order and beyond we anticipate that a significant portion of the higher order results for C 1,2 will still be captured by the moments M κ 1 and M κ −1 , such that C 1 M κ 1 and C 2 M κ −1 . In the subscript θ g ∼ θ g in the second line of Eq. (1.5), the angle θ g , defined below in Eq. (2.4), refers to the largest kinematically allowed groomed jet radius for a given measured jet mass value. The subscript signifies that the cross section is calculated only in the large θ g region where the corrections related to θ g measurement can be treated in the fixed-order perturbation theory, which is the most relevant region for C 2 . There are other regions of the groomed jet radius spectrum that we will discuss, where this is not the case. We will show that carrying out calculations of these moments via the double differential cross section offers us unique insights as to which contributions to the cross section lead to deviations from the LL definitions of the Wilson coefficients in Eq. (1.4). By identifying the relevant effective theories, we will see that these contributions come precisely from the regions of phase space that break the strong-ordering assumption of the partonic radiation. Restricting ourselves to the large θ g region allows us to determine the corrections that are consistent with the geometry which is a key part of the definition of the parameters Ω • • 1κ , Υ κ 1,0 , and Υ κ 1,1 and thus the terms needed to determine higher order corrections to the coefficients C 1 and C 2 .
Following the strategy formulated in Ref. [15], we preferentially work with the cumulant of the R g distribution: (1.7) We will carry out the analysis for the cumulative cross-section dΣ κ (R g ) at next-to-leading-logprime (NLL ) accuracy, where in accordance with the standard convention, the NLL resummation is supplemented by including all O(α s ) terms in the factorization formula. In this work we focus on the case of e + e − collisions, such that in the dijet limit, only the quark initiated jets will be of relevance. Note, however, that for the observable we are considering the leading non-trivial dependence on R g enters at O(α s ). This is because at least one emission off the energetic jet initiating parton is required to set a non-zero value of the groomed jet radius. We will come across scenarios where this R g dependent contribution appears in fixed-order perturbative terms. In order to consistently implement resummation in these cases we will include additional O(α 2 s ) cross terms resulting from combining the aforementioned O(α s ) fixed-order corrections with terms needed for a reasonable definition of NLL resummation. In other cases, where the kinematic hierarchies allow for the R g dependence to be factorized (thus facilitating resummation of logarithms of R g ) the standard prescription for NLL resummation applies.
The outline for the rest of the paper is as follows: In Sec. 2 we discuss the kinematic constraints on the joint R g and m J distribution and the relevant effective theory modes that enter the analysis. The EFT analysis identifies three different regions of resummation. We present the factorization theorems for each of these regions in Sec. 3 and discuss the calculation of the moments M q 1 and M q −1 . We briefly discuss extension of results in this paper to pp collisions leaving an in-depth phenomenological analysis to future work. In Sec. 4 we describe the resummation in each of the three regimes, as well as the resummation of boundary corrections to the cumulative cross section upon shifting the soft drop condition as in Eq. (1.6). This section also compiles the key formulae for the factorized cross sections in the large, intermediate, and small R g regimes, needed for the C κ 1 and C κ 2 calculations. The various factorized cross sections are then combined in Sec. 5 to arrive at the final result for the cumulative cross section that smoothly interpolates across these three regimes. In Sec. 6 we employ the results for matched cross section to calculate and perform a numerical analysis of the moments C q 1 M q 1 and C q 2 M q −1 , and compare the results with parton shower Monte Carlos and earlier coherent branching results. We then conclude in Sec. 7. Various more technical results and derivations are discussed in appendices.

Kinematic constraints on the R g distribution and mode analysis
In this section we discuss how the kinematic constraints imposed by the groomed jet mass measurement and the grooming condition result in the modes relevant for our effective field theory description of the double differential cross section. In this work we focus our analysis on jet measurements in e + e − collisions, while describing generalizations needed for the pp case.
Before discussing the various regimes that result from hierarchies within the kinematic bounds on R g , we first briefly review the soft drop algorithm and set up the notation. Soft drop reclusters the jet using the Cambridge/Aachen algorithm [47,48], and navigates the clustering history backwards, testing at each node the condition 2 where in the pp case R ij = ((η i − η j ) 2 + (φ i − φ j ) 2 ) 1 2 . So long as the tests keep failing, branches with lower energies or p T s are dropped, and the algorithm proceeds to the next node. As soon as one node passes the test, the algorithm stops and the remaining branches are kept in their Figure 1: Kinematic constraints on the range of R g for β = 0, 1, 2 with Q = 1000 GeV and z cut = 0.1. The R g is constrained between θ c (m 2 J ) and θ g (m 2 J , Q cut , β). The dot-dashed vertical green line corresponds to the groomed to ungroomed transition. The shaded brown region demarcates the phase space where hadronization corrections are critical. In the bottom-right panel we show the distribution of the double differential distribution from Pythia for representative values of z cut = 0.1, β = 1 and log 10 m 2 J /E 2 J = −1.5. We also display locations of the NLL kinematic boundaries and a rough estimate of the angle R NP g below which hadronization corrections become significant.
entirety. To keep expressions compact, it is convenient to introduce the quantity
In the following we take R ee 0 = π/2, so Eq. (2.2) reduces to Q cut = 2 β Q z cut . This quantity should not be confused with the (ungroomed) jet radius, that we assume is large, R ∼ 1 with R ≤ π/2.
We now discuss the expected bounds on the angular distribution, and define the different regimes of the analysis. Measuring the groomed jet mass m J limits the possible range of R g to 1 2+β , (2.4) where the expression for θ g holds for both pp and e + e − cases with the corresponding definitions of Q cut in Eq. (2.2). If we had R g < θ c (m 2 J ) then the jet mass would be constrained to be smaller than the desired value m J , and hence kinematically forbidden. While, on the other hand, having R g > θ g implies that the large-angle radiation passing the jet grooming constraint would necessarily yield a value of the jet mass larger than the one measured. These kinematic bounds on R g are displayed for β = 0, 1, 2, Q = 1000 GeV, z cut = 0.1 and R = 1 in Fig. 1. (The brown shaded regions will be discussed below.) In the final panel of Fig. 1 we also show the parton level θ g distribution of the cross section with fixed m J from a Pythia 8.2 simulation. We see that the kinematic end points θ c (m 2 J ) and θ g (m 2 J , Q cut , β) in Eq. (2.4) are consistent with the results from Pythia 8.2 parton shower.
For larger jet masses where θ g (m 2 J , Q cut , β) > R, the jet radius constraint limits R g ≤ R. This happens for jet masses given by with an analogous result for m 0 in the pp case (see Eq.(3.33) of [49]). This region is referred to as the ungroomed region in the context of single differential jet mass cross section, and here power corrections to the jet mass spectrum of O(m 2 J /(QQ cut )) cannot be neglected. Technically, in this region the grooming is still active but here the effects of grooming on the single differential soft drop m 2 J cross section are described via a fixed order treatment. Having an additional measurement of groomed jet radius, however, further necessitates resummation for logarithms of R g , which are also important in the ungroomed region.
In the following we will focus on the SDOE region defined by Eq. (1.1) where the following small angle approximation can be assumed: For the pp case the angles are cut off by R g ≤ R/ cosh(η J ) and hence the condition in Eq. (2.6) is replaced by θ g R/ cosh(η J ). For the purpose of setting up the effective theory description, we will interpret the inequality Eq. (2.6) in a hierarchical sense. With the inclusion of grooming, we therefore have 4 expansion parameters in our EFT namely θ c , θ * g , R g , and z cut . In the SDOE region we also have θ c θ g , which then leads us to consider three possible hierarchies between our expansion parameters.
1. Large groomed jet radius: . Intermediate groomed jet radius: 3. Small groomed jet radius: The main analysis in this work will focus on these scenarios, and we will extrapolate the results to the cases where θ g R or m J ≥ m 0 by encoding a basic description of these regions, but without trying to enforce the same resummation precision that we achieve for regions 1, 2, and 3. Since the perturbative moments in Eq. (1.5) are naturally computed through an integral over the groomed jet radius, their determination for a specific value m J require either combining or considering transitions between contributions from all three regimes.
We now return to discuss the brown shaded region in Fig. 1. With the addition of the R g measurement, the m J dependent condition in Eq. (1.1) for the single-differential spectrum should be further qualified to provide a restriction in the two-dimensional m J -R g plane. For low jet masses, nonperturbative effects become O(1) and the SDOE approximation in Eq. (1.1) no longer holds. This happens for jet masses in the soft drop nonperturbative region (SDNP) defined by For such low jet masses, any emission that stops soft drop is nonperturbative with a virtuality p 2 ∼ Λ 2 QCD . When we include an additional measurement of R g within the bounds in Eq. (2.3), the emission that stops soft drop can become non-perturbative for small R g , irrespective of the value of m J . The corresponding angle is obtained by solving the constraints p 2 ∼ Λ 2 QCD and z = z cut R β g which implies that the emission will be nonperturbative when R g 2 The angle (R g ) NP is the same as the θ ΛCS defined in Ref. [24] in the context of the SDNP region. In Fig. 1 we shade the regions R g R NP g brown for the various values of β in order to highlight that in these regions hadronization corrections are O(1) and cannot be ignored. The impact of the nonperturbative region grows with increasing β, so while the perturbative constraints defined in Eq. (2.3) are applicable for the entire range for β = 0, the NP region already covers a large part of the allowed phase space for β = 2. Thus for the combined R g and m J measurements, the SDNP region is defined by Eq. (2.8) which automatically implies Eq. (2.7), and the SDOE region is defined by which automatically implies Eq. (1.1). Since the partonic moments M κ 1 and M κ −1 defined in Eq. (1.5) involve integration over the entire allowed range of R g for a given jet mass, terms in the perturbative expression can become sensitive to small renormalization scales, which must be frozen at a scale µ 1 GeV to ensure they remain valid. If we studied hadronic versions of the moments M κ 1 and M κ −1 , then in some cases (like β = 2) they would differ substantially from their partonic counterparts (which are related to C κ 1,2 ), due to significant contributions from the brown region of Fig. 1. In the next section, we take detailed look at each of the regimes 1,2,3 and determine the corresponding factorization formulae utilizing a mode analysis in SCET. Figure 2: The mode picture for three regimes of large (left), intermediate (center) and small (right) groomed jet radius, labeling emissions with their energy fraction z and angle θ relative to the emitting particle. The area shaded in yellow corresponds to soft dropped radiation, while the area hatched in magenta is vetoed by the measurement (black lines), where dashed segments do not enforce vetoes. Emissions in the area shaded in grey lie out of the initial (ungroomed) jet.

Factorization
In this section we discuss the modes that arise in the three regimes discussed above. These modes are best depicted in the Lund plane shown in Fig. 2, where we have chosen the variables θ, the angle of a soft/collinear emission relative to the jet axis, and z, the energy fraction relative to the total jet energy E J (or p T in case of pp collisions). Every point on the phase space corresponds to an emission off the jet at a given angle and with given energy fraction. The choice of lnθ −1 -lnz −1 helps us to focus on the soft-collinear emissions that are far away from the origin [6]. In these coordinates distinct points correspond to emissions that are exponentially far apart. Here the jet mass measurement corresponds to a line with slope −2, the equation of the soft drop constraint is a line with slope β, and the equation θ = R g is a vertical line. The three regimes are shown in the three panels of Fig. 2, and differ due to the location of the vertical R g -line. At LL accuracy, in the soft-collinear limit, the probability of an emission in a certain region of the Lund plane is uniform in the logarithmic coordinates. Emissions that are vetoed at LL accuracy correspond to the area hatched in magenta, while the area shaded in yellow identifies soft dropped emissions. The various modes that enter the factorization analysis are located on the boundary of the vetoed region, at the intersections of the various constraint lines, and are indicated by colored circles. Since these modes are far from the origin they are either soft, collinear, or both soft and collinear. It is useful to determine the momentum scaling of the modes represented by colored dots in Fig. 2. We start with the modes that are common to all the scenarios, and in subsequent subsections describe the additional collinear-soft modes that are specific for each of the three regimes. In what follows, we will specialize our analysis to an e + e − collider. In the dijet limit, the differential jet mass cross section cumulative in R g can be expressed as (3.1) This notation holds for both hadron and parton level cross sections. Note that in the dijet limit, at leading power we have E J = Q/2, and hence we drop the differential dΦ J = dE J here and in the expressions below. The normalization factor N e+e − q obtains contribution both from hard modes with p µ ∼ Q, as well as from the hard-collinear mode with θ = R and z ∼ 1 which scales as Here and below we display the momenta scaling in the light-cone coordinates p µ = (p + , p − , p ⊥ ) defined with respect to a light-like reference vector n J = (1,n J ) in the direction of the jet axiŝ n J , such that wheren µ J is an auxiliary light-like four-vector satisfyingn 2 J = 0 andn J · n J = 2. The hard and hard-collinear modes are the same mode for R/2 ∼ 1. We assume that the jet has been identified via an IRC safe jet algorithm and the details of auxiliary measurements outside the jet will be irrelevant to our analysis here. The global-soft (S G ) mode is located at the intersection of the soft drop line and the jet radius constraint at θ = R, and also contributes to N e + e − q . These S G modes account for the radiation that was initially clustered in the jet but failed to pass soft drop. These modes have momentum that scales as We find it helpful to isolate the hard and global-soft contributions as follows: where the effect of grooming is entirely accounted by the global-soft function S q G . Here ⊗ Ω indicates that the angular integrals in the two functions involved cannot be done independently since both the modes have the same angular scaling, i.e. they cannot be distinguished by their angular separation. The sum j in the convolution in Eq. (3.5) sums over the individual Hard partons. This means that there are new logarithms, known as non-global logarithms (NGL's) in the ratio of the virtualities of the two modes that are involved in the angular convolution, and which can only be resummed by doing a more sophisticated resummation or RG evolution. Analytical resummation using the complete factorization formula is difficult and Monte Carlo techniques are usually used for resumming these NGLs.
Following [15], at NLL accuracy, the NGL's can be accounted for by writing where the new functions N q and S q G are obtained by doing independent solid angle integration and the leading effect of the NGL's is accounted for by the function Ξ q G . The global-soft function S q G is a matrix element of soft-Wilson lines with the modes of global-soft scalings, and has been extracted at two loops [6,50]. The hard function N q for dijet production in e + e − collisions is known up to three loops for the hemisphere case where it is determined by the quark form factor [51][52][53][54][55][56], while for calculations with a jet of radius R see Refs. [57][58][59][60][61]. Note that the N q (Q, R, µ) function appears universally for any e + e − dijet observable involving a dijet based e + e − measurement. We note that the NGL's associated with the function S q G do not affect the measurement of our observables directly but only modify the overall normalization. While this is important in determining the relative fractions of quark and gluon jets in pp collisions, in this paper we only consider the e + e − case and hence only deal with quark jets. At the same time, the observables we are interested in are the normalized moments of the double differential cross section which actually do not depend on S q G and N q functions. We will therefore ignore the Ξ q G function in our numerical analyses and discuss only the other angular averaged functions. For consistency we will leave Ξ q G explicit in the formulae in this section.
Next we describe the collinear mode (C) represented by a blue/purple dot in Fig. 2. This is the largest energy mode, z ∼ 1 that contributes to the jet mass. Requiring m 2 J E 2 J results in the mode being at the smallest angle. Thus, the energetic emissions represented by this mode always pass grooming, and have momenta scaling as Note that the corresponding angle for these modes is θ c = 2m J /Q. As shown in Fig. 2, additional soft-collinear modes arise from the interplay between the jet mass measurement, the grooming condition, and the R g constraint, and depend on whether we are in the regime with large, intermediate, or small groomed jet radius. We will discuss these modes and formulate the factorization formulae for the cross section in Eq. (1.7) for each of these regimes in turn. In general the ungroomed fat jet of radius R is typically initially isolated with the anti-kT algorithm and is subsequently reclustered in the soft drop procedure using the C/A algorithm. Thus, in addition to NGL's, sequential clustering algorithms give rise to "clustering logarithms" associated with independent emissions that appear at the same order as the NGL's. They can also be accounted for using a multiplicative function in the manner of Eq. (3.6). We describe them in detail in App. C.2 and discuss them in each region of factorization in the following subsections.

Large groomed jet radius
We start for the case represented on the left of Fig. 2, where the groomed jet radius is comparable with the maximum angle R g θ g (m 2 J , Q cut , β). Here soft drop is stopped by a wide-angle emission with relatively large virtuality, described by a collinear-soft mode (orange), with (+, −, ⊥) momentum scaling If other such emissions pass soft drop, their virtuality is large enough to contribute to the groomed jet mass. In this case there is a single collinear-soft mode in addition to the modes discussed above, and the factorization formula is given by a generalization of the jet mass factorization formula [6] to incorporate a more differential collinear soft function: cut , β, µ , cut , β, µ .
In the second equality we have made use of Eq. (3.6) to perform the angular integrations. In Eq. (3.9) the contributions from collinear regions are encoded in the jet function J q (m 2 J , µ), which is universal across various ungroomed and groomed event shapes and known up to three loops [62][63][64][65], and by the collinear-soft function S q cut , β, µ that we introduce here, together with renormalization group evolution for the scales between these functions. Since these modes contribute additively to the groomed jet mass measurement, they enter with a convolution in Eq. (3.9). For this integration we use + = m 2 J /Q, the small momentum component of collinear(soft) radiation. Extracting the overall factor of Q 1 1+β cut allows for rewriting S c as a combination of only the arguments in brackets. Note that in the following discussion and also in later sections, we will keep the κ index generic when discussing factorization functions, with κ = q being the relevant case for e + e − collisions.
The collinear-soft function S κ c + Q 1 1+β cut , β, µ in Eq. (3.9) is closely related to the collinear-soft function that appears for the single differential jet mass distribution S κ c + Q 1 1+β cut , β, µ , but with an additional constraint from the groomed jet radius. Including the R g measurement at one loop for e + e − collisions yields the following result: where C q = C F = 4/3, C g = C A = 3 are the standard SU(3) fundamental and adjoint quadratic Casmirs. In general, we note that results for e + e − collisions can be expressed in terms of tan Rg 2 , which allows for a smooth transition to R g R region for large jet masses. For the kinematic range we explore, tan cut , which applies for both e + e − and pp collisions, We elaborate on this function further in Sec. 4.3 when we discuss the resummation of large logarithms in this regime.
For the single differential jet mass cross section in the SDOE region, Ref. [6] showed that the Q and z cut dependence enters the collinear-soft function S κ c through the combination + Q 1 1+β cut , as we have displayed in Eqs. (3.9) and (3.11). For the double differential cross section we now show that the additional measurement R g appears in the combination 2 + /R g , or equivalently cut . For a given set of momenta {p µ i } that enter the calculation of the collinear-soft function we perform the rescaling Hence, the angle of these subjets or particles in the new (dimensionless) coordinates {k a i , k b i , k ⊥ i } is given by Similarly the relative angle between any two particles or subjets, θ ij , after this rescaling becomes 2+βθ ij , where the rescaled relative angle is only a function of the new coordinates, θ ij =θ ij (k i , k j ). In addition to the soft drop test, and the jet mass measurement, that already appear for the single differential distribution, we now have an additional comparison of angles θ ij with R g , such that the measurement function now additionally involves Here the i, j ∈ J product is over all subjets i and j which are present when the jet grooming has terminated. Eq. (3.14) together with the arguments for the single differential case from [6] demonstrate that only the variable combinations R g Q 1 1+β cut and + Q 1 1+β cut appear. Finally, we note that the NGL's for the modes with θ ∼ R described by the function Ξ q G (z cut , β, R) only affect the overall normalization, they do not affect the normalized moments that we are interested in for C 1 and C 2 , and hence we ignore Ξ q G for the rest of the paper. NGLs usually appear due to correlation between emissions at the jet boundary. In our case, we also have a boundary associated with the groomed jet and it is natural to ask whether there are NGLs that appear here. A calculation in App. C.1 shows that for the large R g regime, the associated logarithms are not large logarithms, and hence can be treated in fixed order perturbation theory.The same will also apply to the abelian clustering logarithms described in App. C.2.

Intermediate groomed jet radius
We next turn to the case of groomed jet radius measurement in the intermediate regime, whose modes are depicted in the central panel of Fig. 2. Because here R g θ g , the soft drop is now stopped by an emission with smaller angle (marked CS g in yellow). The virtuality of these emissions is too low to contribute to the jet mass measurement; however, the groomed jet mass is affected by emissions with larger virtuality that see the groomed jet boundary (marked CS m , shown in red). We can think of these two modes as a refactorization of the collinear-soft, CS, mode from the large groomed jet radius regime into CS m and CS g modes at the same angle R g with the following (+, −, ⊥) momentum scaling: Given the scaling in these formula, we see that both the resulting collinear-soft functions will give emissions at the same angle while being hierarchically separated in their energy. This is a similar situation to the one described in [41] with the result that the factorization for this process is complicated by the presence of NGLs. Specifically, it leads to a multi-Wilson line structure for the matrix element of the CS g modes, with a Wilson line for each final state CS m emission. We have already encountered for a similar hierarchy of energies between the global soft and the hard function in Eq. (3.5), which accounted for wide angle emissions (θ ∼ 1) leading to a multi-Wilson line structure for the function S q G in Eq. (3.5). Hence, following Eq. (3.6), the factorization formula for this regime is given by The hard, soft, and jet functions are the same as in Eq. (3.9). Consistent with the mode picture, we now find two collinear-soft functions, where the subscripts c g and c m specify the modes that are responsible for stopping the groomer and contributing to the mass measurement, respectively. Here, as in Eq. (3.5), ⊗ Ω indicates that the angular integrals in two functions involved cannot be done independently and the sum j in the convolution sums over contributions from various CS m partons. However, unlike the NGL's in the global soft function, the NGL's associated with the S q cm , S q cg functions do affect the measured quantities R g and m J and hence must be more carefully accounted for. To do this we write where the functions S q cm , S q cg are obtained by doing independent solid angle integration over each function and the NGL's are included in the function Ξ q S . The double differential cross section is given by cut , β, µ (3.18) A detailed calculation of the lowest order result for Ξ q S is presented in App. C.1, yielding while the calculation of abelian clustering logarithms is presented in App. C.2. We will make use of these results in Sec. 6 to demonstrate that the numerical impact of leading NGLs and abelian clustering logarithms is negligible at the order we are working (within our perturbative NLL' uncertainties). For this reason we do not need to bother with resummation of these logarithms. At one loop the results for the two collinear-soft functions in Eq. (3.18) are Here we define the plus functions to integrate to zero over the interval x ∈ [0, 1]. The result in Eq. (3.20) for S κ cg was presented in [15], while the calculation for S κ cm is analogous to that of jet mass (or jet angularities) with a specified jet radius R, here replaced by R g [57]. Expanding away terms that are needed for a large R g but are power corrections for this intermediate R g region, we have the following relation between functions in the large and intermediate R g regions This refactorization can be checked from Eq. (3.20) by making a comparison with with Eq. (3.11).

Small groomed jet radius
Last, we consider the regime θ c R g θ g , where the groomed jet radius is set by the smallest opening angle compatible with the jet mass measurement. Here the collinear-soft radiation responsible for stopping grooming CS g has an angle comparable to the more energetic collinear emissions but does not contribute to the jet mass measurement. The modes CS m and C thus collapse into a single mode that has the following scaling set by the groomed jet radius: The second relation shows that in the small R g region, the scaling of the mode is equivalent to the momentum scaling of the collinear mode in Eq. (3.7). For this case the factorization formula for the cross section reads Following the same argument that lead to Eq. (3.16) for the case of intermediate groomed jet radius, the factorization formula given here also contains NGLs associated with the function pair S cg , C q and the function pair S q,k G , N q,k q which each encode emissions at two different sets of angles, R and R g respectively, but are hierarchically separated in energy. As for the intermediate case, we can write the factorization in terms of the angular averaged functions and a function Ξ q C (z cut , R g , β). The calculation for Ξ q C is very analogous to that of Ξ q S in Eq. (3.17) and is discussed in App. C.1.
In this EFT regime, power corrections in the ratio θ c /R g = 2m J /(QR g ) cannot be neglected, and are captured by the new dimensionless collinear function C q [m 2 J /(Q 2 R 2 g ), QR g , µ] in Eq. (3.24), which describes energetic emissions that set the groomed jet mass. The collinear-soft radiation with similar opening angle, but much lower energy, responsible for stopping grooming, is still described by the same S q cg function as in the regime of intermediate groomed jet radius. Consistency of factorization with the case of intermediate groomed radius requires (3.25) We have computed this new jet function at O(α s ), finding for quarks This result provides an explicit check on the refactorization in Eq. (3.25) when expanding away power corrections in the ratio m J /(QR g ). Note that C q depends on the renormalization scale µ only through the Dirac delta function term in the first line of Eq. (3.26). Therefore, as expected by RG consistency with the remaining m J independent functions in Eq. (3.24), the anomalous dimension for C q is independent of the groomed jet mass. Since in this small R g regime we have 4m 2 J /(Q 2 R 2 g ) = θ 2 c /R 2 g ∼ 1, the terms in C q with more non-trivial m J dependence do not involve potentially large logarithms. We also note that these non δ-function terms in Eq. (3.26) vanish for R g < 2θ c (rather than θ c , as expected from Eq. (2.3)). This is a one-loop accident due to the single-splitting geometry of the event at this order, so the phase-space region θ c < R g < 2θ c will be filled by higher-order corrections and RG evolution. Having discussed the factorization formulae for the double differential cross section dΣ(R g )/dm 2 J in the three regimes we now turn to how we use this cross section to determine the moments M q 1 and M q −1 . According to the definition in Eq. (1.5), the moment M q 1 (related to the coefficient C q 1 of the shift power correction), simply requires taking the normalized first angular moment of the double differential distribution. Explicitly, We can remove the need to explicitly take the R g derivative in Eq. (3.27) by integrating by parts, such that (3.28) Note that despite C q 1 being determined by contributions from the large R g region, we are free to use the cross section which interpolates smoothly between all regions of R g here. This is because when we integrate this moment starting from the lower limit θ min contributions from the small R g or intermediate R g regions are power suppressed. We will describe the resummation of the large logarithms of jet mass and groomed jet radius present in the double differential distribution in three regimes in the next section, and combine the results from each regime in Sec. 5 to obtain a consistent description across the entire m J -R g phase space of interest. This result will then be used to calculate the moment M q 1 (m 2 J ) via Eq. (3.28) and will be interpreted as a result for C q 1 (m 2 J ) at leading power. Next we turn to the moment M q −1 that is related to the coefficient C q 2 (m 2 J ) which is the Wilson coefficient for the power corrections at the boundary of the groomed region and hence is evaluated when the soft drop condition is just satisfied, and was defined above in Eq. (1.5) using the soft drop condition with a small shift in the threshold energy cut by ε in Eq. (1.6). The derivative with respect to ε evaluated at ε = 0 allows us to to implement the δ function in Eq. (1.4) as can be seen by expanding the shifted measurement in Eq. (1.6) to first order in ε. 3 This is the effective soft drop condition that we will implement on our final state. For the purpose of considering the renormalization group consistency, it is simplest to first do the resummation, and only take the derivative with respect to ε after obtaining the RG evolved functions.
involves an an average value of the inverse groomed jet radius it appears to be more sensitive to the small angle regime of the θ g spectrum. However, as part of the definition in Eq. (1.5) we also include the restriction to the large R g region (θ g ∼ θ * g ), since as discussed earlier, the non-perturbative parameters and Wilson coefficient C q 2 encode the geometry associated to this region as part of their intrinsic definitions. Therefore to compute the partonic coefficient C q 2 , we must restrict ourselves solely to the region where large R g regime contributes. We will describe in Sec. 5.2 how this is accomplished. Note that this is also consistent with the fact that the nonperturbative brown shaded region in Fig. 1 should not contribute to the perturbative Wilson coefficient C q 2 . From the mode analysis in Fig. 2 we saw that the large R g and the small R g regimes constitute the upper and lower boundaries of the R g values with intermediate R g filling the bulk, so the large R g region is separated from the NP region shaded in brown. Finally, we will show in Sec. 4.5 that taking LL limit of our calculation of the soft drop boundary cross section in the large R g regime derived in Sec. 4.4 reproduces the LL result for C q 2 from Ref. [24].

Extension to pp collisions
While in this paper our focus lies on the simpler case of e + e − collisions, we, however, note that most of the expressions derived here and in the subsequent sections have straightforward extensions for calculations for pp collisions. We discuss them briefly here and leave a more indepth phenomenological analysis of the double differential measurement on groomed jets in pp collisions to future work. Factorization for jets in pp collisions can be studied either in the context of exclusive N -jet measurements [57,66], or while being inclusive over the external radiation [8,15] and restricting to small jet radii, R 1. In the following discussion, we will focus on the inclusive case where a single jet is being groomed and studied. In hadron collisions, we will have jets initiated by both quarks and gluons in the final state, and thus the cross section will involve a sum over κ = q, g. The equivalent of the normalization factor in Eq. (3.1) will now also account for parton distribution functions and the hard scattering cross section. In particular, the following analog of Eq. (3.1) will now hold true for inclusive jet measurements in the pp collisions: where we sum over the contributions to the cross section for the jet initiating parton κ = q, g.
Here, the large momentum component of the jet is no longer the center of mass energy Q, but now must be related to the jet p T as p − J = 2p T cosh η J , and thus the normalization function N κ now additionally depends on the phase space variables of the jet Φ J ≡ {p T , η J }. We also display its dependence on additional quantities denoted by '. . .' that include any other variables pertinent to the activity outside of the jet, such as veto on additional jets or momenta of any non-colored particles. The function N pp κ also contains the parton distribution functions. In inclusive measurements, the meaning of the global-soft radiation remains the same as above: the global soft mode S G continues to describe the radiation from the jet initiating parton, or any external soft radiation clustered with the jet, that has been groomed away. Thus the decomposition of the normalization factor in Eq. (3.5) and the further factorization of the leading NGL effects at NLL accuracy in Eq. (3.6) continue to hold for the pp case: where we make use of the pp version of the definition of Q cut in Eq. (2.2). Here, the function Ξ κ G (z cut , β, R), as before, describes the leading NGL effects for the jet initiated by parton κ.
We also can observe that the component of the cross section dΣ κ in Eq. (3.29) that is sensitive to the jet mass and groomed jet radius, i.e. the collinear and collinear-soft radiation that is kept after soft drop, will be described via the same factorization functions, such as the jet function J κ in Eq. (B.6) and the collinear function C κ (shown only for the quark jets) in Eq. (3.26), and the collinear soft functions S κ c , S κ cg , and S κ cm in Eqs. (3.11) and (3.20). This also include the functions Ξ κ S and Ξ κ G accounting for leading NGL's in Eqs. (3.17) and (3.24), but with an additional dependence on the jet initiating parton κ. Thus, from Eqs.

Resummation
Having setup the EFT factorization in the three regions, with the formulae in Eqs. (3.9, 3.18, 3.24), we next implement resummation of large logarithms by exploiting these equations. The resummation in SCET is carried out via a UV renormalization in the EFT and a subsequent renormalization group evolution (RGE) of these matrix elements. The RGE evolves these matrix elements from their natural scales µ i where logarithms are minimized to a common final scale µ and results in resummations of large logarithms of ratios of the µ i scales. The resummed cross sections will then allow us to extract the moments M q 1 (m 2 J ) and M q −1 (m 2 J ) introduced in Eq. (1.5). For simplicity, we consider jet production in an e + e − collider, and discuss resummation for the double differential cross section in the three regimes discussed above.
As we mentioned above in Sec. 1 we will carry out the resummation at NLL accuracy, where in addition to the NLL resummation, we will include NLO fixed order corrections to the perturbative functions in the factorization formulae in Eqs. (3.9), (3.18) and (3.24). In Tab. 1 we summarize the orders to which the cusp anomalous dimension Γ cusp [α s ], various non-cusp anomalous dimensions γ i [α s ], and β function are needed to achieve LL, NLL and NLL accuracy. We collect expressions for the evolution kernels, anomalous dimensions, and details of evolution in App. A. In the intermediate R g regime, the R g dependence appears directly through the evolution scales, and there are no further complications. However, in the large R g regime we see from Eq. (3.11) that the R g -dependence is only seen at one loop in the collinear-soft function through a fixed-order term. Thus in the large R g regime we must include additional R g dependent fixed order terms as indicated in the next to last column of Tab. 1. Similarly, in the small R g regime, the measurement of the jet mass appears as a fixed order correction on top of the R g measurement. As we can see from Eq. (3.26), the leading order (LO) distribution for m J > 0 starts at one-loop. Thus, in our NLL treatment we will include additional fixed order terms relevant for these regimes, as indicated by entries in the last two columns of Tab. 1.
Here a generic function F in the factorization theorem is assumed to have a fixed order expansion of the following form: where L F are logarithms (transformed to Laplace space, m 2 J → s, see App. A for details) and the a F mn are coefficients which may contain non-logarithmic dependence on the variables. For the large R g factorization formula, the collinear-soft function S κ cut , β, µ receives a one-loop R g dependent correction that does not involve logarithms of µ, which we indicate by α s a Sc 10 R g /θ g . This one-loop term constitutes the LO distribution of R g in this regime and depends on R g through the dimensionless ratio R g /θ g , where θ g was defined above in Eq. (2.4). At NLL order one would additionally add the R g dependent O(α 2 s ) terms in the large R g regime, in order to include the dependence on R g to one higher order. In our analysis we only include the O(α 2 s ) terms indicated in Tab. 1 which are determined by lower order anomalous dimensions and the one-loop constant terms, which omits the currently unknown O(α 2 s ) R g -dependent term in ∆S κ c (we will estimate the impact of this term in our numerical uncertainty analysis). Finally, as R g → 0, the ratio R g /θ g becomes a power correction leading to factorization in the intermediate regime, where the standard treatment of NLL counting applies.
A similar reasoning applies to the regime of small R g , where the collinear function C q defined in Eq. (3.26) receives one-loop, non-logarithmic corrections depending on m J . As mentioned above, for m J > 0, this is the LO distribution in the small-R g regime, which we denote by α s a C 10 (θ 2 c /R 2 g ), with θ c defined above in Eq. (2.4). We include the additional O(α 2 s ) fixed-order terms in this regime at NLL in the same multiplicative way, as summarized in the table. Thus we note that the treatment of terms is slightly different in the three regimes. We will further justify our specific approach through the smooth matching between regimes as well as parametrize the uncertainty due to the unknown non-logarithmic 2-loop pieces in Sec. 5. Furthermore, we note that, up to NLL accuracy, we do not need to consider O(α 2 s ) corrections to the µ i -dependent SCET functions parametrized by F in Tab. 1 (such as the jet, collinear-soft and global-soft functions) as their inclusion with O(α s ) LO distributions in the small or large R g regimes will only enter at O(α 3 s ).
In the following, we first discuss the simpler cases of resummation in the small and intermediate R g regimes and set up the relevant notation. Then we turn to the large R g region, where the interplay between the LO term depending on R g /θ g and RG evolution requires more complicated calculations with Laplace transforms. Finally, we describe the resummation for the boundary corrections relevant to the calculation of C q 2 (m 2 J ), where we must account for the shifted soft drop condition in Eq. (1.6). Lastly, we note that the double differential cross section is affected by nonglobal logarithms, which require resummation for a complete treatment of the double differential cross section in all regions. In Sec. 3 we have shown how these affect the factorization formulae for our cross section in each of the regimes, appearing in the functions Ξ q G , Ξ q S , and Ξ q C . In our numerical studies in Sec. 6 we will consider the leading order NGLs as well as abelian clustering logarithm that appear in these functions. We find that for C q 1 and C q 2 they only have an impact that is comparable to other missing O(α 2 s ) corrections, and hence resummed results for Ξ q G , Ξ q S , Table 1: Highest order of the anomalous dimensions and fixed order terms that are needed for the given accuracy. The last two columns indicate additional fixed order terms that are added in the large R g and small R g regime respectively.
and Ξ q C are not necessary for the calculation of these coefficients at NLL accuracy. Below in sections 4.1, 4.2, 3.3 we therefore explain how to carry out the resummation of logarithms associated to the scales appearing in the remaining functions in the factorization formulae.

Small R g
Resummation takes on the simplest form in the regime of small R g , where the factorization structure is purely multiplicative. As written in Eq. (3.24), each function depends on a common and arbitrary renormalization scale µ. Any one specific choice of µ induces in the fixed-order ingredients large logarithms of the ratio of widely separated scales (indicated by the function arguments), which can be minimized by choosing instead using a different scale µ i for each function. These scales are then connected to a common scale µ via renormalization group (RG) running, which is how the resummation of large logarithms is carried out in SCET.
Including the appropriate resummation factors, as discussed in App. A.1, the RG-evolved cross section in the small R g regime reads where the dimensionless prefactor N evol q is given by The resummation kernels are defined in Eq. (A.10). For compactness, in the subscripts of scales µ X and kernels K X and ω X , we use the shorthand notation gs, cs g and c for the symbols S κ G , S κ cg and C κ respectively. The functions are RG evolved from their respective scales µ i to a common scale µ, where natural scales are Here the superscript 'can.' denotes the canonical choice of the scale. The function C κ that appears in Eq. (4.2) is a generalization of Eq. (3.26) which includes two-loop terms, given by Following the notation in Tab. 1 we have identified the LO small R g distribution for non-zero jet masses as In the square brackets of Eq. (4.5) we have included the logarithms of µ/(QR g /2) that are generated at two-loops due to the anomalous dimension of the C q function and the running coupling. Furthermore, we have switched to the distribution differential in ln m 2 J , where the Jacobian factor m 2 J is absorbed (in the case of quarks) in the combination turning the distributions in Eq. (3.26) into ordinary functions, and leading to the expression in Eq. (4.6) for the LO distribution. The parameter a C 20 is included in order to parameterize the uncertainty resulting from our lack of knowledge of the two-loop non-logarithmic corrections in the small R g regime. For this uncertainty analysis we are making an approximation that the non-logarithmic corrections have the same functional form as the leading order distribution. At the same time, the modified argument of the function a C 10 in the second line reflects the difference in the end points of the LO and NLO distributions. In the LO distribution in Eq. (3.26) the end point of the R g distribution is given by R min g | LO = 2θ c . In contrast, the calculation of minimization of R g for a given jet mass for three particle configuration at NLO yields an end point at R min g | NLO = (8/5)θ c in the small angle limit, which was computed in appendix B of Ref. [7]. In Sec. 6 we will vary the parameter a C 20 in the range [−2π, 2π] to obtain an estimate of the uncertainty from missing O(α 2 s ) terms.

Intermediate R g
Resummation in this regime is conveniently performed in Laplace space, where convolutions over the groomed jet mass reduce to ordinary products. Taking a Laplace transform in the variable x conjugate to m 2 whereJ q ,S q cm , andΞ q S are the Laplace transforms of J q , S q cm , and Ξ q S respectively. The one-loop results forJ κ andS κ cm are provided in Eqs. (B.8) and (B.18) respectively, and that ofΞ κ S for the leading term in Eq. (C.7). Details of the Laplace transform between momentum and position space are presented in App. A.1.
At NLL accuracy, the simple one-loop expressions of the functions in Eq. (4.7) make it straightforward to take the inverse Laplace transform. To simplify the final formulae after taking inverse Laplace transform, we introduce an alternative notation for the Laplace space expressions, where we notice that the dependence of the functions on the r.h.s. on their first argument is purely logarithmic. Note that in the functionΞ κ S we make use of the physical µ can. csg scale defined in Eq. (4.4) in the argument of the running coupling. In Eq. (4.8) we have made explicit the Laplace space logarithms and the factors of running coupling. In terms of this notation, the cross section (for m 2 with the same prefactor N q evol defined in Eq. (4.3). We have identified the RG invariant single logarithmic kernelω asω (4.10) The first argument involving the derivative ∂ Ω in the Laplace transforms of the jet and the collinear-soft functions in Eq. (4.9) is understood to replace the logarithms that appear in the corresponding Laplace space expressions in the notation described in Eq. (4.8). In the subscript µ csm we introduced the shorthand label cs m for the symbol S κ cm . The collinear scale in Eq. (4.4) is now replaced by the two mass-dependent canonical scales In the small R g limit when R g = θ c = 2m J /Q both of these scales merge such that µ can csm = µ can. J = µ can. c = m J . Indeed, as required by RG consistency, the evolution factors in Eq. (4.9) reduce to the ones in Eq. (4.2) in this limit. Details of the Laplace transform between momentum and position space and formulae for evolution kernels in Eq. (4.9) are presented in App. A.1, and the one-loop results for the anomalous dimensions relevant for NLL running are collected in App. A.2.

Large R g
We now derive a formula for the resummed cross section in the region θ 2 c R g θ g (m 2 J ). The factorization structure and the one-loop ingredients relevant to this regime were discussed in Sec. 3.1. NLL resummation requires introducing a number of new functions, whose notation we describe below, while leaving presentation of the full expressions to App. B.1.
The key ingredient of our analysis is the collinear-soft function displayed at one loop in Eq. (3.11). Since its anomalous dimension is unaffected by the R g measurement, it has the same RG evolution as the collinear-soft function that enters the single differential jet mass distribution. At the same time, the R g dependence only arises at one loop, which makes the O(α s ) correction to the collinear-soft function the leading order contribution to the double differential cross section. This motivates rewriting the collinear-soft function in the following multiplicative manner: cut , β, µ] is singly differential in the jet mass and is responsible for RG evolution, whereas the R g dependent corrections are incorporated via ∆S κ c , which is by construction independent of the renormalization scale. Eq. (4.12) is just a reorganization of the perturbative series. However, including one-loop corrections in both functions on the r.h.s. allows us to supplement our predictions with terms of the form {leading RG logarithms}×{leading R g corrections}, which when taken together first appear at O(α 2 s ). In fact, we define NLL accuracy for the double differential distribution in the large R g region based on this reorganization. However, we stress that such a prescription does not capture all O(α 2 s ) terms in the collinear-soft function and in particular misses the genuine two-loop R g dependent corrections. We obtain an estimate for the uncertainty from these missing two-loop terms via the procedure described below.
To make the ∆S κ c term in Eq. (4.12) explicitly independent of µ, we can further rewrite it as where we have defined the one-loop contributions . This separation is again just a reshuffling of the perturbative series, as the O(α 2 s ) cross terms are predicted from RG consistency of the double differential cross section. The construction to make the R g dependent correction ∆S κ c in Eq. (4.13) explicitly independent of µ can similarly be generalized at higher orders. Finally, we note that expanding Eq. (4.12) to O(α s ) we recover the one-loop expression of the collinear-soft function given above in Eq. (3.11). As a further check, the R g dependent term vanishes when R g = θ g (˜ = 1 in Eq. (4.15)), and the collinear-soft function cumulative in R g of Eq. (4.12) reduces in this limit to the collinear-soft function single differential in the groomed jet mass.
Since our treatment lacks the 2-loop non-logarithmic R g dependent corrections, we have parametrized the resulting uncertainty in Eq. (4.15) via the parameter a Sc 20 which we will vary in our uncertainty analysis. This estimate utilizes an approximation that the two-loop corrections have the same functional form as that of the one-loop term, but with an unknown normalization. In the numerical studies presented below we will consider variation of a Sc 20 in the range [−2π, 2π]. We will return to discussing this when we study the numerical results in Sec. 6.
To carry out resummation in Laplace space the transform of the collinear-soft function is defined with respect to the variable + Q 1 1+β cut : cut , β, µ . (4.16) Eqs. (4.12) and (4.13) in Laplace space read In App. B.2 we show that carrying out RG evolution in Laplace space and transforming back where the Laplace space logarithms in the jet and collinear-soft functions are replaced by the argument in brackets. Here, in analogy with Eq. (4.8), we have defined the alternative notation for the Laplace transform of the collinear-soft functioñ Note that since all the scales and anomalous dimensions for this regime are independent of R g , so is the function of the derivative operator in Eq. (4.21). The effect of the cumulative R g measurement on Eq. (4.20) is accounted for entirely by the kernel Q κ defined as which is µ-independent by construction and vanishes at R g = θ g (m 2 J , Q cut , β). Here ∆S [β 0 ] c is expressed in the notation equivalent to that in Eq. (4.22), see Eq. (B.26). The explicit expression for Q κ is given in in Eq. (B.33). To make a connection with the notation introduced in Tab. 1, we can identify the LO distribution in the large R g regime that multiplies the NLL evolution factors in Eq. (4.21) as:

Boundary Corrections in the Large R g region
Next we discuss boundary corrections to the soft drop cross section in the regime where R g θ g (m 2 J , Q cut , β), relevant for computing the ε derivative in Eq. (1.5). We first compute the bare collinear-soft and global-soft functions with the shifted soft drop condition in Eq. (1.6). This will allow us to identify the relevant corrections to the double differential cross section at first order in the shift parameter ε. We will find that both the function of the derivative operator dσ q (∂ Ω ) and the evolution kernel Q q in Eq. (4.20) receive O(ε) modifications.

Boundary corrections to bare soft matrix elements
The O(ε) corrections to the soft matrix elements lead to logarithmic and non-logarithmic terms proportional to R −β g Qε/Q cut . Again, we keep the κ index general here and will use results for κ = q for the e + e − case. This can be seen from the bare result for the collinear-soft function with the shifted soft drop condition defined in Eq. (1.6) given by where the condition to pass or fail the R g constraint are written respectively as cut , β, µ (4.32) We have used the δ sd to perform the integral over p − , and have isolated the contribution from the soft-drop failing piece, proportional to δ( + ), so as to regulate the integral for p + → 0. The term proportional to δ( + ) in the third line is no longer scaleless as it is bounded from below. Solving the remaining integration, we find for β = 0 and β > 0 respectively: .

(4.33)
We see that for β = 0 the > 0 regulator results in a UV pole, while for β > 0 the result is a power correction and dimensional regularization is not needed. Here, we have defined 5 We can carry out a similar calculation for the O(ε) boundary corrections to the global-soft function. In full analogy with Eq. (4.31), we find that up to O(ε 2 ) terms where the additional argument Θ sd (ε) denotes now the shift in the soft drop failing condition. The result for the O(ε) correction to the bare global-soft function with shifted soft drop reads 36) 5 We note that for a > 0 the function L a 0 (x) is equivalent to 1/x 1−a which has an integrable singularity as x → 0.
where the additional terms not shown for β > 0 are O(1) and not relevant to the discussion here. The details of this derivation and complete results for β > 0 are provided in App. B.3. We see that the additional single logarithmic UV divergence for β = 0 case has the opposite sign to that in Eq. (4.33), as required by RG consistency. Interestingly, we now find that for β = 0 the oneloop collinear-soft and global-soft functions develop a non-zero non-cusp anomalous dimension proportional to ε. With conventions for prefactors in anomalous dimensions, summarized in Eqs. (A.13) and (A.16), we find Hence, for β = 0, the NLL resummed cross section for shifted soft drop will involve additional running between the global-soft and collinear-soft virtualities due to this non-cusp anomalous dimension. We now turn to addressing how the various O(ε) corrections are combined with the other EFT matrix elements so as to arrive at the result for resummed cross section with shifted soft drop condition.

Assembling the ingredients
cut , β . (4.38) Here the collinear soft function for single differential jet mass measurement includes the additional logarithm that is generated for β = 0 in Eq. (4.33), as signaled by the extra argument δ β,0 γ S κ c 0 (εz cut ), and is modified as follows: while the term ∆S κ c,ε is defined analogously to ∆S κ c given in Eq. (4.13): (4.40) The ∆S piece is the same as that defined in Eq. (4.14) and is again included to cancel the µ dependence due to the running coupling at NLL . Finally, from Eq. (4.33) we have Here the Θ(β > 0) signifies that the first term is only present for the β > 0 case, whereas the second term applies for all β ≥ 0. As we shall see later, this second term contributes to the bulk of the boundary corrections, and hence to C 2 (m 2 J ). Note that in defining ∆S κ c,ε in this way we factored out the R −β g power present in Eq. (4.33) for β > 0. Comparing this scaling with the O(ε) correction to the global-soft function in Eq. (4.36), we see that the R −β g power dominates the sin −β (R/2) power in the R g R approximation which is formally still valid in the large R g region. Thus we expect boundary contributions to the global-soft piece (for β > 0) to be numerically insignificant, and we will leave them out from our numerical predictions. Finally, since we lack the knowledge of the non-logarithmic two-loop boundary corrections in S κ c that cannot be predicted by the anomalous dimensions we have included an additional parameter a Sc 20,ε in spirit of Eq. (4.15) in facilitate estimating the uncertainty resulting from these terms. We will come back to this again in Sec. 6 when studying the perturbative uncertainties.
Taking the Laplace transform of Eq. (4.38), performing an expansion up to O(ε), and dropping terms that are O(α 3 s ) we explicitly get cut , β, µ (4.42) c,ε term is given in Eq. (B.50). The first line reproduces the O(ε 0 ) result in Eq. (4.12) and vanishes when taking the ε derivative in Eq. (1.5). The second line describes boundary corrections induced by the additional logarithm present for β = 0, while the remaining lines encode boundary corrections that are R g dependent. The product of the two ∆S [Cκ] c functions in the last line is important to ensure a smooth transition to the regime R g θ g , which justifies including the boundary corrections ∆S κ c,ε multiplicatively in Eq. (4.38). Likewise, we express the global-soft function with O(ε) corrections as with the result for µ-independent ∆S κ G,ε piece given in Eq. (B.43).

Cumulative boundary cross section
With all the ingredients for the shifted soft drop condition at hand, we are finally in the position to write resummed results for the cumulative cross section. Eq. (4.19) is still formally valid, but the global soft and collinear-soft functions that appear there now include boundary corrections according to Eqs. (4.38) and (4.43). As a consequence, the equivalent of Eq. (4.20) has two kind of terms: They describe boundary corrections to the NLL RG evolution and the R g dependence respectively, and we now examine them in turn. The first term reduces to the O(ε 0 ) cross section for positive β but, for β = 0, includes the additional one-loop fixed order logarithms and running due to the O(ε) one-loop non-cusp anomalous dimension displayed in Eq. (4.37). Since these corrections do not depend on the jet mass, they can be absorbed in a generalization of the operator dσ κ (∂ Ω ) in Eq. (4.21) yielding while the evolution kernel in brackets is unchanged. The operator is obtained from the following O(ε) modifications to its analogue in Eq. 3. The collinear-soft function has additional running: K cs µ, µ cs , δ β,0 γ S q c 0 (ε, z cut ) = K cs (µ, µ cs ) − δ β,0 η γ 0 (ε, z cut ), µ, µ 0 , where the LL single logarithmic evolution kernel η(Γ, µ, µ 0 ) is given in Eq. (A.11) from the one-loop non-cusp anomalous dimension.
4. The normalization factor N q evol has additional running due to the global soft function: K gs µ, µ gs , δ β,0 γ S q G 0 (ε, z cut ) = K gs (µ, µ gs ) + δ β,0 η γ 0 (ε, z cut ), µ, µ 0 . (4.47) The second term in Eq. (4.44) consists of the residual O(ε) non-logarithmic, R g -dependent corrections to the collinear-soft and global-soft functions in Eqs. (4.40) and (4.43), which apply for every β ≥ 0. Since these terms do not alter evolution, but have a nontrivial dependence on the jet mass and R g , they require a different kernel than Eq. (4.25): c,ε υ, β, µ cs 1 + ∆S [Cκ] c υ, β, µ cs , but use the same function of the derivative operator, dσ κ (∂ Ω ) defined in Eq. (4.21), valid for ε = 0: The explicit expression of the evolution kernel Q κ ε between jet and soft scales for shifted soft drop is given in Eq. (B.51). Note that the term ∆S κ G,ε for β > 0 at NLL accuracy contributes only to the boundary correction in the single differential jet mass cross section, and is thus independent of R g . In the following, we will simply ignore this term as this is also power suppressed compared to the term in the first line in Eq. (4.49).
At this point, we have the soft drop boundary cross section which is cumulant in the groomed jet radius. In order to evaluate C 2 , following Eq. (1.5), we need to switch to the cross section differential in the groomed jet radius and further take the ε derivative. We perform this explicitly in App. B.4, where we compile the results for the boundary cross section differential in R g . The large R g regime discussed so far is the only relevant geometry in the definition of C q 2 (m 2 J ). For completeness, we discuss in App. B.3 the boundary corrections in the intermediate and small R g regimes and show consistency with the results in the large R g regime as R g /θ g (m 2 J ) → 0.

Recovering the Leading Logarithmic Coherent branching result
In this section we make use of the results derived above for the cumulative cross section and the boundary corrections to demonstrate that the EFT that we have developed reproduce the computation of C 1 and C 2 at LL accuracy within the coherent branching formulation from [24]. This also makes explicit how the EFT encodes the strong angular ordering limit used there. The LL results for C 1 and C 2 were derived by considering a chain of emissions that are ordered in their respective contribution to the jet mass. 6 The emission at the largest angle is required to satisfy soft drop and set the jet mass. This leads to the following expressions for the two moments: where C q 0 (m 2 J )/m 2 J is the single differential jet mass distribution and is given at this order by, where the result for C 0 (m 2 J ) is valid up to LL. In Eqs. (4.50) and (4.51) we have rewritten the soft drop condition in Eq. (2.1) terms of Q cut . In the expression for C q 1 and C q 0 , the condition is imposed on the last emission with energy fraction z and angle θ. In case of C q 2 we see, however, that the soft drop Θ function is replaced by δ function in order to capture the corrections at the soft drop boundary. The radiator, R q (θ 2 c ), results from all the emissions at larger angles that are either vetoed by the jet mass measurement or that fail soft drop (and the virtual emissions). The radiator resums double logarithms of θ g /θ c , which however cancel out in the ratio in Eqs. (4.53) and (4.51). The splitting kernel p gq is given by Terms that are power corrections in the z → 0 limit are formally beyond LL accuracy, and we will drop them in the following; in the EFT formulation, these are captured by the NLO corrections to the jet function. Simplifying Eq. (4.50) and using the LL limit of p gq (z) in Eq. (4.52) we get where we have expressed the results in terms of the angular bounds θ c and θ g defined in Eq. (2.4).
To establish a connection with our results, we must compare with the large R g regime of Sec. 3.1, since due to the associated geometry for power corrections in the SDOE region of the single differential jet mass distribution, this is the primary regime that is relevant for defining the coefficients C 1 and C 2 . 7 Thus, the EFT prediction for the LL double differential cross section in this regime is obtained from the cumulant in Eq. (4.20), (4.54) In the second line we have used that the evolution in this regime is independent of the groomed jet radius and identified the LO differential operator with the exponential of the radiator R q . This follows by setting the Laplace space fixed order ingredients in Eq. (4.21) to unity. We now notice that the fixed-order corrections in Eqs. (4.53) and (4.51) are described by a single splitting function, in contrast with the convolution in Eq. (3.9): this implies in particular that the coherent branching result ignores corrections induced to the kernel Q q by evolution in Laplace space. We can therefore take the limit of no evolution by setting Ω → 0 in Eq. (B.33), obtaining where we have used that under the integral Furthermore taking Ω → 0 limit corresponds to using the same scale in both jet and the collinearsoft function. This is analogous to using a single splitting function p gq in the coherent branching formalism that is not factorized into jet and collinear-soft function. In Sec. 6 we will observe this difference to have some impact on the results. Next, we choose the natural scale choice for µ cs in Eq. (4.24) noting that This corresponds to setting θ = θ g in Eq. (4.53), with the difference between the two expressions being beyond LL accuracy. Substituting Eq. (4.55) in Eq. (4.54) gives the integrand of Eq. (4.53). 7 We note that the coherent branching description is based on ordering in a single variable, the jet mass contribution of each emission, ρi = ziθ 2 i . This is in correspondence with the mode picture in the large Rg regime where the measurement of jet mass provides a distinction between the various EFT modes (see Fig. 2, left). This is, however, not the case with small and intermediate Rg regimes: in the case of small Rg regime, the primary measurement is the groomed jet radius which distinguishes the modes based on their angular separation. In the intermediate Rg regime an analogous treatment in coherent branching would presumably involve ordering in two different variables.
The LL coherent branching result for C 1 then trivially follows by including the additional factor θ/2 in the angular integration.
To obtain the LL coherent branching result for C 2 from the EFT formalism we can follow the same steps. We first notice that the anomalous evolution induced by boundary corrections at β = 0 affects the cross section beyond LL, and we can thus ignore this correction. The relevant cumulant cross section for large R g is then given by Eq. (4.49). Differentiating in R g we get up to O(ε 2 ) corrections The first line is obtained from the same manipulations as in Eq. (4.54), while the explicit expression in the second line follows from differentiating the kernel Q q ε given in Eq. (4.48), which reduces to Eq. (B.51) when discarding NLO corrections. In order to compare with the coherent branching result, we again expand out this result in the limit Ω → 0. Using the expansion in Eq. (4.56), (4.60) Including the factor of m 2 J /Q 2 , according to our definition of C κ 2 from Eq. (1.5), this yields the LL coherent branching result of Eq. (4.53). Finally setting θ = θ g using the δ function and noting that we find that the arguments of the running coupling in the two expressions also match.

Matching the three regimes
Here we describe our strategy for consistently matching the results in various regimes so as to obtain a smooth interpolation across the entire phase space of m 2 J -R g measurement. Which of the three EFTs come into play will crucially depend on the range over which the groomed jet radius varies for given values of z cut , β and the jet mass. In the following we first describe the basic strategy for implementation of the matched cross section in Sec. 5.1 using profile functions for scale choices [67,68] and weight functions which control the contribution from different EFTs [27]. In Sec. 5.2 we generalize the weight function construction to depend on two variables that track where we are in our two dimensional phase space so that the appropriate factorized formula can be used. Next in Sec. 5.3 we describe the profile functions for each of the large, intermediate, and small R g regimes, and then in Sec. 5.4 how we vary them to obtain an estimate for the perturbative uncertainty.

Implementation of matched cross section
Here we describe our implementation of the complete cumulative-R g and differential jet mass soft drop cross section. This will then set up the groundwork for evaluation of the moments C κ 1,2 (m 2 J ). An important ingredient are the profile functions for our different EFT regions, which are the functions used for the renormalization µ i scales in the factorization formulae, and which depend on both m J and R g . First we recall the set of scales that are relevant for the various kinematic regions: Large R g region : Plain jet mass region : In this analysis we will focus exclusively on the groomed region described by our EFT, although we will also extrapolate our results into the ungroomed region. In the groomed region, we construct the interpolation between different EFT regimes via the following formula: where dΣ i (m 2 J ) for i ={large, int, small} stand for the cumulative-R g factorized cross section formulae in the large, intermediate and small R g regimes respectively. We have also indicated the choice of profile scales used for jet, global-soft and collinear-soft scales used in these factorization formulae utilizing the notation introduced in Eq. (5.1). In Eq. (5.2) the cross sections in the large and small R g regimes are evaluated with their respective profile scales, whereas the intermediate R g uses either the profile scales in Eq. (5.1) (in the two subtraction terms), or a set of hybrid profiles denoted by µ hyb . The hybrid µ hyb profiles are designed such that they interpolate between the µ small , µ int and µ large sets as follows: where the exponents 0 ≤ w i ≡ w i (m J , R g ) ≤ 1 should turn on/off in the appropriate regions and add up to one [27]. We give an explicit construction of these weight functions below for our case. For now, we note that in the large R g region we will have w large = 1 (and hence w small = w int = 0) such that the µ hyb profiles are the same as the set µ large . The term dΣ int | µ large in Eq. (5.2) corresponds to evaluating the intermediate R g cross section with large R g profile scale settings, i.e. by setting µ csg , µ csm → µ cs . This leads to exact cancellation of the terms dΣ int | µ hyb and dΣ int | µ large in the R g θ g region where w large = 1, and hence the matched cross section is given by The term dΣ int | µ small here denotes the intermediate R g cross section evaluated with the µ small profiles, i.e. by replacing µ csm , µ J → µ c . Thus the difference dΣ small − dΣ int | µ small supplements Eq. (5.4) with the power corrections that are missing in both the intermediate and large R g factorizations in the θ c R g region. Similarly, in the small R g region we will have w small = 1 (and thus w large = w int = 0). As a result, the µ hyb profiles will coincide with µ small and the terms dΣ int | µ hyb and dΣ int | µ small in Eq. (5.2) will cancel exactly. Hence, the dΣ matched in this region will be given by By using µ large profiles in the intermediate R g cross section we turn off the additional resummation relative to the large R g regime, so that the second term dΣ large − dΣ int | µ large captures the power corrections that are missing in the intermediate and small R g factorizations in the R g θ g region. Finally, if there is a contribution purely from the intermediate region where w int = 1, the µ hyb profiles will be the same as the bulk intermediate R g profiles in Eq. (5.1). Here the term dΣ int | µ hyb in Eq. (5.2) provides for the essential intermediate R g resummation, whereas the other two terms contribute the power corrections from the large and small R g regions. We will find that the intermediate regime only has a dominant weight for larger values of β, and even then in a rather small range of phase space. When the intermediate regime is not relevant the formalism we develop simply implements a smooth transition between the small and large R g EFTs with w int = 0 throughout the entire range of R g (for a fixed m J ). The smooth variation of w i (m 2 J , R g ) will ensure consistent interpolation between all these regions.

Weight functions
Here we describe the implementation of weight functions w i (m 2 J , R g ) that appear in Eq. (5.3). We first recall that the transition from the large R g EFT, R g θ g (m 2 J , Q cut , β), to the intermediate EFT corresponds to the canonical CS g and CS m scales becoming well separated. On the other hand, transitioning from the small R g EFT, θ c (m 2 J ) R g , to the intermediate EFT happens when R g θ c (m 2 J ), resulting from separation of the canonical CS m and the collinear scales. We can therefore characterize the regions with two power counting parameters defined as C κ factorizes into J κ and S κ cm : S κ c factorizes into S κ cg and S κ cm : Here the µ scales stand for the virtuality of the modes determined by {m 2 J , R g , Q, Q cut , β} such that which is consistent with Eqs. (3.25) and (3.22). 8 Note that λ min grows with R g while λ max decreases with R g , and in particular they become equal when R g = θ t , where . (5.8) Whether or not the additional resummation in the intermediate R g region is essential is given by the condition λ min (θ t ) = λ max (θ t ) 1. In our numerical analysis below we replace the strong inequality by a factor of 3, such that for intermediate EFT resummation to be significant we require When this condition is not satisfied, we directly match the cross section in region R g θ g to θ c R g . We refer to this as the "two-EFT scheme" since it does not include the intermediate R g regime. In the opposite situation, when the condition in Eq. (5.9) is valid, then the additional resummation in the intermediate R g region is needed, thus defining a "three-EFT" scheme. For a fixed value of {Q, Q cut , β}, this happens at smaller jet masses. In the three-EFT scheme we can further identify the two angles where respectively the small-R g and large-R g power corrections are equal to 1/3, This suggests that intermediate R g resummation is important for θ t,min ≤ R g ≤ θ t,max . Finally, we introduce a transition function which smoothly interpolates between 0 and 1, where r t controls the rate of the transition and is normalized with respect to the range of integration. In our numerical analysis in Sec. 6 we use r t = 20. This transition function is used to construct the weights w int and w large for the two-and three-EFT schemes: 2-EFT scheme: w large = X(R g , θ t ) , w int = 0 , (5.12) 3-EFT scheme: w large = X(R g , max{θ t,max , θ t }) , w int = 1 − X(R g , θ t,max ) X(R g , θ t,min ) , and w small = 1 − w large − w int in either of these cases. In the 2-EFT scheme, the large R g EFT is employed for R g > θ t and small R g for R g < θ t . Here the intermediate R g resummation is not present and the piece dΣ int | µ hyb in Eq. (5.2) serves mainly to provide an interpolation between these two regimes. In the 3-EFT scheme, the intermediate R g resummation is effective in the range θ t,min < R g < θ t,max . Note that w large in this scheme reduces to the corresponding weight in the 2-EFT scheme when Eq. (5.9) is valid with equality sign. The transition function X(R g , θ t,min ) in w int turns off the intermediate R g resummation for R g θ t,min . Furthermore, in the 2-EFT scheme, we have θ t,min > θ t,max which automatically drives w int to 0.

Profile functions
With the exception of the hard and global-soft functions, the logarithms that appear in the various matrix elements in the soft drop cross sections depend on the jet mass and groomed jet radius, and thus vary along the spectrum. Hence, to minimize these logarithms, the renormalization scales must also vary accordingly in the resummation region as discussed in Sec. 4.2. Additionally, for larger jet masses where R g R we enter the ungroomed region where the soft drop resummation must be turned off. For smaller jet masses, when we approach the nonperturbative region the soft scale must be frozen to an O(1) GeV value. This requires some care since the boundary of the nonperturbative region is a non-trivial function of both m J and R g measurements. Hence, in order to consistently describe the entire jet mass spectrum, we make use of profile functions [67,68] which depend on the jet mass and groomed jet radius. In this work we limit ourselves to the profile functions for e + e − collisions and leave the extension to the pp case to future work.

General strategy and canonical relations
In this section we first review the canonical relations between the scales that form the basis for implementation of the profile functions. To set up the notation we first consider a relatively simple scenario of plain jet mass resummation for a jet with radius R. The relevant canonical scales are the hard, jet, and the (ungroomed) soft scales: such that the endpoint of the spectrum is at m J = m ee max = Q tan(R/2), where the jet and the soft scales merge with the hard scale. Thus, the seesaw relation between the three scales is at the heart of the soft-collinear factorization in what is often called SCET I . When we include soft drop grooming, the canonical hard and the jet scales remain the same, whereas only the soft sector is modified, splitting into the global-soft and collinear-soft scales: We see that µ cs is independent of the jet radius R. These scales are relevant for the single differential jet mass distribution as well as in the large R g region. For the single differential jet mass case the canonical scales were first discussed in Ref. [6], with the generalization to include R dependence given in Refs. [8,49]. We note that the collinear-soft and global-soft scales are related by the following condition: This relation can be understood from RG consistency, noting that the z cut and β dependence must cancel within the soft sector, as jet and hard functions do not depend on these parameters. The relation in the second line is obtained using Eq. (5.14) and we will use it below in implementing the profile for our jet scale.
In the intermediate R g region, the collinear-soft sector further splits into CS g and CS m modes, with the canonical values of the corresponding scales being (5.17) They satisfy the following see-saw relation: which holds for all values of R g . This relation will be used below to define the µ csm profile. It is analogous to the relation in Eq. (5.16) as can be seen by setting R g = R. We will elaborate on this more below. Lastly, we consider the canonical scales in the small R g region. In addition to the µ gs and µ csg scales given above, the jet and the CS m scale in the intermediate region merge into a single collinear scale. As we saw in Sec. 3.3, the logarithms in the collinear function C κ are minimized for µ ∼ QR g /2, hence the canonical scale is given by 19) and the canonical relation between three scales in the small R g region reads (5.20) In deriving this relation we have replaced R g /2 → tan(R g /2), capturing some of the power corrections in R g /2 that are important in this region. Again, we will use the see-saw relation in Eq. (5.20) to fix the µ c profile. Furthermore, given the consistency between the small and intermediate regime, the following canonical relation also holds: This suggests that, unlike for groomed to ungroomed transition for the single differential jet mass distribution, there is a smooth transition between the intermediate and the ungroomed region owing to the additional R g measurement. Here the soft drop resummation is automatically turned off as µ csg approaches µ gs , whereas µ csm → µ s provides for the essential plain jet mass resummation in the ungroomed region. We stress, however, that the intermediate R g regime still lacks the power corrections of O(m 2 J /(QQ cut )) in the ungroomed region. Next, we note that upon decreasing R g the µ csm scale meets the jet scale as R g approaches the lower limit of the angle of the collinear mode: On the other hand, from Eq. (5.17) we see that µ csg keeps decreasing on reducing R g . Being the smallest of all the scales, for a given jet mass it becomes nonperturbative first for R g ∼ (R g ) NP defined in Eq. (2.8). Thus we will first implement the µ csg profile so as to consistently freeze it to an O(1GeV) value in the nonperturbative region.

Rescaled variables
To simplify the expressions we will work with the following unit-interval variables corresponding to jet mass and R g : Here we use µ can. N = Q tan(R/2) to keep it distinct from µ N which we will vary to study perturbative uncertainty. From Eq. (2.5) we see that for m 2 J > m 2 0 we enter the ungroomed resummation region. Here we turn off the soft-drop dependent resummation in the soft sector by merging the collinear-soft and global-soft scales to the ungroomed soft scale. In terms of ξ, the transition point is given byξ We note that as a result of the cos 2+β (R/2) the transition to ungroomed region happens slightly earlier than the value of ξ = ξ 0 for which the global-soft and collinear-soft scales merge: In the following we will, however, ignore the difference betweenξ 0 and ξ 0 and take ξ 0 in Eq. (5.28) to correspond to the groomed and ungroomed transition point. From Eq. (2.3) the phase space in the variables ξ and ψ is then given by 29) and is shown in Fig. 3. Here ψ (ξ) corresponds to ψ for R g = θ g (m 2 J , Q cut , β), and is given by For ξ > ξ 0 , i.e. in the ungroomed region, R g is bounded above by the jet radius R, which limits ψ ≤ 1.
We now summarize the canonical scales in terms of ξ and ψ variables: 1. Canonical scales for plain jet mass: 2. Canonical scales for large R g region (with µ N and µ J same as above): 3. Canonical scales for intermediate R g region (with µ N , µ J and µ gs same as above): 4. Canonical scales for the small R g region (with µ N and µ gs same as above): We will find the expression of µ cs in terms of ψ in Eq. (5.32) useful to construct the µ cs profile below.
Lastly, we note that in the small angle limit, we have ψ tan(R/2) = tan(R g /2) ∼ R g /2. Since using the variable ψ allows us to smoothly connect all the profiles from small to large values of R g , we also additionally replace all the appearances of R g /2 in the factorization formulae above by ψ tan(R/2). The difference between the two is a power correction, and when R g is O(1) it automatically results in the correct expressions. This includes all the fixed-order results for factorization functions in App. B, the appearances of R g 's in the logarithms (which are simply the canonical scales displayed above), the pre-factors in soft drop boundary cross section (see Eq. (4.49)), in the ratios υ = (R g /θ g ) 2+β that shows up in the resummation kernels in the large R g region (see Eqs. (4.25) and (4.48)), in the weight functions in Eq. (5.12), and as well as any non-logarithmic dependences on R g (such as in the collinear function in the small R g region in Eq. (4.5)).

Profile functions for double differential cross section
We now turn to implementing the profile functions based on these canonical scale choices and relations. Our goal is twofold: accounting for freezing of any perturbative scale that approaches Λ QCD and providing a way to estimate perturbative uncertainty though scale variation. We consider profile variations in the next section. We base the freezing behavior of all the scales on the the µ csg profile, defined as follows: In the resummation region µ csg (ψ) ∼ µ gs ψ 1+β , consistent with Eq. (5.17). The transition of the µ csg scale into the NP region is determined by where n 0 is an O(1) number that controls the value for the perturbative scale we freeze to which is > Λ QCD . Recall that it is mandatory for the scales to be frozen in this manner, since the entire formalism is based on perturbative expansions of things like the anomaous dimensions. In the numerical analysis presented in Sec. 6 we use n 0 = 0.75 which from Eq. (5.36) corresponds to onset of the freezing of the coupling at µ = 1.5 GeV. Note that µ csg has the lowest virtuality, so it is always the first scale to enter the NP region when the jet mass grows smaller. To freeze in turn the scales with larger virtuality as we venture into the NP region while maintaining the canonical scaling relations, we now simply insert Eq. (5.35) in the canonical relations derived above. First, as seen from Eq. (5.32), the collinear-soft scale can be derived using µ csg (ψ) by setting R g = θ g (m 2 J , Q cut , β). Thus we have We note that, compared to the plain jet-mass soft scale in Eq. (5.13), the collinear-soft scale should become nonperturbative at yet smaller jet masses in the SDNP region defined by Eq. (2.7). This which is consistent with Eq. (2.7).
Having defined the µ cs (ξ) and µ csg (ψ) profiles we can now use Eq. (5.18) to implement the µ csm profile as follows: while using Eqs. (5.14) and (5.16) we implement the profile for the jet scale as follows: Finally, using Eq. (5.20) the µ c (ψ) profile is given by We have so far defined the profile functions only in the groomed region for ξ < ξ 0 . A complete analysis of the distribution in the ungroomed region is beyond the scope of this work. We will, however, extend our results into the ungroomed region making use of the following extensions of the collinear-soft profile scales: Intermediate/Small R g : µ csg = µ gs , µ csm = µ s (ξ) = µ N ξ (ξ > ξ 0 ) . The large R g transition occurs when we cross the dot-dashed green line in Fig. 3 for values of ψ that place us near the upper solid magenta line, whereas the intermediate/small R g transitions occur when we cross the dot-dashed green line when we are far from all solid lines, or close to the lower blue line, respectively. Note that in the large R g regime the resummation related to jet-grooming must be immediately turned off as ξ > ξ 0 . This is the same as the groomed to ungroomed profile transition for the single differential jet mass case, which was discussed in Ref. [49]. In case of intermediate R g regime, however, this happens more smoothly owing to its more factorized nature. Finally, the jet and the collinear scales remain the same in the ungroomed region.
In Fig. 4 we show the profile scales resulting from the above formulas picking β = 1. In the left panel we show them as a function of the groomed jet radius for a fixed jet mass. The scales satisfy the joining condition for the end points of the R g spectrum: µ cs , µ csg , and µ csm merge at the end point R g = θ g (m 2 J , Q cut , β), whereas µ J and µ csm merge into µ c at the lower end point R g = θ c (m 2 J ). In the right panel of Fig. 4 we fix the R g to a representative value of R g = 0.3 and plot the scales as a function of the jet mass. Our choice of R g = 0.3 < R implies that the entire allowed mass range lies within the groomed region. Once again we see the expected behavior when the scales join. The profile functions for other values of β look qualitatively similar.
As we discussed above in Sec. 5.1, the implementation of the matched cumulative R g cross section involves constructing hybrid profiles µ hyb , that select the appropriate set of scales from each regime following Eq. (5.3). These involve the weight functions w i (m 2 J , R g ) that we constructed in Sec. 5.2. In Fig. 5 we show the hybrid profiles for β = 0 and β = 2 for a given jet mass value in the upper left and right panels respectively. In both these plots we show the scales only in the region where they are valid (except for the global-soft scale that is common to all). The µ c scale splits off into µ J and µ csm as we transition out of the small R g regime, and the µ csm and µ csg scales merge into the single collinear soft scale µ cs as we enter the large R g regime. In the bottom panels we show the weight functions that correspond to that specific choice of jet mass and and soft drop parameters. In the case of β = 0 we transition directly between the small and large R g regimes, thus providing an example of the 2-EFT scheme. In this case the scale µ csm only serves to provide an interpolation into the large R g region. On the other hand, we see that for β = 2 there is a significant region where w int > 0 and the intermediate R g resummation is active.

Perturbative uncertainty
We now briefly discuss our prescription to assess the perturbative uncertainty of the EFT results by means of scale variations and varying the unknown O(α 2 s ) non-logarithmic terms necessary for NLL resummation in the large and small R g regimes through the parameters a C 20 , a Sc 20 , and a Sc 20,ε that appear in Eqs. (4.5), (4.15) and (4.41) respectively. Note that we will not consider variations in the ungroomed region, leaving a more detailed analysis of this transition to future work, and hence we often cutoff our plots with uncertainty bands prior to hitting this transition.

Scale variations
In this analysis we will consider the following three types of scale variations: 1. Vary the hard and global-soft scales by a multiplicative factor, and let this variation propagate to other scales through the canonical relations: Here we will vary e gs from the default value of e gs = 1 to e gs = 1/2 and 2. Note that we vary the µ N and µ gs scales simultaneously. This is to ensure that the variation does not affect the groomed to ungroomed transition point ξ 0 defined in Eq. (5.28). This variation will affect all the scales in Eqs. . This is done such that the variations are frozen to the default value at the end point. More specifically, we first define the variation of the µ csg scale as µ vary csg (ψ, α) ≡ µ gs f vary ψ 1+β α f run ψ 1+β (5.46) and use it to derive variations of other scales. Here we make use of the following trumpet function f vary (ψ) to turn off the variation at the end point ψ = 1: The choice α = 0 returns the default profile choice, and varying α = ±1 allows for variation in the resummation region up to a factor of 2. In accordance with Eq. (5.38), the variation of µ cs scale is given by Having defined α-variations of µ csg and µ cs , the corresponding variations of µ J , µ csm and µ c scales are derived using Eqs. (5.41), (5.40) and (5.42) by replacing µ csg,cs → µ vary csg,cs .
3. A variation that relaxes the canonical see-saw relations that we used to derive the µ csm , µ J , and µ c scales. The two sets of canonical relations, Eqs. (5.16) and (5.18), and Eqs. (5.14) and (5.21), are a consequence of soft drop and cumulative-R g constraints respectively. The two sets are equivalent when R g → R, where both constraints are lifted. For convenience, we define an auxiliary plain-jet mass soft scale µ s (ξ, α) using Eq. (5.16): We now parameterize the deviations in the µ csm scale upon modifying the soft-drop see-saw relations in Eqs. (5.16) and (5.18) with a small exponent ρ: Solving for µ csm yields the variation Next, we apply the same strategy to the µ J and µ csm scales by entering deviations in Eqs. (5.14) and (5.21) through a small exponent γ: Note that we have set ρ = 0 in Eq. (5.50) as these two variations are independent. Solving for the scales µ J and µ c yields the variations

Parameterizing variation of the O(α 2
s ) non-logarithmic terms We now consider the variations of the two-loop parameters a C 20 , a Sc 20 , and a Sc 20,ε introduced in Eqs. (4.5), (4.15) and (4.41) respectively in order to assess uncertainty of our lack of knowledge of these terms. We first note that in combining the cross section through the recipe in Eq. (5.2) we are relying on the fact that the intermediate R g cross section reproduces the EFT result on either end of the R g spectrum up to power corrections. This is guaranteed to hold from Eqs. (3.25) and (3.22) when all the functions are expanded to a given order. Given our complete knowledge of O(α s ) pieces, we can use the construction in Eq. (5.2) to NLL accuracy. However, when we include O(α 2 s ) pieces in the small and large R g regimes to account for NLL resummation, the missing non-logarithmic pieces parameterized by a i 20 spoil this delicate cancellation when they extrapolated into the intermediate regime. We solve this issue by explicitly multiplying the O(α 2 s ) fixed-order pieces in the large and small R g regimes by the weight factors w i (m 2 J , R g ) so that they are set to zero when extrapolated into other regimes: dΣ small ≡ dΣ [1] small (a C 10 ) + w small (m 2 J , R g ) dΣ [2] small (a C 20 ) , (m J > 0) , large + dΣ [1] large (a Sc 10 ) + w large (m 2 J , R g ) dΣ [2] large (a Sc 20 ) , large term in the large R g regime is simply the single differential jet mass cross section, and does not contribute to the differential R g measurement. By demanding m J > 0, the O(α 0 s ) δ-function piece in the small R g regime does not contribute. The O(α s ) terms depend on the LO expressions in Eqs. (4.6) and (4.27) in the large and small R g regimes respectively. The O(α 2 s ) terms are the additional pieces required to achieve NLL accuracy which include the cross terms shown in Tab. 1, as well as the parameters a C 20 and a Sc 20 that are included for uncertainty estimation. Thus, we achieve the reasonable outcome that the variations induced by these parameters will be limited to the region where the corresponding weight function Having constructed the matched cumulative R g cross section, we will use it directly to evaluate C κ 1 (m 2 J ). The uncertainties resulting from scale variations and 2-loop parameter variations will then be directly translated to uncertainties in predictions for C κ 1 (m 2 J ). For computing C κ 2 (m 2 J ) we will, however, only make use of the large R g boundary cross section at the differential level. Since we have no contribution from the small and intermediate R g regime, we adopt the following prescription following Eq. (5.53): dσ [1] (ε) dm 2 J dθ g θg∼θ g + dσ [2] (ε) dm 2 J dθ g θg∼θ g . (5.54) From this equation we see explicitly that the cutoff on the large R g cross section appearing in the determination of C q 2 is directly determined by w large .

Numerical results and comparison with Monte Carlo simulations
Having set up the formalism, provided the method to match across different regimes and described our prescription for quantifying the perturbative uncertainty, we are now in position to present the numerical results for the matched cumulative-R g cross section and the moments C q 1 (m 2 J ) and C q 2 (m 2 J ). We build our c++ code based on the core of the SCETlib package [69], extending it for soft drop observables, and cross checked the numerical results using Mathematica. We will first present the results for the matched double differential cross section in Sec. 6.1, and then use it to calculate C q 1 (m 2 J ) in Sec. 6.2 and C q 2 (m 2 J ) in Sec. 6.3 in the manner described in earlier sections. For definiteness, we fix z cut = 0.1, E J = 500 GeV and R = 1 throughout the whole analysis, but we sample different values of the soft drop angular weight β. Along with the numerical results for the moments we will also present a comparison with parton shower Monte Carlo predictions. Figure 6: The matching prescription smoothly interpolating between the three EFT regimes is shown for the distribution cumulative in R g , at NLL accuracy, for different values of β and a fixed value of m J . The large R g (red dashed) and small R g (green dashed) predictions include fixed-order corrections that are relevant in the respective plotted range; the profile-improved intermediate prediction (blue, dot-dashed) does not yet include these corrections while resumming large logarithms in the central region. We obtain the matched result (solid black) by weighting the three curves with the exponents w i (R g ) shown in the bottom panels, and subtracting the overlap (grey) as described in the text.

Results for the matched cross section
Following the strategy outlined in Sec. 5.1, we now utilize the profile and transition functions described in Sec. 5.2 and Sec. 5.3 to smoothly match the cross section across all our EFT regimes. Here we consider explicitly the cross section cumulative in R g which determines C q 1 as well as the normalization factor C q 0 it depends on, whereas C q 2 will be calculated at the differential level using the prescription in Eq. (5.54). Since our focus is on making predictions for C q 1,2 , where NGL effects turn out to be small, we will not include the effects of the NGL functions Ξ q G , Ξ q S , and Ξ q C in our results for the double differential cross sections (despite the fact that the resummation of NGL effects may be important). We leave inclusion of NGL effects in the double differential cross section to future work. In Fig. 6 we show the various components that enter the calculation of the matched cross section for β = 0, 1, and 2 plotted as a function of R g for a representative jet mass value of log 10 (m 2 J /E 2 J ) = −1.7. In order to explain the legend and the meaning of each curve, consider the first panel with β = 0. The solid black line is the matched cross section that is evaluated after combining the cross sections in the three regimes according to Eq. (5.2). The results for cross section in the large R g EFT, dΣ large | µ large , and in the small R g EFT, dΣ small | µ small , evaluated with their corresponding profiles are shown by the red dashed and green dashed curves respectively. The overlap of the intermediate R g with either of these regimes are shown in gray-solid for dΣ int | µ large (labeled as "Large R g w/o O (R g /θ g ) 2+β "), and gray-dotted for dΣ int | µ small (labeled as "Small R g w/o O θ c /R g )"). We see that the dΣ int | µ large term matches with the dΣ large | µ large term as R g → 0, such that the power correction O (R g /θ g ) 2+β becomes small. Likewise, the dΣ int | µ small cross section merges with the small R g result as R g → θ g , up to the O θ 2 c /R 2 g ) power correction. Finally, the blue-dot-dashed curve that spans the entire range of R g values is the interpolation Σ int | µ hyb that make use of the hybrid profiles described in Eq. (5.3). The hybrid profiles are constructed to reproduce the large-R g and small-R g profiles in their respective regions. The weight functions w i (m 2 J , R g ) given in Eq. (5.12) demarcate the boundary of each region and govern the behavior of the hybrid profiles. Thus, we see that Σ int | µ hyb matches exactly with the dΣ int | µ large (gray-solid) in the large R g regime and dΣ int | µ small (gray-dotted) in the small-R g regime.
As mentioned in Sec. 5.2, for the cases β = 0, 1 the intermediate R g resummation is not present (for the jet mass value shown in Fig. 6) as there is simply no room for this EFT to be valid according to Eq. (5.9), whereas for β = 2 there is a significant region where the intermediate R g EFT is valid. For each of the three β values that we show, we see that our choice of the weight functions corresponds with the size of the power corrections as determined by difference of the large and small R g cross sections with the corresponding intermediate R g overlaps.
Next we show in Fig. 7 our estimate of perturbative uncertainty on the NLL cumulant cross section, again for β = 0, 1, and 2. The uncertainty bands in the top left plot are obtained by varying the profile scales as described in Sec. 5.4 (magenta) and by varying the quantities a C 20 and a Sc 20 that parametrize our ignorance of the two-loop constants terms in respectively the small R g and large R g EFTs (green). Specifically, we vary both a Sc 20 and a C 20 in the range [−2π, 2π], and take the envelope of the two independent variations. We notice that scale variation is the dominant source of uncertainty at large R g , while at small R g the perturbative uncertainty is mainly set by the a C 20 term. As discussed below Eq. (4.6), the lower endpoint of the distribution receives perturbative corrections starting at NLO. It is clear from Eq. (4.5) that this effect is captured by variations of a C 20 , resulting in a non-vanishing error band below the endpoint of the central value. This effect has no equivalent in the scale variations, where by construction we set a C 20 = 0 and use the same LO endpoint as for the central value.
The top right plot in Fig. 7 shows for β = 0 the breakdown of the scale uncertainty into the three independent variations listed in Sec. 5.4: overall scale uncertainty obtained by simultaneously varying the hard and soft scales by a factor 2 (red band); variations of collinear-soft scales by a trumpet factor, setting α = ±1 in Eqs. is obtained with the matching prescription described in the text, while the colored bands probe respectively scale variation (magenta), and variation of the fixed order two-loop constants a Sc 20 and a C 20 , collectively denoted as a i 20 variation (green). In the top right plot, the envelope of scale variations for β = 0 is broken down into various components: overall variation by a common factor (red), variation by a trumpet factor preserving canonical relations (yellow), breaking of canonical relations (green). and (5.51) (green band). We see that the three sources of uncertainty have similar size, with the trumpet factor playing the most relevant role. A similar qualitative behavior occurs for β = 1, 2, so for simplicity we do not plot the breakdown for these cases, and only show the envelope of scale variations and the fixed-order uncertainty due to the parameters {a i 20 } in the lower panels of Fig. 7.
Before moving on to the numerical results for C q 1 and C q 2 , we anticipate that the theoretical uncertainty for these two observables will look rather different than the one just shown for the double differential cumulative cross section. In particular, the part of uncertainty that only depends on the jet mass but not on R g will cancel in the ratio of the angular moments in Eq. (1.5) and the cumulant distribution, since effectively this means they are normalized point by point in the m 2 J spectrum. We anticipate that this will lead to smaller uncertainties for both the C q 1 and C q 2 coefficients. We now show predictions for C q 1 , based on the NLL expressions for M 1 obtained by evaluating Eq. (1.4) on the resummed cross section described in Sec. 4 and shown in Sec. 6.1. We then compare these EFT results against Monte Carlo simulations and numerical results from coherent branching. We explicitly present predictions for quark jets, expecting gluon jets to follow the same qualitative behavior. We set again E J = 500 GeV, R = 1, and z cut = 0.1, and explore the three cases β = 0, 1, 2 across the mass range −3.5 < log 10 (m 2 J /E 2 J ) < −0.5. At low jet masses, the double differential distribution on which we base our predictions itself becomes nonperturbative. However, these NP effect should not be included in C q 1 , which is by definition a purely perturbative coefficient. What portion of the jet mass range is significantly affected by NP corrections crucially depends on the value of β. Our method of calculating C q 1 automatically deals with this as discussed in Sec. 5.3.
In Fig. 8 we present central values for our EFT results for C q 1 at both LL and NLL orders. These results are compared to results from coherent branching and from Monte Carlo simulations. For the latter, we consider predictions from Pythia 8.2 [70], Vincia 2.2 [71] and Herwig 7.1 [72], where the jet reconstruction and soft drop are performed using FastJet 3.3 [73]. We signal the region where NP effects are relevant by showing the transition from SDNP to SDOE region, as Figure 9: Error estimate for the NLL and LL effective theory predictions for the perturbative coefficient C 1 . We do not estimate LL uncertainty in the ungroomed region, where the EFT description loses validity. NLL error bands are tiny and examined in more detail in Fig. 10. defined in Eq. (1.1). Note that Pythia, Herwig and Vincia return very similar results for C q 1 [24]. At intermediate and large masses the NLL predictions differ significantly from the ones from coherent branching, in better agreement with the Monte Carlo, while the LL EFT prediction lies in-between these two results. In the SDNP region the NLL and LL coherent branching curves agree very well with each other, while the LL EFT result shows some deviation. The fact that our LL EFT predictions differ from the LL coherent branching derived in [24] is not too surprising: the EFT description is more refined in that it involves matching the three regimes via profile scales and including additional sources of resummation, but (as we discussed in Sec. 4.5) the LL-EFT result lacks some of the corrections that are present in the splitting kernel p gq (z). These power corrections have a small effect for C q 2 , thus the difference is driven by the additional resummation included in the LL EFT results.
In Fig. 9 we reconsider the EFT predictions, now including also our estimate of their perturbative uncertainties (repeating the Monte Carlo simulation results for reference). While the LL uncertainty band (blue) is obtained from scale variation, the NLL uncertainty band (orange) is obtained from the envelope of scale variation and variation of the two-loop fixed-order constants {a i 20 }. Since our estimates for the perturbative uncertainty does not cover the ungroomed region, we exclude values m J > m ee 0 in the plotted range. We observe good convergence between the two orders, with only the β = 2 results being slightly outside the LL uncertainty band. The most striking feature is the very small size of the NLL uncertainty band, whose width becomes notice- able only for β = 2 (low panel). Given that for the NLL cumulant cross section the variations are rather large (see Fig. 7), the tiny uncertainty bands for C q 1 are due to cancellations between numerator and denominator of Eq. (3.27). Indeed, we know that the leading double logarithms cancel in the ratio to all orders in perturbation theory.
To further study the small NLL theoretical uncertainty, in Fig. 10 we show the relative difference of the variations to the central curve. As we did for the cumulant in Fig. 7, we separate scale variation (magenta) and the variation of the two-loop constants (green). Again, we do not show the perturbative uncertainty beyond the groomed region. We first notice that in the SDOE region the uncertainty, of a few percent, is dominated by the unknown two-loop term. Scale variation gradually takes over in the SDNP region, which becomes larger as β increases (for β = 2, scale variation grows up to ∼ ±15% at the lowest jet mass we consider, but for clarity we limit the range shown on the vertical axis to ±10%). For β = 1, 2 the enlarged uncertainty at high jet mass comes from varying the two-loop constant a C 20 in the small R g region, whereas variations of the large R g term a Sc 20 are roughly constant across the mass range, and effectively negligible. The situation is different for β = 0, where varying a Sc 20 gives the dominant contribution. That the variations from the small R g region are less relevant for this β was already clear from Fig. 10. In the top right plot of Fig. 7 we further break down the scale variations for the β = 0 case into the three different sources, where the same descriptions made for Fig. 7 apply.
In Fig. 10 we also show the shift to the central curve obtained when including the leading O(α 2 s ) NGLs and abelian clustering logarithms (black dashed). We observe that the leading NGLs always lie within the small error bands, and in particular are well within our NLL uncertainty bands in the perturbative SDOE region. This justifies our choice of not resumming these classes of logarithms, treating them as the same order as other (missing) two-loop contributions. For this reason we have neglected them from our central predictions for C q 1 . Finally, we tested the dependence of our results on the choice of weight functions discussed in Sec. 5.2. Specifically, we modified the inequality in Eq. (5.9) to ≥ 1 2 and ≥ 1 4 (thus varying the range where the intermediate R g EFT is included) and changed the transition rate r t between EFTs in Eq. (5.11) by a factor 2. We find that in both cases the variations lie within our uncertainty bands.

Results for C q 2
Based on the definition in Eq. (1.5), we obtain numerical results for C q 2 (m 2 J ) by using only the large R g EFT, setting up the calculation at the differential level, and weighting the large R g cross section according to Eq. (5.54). We consider the same kinematical values as for C q 1 , and plot the comparison between our LL and NLL EFT results, the LL coherent branching prediction, and MC data in Fig. 11. For this observable the spread between different MC results is much larger, with Herwig predictions being systematically lower than Pythia and Vincia. We see that our calculations are in good agreement with the MC in the region of validity of the EFT, with a preference for the Herwig results for larger β. Unlike the treatment in [24], here we do not make the small jet radius approximation in the coherent branching LL predictions, i.e. we distinguish tan(R g /2) = R g /2. This ensures that the coherent branching and EFT calculations both have their transitions to the ungroomed region at the same value of the jet mass. Although the impact of these corrections for C q 1 were minimal, here they induce considerable corrections to C q 2 at large mass, shifting the peak of the distribution towards larger values and improving the agreement with MC data. This happens because the boundary cross section, on which the calculation of C q 2 is based, is highly peaked at larger R g , where the R g /2 1 approximation is less justified. We stress, however, that a more accurate description of this region would involve a an extension of our EFT to account for power corrections O(m 2 J /(QQ cut )) that become relevant at the groomed-ungroomed transition point, m (ee) 0 defined in Eq. (2.5), which we do not attempt here. We further notice in Fig. 11 that the LL coherent branching result is systematically larger than the LL EFT result. Compared to NLL , the LL coherent branching result is similar for β = 0 but systematically larger for β = 1, 2.
As noted in Sec. 4.5, our LL EFT predictions include extra resummation in Laplace space, which is absent in the LL coherent branching results, while the latter include some additional power corrections through the splitting function p gq . We find that at the large jet masses (log 10 (m 2 J /E 2 J ) −1.5) the disagreement is almost entirely due to the power corrections, while at smaller jet masses the additional evolution in the EFT becomes more important. Finally, in the bottom right plot we show the effect of including the extra anomalous dimension γ 0 (ε, z cut ) from Eq. (4.37), which is present only in the case β = 0. The NLL result is shown before (dashed magenta) and after (red) including the fixed-order logarithms and additional running induced by γ 0 . For more clarity, we zoom on the range [−3.5, −1.5] in the logarithmic jet mass variable, leaving out the ungroomed region. We observe that this contribution is important across the whole plotted range, and brings the calculations much closer to the MC predictions.
We also show the LL and NLL EFT results with theoretical uncertainties in figure Fig. 12, again restricting to the region m J < m ee 0 where grooming is active. The NLL errors are estimated again by combining scale variation with variations of the two-loop non-logarithmic term via a Sc,ε 20 introduced in Eq. (4.41). We see that the error bands shrink from LL to NLL , but unlike C q 1 , their size is still significant at NLL . For this reason we do not show relative variations here. For β = 0 the NLL curves have excellent agreement with the MC data, while β = 1 and β = 2 show some tension, especially at large jet masses, which can likely be attributed to power corrections that we have not added to the EFT results, as inferred from the comparison with the LL coherent branching results for β = 1, 2 in Fig. 11. Note that since the C q 2 calculation is carried out with the large R g cross section, that the NGL and clustering effects present for the corresponding double differential cross section only modify the normalization, and hence cancel out in the ratio of moments that determine C q 2 . While the NLL uncertainty in Fig. 12 has roughly the same size for all β, the LL band becomes larger as β decreases. To investigate this further, we show the breakdown of variations for β = 0 in the top right panel of Fig. 12. The LL uncertainty, entirely determined by scale variation, is once more decomposed into the variation of the collinear-soft scale by a trumpet factor (yellow) and the breaking of the canonical relations (green). We find that the overall variation of the hard and soft scales is irrelevant for C q 2 , and hence have not displayed it here. In the same plot, the NLL band is decomposed into scale variation (red) and variation of the two-loop constant a Sc,ε 20 (blue), where we vary a Sc,ε 20 ∈ [−π, π]. We see that the variation of the two-loop constant dominates over scale variation across the whole range in the jet mass. This is because the denominator in Eq. (1.5) does not involve a Sc,ε 20 (which parametrizes uncertainty in the boundary cross section), and hence unlike other variations does not have a cancellation in the ratio of moments used to compute C q 2 . We also checked that variations of the weight functions through the parameters in Eqs. (5.9) and (5.11) have a small impact compared to the theoretical error. (Note that even though C q 2 uses only the large R g regime, it still involves the weight function as in Eq. (5.54).)

Conclusion
In this paper we have set up an effective field theory description of processes involving jets which are groomed using soft drop, where the jet mass and groomed jet radius are simultaneously measured. Focusing on dijet production in e + e − collisions, we presented factorization formulae formulated using SCET in order to resum large logarithms in the cross section which is differential in these two variables. We also combined the three EFT regimes that span the phase space of the double differential measurement to obtain a coherent description of the double differential cross section across the whole range of kinematics where grooming is active.
As a first application of this formalism, we computed two angular moments of the double differential cross section: the average value of the groomed jet radius, M q 1 (m 2 J ), and the average value of the inverse groomed jet radius for emissions at the soft drop boundary, M κ −1 (m 2 J ). Apart from being interesting observables in their own right, they are closely related to the Wilson coefficients C q 1 (m 2 J ) and C q 2 (m 2 J ) which describe the leading nonperturbative corrections to the soft drop jet mass cross section [24]. By making the correspondence C q 1 M q 1 and C q 2 M q −1 , we use our NLL calculation of the moments to determine these Wilson coefficients with higher precision than previously available. This leads in turn to a more accurate description of the leading nonperturbative corrections to the groomed jet mass.
We determined the regions of validity of three different EFT regimes based on their value of the groomed jet radius R g , combined them using profile scales, and estimated the theoretical uncertainties on the double differential cross section and its angular moments via scale variation. The NLL improved computation involved, in certain regimes, fixed order ingredients at O(α 2 s ) accuracy. We parameterized unknown 2-loop terms that are not predicted by the renormalization group equations, and included them in our estimate of perturbative uncertainty. We also compared our predictions with parton shower Monte Carlo simulations and found that the two agree within uncertainty in the region of validity of the formalism. Although we did not resum the nonglobal logarithms in our factorization formulae, we showed that they cancel in C q 2 . We computed the leading effects on C q 1 finding them to be not larger than other missing O(α 2 s ) effects included in our uncertainty estimate.
While we have focused on the application of simultaneous measurement of the groomed jet radius and the jet mass, we expect the EFTs explored in this work to have a much wider applicability, and our work to pave the way to analogous computations for other groomed jet observables. Single differential jet based observables, such as the jet mass and jet angularities, have exciting prospects for measurements of standard model parameters such as the strong coupling constant α s and the top mass m t , with grooming playing a key role in reducing contamination in jets. Having a double differential cross section prediction at hand offers complementary ways of extracting these parameters.
The double differential cross section includes a resummation of large logarithms of the groomed jet radius, and this resummation remains relevant even when m J is taken to be in the ungroomed region. The study of power corrections for the multi-differential cross section in this region are also rather interesting. In the future work, we envision extending our double differential predictions to the ungroomed region, to achieve a more complete description of the double differential cross section.
One of the greatest advantages of jet grooming is the decrease to the impact of nonperturbative corrections compared to those in ungroomed jets. Grooming also significantly reduces contamination from multi-parton interactions, and is thus an invaluable tool for carrying out phenomenology with jets, and for analyzing heavy ion collisions. In hadronic collisions, measuring observables on groomed jets also allows for much cleaner access to the initial state hadron struc-ture, which is one of the principal goals of nuclear physics. By bringing under control the small nonperturbative effects from final state radiation, we are one step closer to the accurate extraction of hadron structure functions such as the transvere momentum dependent parton distributions. Additional natural areas for future explorations of multi-differential groomed jets include for transverse momentum observables and for the phenomenology of jets in heavy ion collisions.

A.1 Laplace transforms and RG Evolution
Here we present the details of the resummation of large logarithms via RG evolution. Our notation follows that of Ref. [74]. We consider a generic MS renormalized SCET matrix element F(q, µ), with j F the mass dimensions of q. When F(q, µ) enters factorization theorems in a Laplace convolution over the momentum space variable q, its renormalization group equation (RGE) reads where γ F is the anomalous dimension in momentum space. In SCET such functions are ubiquitous and result from the fact that UV renormalization of soft and collinear matrix elements accounts for double poles in the dimensional regulator, that are in turn reminiscent of the double IR poles of the corresponding QCD matrix elements. The presence of double poles in counterterms, and the convolution nature of the factorization theorem, lead to nontrivial running such as the one shown in Eq. (A.1). To solve of Eq. (A.1), we consider the Laplace transform of F given by such that x has mass dimensions −j F . The formula for inverse Laplace transform reads where the real number c determines a vertical contour in the complex plane that leaves any pole of F to its left. In position space, the RGE is simply multiplicative: where the anomalous dimension 9 γ F (x, µ) consists of a cusp piece Γ F [α s ], responsible for resumming double logarithms, and a non-cusp part γ F [α s ]. The cusp part, Γ F [α s ], is proportional to the cusp-anomalous dimension such that has a formal solution in terms of the evolution operator The evolution kernels K F and ω F are obtained from integrals of the anomalous dimensions, involving the QCD beta function β[α s ], whose perturbative expansion we give in Eq. (A.13). We also state the results for functions F (µ, µ F ) that enter the factorized cross section as multiplicative factors (and not through convolutions) with µ F being the energy scale specific to F . Here the RGE is directly given by with the solution for evolution from scale µ 0 to µ being where the evolution kernels K F and ω F are again given by Eq. (A.7). The kernels in Eq. (A.7) depend on the specific function F (or F ) only through the anomalous dimensions: we can thus rewrite them in terms of two universal kernels whose explicit NLL expressions are where r = α s (µ)/α s (µ 0 ). Here the perturbative expansion of a generic anomalous dimension Γ[α s ] and β[α s ] is written as and the terms relevant for NLL resummation read

A.2 One-loop anomalous dimensions
We state here one-loop results for the anomalous dimensions of the various functions appearing in Eq.
cut , β, µ) = and the inverse Laplace transforms for J κ , S κ cm and S κ c are defined as follows, such that the mass dimensions of x, y, and z are respectively The cusp pieces in Eq. (A.16) are multiples of the universal cusp anomalous dimension Γ cusp [α s ] and are given by 19) and the one-loop noncusp pieces are We can now check the RG consistency relation for the factorized cross section for doubly differential m 2 J and R g measurements. Written in position space, the dijet cross section for jet masses is given in the intermediate regime by Eq. (4.7). By differentiating with respect to the renormalization scale, it is immediate to show that the sum of the cusp terms must vanish, and that the noncusp pieces must also add to zero from RG consistency, as can be easily verified at one loop from Eq. (A.20). This puts nontrivial constraint on the β dependence of the global-soft and collinear-soft non cusp anomalous dimensions. As a consequence of the refactorization in Eq. (3.22) in the regime of large groomed jet radius, the relation is also enforced by RG consistency.

B One-loop results
In this appendix we collect the one-loop results for the ingredients appearing in the double differential cross section. We start from the functions that enter the factorized expression for the cross section cumulant in R g and describe their evolution, then we move on to describing the ingredients of the boundary cross section that enters the definition of C q 2 .
B.1 Results for the cumulant distribution The global soft function accounts for the soft radiation which is groomed away, namely lies within the original jet but fails the soft drop test. The one loop correction for e + e − collisions requires computing where C q = C F and C g = C A , the condition for failing soft drop is with sin 2 (θ p /2) = p + /(p + + p − ), and Θ R constrains radiation within a jet of radius R, Solving the integration in Eq. (B.1) we get Expanding for small jet radius, R 1, one retrieves the familiar expression

B.1.2 Jet function
The one-loop result for the jet function is well known [62,63] and given by In the notation of Eq. cut , β, µ = α s C κ π (µ 2 e γ E ) Γ(1 − ) dp + dp − (p + p − ) 1+ Θ yielding the result presented in Eq. (3.11). Given our multiplicative treatment of this function, see Eq. (4.12), we need to separate the terms that contribute to RG evolution from the ones that depend on R g . The R g independent function coincides with the collinear-soft function for single differential jet mass distribution. Its one-loop result in the MS scheme is Likewise, Laplace transforms of other pieces appearing in Eq. (4.17) for NLL resummation in the large R g region are also defined with respect to the momentum variable + following Eq. (4. 16), such that at one-loop we find where the incomplete Gamma function is defined by In our alternative Laplace space notation, we havẽ S κ c L cs , β, α(µ cs ) = 1 + α s (µ cs )C κ π − 1 + β 2 + β L 2 cs − , β, α s (µ cs ) (B.33) where we can equivalently write

B.3 Results for the boundary cross section
We now provide results and derivations for boundary corrections to soft matrix elements. We first consider the case of global-soft function. We then present results for the boundary corrections to the S cg function in, that would be needed to compute the boundary cross section in the small and intermediate R g regions. Finally, we give an explicit result for the kernel Q κ ε defined in Eq. (4.48).

B.3.1 Global-soft function with shifted soft drop
Here we derive the one-loop expression of the global-soft function with shifted soft drop condition, following the analogous derivation preformed for the collinear-soft function in Sec. 4.4. In analogy with Eq. (4.31) we can write S κ [1],bare G Q cut , R, β, Θ sd (ε), µ = α s C κ π (µ 2 e γ E ) Γ(1 − ) dp + dp − (p + p − ) 1+ Θ sd (ε)Θ R = S We first consider the case of β = 0 for which the result is given by S κ [1],bare G,ε Q cut , R, β = 0, µ + virtual term = α s C κ π 1 UV + 2 log µ Q cut tan R 2 . (B.39) Note that we only considered the real emission term in the measurement function in Eq. (B.35). We then demand that the virtual term, which we have not explicitly included, appropriately turns the IR pole in Eq. (B.39) into a UV pole.
Next we consider the case of β > 0. Here, we note that the integral is power law IR divergent and will be set to zero in dimensional regularization. This results from the fact that we are trying to factorize the soft function at subleading power but without considering the full set of operators that appear at this order. Hence, in order to resolve this ambiguity it is helpful to consider the combined measurement resulting from the shifted soft drop condition in Eq. (B.36), the R g constraint in Eq. (4.29) (without small-angle approximation), the jet radius constraint in Eq. (B.3), along with p + measurement applied to the matrix element defined by a single soft mode (i.e. not factorizing it into a global-soft and a collinear-soft mode). In this case, the measurement function is given by δ full p + , + , R g ≡ Θ R Θ sd (ε)Θ Rg δ( + − p + ) + Θ sd (ε)δ( + ) + Θ R δ( + ) − δ( + ) = Θ R Θ sd (ε) Θ Rg δ( + − p + ) − δ( + ) , (B.40) where the three terms in the first line correspond to emissions that pass the jet radius constraint, that fail and the virtual term. The ones that are clustered in the jet radius are then tested for (shifted) soft drop condition and the R g constraint as in Eq. (4.28). Writing Θ sd (ε) = Θ sd + εδ sd and using the δ function to fix p − to Eq. (B.37), we find ∞ 0 dp − Θ R Θ Rg δ( The first term in the second line when expanded in small-angle approximation matches with the corresponding term in the measurement function that we used above for the + -differential collinear-soft function S c in Eq. (4.28). The second term proportional to δ( + ) is of interest to us and when integrated over p + leads precisely to the bare O(ε) correction in Eq. (B.38) but now with p + integral bounded below due to R g measurement, which then renders it IR finite. This R g -dependent term in the factorized cross section is captured in the O(ε) modification to the collinear-soft functions: the + -differential S c , as seen in Eq. (4.32), or equivalently in the cumulative-R g S κ cg function if we are considering intermediate R g regime as we show below. Thus, we can consider the integration in Eq. (B.38) including the lower bound on p + due to R g measurement, and identify the pieces depending on R and R g as parts of the global-soft and the collinear-soft functions respectively. We then arrive at the following result for β > 0: The result for S where as in the case of R g -dependent collinear-soft function in Eq. (4.13) we have included a β 0 term to cancel the µ dependence due to the running coupling.

B.3.2 Boundary Corrections for intermediate and small R g
Here we discuss the soft drop boundary cross section in the intermediate and small R g regimes.
In both of the regimes where R g θ g it is the S κ cg function that receives O(ε) corrections upon shift to the soft drop condition. This allows to treat the two cases where on equal footing and results in considerable simplifications, since none of the boundary corrections now depend on the jet mass. In this work we only provide the details of the boundary corrections in the soft sector in these regimes and leave further studies of the corresponding resummation to future work.
As commented below Eq. (B.42), the results for boundary corrections to S κ cg are closely related to the corresponding corrections to the global-soft function. Following the reasoning described there, we find , Again, as required by consistency with RGE and the large R g case, for β = 0 we find a UV pole, which cancels against the global-soft function in Eq. (4.36). Hence, we have an additional running between the scales µ gs and µ csg resulting from the non-cusp anomalous dimension cut , Θ sd (ε), β, µ = S κ cg R g Q 1 1+β cut , β, µ − 2δ β,0 Qε Q cut α s (µ)C κ π L csg (B.46) cut , β, µ ∆S κ cg,ε R g Q The factor 4 9 is the final coefficient taking into account clustering effects to the NGLs for the single differential distribution, and we find that it is unaltered by the jet mass measurement.
The leading NGLs in the other two regions can be obtained from the intermediate R g expression in Eq. (C.3) by taking appropriate limits. For the large groomed jet radius regime R g = θ g and the NGL contribution goes to 0 as was expected. For the small groomed jet radius regime we can formally set θ c = R g and the dominant contribution in this EFT is which is larger than the contribution for the case of intermediate groomed jet radius. Finally, we can switch back from Eq. (C.3) to distribution space, finding This function enters the double differential cross section in a convolution with the other jet mass dependent ingredients in Eq. (3.17). The Laplace transform of Ξ q S (r, β) is defined as follows: Ξ q S (q, β) ≡ To get a rough estimate of the impact of these leading NGLs on the double differential distribution, we find it more convenient to replace the convolution with the cumulative value in Eq. (C.4), such that which gives an upper bound to the absolute size of the O(α s ) contribution.

C.2 Abelian clustering logarithms
In this section we account for the leading logarithms arising from clustering effects of the C/A algorithm while implementing soft drop for independent emissions. To this aim, we first consider a naive implementation of the soft drop algorithm where clustering logarithms are absent and see how clustering effects modify this result. The naive implementation simply tests for the SD condition for each particle independently, i.e. the allowed phase space for each emission is Θ(z i − z cut θ β i )Θ(R g − θ i ) + Θ(z cut θ β i − z i ) . (C.10) We also note that in this calculation we have not considered the scenario when both partons naively fail the SD condition but pass if clustered together. This is evaluated in [6] and would lead to a sub-leading single logarithmic contribution α 2 s log X compared to the correction explicitly considered here.