Joint resummation of two angularities at next-to-next-to-leading logarithmic order

Multivariate analyses are emerging as important tools to understand properties of hadronic jets, which play a key role in the LHC experimental program. We take a first step towards precise and differential theory predictions, by calculating the cross section for e+e− → 2 jets differential in the angularities eα and eβ. The logarithms of eα and eβ in the cross section are jointly resummed to next-to-next-to-leading logarithmic accuracy, using the SCET+ framework we developed, and are matched to the next-to-leading order cross section. We perform analytic one-loop calculations that serve as input for our numerical analysis, provide controlled theory uncertainties, and compare our results to Pythia. We also obtain predictions for the cross section differential in the ratio eα/eβ , which cannot be determined from a fixed-order calculation. The effect of nonperturbative corrections is also investigated. Using Event2, we validate the logarithmic structure of the single angularity cross section predicted by factorization theorems at Oαs2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s^2\right) $$\end{document}, highlighting the importance of recoil for specific angularities when using the thrust axis as compared to the winner-take-all axis.


Introduction
Hadronic jets play a central role in collider physics as proxies of the hard quarks and gluons produced in short-distance interactions.A precise theoretical understanding of jet properties is often key to establishing measurements which can either further confirm the predictions of the Standard Model at higher accuracy, or identify deviations that could hint at New Physics.Progress in QCD calculations involving jets has been impressive in the last years, boosted by the demands of the LHC experimental program.Jet properties are typically studied using two complementary tools: analytic resummation and Monte Carlo parton showers.The latter offer a fully exclusive description of the final state, enabling the user to perform any measurement, but their formal accuracy is currently limited to leading-logarithmic order. 1 Corrections to the traditional parton shower method have been considered lately, e.g. by incorporating additional information about interference effects [2][3][4] and higher-order splitting functions [5][6][7][8].On the other hand, with an analytic approach one can often achieve higher logarithmic resummations and obtain uncertainty estimates that can be validated by comparing different orders.However, this approach is traditionally limited to single differential measurements.
Inspired by ref. [9], we recently took the first step towards a precise and more differential characterization of jets by constructing an effective field theory (called SCET + ) and using it to derive a factorization formula [10], which enables the simultaneous resummation of two independent observables to higher logarithmic accuracy.This opens up the possibility of performing multivariate analyses, including correlations with controlled theory uncertainties.Applications of our framework are particularly relevant in the context of jet substructure studies (see e.g.ref. [11] for a recent review), where a more detailed characterization of the QCD radiation pattern within a jet is exploited to obtain crucial information about the hard scattering process, thereby providing innovative ways to search for New Physics.This generally involves multi-differential cross sections, with several independent measurements performed on a single jet and the possibility to exploit shared information content among these observables.Furthermore, several of the most powerful discriminants of quark-vs.gluoninitiated jets or of QCD jets vs. boosted hadronically decaying heavy particles are formed by taking ratios of two observables, as is done for N -subjettinesses [12,13], energy-energycorrelation functions [14,15] and planar flow [16].These are typically not infrared-and collinear-safe [17] but still Sudakov safe, meaning that they can nevertheless be properly defined and calculated by marginalizing the corresponding resummed double differential cross section [18].
As a first step in understanding multi-differential cross sections beyond next-to-leading logarithmic (NLL) accuracy, we demonstrate in this paper how to exploit our theoretical framework for the case of the simultaneous measurement of two event shapes for e + e − collisions in the dijet limit, where all-order resummations are essential to obtain reliable theory predictions.In order to avoid complications related to non-global logarithms [19], we restrict ourselves to event-based observables and postpone the case of jet-based measurements at NNLL to a future publication.Our focus here is on the family of infrared and collinear safe angularities [20], which generalize the classic event-shape variables thrust and broadening, and characterize the energy distribution of final-state particles as function of the angle with respect to some axis.Calculations of single angularities with respect to the thrust axis were carried out at NLL accuracy in refs.[20][21][22][23].Recently they have been analyzed up to NNLL including next-to-next-leading fixed-order (NNLO) corrections [24] for the purpose of a precision determination of α s (m Z ) from LEP data which will provide complementary information to analogous precision fits based on event shapes like thrust [25,26] and C-parameter [27].Furthermore, (generalized) angularities measured on individual jets [16,28] are useful tools to investigate jet substructure [29].
We calculate the jointly resummed cross section of two angularities up to NNLL using SCET + and obtain all one-loop ingredients entering the factorization formulae for the various kinematic regimes characterizing the double differential spectrum.We match the resummed distribution to the corresponding fixed-order QCD result at NLO to obtain a reliable description of the cross section beyond the dijet limit.Theoretical uncertainties are provided by suitable "profile functions" designed to produce scale variations that smoothly interpolate between the distinct kinematic regions where resummations must be handled differently.We also investigate nonperturbative corrections, and compare the results of our numerical analysis to the parton shower of Pythia 8.2 [30].
By projecting our double differential cross section, we also obtain predictions for the cross section differential in the ratio of two angularities, which cannot be determined from a fixedorder calculation.Furthermore, we analyze single angularity distributions up to NNLL+NLO and compare the fixed-order expansions from our resummed distributions against numerical results from Event2, for angularities calculated with respect to the thrust axis and winnertake-all (WTA) axis [31,32].Our analysis demonstrates that for the WTA axis the same factorization formulae can be used for the whole range of angularities, even for those measurements that would be sensitive to recoil by soft radiation if the thrust axis was used.
The paper is organized as follows: sec. 2 describes in detail our theoretical framework and the analytic input for our numerical analysis.After reviewing the factorization for the double differential angularity distribution, we collect all relevant fixed-order ingredients (correcting typos in the literature) and anomalous dimensions.We describe in detail our scale choices and procedure to estimate the perturbative uncertainty.In sec. 3 we show numerical results for single and double differential angularity distributions at NNLL+NLO accuracy, as well as for the ratio of two angularities.We conclude in sec.4.

Angularities
The angularities e α are a one-parameter family of global e + e − event shapes, defined as [20] where Q is the center-of-mass energy.The sum runs over all particles i, where E i denotes the energy and θ i the angle of the particle with respect to an appropriately chosen axis.Smaller values of angularities correspond to more collimated radiation, where the parameter α determines the weight of the angle.Our convention for α is such that for small angles, i.e. α = 2 corresponds to thrust [33] and α = 1 to (total) broadening [34,35].
For angularities with α 2, the direction of the thrust axis is insensitive to (recoil by) soft radiation, but as α → 1, and certainly for α ≤ 1, this effect cannot be ignored [36].Thus we find it convenient to use an axis that is recoil-insensitive [32].This is accomplished by clustering the event with exclusive k T [37], which splits the event into two jets, using the WTA recombination scheme [31,38]. 2 The angle θ i in eq. ( 2.2) will be taken with respect to the axis of the jet the particle belongs to, so there is no global axis for the event.
In our previous publication [10], we focussed on jet-based angularities [16,28].However, since the correlation between soft radiation inside and outside the jet makes these observables theoretically more complicated, introducing non-global logarithms [19], we shall limit ourselves here to event-based angularities.

Power counting and modes for double angularity measurements
We calculate the e + e − → 2 jets cross section differential in two angularities e α and e β taking into account the fact that the phase space is characterized by three different regions (fig. 1 each one with its own factorization theorem that enables the resummation of logarithms of e α and e β .Regime 1 and 3 correspond to the boundaries and were discussed in ref. [9], while we obtained the factorization theorem for regime 2 describing the bulk of the phase space in ref. [10].We will briefly review how the regimes in eq.(2.3) arise, and present the factorization theorems in the next section.The relevant modes (degrees of freedom) in the framework of Soft-Collinear Effective Theory (SCET) [40][41][42][43] are summarized in table 1.In SCET, the real radiation in the two-jet region is either collinear or soft.The corresponding momenta have the following parametric scaling in terms of light-cone coordinates  Here n µ and nµ are light-like vectors along the axes used to define the angularities in eq.(2.2), and The scaling of λ c and λ s in eq.(2.4) is fixed by the measurement: the parametric size of the contribution of collinear or soft radiation to the angularities simplifies to Assuming α > β for definiteness, this implies λ c ∼ e 1/β β and λ s ∼ e α .A consistent theory is obtained if λ α c ∼ λ s or λ β c ∼ λ s , which correspond to regime 3 and 1 in eq.(2.3).In regime 2 there is an additional collinear-soft mode, whose power counting is uniquely fixed by requiring that it contributes to both e α and e β , This leads to and λ − cs ↔ λ + cs for the collinear-soft mode in the other direction.These extensions of SCET have been named SCET + [10,[44][45][46].As one approaches regime 1 and 3 from regime 2, the collinear-soft mode merges with the collinear mode or soft mode, respectively.

Factorization
Before presenting the factorization theorems for the various regimes, we want to point out that these all describe the full cross section up to power corrections, Regime 2 resums the most logarithms but also involves two expansions.Starting from regime 2 and approaching either of the phase space boundaries, one of the power corrections becomes of order one and the other smoothly matches onto the power correction for regime 1 or 3, respectively.We will discuss how to combine these formulae to obtain predictions throughout phase space in sec.2.8.
In regime 1, the power counting in eqs.(2.3) and (2.6) implies that collinear and soft radiation both contribute to e β but only soft radiation contributes to e α .This leads to the following factorization theorem [9], Here, σ0 denotes the Born cross section, the hard function H contains hard virtual corrections, the jet functions J describe the contribution of collinear radiation to e β , and the soft function S accounts for the contribution of soft radiation to e α and e β .Their expressions at one loop are collected in sec.2.4.The delta functions simply sum the various contributions, since angularities are additive.This is basically the factorization theorem for a single angularity e β [22,47], with a soft function that is differential in e α too.
Similarly, in regime 3 only collinear radiation contributes to e β , but collinear and soft radiation contributes to e α , leading to [9]3  (2.12)This is the factorization theorem for e α but with double differential jet functions.Finally, the factorization theorem for regime 2 is given by [10] where the collinear-soft function S accounts for the contribution of collinear-soft radiation to e α and e β .The jet functions J are the same as in eq.(2.11) and the soft function S is the same as in eq.(2.12).
Expanding σ 1 and σ 3 in the SCET + regime described by σ 2 (when resummation is turned off), we obtain the following consistency relations4  We have verified these relation at one-loop order using the expressions in sec.2.4.In the second case, the power corrections turn out to vanish at this order, so we could not determine the exponent p > 0.

Fixed-order ingredients
In this section we collect all fixed-order ingredients needed for our numerical analysis, some of which we calculated ourselves.We use the perturbative expansion where F = H, J, S, S , and give F (0) and F (1) .The following shorthand notation for plus distributions is used (2.16)

Hard function
The hard function entering all aforementioned factorization theorems encodes virtual corrections in the q q-production at the hard scale Q, and is given by the square of the Wilson coefficient in the matching of QCD onto SCET currents [49,50], (2.17)

Jet functions
The single differential jet function in eqs.(2.11) and (2.13) is [32] Note that the constant terms differ from those obtained in ref. [22] because we employ the WTA axis.For angularities with α 2, this is the only difference between using the thrust axis or the WTA axis in the factorization formula. 5he double differential jet function in eq.(2.12) is given by [9] 6 (2. 19) In principle one can perform the derivatives, but we find it more convenient to work the cumulative distributions to avoid complicated plus distributions.Note that we can perform the necessary convolutions using cumulative distributions, as discussed in sec.2.5.3.
Integrating the double differential jet function over e β yields the single differential jet function of e α [48] d(Q β e β ) J (1) (2.20) This is obvious from comparing eq.(2.11) with the factorization theorem for a single angularity e α , as the only difference is that the double differential jet function is replaced by this single differential jet function.

Soft functions
The soft function encodes the effects of soft radiation, which is described by a matrix element of eikonal Wilson lines (along the two outgoing quarks) on which the appropriate measurement is performed.The soft function for a single angularity in eqs.(2.12) and (2.12) is given by [22] 7 For the double differential soft function in eq.(2.11) one obtains [9] 6 Note that this expression is more complicated because we chose to write it in terms of Q β e β instead of Qe β .In particular, the first two lines of the one-loop expression correspond directly to the single differential soft function in eq.(2.21), but for e β .
Integrating the double differential soft function over e α produces the single differential soft function of e β [10], This is obvious from comparing eq. ( 2.11) with the factorization theorem for a single angularity e β , as the only difference is that the double differential soft function is replaced by this single differential soft function.

Collinear-soft function
Finally, the collinear-soft function that enters in eq.(2.13) is given by [51] This is the simplest double differential function, as it contains pure logarithms.To see that it involves a single scale, it is the easiest to consider the double cumulative distribution, which only involves logarithms of e β−1 α e 1−α

Resummation
The factorization theorems enable the resummation of large logarithms of e α and e β through renormalization group (RG) evolution, since each ingredient is only sensitive to a single scale.By evaluating the ingredients at their natural scale, where they contain no large logarithms, and evolving them to a common scale, these logarithms get exponentiated.In this section we give the form of anomalous dimensions and evolution kernels, with explicit expressions provided in app. A. We also discuss how to perform convolutions with cumulative distributions.

Anomalous dimensions
The RG equation of the hard function is Here Γ cusp is the cusp anomalous dimension [52], and γ H (α s ) the non-cusp contribution.
Similarly, for the jet functions for the soft functions and for the collinear-soft function (2.28) Using these expressions, one can verify that the cross sections in eq.(2.11), (2.12) and (2.13) are µ-independent up to the order that we are working, if the following relations hold (2.29) We have checked this equation at one-loop order, and use it to extract the two-loop non-cusp anomalous dimensions, taking the known results for the hard function and soft function to fix all the others.We stress that an essential ingredient to achieve NNLL accuracy in our analysis is provided by the novel calculation of the two-loop soft anomalous dimension in [53].
Cusp and non-cusp contributions to the anomalous dimensions are collected in app. A.

Evolution equations
For the hard function, the solution to RG equation in eq.(2.25) is given by . (2.30) Here K H and ω H are given by where γ H in the subscript denotes the non-cusp anomalous dimension and These integrals can be performed analytically in a perturbative expansion, see app. A. Similarly for the jet and soft function (F = J, S), where t J = Q β e β , t S = Qe α .The evolution kernel U F is given by (2.34) where j J = β and j S = 1.We use the plus distribution and We do not need the evolution kernel for the collinear-soft function, as we choose the collinearsoft scale as the endpoint of our evolution.

Convolutions with cumulative distributions
The most complicated step in our numerical evaluations is the convolution of the evolution kernel in eq.(2.34) with the one-loop double differential jet and soft functions in eqs.(2.19) and (2.22).To avoid subtleties with plus functions, we perform these convolutions using cumulants, as follows.For the cumulative distributions F and G, if we want to perform the convolution of F and G and take the cumulant of the result, we can rewrite Note that since G(y c − x) − G(y c ) vanishes for x → 0, the final integral does not require a plus prescription for F (x).In our case it is convenient to take G to be the cumulant of the double differential jet or soft function, so its derivative is never needed.

The next-to-leading order cross section
In this section we first present the calculation of the double angularity cross section at NLO.We subsequently decompose this result into a singular and nonsingular component.By adding the latter to our resummed cross section, the matching at NLO is achieved.We also give the nonsingular contribution for single angularity measurements.

Calculation
Since the virtual corrections are already included in our factorization theorems, we only need the real contribution to calculate the double differential cross section at O(α s ).The final state consists of three massless partons, which can be characterized by their energy fractions x i in the center of mass frame, normalized to partons 2 and 3 get clustered together into one jet by exclusive k T algorithm [54], because the angle θ 23 is smaller than θ 12 and θ 13 .The other jet then only consists of parton 1. Due to the WTA recombination scheme, the jet axes are along the momenta of particles 1 and 2, so that the angularity is determined by the energy fraction x 3 and the angle θ 23 , The cross section is then calculated numerically using for x q , xq < 1.Specifically, we sample logarithmically in 1 − x q and 1 − x q, using a cutoff that is outside our plot ranges.
Our result is shown in fig. 3 for three pairs of angularities with exponents (α, β).From the double differential jet and soft function in eqs.(2.19) and (2.22), we see that in the resummation regime the phase-space boundaries are at one-loop order.Note that the lower boundary is slightly shifted compared to the canonical expression in eq. ( 2.3).The upper boundary e β ≥ e α corresponds to cos θ 23 = 0 and not to one of the phase-space boundaries ( At NNLO, the phase-space boundaries in the resummation regime are β . (2.41) The lower boundary follows from considering two one-loop jet functions, whose contribution e α i , e β i to the respective hemispheres each satisfies eq.(2.40).8

Fixed-order nonsingular
We start by showing the nonsingular cross section for a single angularity in fig.2, for four representative angularity exponents β = 0.5, 1.2, 2 and 3.Here Q = 1000 GeV. 9 The full NLO cross section is normalized to 1.The singular and nonsingular cross sections are rescaled by the same amount as the full NLO.For small values of the angularities, the singular contribution to the cross section dominates and the nonsingular cross section is power suppressed, demonstrating the validity of the factorization theorem at this order. 10On the other hand, at large values of the angularity, the singular and nonsingular contributions become equal in size, and matching to the NLO is important to correctly describe the cross section in this region.Indeed beyond the endpoint of the distribution, where the cross section vanishes, the singular and nonsingular are exactly equal and opposite in sign.Thus the resummation must be turned off in this region through an appropriate scale choice (see sec.2.7.1) to maintain this cancellation.The bump in the cross section in the fixed order region arises because we use the WTA axis rather than the thrust axis.
Moving on to two angularities, fig. 3 shows the total NLO cross sections and the corresponding NLO nonsingular for (α, β) = (2, 0.5), (2, 1.2) and (3,2).The singular cross section is not shown separately as its shape can be derived by comparing the total NLO with the NLO nonsingular cross section.It happens to be constant in the SCET + region in this double-logarithmic plot.As expected, the nonsingular is relevant in the fixed-order region of phase-space, where the angularities are large.The feature which we observe in the fixed-order region is the two-dimensional analogue of the bump we saw for one angularity.It occurs for the WTA axis, where regions in phase space in which two particles carry the same energy fractions (x 1 = x 2 or x 2 = x 3 ) lead to those sharp edges.

Scales and uncertainties
In this section we specify the central scale choices used to achieve resummation, and describe in detail the scale variations used to estimate the perturbative uncertainty.Since the spectra we are dealing with have distinct kinematic regions where resummations must be handled differently, we use angularity-dependent soft and jet scales given by "profile functions", a method previously applied to e.g. the thrust event shape [26] and the B → X s γ spectrum [55].We start with the single angularity distribution and then extend our discussion to the case of two angularities, introducing profile functions depending simultaneously on both e α and e β .The hard scale is µ H = Q and will not be varied to estimate the resummation uncertainties.

Single angularity
The canonical scales for the single-e β resummation are with L β ≡ log 10 e β .Given these expressions, we find it convenient -in particular in view of the case of two angularities discussed later -to construct profile functions in terms of the logarithms of the angularities, rather than the angularities themselves.
For the central value of our predictions we take The function (2.44) smoothly connects the canonical region t ≤ t 1 to the fixed-order region t ≥ t 3 using a cubic polynomial.
For t 1 we take the value of e β where the NLO nonsingular cross section is 10% of the NLO singular cross section (see fig. 2).This leads to (β, t 1 ) = (0.5, −0.795) , (1.2, −0.82) , (2, −0.90) , (3, −0.98) . (2.45) For the t 3 parameter, we take the point where the NLO singular vanishes, which is for all angularities.To simplify our scale choices in the double differential case, we do not introduce a profile to handle the transition to the nonperturbative regime but instead freeze α s below 2 GeV to avoid the Landau pole. 11e now consider a range of scale variations to estimate the perturbative uncertainty.
(i) Fixed-order uncertainty: We simultaneously vary all scales µ i , including µ H , by a factor of 2 or 1/2.
(ii) Resummation uncertainty: Following refs.[56,57], we vary the jet and soft scale according to log 10 (µ vary where we take the following values for the parameters a and b and This form of f vary corresponds to a factor 2 variation in the canonical region but no variation in the fixed-order region, where we are not allowed to vary µ H , µ J and µ S independently of each other, and a smooth transition in between. For a = b = 0, eq.(2.47) reproduces our central scale choice.Scale variations with a = 0 and b = 0 preserve the canonical relation We choose the size of these variations in eq. ( 2.48) such that the smallest scale varies by a factor of 2 or 1/2 in the canonical region. 12Setting a = 0 and b = 0 does not preserve eq.(2.50).We also impose that these variations vanish for β → 1, since µ S and µ J coincide in this limit, and agree with the choice for thrust in ref. [56].Furthermore, the deviations from eq. (2.50) for β > 1 are required to be of the same size as the deviations from (iii) Variations of the transition points: Since there is also a certain amount of arbitrariness in choosing the transition points t 1 and t 3 , these get varied as well (but only one at a time).For t 1 we consider the following alternatives to eq. (2.
This corresponds to the point where the total NLO cross section vanishes.
Our final uncertainty band is obtained by adding the fixed-order uncertainty in quadrature to the resummation uncertainty, which is obtained by taking the envelope of the a and b variations above, and the variations of the transition points.The resummation uncertainty is dominated by the a and b variations, whereas the uncertainty from the variations of the transition points is rather small (but still contributes to the envelope in parts of the transition region).
Table 2. Canonical scales for the measurement of two angularities in the three regions, with L α ≡ log 10 e α and L β ≡ log 10 e β .

Two angularities
Next we consider the case where two angularities are measured.We will first construct running scales for the central value of our prediction, using the canonical values of the scales in the three regimes listed in table 2 as a starting point.Although in regime 1 µ S ∼ e α Q ∼ e β Q, we choose µ S ∼ e β Q as our canonical scale because the resummation in this regime is governed by e β (and conversely the jet scale in regime 3 involves e α ).The collinear-soft scale merges with the soft scale in region 1 and with the jet scale in region 3.
We take for the jet and soft scale and we fix the collinear-soft scale using the canonical relation .55)The transition between the SCET regions is handled by the function g, (2.56) In the intermediate region, x 1 ≤ x ≤ x 2 , g is given by the cubic polynomial that is continuous and has a continuous derivative.The first three arguments of g in eq.(2.54) directly follow from the canonical scales in table 2. The transition points x 1 were chosen as a fraction t R1/R3 of the total distance (in logarithmic space) between the two phase-space boundaries in eq.(2.40).For the central profiles we choose t R1 = 0.8 and t R3 = 0.95, which corresponds roughly to the region where the nonsingular terms (from the boundary regimes) are 10% of the singular one.The transition points x 2 were chosen at the phase-space boundary.For example, for µ S we start the transition at ln e α − log 10 2 (α−β)/β e α/β β = tR1 log 10 e β − log 10 2 (α−β)/β e α/β β , (2.57) and end at e α = e β .With the definition above, the profile scales remain constant at their canonical regime 3 values, beyond the NLO phase-space boundary, such that the e α resummation is turned off here.This implies in particular that also the NLL results know about the NLO phase-space boundary.We will add a comment below, how our results change if the canonical (instead of NLO) phase-space boundary would have been used in the profile scales.
The transition to the fixed-order region is controlled through the function h in eq.(2.44).Due to the argument of h in eq.(2.54), the fixed-order region min(L α , L β ) ≥ t3 (2.58) has a square shape.We have checked that other choices have minimal impact on the result.The transition points are taken from the single angularity case: for t1 we take the minimum (which corresponds to a larger transition region) of t 1 for the single angularities e α and e β from eq. (2.45), and t3 = t 3 = −0.33.
As for the single angularity spectrum, several scale variations are taken into account.
(i) Fixed-order uncertainty: We simultaneously vary all scales µ i by a factor of 2 or 1/2.
(ii) Resummation uncertainty: Extending the one-dimensional case, jet, collinear-soft and soft scales are varied according to with a(α) = ± min(α, 1), b(α) = ±(α − 1)/(3α), as in the single angularity case in eq.(2.48).The transitions are governed by the functions which have the property that h J = 1 in regime 1 and 2, and 0 in regime 3, and h S = 1 in regime 2 and 3, and 0 in regime 1.Thus the scale variations in regime 1 and 3 in eq.(2.59) are the usual single angularity ones.In these regimes the collinear-soft scale is not independent and thus needs to be varied in tandem with the soft or jet scale it has merged with.In the intermediate regime we have used eq.(2.55) to determine the collinear-soft scale, but setting b = 0 there.This is necessary, because otherwise µ S is varied by much more than a factor of 2 when the angularity exponents α and β are close to each other.
(iii) Variations of the transition points: We vary t1 , using the maximal and the minimal value of the t 1 variations of the two single angularities considered in sec.2.7.1.Similarly, for t3 we use the variation of each of the single angularities.To vary the transition between the boundary theories, we vary tR1/R3 by taking tR1 = 0.7, 0.9 or tR3 = 0.9.These values are motivated by looking at the contour where the nonsingular terms are 10% of the singular one (focusing on the resummation region).
As in the one-dimensional case, the total uncertainty is obtained by adding the fixed-order uncertainty in quadrature to the envelope of the resummation variations and the variations of the transition points.

Differential vs. cumulant scale setting
We implement our choice of scales at the differential level, i.e. we calculate d 2 σ/(de α de β ) using scales evaluated at e α and e β .An alternative is to use cumulant scale setting, evaluating the scales at e c α and e c β .Differentiating this introduces derivatives of the scales, For the unprimed orders in table 3, such as NNLL, differential scale setting does not capture all the logarithms, as discussed in detail in e.g.ref. [23].However, our scales in secs.3. Perturbative ingredients needed at different orders in resummed perturbation theory.The columns correspond to the loop order of the fixed-order ingredients, the non-cusp and cusp anomalous dimensions, and the QCD beta function.
2.7.2 undergo fairly rapid changes in transition regions, leading to artefacts from the terms involving the derivatives of scales, when using the cumulant scale setting.We investigate this issue by supplementing our cross section with differential scale setting with the additional terms on the right-hand side of eq.(2.62).By using the canonical scales to determine the scale derivates d ln µ i,j /de α,β in these terms, we maintain the required formal accuracy while avoiding artefacts from derivatives of our profile scales encountered with cumulative scale setting.For example, in region 2 In fig. 4 we compare the standard differential scale setting (left panel) and cumulative scale setting in eq.(2.61) (right panel) to the alternative procedure we just described (middle panel).The cumulative scale setting leads to clear artefacts in the transition to fixed-order and to the boundaries of phase space, which are due to the derivatives of profiles scales in eq.(2.62), which undergo a rapid transition.For example, the boundary of the box in eq.(2.58) is clearly visible.Our alternative approach avoid these artefacts, by using canonical scales in the derivatives of scales.However, the alternative approach has a major disadvantage: Since the canonical scales do not turn off properly in the fixed-order region, the singular-nonsingular cancellation is spoiled there.Thus we are left with using standard differential scale setting in the results presented in sec.3, even though not all logarithms are captured.

Matching
Given that we have a different factorization theorem for each of the regions of phase space, we would like to obtain an expression for the cross section which is valid everywhere.This is achieved by matching the cross section predictions from the various regions [58] where each of the cross sections is differential in e α and e β .The first line describes the cross section for regime 2, using the scales which are appropriate for this regime.The second line ensures that in regime 1 we reproduce the correct cross section.This is achieved by including the nonsingular contribution obtained by adding the cross section of regime 1 and subtracting the one of regime 2 evaluated at the scales of regime 1 (i.e. the overlap).Note that in regime 1, the R2 scales merge into the R1 scales, such that the second term on the second line cancels against the first line.This procedure is similar to the construction of the fixed-order nonsingular when a single type of logarithm is resummed.Similarly, the third line describes the nonsingular correction from the regime 3 boundary of phase space.The last line corresponds to the fixed-order nonsingular, shown in fig. 3 above.A smooth transition between the regimes is achieved by the profile scales discussed in sec.2.7.
In fig. 5 we show the contributions (besides the fixed-order nonsingular, already discussed in sec.2.6.2) which make up the total NNLL cross section. 13As shown in the left panel, in the bulk of the phase space the cross section is already captured by the regime 2 cross section.The nonsingulars from regimes 1 and 3 correct the regime 2 cross section close to the phase space boundaries and cause the cross section to vanish outside the boundaries.

Nonperturbative effects
We have also studied the effect of nonperturbative corrections, which we first discuss for single angularities before extending to the double angularity cross section.We restrict to e α with α > 1 because only in this case does the soft function capture the (dominant) nonperturbative corrections.We can factorize the soft function [55,59,60] into its perturbative contribution S pert and nonperturbative contribution F .F is dominated by momenta of the order Qe α Λ QCD , and its integral must be one, since nonperturbative effects do not change the total cross section.Expanding eq.(2.65) for Qe α Λ QCD , where the leading nonperturbative correction is characterized by the parameter Ω α , with a calculable dependence on α [61][62][63][64] 14 We take Ω = 0.323 GeV with 16% uncertainty [26].Since a shift is rather crude, we implement nonperturbative effects in our analysis using the following functional form for the nonperturbative contribution F in eq.(2.65)15 which is normalized and has the first moment required by eq.(2.66).
We have tested eq.(2.67) using Pythia, applying two methods to extract Ω α .Method 1: Ω α is obtained by taking the difference of the first moment of Pythia cross sections at hadron-and parton-level, This follows directly from eq. ( 2.65) and the definition of Ω α in eq.(2.66), but it assumes that the convolution in eq.(2.65) is also valid for the nonsingular cross section.The resulting distribution for Ω α , shown in fig.6, approximately exhibits the α-dependence of eq.(2.67).However, it clearly breaks down for large values of α, where Ω α becomes negative.One possible explanation is that the above assumption on the nonsingular cross section is not justified.
We have therefore attempted to extract Ω α in a second way.Method 2: We performed the convolution of the parton-level Pythia prediction with eq.(2.68) and determined Ω α by minimizing the distance between the resulting distribution and Pythia's prediction at hadron level 16 .This also approximately exhibits the α-dependence in eq.(2.67), but now overshoots it for large values of α.Note that this method also relies on the assumption that eq.(2.65) extends to the nonsingular cross section.Moving on to two angularities, in regime 2 and 3 the nonperturbative effects for e α arise from the soft function S(Qe α , µ) discussed above.In regime 1, we encounter the double differential soft function, which can be factorized in a way similar to eq. (2.65) (2.70) The leading nonperturbative corrections take on a particularly simple form since nonperturbative correlations vanish at this order.In our numerical analysis we set where the effect of nonperturbative correlations are encoded in (2.73) F cor was imposed to be a polynomial of degree 2 in k 1 and k 2 that introduces correlations such that F in eq. ( 2.72) remains normalized and produces the first moments required by eq.(2.71).
We explored correlations by varying the size of the correlation parameter c ∼ 1/Λ 2 QCD .In regime 2 and 3 the nonperturbative effects involving e β are suppressed because Qe β Qe α .We will nevertheless use eq.( 2.72) in these regimes as well.For the leading nonperturbative correction, this seems reasonable from the point of view of continuity.The cross section for the ratio of angularities is particularly interesting, because nonperturbative corrections contribute to any value of the ratio, since this integrates over a line that goes through (e α , e β ) = (0, 0).

Single angularity
We start by presenting results for the cross section of a single angularity e β in fig. 7. Shown are our predictions at NLL, NNLL and NNLL+NLO order (defined in table 3) for angularity exponents β = 0.5, 1.2, 2, 3 with Q = 1000 GeV.The bands show the perturbative uncertainty estimated by varying the profile scales, as described in sec.2.7.1.Our predictions for the central curves are normalized to 1. 17 The variations are not normalized to 1, but rescaled by the same amount as the corresponding central curve.As expected, the uncertainty bands reduce at higher orders, and overlap between the different orders over most of the range.The one exception is the NNLL vs. NNLL+NLO in the fixed-order region.This is not surprising, because in this region the matching with NLO cannot be neglected.The bump in the spectrum at large values of e β comes from the matching with NLO and is due to our choice of using the WTA axis, as observed already in sec.2.6.2.Interestingly, for smaller values β the range of e β values gets squeezed, such that there is a fairly rapid transition from nonperturbative region at small e β to the fixed-order region at large e β .
In fig.8 we show our NNLL+NLO result for β = 1.2 and 3 with and without nonperturbative corrections, included using the procedure described in sec.2.9.We compare these to Pythia with and without hadronization.We also show the ratio with NNLL+NLO to make it easier to distinguish these curves.The effect of hadronization on our perturbative prediction is very similar to the difference between Pythia at the parton and hadron level.The curves are not the same, but this difference is already present before including nonperturbative effects.Interestingly, although Pythia does not have a clear bump in the fixed-order region, there seems to be some feature at the same spot.
The corresponding plot for Q = 91.2GeV is shown in fig.9.Here we added also a nonperturbative uncertainty band, which was obtained by varying Ω within its uncertainty [26] and adding in quadrature the envelope of the variations obtained by considering with a = 1, 2, 3 and 4. F (Qe α , 1) coincides with F (Qe α ) in eq.(2.68).These alternative functional forms F are all normalized and have the same leading nonperturbative correction, thus probing the effect of subleading nonperturbative effects.Indeed, the uncertainty band in fig. 9 grows significantly at small values of the angularity, because of the sensitivity to the shape of F and not just its first moment.For Q = 1000 GeV this uncertainty is very small, which is why we do not show the corresponding plot.

Single angularity distributions from Event2
To test the factorization framework, especially for the WTA axis choice, we compared the fixed-order expansions from our resummed single differential cross sections against numerical results from the Event2 generator [66].For this purpose, we ran Event2 with n f = 5 and an infrared cutoff ρ = 10 −10 and generated one trillion events.To be explicit about what is compared here, we write the expansion of the cross section as In fig. 10, we plot the difference between the Event2 output and our singular contributions to the NLO and NNLO coefficients A and B for angularity exponents β = 1.2 and 2.Here we consider both the thrust axis and the WTA axis.Assuming that recoil effects can be ignored, the only difference at this order can be traced back to the constant in the (cumulative) oneloop jet function, which for the thrust axis was calculated in ref. [22].In app.B we collect our (N)NLO singular results for A and B for several angularity exponents and both axis choices.For the WTA case, the difference between Event2 and our singular cross section goes to zero at small values of the angularity (within statistical uncertainty) for both β = 1.2 and 2, while for the thrust axis, this does not happen for β = 1.2.In this case, instead, the difference approaches a straight line, which corresponds to the single logarithms one expects from recoil.Interestingly, the turn off of the nonsingular contribution takes place substantially faster for the WTA axis.
At very small values of e β , the comparison breaks down due to infrared cutoff effects in Event2.More specifically, Event2 regulates infrared divergences by cutting on the invariant mass of pairs of partons, (p i + p j ) 2 > ρ Q 2 .By applying this prescription to the SCET modes for the single angularity distribution, we conclude that Event2 is expected to deliver reliable results for values of e β down to about ρ min(β/2,1) at NLO, and about ρ min(β/2,1/2) at NNLO.The further restriction at NNLO stems from the fact that at this order two soft emissions arise.We stress that this is simply an order-of-magnitude estimate, and judging from our numerical results, the true cutoff seems to be somewhat higher.

Two angularities
In fig.11 we show our results for the normalized cross section differential in the angularities e α and e β at NLL and NNLL+NLO order, compared to Pythia (parton level), where we again take Q = 1000 GeV.The difference between the NLL and NNLL+NLO is not very large, except in the fixed-order region.However, as is clear from our one dimensional plots, the uncertainties at NLL are pretty large.Indeed, the only reason that the cross section vanishes at NLL at the NLO phase-space boundaries is simply due to our choice of profile scales.If we would have turned off our profile scales at the canonical boundary instead, the peak region of the NLL cross section would be broader and extend (slightly) over the NLO phase-space boundary (dashed line).As discussed in sec.2.6.2, the sharp feature that the NNLL+NLO cross section exhibits in the fixed-order region is analogous to the bump of the single differential distributions, and is due to our choice of using the WTA axis.For (α, β) = (2, 0.5) the peak of the distribution is close to the phase-space boundary corresponding to regime 3, while for the other angularity combinations it sits more in the middle between regime 1 and 3.
Comparing our results to Pythia, we see that the Pythia cross section is closer to our NNLL+NLO than NLL cross section.There are however notable differences: In Pythia there is no sharp feature in the fixed-order region.Although it is expected that this would be somewhat washed out in Pythia, it is surprising that there is no visible remnant (the corresponding bump for the one-dimensional distribution is still noticeable for Pythia in fig.7).The largest difference between the Pythia and the NNLL+NLO distribution is for (α, β) = (2, 0.5).Surprisingly, the peak for Pythia extends outside the NLO phase-space boundary in eq.(2.40).That there are large differences in this case is not so surprising, because we have already seen that for β = 0.5 the resummation region gets squeezed such that there is a quick transition between the fixed-order region and the nonperturbative region.
The uncertainties are taken from the scale variations for the two-dimensional distributions using the procedure outlined in sec.2.7.2.As for the single angularity distributions, we have normalized the central curve, and rescaled the scale variations with the same factor.However, note that unlike for the single angularity case, the resummation region can contribute to all values of r.It is reassuring to see that the uncertainties decrease at higher orders, and the uncertainty bands overlap.The reason that the uncertainty is so large (α, β) = (2, 0.5), is that the peak is close to the region 3 boundary, so almost all r values are affected by large resummation uncertainties (see fig. 11).Given how close the central curves are, compared to the size of the uncertainty bands, our estimate is probably quite conservative.The nonperturbative effect on the cross section differential in r is shown in fig.13 for (α, β) = (3, 2).Referring to sec.2.9 for the notation, the uncertainty band here includes both the c-variation (within -2 to 2), 19 the variation of Ω within its uncertainty of 16% and of F (with a = 1, 2, 3, 4) as described for the single angularity case.We constructed a separate envelope for each of them and added these three uncertainties in quadrature.For Q = 1000 GeV, the correlations probed by c dominate the uncertainty.For Q = 91.2GeV, the subleading nonperturbative corrections estimated by varying a are the largest instead, for

Conclusions
In this paper we presented our calculation of the cross section for e + e − → hadrons differential in two angularities.We simultaneously resummed the logarithms of each angularity, employing the SCET + framework we developed in ref. [10].The resummation was performed at NNLL accuracy and matched to NLO, thereby obtaining a prediction that is valid throughout the phase space.By using exclusive k T clustering with the WTA recombination scheme, we could ignore the issue of recoil.We performed a detailed numerical study, assessed the perturbative uncertainties through variations of each of the various scales entering factorization, and studied the impact of the leading nonperturbative corrections.
The one-loop matching with the full QCD calculation shows that our SCET + factorization correctly captures the singular limit at this order.We extended this check for the factorization theorem of a single angularity to O(α 2 s ) by using Event2, and found agreement.We also showed that the effect of recoil can not be ignored for the angularity exponent β = 1.2, highlighting the advantage of the WTA axis.In the fixed-order region the cross section has some sharp features.These arise because the position of the WTA axis can change abruptly depending on the precise momentum configuration, as we are no longer in the dijet region.
We have tested the perturbative convergence of our resummed calculation, finding that the uncertainty bands at higher orders become smaller and (mostly) overlap with those at lower orders.For the double angularity distribution, Pythia seems closer to our NNLL+NLO prediction than our NLL prediction, though the sharp feature in our predictions that arises in the fixed-order region is washed out.Of course a benefit of our calculation is that it provides an estimate of the perturbative uncertainty, and is systematically improvable.
We also considered the cross section differential in the ratio of two angularities, which is not infrared safe but still Sudakov safe.This is interesting to investigate because many further jet substructure observables fall in this category, and calculations have been performed at most up to NLL accuracy (with the exception of the ratio τ 2,1 of 2-to 1-subjettiness with angular exponent 2 [67]).Since the resummation region contributes to most of the plot range for the angularity ratio, the uncertainty on the cross section is larger than for the single angularity measurement, but still reasonable.As may be expected, nonperturbative corrections similarly play a more important role.
In this paper we restricted ourselves to two-jet production in e + e − collisions to have a clean theoretical setup, but it is our goal to extend this analysis to the measurement of (multiple) angularities of jets in LHC collisions at NNLL+NLO.At this accuracy, only very few jet substructure observables have been calculated so far.Our framework allows us to reliably account for correlations between jet observables, and demonstrates the feasibility of performing higher-order resummation for more differential measurements.
The coefficients of the cusp anomalous dimension that enter in eq.(A.1) are [68] Γ 0 = 4C F , The coefficients for the non-cusp anomalous dimension for the hard function are [49,50] γ (A.6) The other non-cusp anomalous dimensions follow from eq. (2.29).

B NLO and NNLO singular terms in the single angularity distribution
The fixed-order single angularity distribution can be written as In our comparison against Event2, we also analyzed angularities with respect to the thrust axis, with exponents β = 1.2, 2, 3.The corresponding NLO coefficients A thr sing (L β ) coincide with the ones calculated with respect to the WTA axis.Differences first appear at NNLO and are due to the non-logarithmic terms in the one-loop cumulative jet function, which for the thrust axis can be obtained from ref. [22].Thus in the NNLO coefficients B thr sing (L β ) the only changes are For the thrust case, the coefficients A thr sing (L 2 ) and B thr sing (L 2 ) agree with the well-known results from the literature (see e.g.ref. [25]).

Figure 4 .
Figure 4. NNLL cross section with differential scale setting (left), with differential scale setting plus extra NNLL terms in eq.(2.62) evaluated using canonical scales (center) and with cumulative scale setting (right).

Figure 5 .
Figure 5. NNLL cross section in regime 2 (left) and nonsingular corrections from regimes 1 (center) and 3 (right), corresponding to the second and third line of eq.(2.64).

Figure 8 .
Figure 8. NNLL+NLO without and with nonperturbative effects, compared to Pythia at parton and hadron level, for β = 1.2 (left) and β = 3 (right).The bottom row shows the ratio with the NNLL+NLO cross section.

Figure 9 .
Figure 9. NNLL+NLO without and with nonperturbative effect, compared to Pythia at parton and hadron level for Q = 91.2GeV and β = 3.The band provides an estimate of the nonperturbative uncertainty, as described in the text, which is sizable due to the smaller value of Q.

Figure 10 .
Figure 10.Difference between the NLO and NNLO terms for the single angularity cross section calculated by Event2 and our singular results.For β = 1.2 with the thrust axis, the absence of a plateau is due to recoil effects.

Figure 13 .
Figure 13.Distribution for the ratio of two angularities from the NNLL+NLO cross section with and without nonperturbative effects, for (α, β) = (3, 2), compared to Pythia at parton and hadron level for Q = 1000 GeV (left) and Q = 91.2GeV (right).The band indicates the uncertainty from nonperturbative effects, as described in the text.
9 T F n f ,