Two- and Three-Loop Data for Groomed Jet Mass

We discuss the status of resummation of large logarithmic contributions to groomed event shapes of hadronic final states in electron-positron annihilation. We identify the missing ingredients needed for next-to-next-to-next-to-leading logarithmic (NNNLL) resummation of the mMDT groomed jet mass in $e^+e^-$ collisions: the low-scale collinear-soft constants at two-loop accuracy, $c_{S_c}^{(2)}$, and the three-loop non-cusp anomalous dimension of the global soft function, $\gamma_S^{(2)}$. We present a method for extracting those constants using fixed-order codes: the EVENT2 program to obtain the color coefficients of $c_{S_c}^{(2)}$, and MCCSM for extracting $\gamma_S^{(2)}$. We present all necessary formulae for resummation of the mMDT groomed heavy jet mass distribution at NNNLL accuracy.

(2) Sc , and the three-loop non-cusp anomalous dimension of the global soft function, γ (2) S . We present a method for extracting those constants using fixed-order codes: the EVENT2 program to obtain the color coefficients of c (2) Sc , and MCCSM for extracting γ (2) S . We present all necessary formulae for resummation of the mMDT groomed heavy jet mass distribution at NNNLL accuracy. At the Large Hadron Collider (LHC), quantum chromodynamics (QCD) has been studied at an unprecedented precision in both experiment and theory. Major advances have especially recently been made in the field of jet substructure in which the internal dynamics of jets are exploited to learn information about the short-distance physics; see Refs. [1,2] for recent reviews. Because jets are identified as restricted regions of phase space, particles that land in a jet can come from numerous sources from outside of the jet that contaminate a simple theoretical picture of the jet's dynamics. To mitigate this contamination radiation, which may come from initial-state radiation, secondary parton scatterings, or secondary proton scatterings (pile-up), analysis techniques broadly referred to as jet grooming have been introduced. The jet groomers systematically identify emissions in the jet that likely came from such contamination and remove them from consideration. The best theoretically understood jet groomers are the modified mass mass drop tagger (mMDT) [3,4] and soft drop [5] algorithms, due to their unique feature of elimination of non-global logarithms (NGLs) [6], leading correlations between in-jet and out-of-jet scales.
The elimination of NGLs of mMDT and soft drop groomers means that calculations of observables can be meaningfully extended beyond leading-logarithmic accuracy. Refs. [7,8] derived an all-orders factorization theorem for observables like the jet mass measured on mMDT/soft drop groomed jets, in the limit that the jet mass is small with respect to the jet energy. This factorization theorem enabled calculation to next-to-next-to-leading logarithmic (NNLL) accuracy for jet observables, resulting in predictions with controlled theoretical uncertainties over a wide, and experimentally-relevant, region of phase space.
With the high theoretical precision achievable with mMDT/soft drop and the continual experimental measurements, a realistic goal of this program is a first extraction of underlying physical quantities from jet substructure. First studies have demonstrated that distributions of observables on these groomed jets are very sensitive to the value of α s [30]. Further, grooming promises to provide an unambiguous theoretical definition of the top quark mass, and a universal robustness to its value [31]. To reach the necessary theoretical precision for these measurements requires not only high logarithmic accuracy, but also high fixed-order accuracy and well-understood non-perturbative corrections. Next-to-next-to-leading fixed order (NNLO) predictions for mMDT and soft drop groomed events have been produced using the MCCSM (Monte Carlo for the CoLoRFulNNLO Subtraction Method) code [32,33] for jet production in e + e − collisions, and a detailed analysis of non-perturbative corrections within the context of the factorization theorem for groomed jets has been presented [34].
For precision extraction of the strong coupling, for example, logarithmic accuracy should extend to next-to-next-to-next-to-leading logarithmic order (NNNLL), the first order at which simple additive matching to NNLO fixed-order is possible. NNNLL is the order to which thrust [45,46] and C-parameter [47,48] are resummed for precision extraction of α s from LEP data. Such high precision is possible because improving logarithmic accuracy requires calculation of universal anomalous dimensions and constants for the functions that appear in the factorization theorem for these observables. In the case of mMDT/soft drop groomed jet mass, the factorization theorem for the process e + e − → hadrons reads where σ 0 is the leading-order e + e − → qq cross section, z cut is a parameter of the mMDT/soft drop groomer, and τ L , τ R are the left-and right-hemisphere masses, respectively. Extending this factorization theorem beyond NNLL requires the two-loop constants for each of the functions on the right, as well as their three-loop anomalous dimensions. The hard H(Q 2 ) and jet functions J(τ ) are known completely through three loops, and the global soft function S(z cut ) is known completely through two loops [49,50]. Consistency of the factorization theorem requires the sum of anomalous dimensions of the functions to vanish so that the cross section is renormalization-group invariant. Therefore, for NNNLL resummation, we only need to calculate the two-loop constants of the collinear-soft function S c (τ, z cut ) and the three-loop global soft function anomalous dimension.
In this paper, we use numerical fixed-order event generators to extract these necessary ingredients. We will restrict our analysis to jet masses groomed with mMDT for simplicity and compactness. In this case, we use EVENT2 [51] to find the two-loop constants c mMDT Sc of the collinear-soft function (in Laplace conjugate space) to be separated into distinct color channels where C F = 4/3, C A = 3, and T R = 1/2 in QCD, and n f is the number of active quark flavors. Using MCCSM, we find the three-loop anomalous dimension of the global soft function γ mMDT S to be where we fix the number of active quarks to n f = 5. These results enable resummation to NNNLL accuracy for jet substructure observables for the first time, which we present in a companion paper [52].
The outline of this paper is as follows. In Sec. II, we review the mMDT grooming algorithm and its all-orders factorization theorem. In Sec. III, we discuss our procedure for extracting the two-loop constants of the collinear-soft function. In Sec. IV, we extract the three-loop anomalous dimension of the global soft function using numerical results from MCCSM. We conclude in Sec. V. Several appendices provide the analytical and numerical details of the results provided in the main text.

II. MMDT GROOMING AND FACTORIZATION THEOREM
In this paper, we restrict our analysis to hadronic final states in electron-positron annihilation where high fixed-order results are available and can be used to extract currentlyunknown data for resummation. As such, we will only discuss the grooming algorithm and observables for e + e − collisions, and we point readers interested in other applications to the cited literature. In this section, we will review the definition of the grooming algorithm, as well as the all-orders factorization theorem that we exploit.
The mMDT groomer [3], or soft drop with angular exponent β = 0 [5], proceeds as follows: 1. Divide the final state of an e + e − → hadrons event into two hemispheres in any infrared and collinear safe way. This could include separation with thrust [53,54], or in a recoilfree manner using broadening axes [55,56].
2. In each hemisphere, run the Cambridge/Aachen jet algorithm [57,58] to produce an angular-ordered pairwise clustering history of particles. The clustering metric appropriate for e + e − collisions is where θ ij is the angle between particles i and j in the same hemisphere.
3. Starting with the left hemisphere and at widest angle, step through the Cambridge/Aachen particle branching tree. At each branching in the tree, test if where i and j are the daughter particles at that branching and z cut is some fixed numerical value where 0 ≤ z cut < 0.5. If this is true, stop and return all particles that remain in the left hemisphere. If it is false, remove the lower energy branch, and continue to the next branching at smaller angle. Repeat for the right hemisphere.
4. Once the groomer has terminated, any observable can be measured on the particles that remain in the two hemispheres.
We will refer to this procedure as the mMDT groomer or simply mMDT throughout this paper.
From this procedure, Ref. [8] proved an all-orders factorization theorem for the cross section differential in the groomed hemisphere masses τ L and τ R , where for mass m i and energy E i of hemisphere i. For τ L , τ R z cut 1, the cross section factorizes into: where σ 0 is the leading-order cross section for e + e − → qq, H(Q 2 ) is the hard function for quark-anti-quark production in e + e − collisions, S(z cut ) is the global soft function for mMDT grooming, J(τ i ) is the quark jet function for hemisphere mass τ i , S c (τ i , z cut ) is the collinearsoft function for hemisphere mass τ i with mMDT grooming, and we suppress the dependence on the renormalization scale µ in all functions. The symbol ⊗ denotes convolution over the hemisphere mass τ i . In most of the expressions that follow, for simplicity we will set σ 0 = 1 GeV −2 . While we will calculate data relevant for the functions in this factorization theorem in this paper, we will do so through the single-differential cross section of the groomed heavy hemisphere mass, defined as The subscript g on the cross section denotes that this is groomed. Additionally, we note that this definition of the heavy hemisphere mass differs from the standard definition in the ungroomed case when the heavy hemisphere mass is defined as: with Q being the center-of-mass energy. When hemispheres are groomed, the grooming eliminates their dominant correlations, and so it is more natural to define the groomed mass with respect to the hemisphere energy, and not the center-of-mass energy.
The factorization theorem is simplest to analyze in Laplace space, where we Laplace transform in both τ L and τ R to eliminate the convolutions. In Laplace space, the cross section becomes a simple product: where ν L (ν R ) is the Laplace conjugate of τ L (τ R ). In this product form, each function in the factorization theorem satisfies a simple renormalization group equation (RGE), where d F is a constant, µ F is the canonical scale, and γ F is the non-cusp anomalous dimension particular toF . Γ cusp is the cusp anomalous dimension for back-to-back light-like Wilson lines in the fundamental representation of color SU(3). Large logarithms of hemisphere masses can be resummed to all orders in α s through this renormalization group equation. It can be solved exactly, and we present its explicit solution through α 3 s in App. B. The order to which logarithms can be resummed through this renormalization group equation depends on the accuracy to which its components are calculated. For the canonical definition of logarithmic accuracy [59], Table I shows the order in α s to which the components of the renomalization group equation are needed. In this table, "LL" denotes leading logarithmic accuracy, "NLL" is next-to-leading logarithmic accuracy, etc. Ref. [8] resummed the mMDT groomed mass distribution to NNLL accuracy, 1 requiring the new calculation of 1 The logarithmic counting of mMDT groomed mass is actually modified from the traditional counting.
mMDT grooming completely eliminates double logarithms of the mass in the cross section, so what is labeled NLL in Table I is actually LL for mMDT groomed mass, and similar for higher accuracy. However, we keep with the counting as in Table I   one-loop constants and two-loop non-cusp anomalous dimensions. The goal of this paper is to extend the order of the known components to accomplish NNNLL resummation.
The cusp anomalous dimension is now known to four-loop order [60][61][62][63][64][65][66][67][68][69][70][71][72][73], and the QCD β-function is as well [74]. The hard and jet functions are known through three-loop order [45,[74][75][76][77][78][79][80], as they are relevant for resummation of a broad class of observables. The global mMDT soft function is now completely known to two-loops [8,49,50], while only the two-loop non-cusp anomalous dimension of the collinear-soft function is known. Thus, to achieve NNNLL accurate resummation of mMDT groomed mass, we need the two-loop collinear-soft function constants, and the three-loop non-cusp anomalous dimensions of the soft and collinear-soft functions. Actually, because the cross section is renormalization-group invariant, the sum of non-cusp anomalous dimensions of the functions in the factorization theorem must vanish: Hence, only the two-loop constants of the collinear-soft function and the three-loop non-cusp anomalous dimension of the global soft function are needed to accomplish resummation at NNNLL accuracy. In the following sections, we numerically extract these values from fixedorder codes. All results needed for solving the renomalization group equations to O(α 3 s ) are provided in the appendices.
Before moving to the numerical extraction, it is useful to comment on the region of phase space in which mMDT grooming acts non-trivially. Assuming that the grooming parameter is small, z cut 1 formally, then only low-energy particles in each hemisphere are groomed away. To leading power in a soft particle's energy E s , the mMDT constraint is where E H = Q/2 is the hemisphere energy, which is half of the center-of-mass collision energy in the soft limit. For such a soft emission that just passes mMDT, its contribution to the hemisphere mass is where the inequality corresponds to a soft emission right at the hemisphere boundary. Taking the upper bound as the parametric scaling of the soft particle's energy, we have Then, using the mMDT constraint, the value of the hemisphere mass for such an emission and mMDT strongly affects the cross section for values of ρ smaller than this.

III. EXTRACTION OF TWO-LOOP CONSTANTS
In this section, we will present a method for numerical extraction of the two-loop constant terms of the collinear-soft function c (2) Sc for the factorization theorem formulated in Laplace space, Eq. (10). Our procedure for doing so will be to relate the cross section differential in the groomed hemisphere mass ρ to the total cross section. We can express the leading-power (LP) differential cross section for ρ → 0 as a sum of terms with support exclusively at ρ = 0 and terms with support away from 0: The terms that have support away from ρ = 0 are defined to integrate to 0 on ρ ∈ [0, 1].
The total cross section can then be written as: where dσ g /dρ is the full differential cross section of the groomed hemisphere mass for all The desired unknown two-loop constants are contained in the δ-function coefficient, D δ,g .
One approach to determining them is to numerically evaluate the integral that remains using the known total cross section for e + e − → hadrons and a fixed-order code. Performing this particular numerical integral is somewhat complicated by the fact that grooming introduces multiple regions and so the integrand exhibits cusps and has support all the way out to the maximal value of ρ. We will significantly rewrite this expression into a more directly useful form so that the numerical integrals we must perform are restricted to the region where grooming turns on, ρ ∼ z cut . First, we separate the integral into two parts, Below 2z cut mMDT grooming does, while above 2z cut it does not, dominantly enforce constraints on the radiation. Then we add and subtract the singular differential cross section for ungroomed heavy hemisphere mass to the second integral on the right of Eq. (19), : For ρ ∼ 1 z cut , mMDT grooming has no effect, so the first integral in this expression is dominated by the region near ρ z cut , and we can drop the step function Θ. Integrals on the right extend to 4, rather than 1 as simple consequence of the factor of 4 difference in definition of the groomed versus ungroomed heavy hemisphere mass, as discussed in Sec. II.
For the integral exclusively of the ungroomed heavy hemisphere mass cross section, we can write it as Assuming z cut 1, the lower bound of the first integral can be set to 0, neglecting power corrections that vanish as z cut → 0.
After implementing all these identities, the total cross section can be written as Further progress can be made by replacing the total cross section σ tot with its expression as an integral over the ungroomed hemisphere mass cross section. Writing the ungroomed hemisphere mass cross section at LP as the total cross section can also be expressed as Using Eqs. (22) and (24), we find the relationship Through two-loop order, the only unknowns in this equation are the collinear-soft function constants c Sc within the D δ,g term and the first two integrals. However, those integrals can be numerically evaluated with a fixed-order code. We will first validate this relationship at one-loop, and then use it at two-loops to determine c (2) Sc .

A. Test at One-Loop
At one-loop order, Eq. (25) simplifies significantly. For ρ > 2z cut , the groomed and ungroomed differential cross sections are identical, and so the third term in Eq. (25) is identically zero. Hence, the one-loop relationship between the groomed and ungroomed cross sections can be expressed as This equation has been written so that everything on the left is known, and the right side can be evaluated numerically. Using the formulas from App. C, the left side evaluates to

2.32896
To evaluate the right side of Eq. (26), we generate about 10 13 e + e − → qqg events in As z cut decreases, the integral values are seen to converge to the exact expected value computed in Eq. (27).

B. Fit at Two-Loops
Having validated the procedure at one-loop, we move on to using it at two-loops to determine the constant terms of the collinear-soft function, c Sc . Unlike at one-loop, all terms in Eq. (25) are generically non-zero, so we have the equality: The right side of the equality is unknown and must be evaluated numerically. The left side is completely known, except for the two-loop constant of the collinear-soft function, c Sc . Using the results of App. C and the singular ungroomed heavy hemisphere mass distribution calculated in Refs. [81][82][83], the left side can be expressed as: In this expression, ζ(3,1) is a multiple zeta value, where The coefficients of the two-loop anomalous dimension of the soft function γ S were calculated in Refs. [8,49,50] and are: The two-loop constants of the soft function c (2) S were recently calculated by the SoftServe collaboration [49,50]: The only unknowns in this expression are therefore the two-loop constants of the collinearsoft function. As with the soft function's anomalous dimension and constants, we break up the collinear-soft constant into distinct color channel coefficients: Sc,n f .
To evaluate the right side of Eq.  29), we also measure the ungroomed heavy hemisphere mass over the same range. We then construct a smooth interpolating function for both integrands in Eq. (29) which we then numerically integrate. Because of the finite cutoffs inherent in EVENT2, we cannot practically extend the first integral all the way down to 0. Instead, we integrate over the reduced range of ρ ∈ 2z cut 400 , 2z cut .
We have additionally verified that increasing or decreasing the lower bound by a factor of 2 only affects the results within quoted uncertainties. With the integrals thus calculated, we then solve Eq. (29) for each color channel's collinear-soft function two-loop constant coefficient.
The results of this procedure are plotted in Fig. 2. Here, we show the result of the twoloop constant extraction at each value of z cut considered. We observe clear convergence as z cut → 0 for the C A and T R n f channel coefficients, with the extracted value denoted by the dotted line, and the shaded band is our claimed uncertainty. For the C F channel, the convergence is less clear, with the extracted value of the two-loop constant appearing to diverge as z cut → 0. In other contexts, it has been noted that the extraction of anomalous dimensions and constants in the C F channel with EVENT2 is problematic [8,81], unless a significantly larger number of events are generated. Thus, we take as our extracted value of the C F channel constant the lowest value determined through this procedure. Further, if the approach to the true value were exponential in log z cut , as would be expected, then the difference between the extracted constants at different z cut values that differ by a factor should uniformly decrease. Thus, we take our uncertainty to be twice the difference between the extracted value at z cut = 0.001 and z cut = 0.00316.
Finally, we obtain for the coefficients in the color decomposition of the two-loop constant of the collinear-soft function as

IV. EXTRACTION OF THREE-LOOP ANOMALOUS DIMENSIONS
We now turn to extraction of the three-loop non-cusp anomalous dimension of the mMDT global soft function, γ S . With the results quoted in the appendices and the extracted value of the two-loop constant of the collinear-soft function, we are able to solve the renormal-ization group equation through O(α 3 s ) order, and construct the differential cross section by inverse Laplace transformation. For ρ > 0, the only unknown piece of the cross section is γ (2) S and so its value can be extracted through comparison of the analytic prediction in the regime where ρ z cut 1 to a next-to-next-to-leading fixed-order e + e − → hadrons event generator. For some fixed values of z cut , we can then fit the analytic cross section to match the predictions of the NNLO code for the value of the anomalous dimension and extrapolate to z cut → 0 to determine the γ S . The first step of this procedure is to solve the renormalization group equations to O(α 3 s ) to identify the cross section predicted by the factorization theorem. The O(α 3 s ) coefficient of the expansion is given in Eq. (C4) where we substituted explicitly the color factors of QCD (C F = 4/3, C A = 3, T R = 1/2), and set the number of active quarks to n f = 5.
We have also just used the central values of the collinear-soft constants, with no accounting for uncertainties. The uncertainties that we will find for γ (2) S will more than account for uncertainties in the collinear-soft function constants, so we can safely ignore them. In this logarithmic form of the cross section, the anomalous dimension γ (2) S appears as a constant offset of the cross section. This suggests a procedure for its numerical extraction.
As the mMDT grooming removes double logarithms in ρ to all orders, the singular term at O(α 3 s ), D C,g , is only quadratic in log ρ. Thus, our task of determining γ S is reduced to finding the constant term of a parabola. A parabola exhibits greatest sensitivity to the constant where its slope is smallest; that is, near its extremum, in our case a minimum. Therefore, we would like to fit for γ (2) S near the minimum of the logarithmic parabolic cross section. However, the factorization theorem is only valid for ρ z cut 1, so the fit is only sensible if the minimum is compatible with this singular limit. These two competing requirements, fitting near the minimum and existing in the factorization regime, will greatly restrict the range of z cut that we can consider.
Taking the derivative of the function D C,g given in Eq. (C4) to determine the point ρ min at which the cross section is minimized, we can then compare ρ min to z cut . We plot the result of this minimization and comparison to z cut in Fig. 3 where the solid curve is the ratio z cut /ρ min , as a function of z cut . The fixed ratios of 1 and 10 are also shown to guide the eye.
We find that only for z cut 0.01 is ρ min smaller than z cut , and only for z cut 0.03 is there an order of magnitude between them. As our factorization theorem requires a parametric separation between ρ and z cut , we will only consider values of z cut larger than 0.03 in our fit for the anomalous dimension. Further, demanding that z cut is relatively small, we will restrict to an upper bound on z cut of 0.1. While a limited range in z cut , and relatively large values, we will see that it is sufficient for an extraction of the anomalous dimension.
For numerical generation of the cross section for the process e + e − → hadrons at O(α 3 s ) we use the MCCSM (Monte Carlo for the CoLoRFulNNLO Subtraction Method) code [32].
MCCSM has previously been used to calculate distributions of standard event shapes [84][85][86] as well as mMDT and soft drop groomed event shapes [33] in e + e − collisions. To validate the implementation of the mMDT grooming we compared the resulting distributions at LO and NLO accuracy from the independent implementations in MCCSM and EVENT2. We found perfect agreement between the predictions of the two codes.
For extraction of the three-loop soft function anomalous dimension, we computed the distribution of the mMDT groomed heavy hemisphere mass for four values of z cut that satisfies the criteria established earlier: z cut = 0.04, 0.06, 0.08, 0.1. Any Monte Carlo integrator must contain a lower cut y min on the phase space that any scaled two-particle invariant of the event must exceed. Double-precision arithmetic limits the lowest value of y min to the range 10 −8±1 . As for the extraction of γ (2) S we need stable predictions for ρ z cut we have to investigate carefully the lower limit on ρ where our predictions are still reliable at a given y min . We found that we can use y min = 10 −8 for computing the phase-space integrals of the real-virtual contributions and y min = 10 −7 for the double-real contributions to obtain realiable distributions for ρ ≥ 5 · 10 −5 [87]. We used these values of the technical cut to generate The events were binned logarithmically, with 10 equal bins in the range −10 < log(ρ) < 0, i.e., we computed the coefficient of the distribution ρ dσg dρ directly. First we performed a flat average of the results of the batches of events. The resulting distributions were fairly stable and smooth except for the first two bins with smallest values of ρ due to occasional mis-binning of the squared matrix element and subtraction terms. We filtered the results of the individual batches of events by selecting the ten batches with largest uncertainty in the second bin. Removing the batches from the average one-by-one, we found compatible results for the distributions, but gradually decreasing uncertainty in the first two bins. We produced our final distributions without removing any of the events, but employing a weighted (by the inverse of the uncertainties) average of the individual results of batches. We have verified that the resulting distributions are consistent with the unfiltered and unweighted events, though with significantly smaller uncertainties in each histogram bin. The resulting distributions are provided in App. D for each value of z cut that was generated, and we also show it for z cut = 0.04 in Fig. 4 (left). S . This is illustrated more directly in the right panel of Fig. 4. Here, we have set the analytic prediction equal to the MCCSM data and solved for γ (2) S at each value of ρ. Below about ρ 0.01, the value of the anomalous dimension at each data point is consistent within uncertainties of one another, demonstrating independence of ρ, as γ (2) S must. For each value of z cut that we generated, we determined the anomalous dimension through the following procedure. Taking the plot on the right of Fig. 4 (and the corresponding plots for other values of z cut ), we fit for the anomalous dimension from data points in the range ρ ∈ [0.0002, 0.01]. 2 This range ensures that the largest value of ρ is still significantly smaller 2 We have verified that the extracted central value of γ (2) S is consistent within uncertainties using a fit range than z cut . The contribution of a data point to the anomalous dimension and its uncertainty are weighted by the quoted variance of that data point. So, the points at larger ρ with smaller uncertainty dominate the fit. The result of this fitting procedure is shown in Fig. 5.
As z cut decreases, the value of the extracted anomalous dimension also decreases, and the central values lie almost perfectly on a line. Extrapolating the points in the z cut → 0 limit produces the desired three-loop anomalous dimension. Performing a weighted fit of the points in Fig. 5 to a line, we find the intercept, that is, γ S , to be: We expect this claimed uncertainty to be a conservative overestimate due to the extremely linear behavior of the central values. 3

V. CONCLUSIONS AND OUTLOOK
Using numerical predictions at NLO and NNLO accuracies, we were able to extract two-loop constants and three-loop anomalous dimensions for the mMDT groomed mass factorization theorem in e + e − → hadrons events. As a result, we are able to resum the groomed mass distribution to NNNLL accuracy, extending the predictions of Ref.
[8] to one higher logarithmic order. We present NNNLL+NNLO resummed and matched predictions for the heavy hemisphere mMDT groomed jet mass in a companion paper [52].
The results presented here also enable higher-precision calculations for groomed mass distributions for jets produced in hadron collisions. The cross section for groomed mass ρ of jets at a hadron collider can be expressed as At a hadron collider, both quark and gluon jets can be produced, so we have to sum over both possibilities. In the factorization limit ρ z cut 1, all radiation that contributes to ρ must be collinear, so the only dependence on the groomed jet mass lives in the universal jet and collinear-soft functions, appropriate for quarks or gluons. The normalization factor of ρ ∈ [0.000075, 0.004], though with larger uncertainties. 3 There are additional uncertainties in MCCSM due to a numerical fit of part of the NNLO matrix element.
We find that the uncertainty stated here in the value of the anomalous dimension is representative of uncertainties that come from this fit in MCCSM. N i (p T , y, z cut , R) encodes the relative fraction of jets in the sample that are quark or gluon flavor, and can in principle depend on all other scales in the event, except for ρ: the jet's transverse momentum p T , rapidity y, the value of z cut , and the jet radius R. The results presented here enable NNNLL resummation of all ρ dependence at a hadron collider.
Extension of the mMDT groomed jet mass factorization theorem to three-loop order is one necessary step for complete perturbative control over the groomed jet mass distribution. and finite z cut contributions should be included. This was studied in the original mMDT paper [3], but a factorization theorem for this regime does not yet exist.
Additionally, our results in this paper are restricted to the mMDT groomer, or soft drop with β = 0. An analysis for general soft drop angular exponent β is much more challenging at three-loops. Our extraction of the three-loop anomalous dimension for mMDT grooming was simplified because mMDT removes double logarithms of the mass, but for β > 0, double logarithms generically exist. What was a problem of finding the constant term of a parabola here turns into finding the constant term of a fifth-degree polynomial for general β. A significant amount of data are known for general soft drop grooming at two-loop accuracy [49,50], but complete resummation to NNNLL accuracy will be challenging. Soft drop grooming with non-zero β can provide a handle on quark vs. gluon fractions in a jet sample at a hadron collider [30] as well as provide just the right amount of grooming for boosted top quarks [31]. Three-loop data for soft drop grooming is therefore highly desired. We hope that the work presented in this paper inspires further precision predictions for jet physics. The QCD β-function is defined to be For expansion of the cross section to three-loop order, we need the β-function to two-loop order [88][89][90]. The first two coefficients are Correspondingly, we need the cusp anomalous dimension to three-loop order. The first three coefficients of the cusp anomalous dimension are [91][92][93][94]: T R n f ,

Non-Cusp Anomalous Dimensions and Constants
The non-cusp anomalous dimensions γ and constants c can be expanded in similar series in α s as the β-function and cusp anomalous dimensions. We have To O(α 3 s ) accuracy in terms with support away from 0, we need the non-cusp anomalous dimension to α 3 s order and the constants to α 2 s order.

a. Hard Function
The first three non-cusp anomalous dimensions of the hard function (quark form factor) are [75,77,95,96]: The first two constants of the hard function are [75,76]:

b. Jet Function
The first three non-cusp anomalous dimensions of the quark jet function are [74,75,77]: The first two constants of the jet function are [45]:

c. mMDT Global Soft Function
The first two coefficients of the non-cusp anomalous dimensions of the mMDT global soft function were calculated previously [8,50] and the third coefficient is new in this work: S = −11600 ± 2000 (n f = 5) .
In γ S , the coefficient of the C 2 F term is known analytically: where Cl 2 (x) is the Clausen function: The C F C A and C F T R n f coefficients in γ (1) S are not known analytically, so we present their decimal expansion to 6 significant figures. The value of γ (2) S is evaluated with five active quarks, n f = 5.
The first constant of the mMDT soft function was calculated in [8] and the second constant was calculated by the SoftServe collaboration [49,50]: Uncertainties in the values of these numerical expressions occur in the last digit, which we do not quote because they are vastly smaller than other relevant uncertainties.

d. mMDT Collinear-Soft Function
By renormalization group invariance, the anomalous dimension of the collinear-soft function can be found to all-orders through its relationship with the other anomalous dimensions in the factorization theorem: As such, we will not present the fixed-order expansion of the collinear-soft function's anomalous dimension here. The first coefficient of the collinear-soft function's constant terms was calculated in Ref. [8] and the second coefficient is new to this work:

Appendix B: Fixed-Order Expansion
The hemisphere mMDT groomed jet mass factorization theorem in e + e − → hadrons collisions is where H(Q 2 ) is the hard function for quark-anti-quark production in e + e − collisions, S(z cut ) is the global soft function for mMDT grooming, J(τ i ) is the quark jet function for hemisphere mass τ i , and S c (τ i , z cut ) is the collinear-soft function for hemisphere mass τ i with mMDT grooming. The symbol ⊗ denotes convolution over the hemisphere mass τ i . The convolution can be transformed into a simple product by Laplace transforming the jet and collinear-soft functions appropriately. The Laplace transform of a function F in the factorization theorem and the Laplace-transformed factorization theorem is All functions in this factorization theorem satisfy a simple renormalization group equation: where µ is the renormalization scale. The expression in parentheses isF 's anomalous dimension where Γ cusp is the cusp anomalous dimension, d F is a coefficient particular to function F , µ F is the canonical scale of F , and γ F is the non-cusp anomalous dimension ofF .
The anomalous dimension equation can be solved iteratively in powers of α s . Through O(α 3 s ), the solution is For the functions in the factorization theorem, the d F coefficients are: The canonical scales µ 2 F of the Laplace-space functions are: Then, to find the cross section in real space through O(α 3 s ), we multiply the Laplacespace functions together, expand in α s , and then inverse Laplace transform. Because mMDT grooming removes double logarithms, the inverse Laplace transformation is relatively easy, and one only needs three transformations: L −1 (log 2 ν) = 2 log τ τ + − π 2 6 δ(τ ) , L −1 (log 3 ν) = −3 log 2 τ τ + + π 2 2 1 τ + − 2ζ 3 δ(τ ) .
The coefficient functions are D δ,g = 1 + α s 4π C F −2 + 12 log 2 + 16 log 2 log z cut + 4 log 2 z cut (C2) + α s 4π 2 C 2 F 4 + 4 3 π 2 − 7 90 π 4 − 72ζ 3 − 18 log 2 + 72 log 2 2 + 96ζ 3 log 2 +c (2) S,C F + 2c (2) Sc,C F − 4π 2 log(4z cut ) − γ (1) S,C F log(4z cut ) − 8 log z cut log(16z cut ) − 8 3 π 2 log 2 z cut + 192 log 2 2 log z cut + 48 log 2 log 2 z cut + 128 log 2 2 log 2 z cut + 64 log 2 log 3 z cut + 8 log 4 z cut The D C,g coefficient is extremely unieldy, so we only quote here the approximate numerical function: z cut = 0.04, 0.06, 0.08, 0.1. For the differential cross section of the mMDT groomed heavy hemisphere mass, MCCSM calculates the ρ-dependent factors A, B, and C at leading, next-toleading, and next-to-next-to-leading fixed order, respectively. These can then be combined to evaluate the cross section through O(α 3 s ) as: Here, α s = α s (µ), the strong coupling evaluated at the renormalization scale µ, β 0 and β 1 are the QCD β-function coefficients, and Q is the center-of-mass collision energy. When µ = Q, the terms proportional to the β-function vanish, and the expression simplifies. To extract γ S , we need the distribution of C as a function of ρ. The value of C for log ρ ∈ [−9.5, −0.5] in steps of 1 and all z cut values we consider are listed in Table II. The entries in this table were determined using the filtering and weighting procedure described in Sec. IV.