UvA-DARE (Digital Academic Repository) How much joint resummation do we need?

,


Introduction
Measurements at colliders often impose restrictions on QCD radiation through e.g.jet vetoes, transverse momentum measurements, or production at threshold.One particular example, that will play a prominent role in this paper, is the family of event shapes for e + e − collisions called angularities [1].These are defined as where Q is the center-of-mass energy, and the sum runs over all final-state particles i with energy E i and angle θ i with respect to a chosen axis.The cross section integrated over e α up to some cut-off e cut α , known as the cumulative cross section, contains a series of logarithms L = log 10 e cut α at each order in perturbation theory, schematically

JHEP10(2019)130
Rows correspond to different orders of fixed-order perturbation theory, denoted by leading order (LO), next-to-leading order (NLO), etc.For e cut α 1, the α s expansion breaks down, because we cannot treat L ∼ 1.In this case we want to sum columns, which correspond to different orders of resummed perturbation theory, called leading logarithmic (LL) order, next-to-leading logarithmic (NLL) order, etc.To be precise, we will treat α s L ∼ 1, which corresponds to resumming logarithms "in the exponent".
This resummation can be carried out either by using a Monte Carlo parton shower, such as Herwig [2][3][4] or Pythia [5], or through analytical methods.The advantage of parton showers is that they allow for any measurement on the fully exclusive final state.However, their formal accuracy is limited to LL order in the large N c limit, and it is challenging to systematically go beyond this order.For a recent discussion of the logarithmic accuracy of the parton shower, see refs.[6][7][8].Some recent improvements in parton showers are the inclusion of higher-order splitting functions [9][10][11][12], corrections to the large N c limit [13][14][15], spin correlations [16], and the simultaneous treatment of small x and collinear and soft logarithms [17].
In principle, one could imagine the simultaneous resummation of ever more observables, which would lead to an increasingly precise parton shower.This is important for the burgeoning field of Machine Learning in jet substructure, see ref. [52] for a review.Here, samples from Monte Carlo parton showers are often employed, thus raising the question to what extent discrimination is based on features of physics or of the Monte Carlo.There has also been some work on approaches that do not require labeled samples though, see e.g.refs.[53][54][55][56].We were inspired by ref. [57], which investigates how sensitive Machine Learning is to details of the final state, studying the discrimination of jets from boosted, hadronic decays of Z bosons from jets initiated by QCD processes.By using a complete basis of observables that probe the N -body final state and increasing N , they find that the discrimination saturates at N = 4.
We are interested in asking how much resummation is needed to reliably describe the jet, focussing on QCD-initiated jets.For simplicity, we restrict ourselves to observables that are azimuthally symmetric, for which the angularities form a basis.Specifically, we consider a set of N angularities e α i , with i = 1, . . ., N , derive the resummed prediction for the cross section differential in a subset of angularities denoted by I = {i 1 , . . ., i n } at -2 -

JHEP10(2019)130
NLL accuracy, and investigate the degree to which any of the other angularities e α j , with j / ∈ I, can be predicted.To this end, we generate events from flat 1 k-body phase space, using Rambo [58] "on diet" [59].We reweigh these events to match the multi-differential cross sections from the input, and then calculate the distributions for other angularities from these reweighed events.As an alternative to analytic resummation, we also use the cross section differential in n angularities from Herwig or Pythia as basis for reweighing.The dependence of this procedure on the number n of angularities that have been jointly resummed is investigated and the optimal subset I of angularities to be used as input for any given n is determined.
Our main conclusion is that there is an order of magnitude improvement for reweighing with n = 2 angularities over n = 1, but that the advantage quickly diminishes for larger values of n.We have investigated the dependence on the input (Herwig, Pythia, analytic), number of particles k of the flat phase space, center-of-mass energy Q, and set of angularities under consideration.None of these lead to a qualitatively different behavior.
The outline of this paper is as follows: we start in section 2 by describing our setup and discussing in detail how we perform the reweighing, determine the optimal input set of angularities I, and estimate the statistical uncertainty.In section 3, we present our analytic calculation of the cross section differential in n angularities, limiting most of the discussion to NLL accuracy.Our main results are presented in section 4, with additional plots relegated to appendix A. We conclude in section 5.

Optimal reweighing procedure
We start this section by presenting our setup, discussing the observables in more detail.We then describe the reweighing procedure of flat k-body phase space with resummed predictions for n angularities, as well as the determination of the optimal set of angularities to be used as input, including the treatment of statistical uncertainties.

Setup
We use the process e + e − → dijets as a case-study for our procedure.The final state of the collision is clustered into two jets using the exclusive k T [60] jet algorithm with the Winner-Take-All (WTA) recombination scheme [61,62].Following ref. [63], we modify the original definition of the angularities in ref. [1] for large angles θ i ∼ 1, as described by eq.(1.1).The angle θ i of each particle contributing to an angularity is measured with respect to the WTA axis of the corresponding jet, and the angularities e α that we consider are the sum of the angularities e Jet J α of the individual jets, i.e. e α = e Jet 1 α + e Jet 2 α .The variable α i ≡ log 10 e α i is used instead of the angularities themselves, which is more natural since the angularity distributions are peaked at small values of e α i . 2  The various distributions of α i are constructed in 16 bins on the interval α i ∈ [−4, 0].We have verified that the chosen binning does not alter the conclusions of our study.
1 "Flat" here implies that there are no preferred points in phase space, so that, up to four-momentum conservation, each point is assigned the same probability.
2 For brevity, we will also refer to the observables α i as "angularities".

JHEP10(2019)130
The binned distributions are constructed either by using analytic predictions (described in section 3) or through 10 6 parton-level events generated by the Herwig 7 or Pythia 8 general-purpose Monte Carlos (described in section 2.2).

Optimal reweighing
To obtain the binned distributions for flat massless k-body phase space, with Q µ = (Q, 0, 0, 0) in the center-of-mass frame, we employ the Rambo technique of refs.[58,59], with a slight modification that improves sampling in the collinear and soft regions through weighted events.In particular, we perform a transformation that distributes the first random number of "Algorithm 1" in ref. [59] logarithmically, with the weight of this event given by the Jacobian.This ensures that the phase space is sampled sufficiently to obtain statistically reliable predictions for small values of the angularities that we consider.We start by reweighing the flat phase space result for n + 1 angularities by the cross section differential in n angularities α i with i ∈ I = {i 1 , . . ., i n }.By integration3 over the n angularities α i , a reweighed cross section dσ reweigh /d α j differential in the remaining angularity α j with j / ∈ I is determined, i.e.
appropriately applied to the binned distributions.Comparisons between the resulting reweighed cross sections and direct determinations of equivalent cross sections can be found in figures 5 and 6.
We define a goodness-of-fit measure for the reweighed angularity distribution by comparing to the resummed distribution for α j To find the optimal set I when reweighing with n angularities, we introduce a global goodness-of-fit variable, defined as the average χ 2 α j of the angularities j / ∈ I, Here N − n is the number of angularities not used as input.
Up to n = 3, we search for the global minimum of χ 2 as a function of the set of reweighed angularities I, denoted by χ 2 min .We will refer to this as the optimal set of angularities in each case, and the corresponding value for χ 2 α i will be denoted by χ 2 α i ,min .

JHEP10(2019)130
For n = 4 and n = 5, we start from the optimal I for n − 1 angularities and iteratively determine the optimal additional angularity to add to I. We have verified that the result of this iterative method is already very close to the global minimum of χ 2 for n = 3.The value of the global minimum χ 2 min is rather sensitive to statistical fluctuations, owing to the finite size of the Monte Carlo samples distributed over 16 n+1 bins. 4Furthermore, it is not expected to follow a Gaussian distribution.To obtain an estimate of the statistical uncertainty, we perform the reweighing procedure over 11 replicas of the event samples.The median of χ 2 min is taken as our central prediction, and the spread of the 7 most central replicas as a reasonable approximation for the spread corresponding to roughly one standard deviation.Results for χ 2 min as a function of n are depicted in figure 7 and in appendix A.1.

Joint resummation of n angularities at NLL
In this section we present our framework for performing the joint resummation of n angularities.We start by drawing the Lund planes in section 3.1, from which the leadinglogarithmic resummation immediately follows.These diagrams also allow us to identify the modes in SCET, for which the factorization formulae are presented in section 3.2.From the renormalization group equations for these factorization formulae we obtain the resummed cross section at NLL accuracy in section 3.3.We describe the matching of the different factorization formulae for different regions of phase space in section 3.5, for which we are guided by the size of the power corrections, estimated in section 3.4.

Lund diagrams and phase-space boundaries
In the collinear limit, the probability P of a particle emitting real radiation can be characterized by the momentum fraction z and angle θ of the radiated particle j with respect to its emitter i, At LL accuracy, only the ∼ 1/z term of the splitting function P i→j (z) is kept, so these emissions are uniformly distributed in x ≡ log 10 (1/θ) and y ≡ log 10 (1/z).The Lund plane [64] spanned by these variables is shown in figure 1, with some emissions indicated by crosses.By identifying 2E i /Q → z i and θ i /2 → θ (the factor of 1/2 is purely for convenience and doesn't affect the leading logarithms), the angularity in eq.(1.1) for a single emission with small z and θ can be written as corresponding to a straight line in the Lund plane with slope −α.Due to their uniform distribution in the (logarithmically-spaced) Lund plane, a single emission will dominate the measurement at this accuracy, as indicated by the red line in figure 1.All emissions JHEP10(2019)130 There can be no emissions in the shaded area, while emissions above this only contribute at higher order.The green and orange dot denote the collinear and soft mode respectively.
above and to the right of this line are more soft or collinear and only enter beyond LL accuracy.There are no emissions below the line, otherwise these would be dominant.The shaded area under the line corresponds to the Sudakov factor describing the no-emission probability, and can be used to calculate the cumulative cross section, where σ0 is the Born cross section and C i = C F (C A ) for quark (gluon) jets.Interestingly, the relevant degrees of freedom in SCET correspond to the points that describe the edges of the shaded region.For the measurement of a single angularity these are indicated by the orange and green dot in figure 1, and correspond to soft and collinear modes respectively.The parametric scaling of the momenta of these modes is most conveniently expressed in terms of lightcone coordinates defined through Indicating the momenta of the soft, n-collinear and n-collinear modes by p µ s , p µ n and p µ n respectively, their scaling is found to be [65] When the simultaneous measurement of two angularities α 1 and α 2 is considered, the straight lines describing the variables in the Lund plane have to cross one another at some point to ensure that the cross section depends on both measurements.Assuming the JHEP10(2019)130 hierarchy α 1 > α 2 for definiteness, three distinct cases can be distinguished, as shown in figure 2. The boundaries of these three regions of phase space for two angularities are Regime 1: which agree with the regions of phase space identified in refs.[40,44,45].In all three cases there are soft (orange) and collinear (green) degrees of freedom.The intermediate regime 2 has an additional collinear-soft mode (blue), which contributes to both measurements since it lies on the intersection of both lines.This method of finding all relevant regions of phase space can be generalized to the simultaneous measurement of an arbitrary number of angularities.There is only one additional subtlety that has to be taken into account when more than two angularities are considered, which we illustrate in figure 3 for three angularities with parameters α 1 > α 2 > α 3 .If the line corresponding to α 3 were to be placed above the position indicated by the dotted line, the angularity α 2 would no longer be connected to the boundary of the region in which emissions are forbidden and hence not affect the cross section.The point at which the dotted line crosses the y-axis is given by The phase space of a cross section involving an arbitrary number of angularities n can be divided into various regimes, as listed explicitly in eq.(3.6) for n = 2.The regime in which the largest number of independent logarithms occur, is the one for which the edge of the forbidden (gray) region in the Lund plane involves every line corresponding to an individual angularity.For n = 2, this corresponds to the center panel in figure 2 and for n = 3, this situation is depicted in figure 3.This region will be denoted by R n (α 1 , . . ., α n ), JHEP10(2019)130 , which shows up in the boundaries of the region of phase space.and its boundaries in phase space are given by y-conditions: The first line consists of the n − 1 conditions on the hierarchy between the points at which each line in the Lund plane crosses the y-axis.The second line contains the n−1 conditions on the hierarchy between the points at which the lines cross the x-axis.As the conditions consist solely of inequalities, this region in phase space is n-dimensional and will be called the "bulk".
Regions that involve fewer logarithms can be obtained by raising or lowering the point where an angularity crosses either axis in the Lund plane, such that two modes (the colored dots) overlap.This can be seen explicitly in figure 2 by starting from the center panel and raising − α 2 until it reaches − α 2 = − α 1 in the right panel, sliding the mode indicated by the blue dot up to the orange dot in the process.In full generality, the boundaries of a region in phase space involving a subset of angularities β 1 , . . ., βm with m < n and {β 1 , . . ., β m } ⊂ {α 1 , . . ., α n } are found to be5 y-conditions:

JHEP10(2019)130
where the B-conditions (boundary-conditions) contain all restrictions on the angularities that are only connected to the boundary of the shaded area in the Lund plane through a single point, i.e. the angularities not involved in the region.As any such regime is characterized by n − m equalities, it represents an m-dimensional region in the n-dimensional phase space, denoted by R m (β 1 , . . ., β m ).By considering all possible combinations of angularities it then follows that there are n n−m distinct regions of dimension m in the phase space of n angularities.

Factorization formulas
The analytical resummation will be performed by making use of the Soft-Collinear Effective Theory (SCET) [26][27][28][29], which describes the infrared limit of QCD.The relevant degrees of freedom are determined by the process and measurements under consideration.The version of SCET that involves the collinear and soft modes in eq.(3.5), known as SCET I , correctly describes the regions of phase space dominated by a single angularity, e.g. the left and right panels in figure 2. Regions of phase space involving multiple angularities (such as the middle panel in figure 2) contain additional collinear-soft modes and are correctly described by SCET + [40,42,43,66].
As the various modes in SCET are decoupled at the level of the Lagrangian [29], cross sections may be factorized into products or convolutions of perturbative functions as long as the contributions of the various modes to the measurements can be shown to factorize as well.In general, each of these functions contains logarithms of the ratio of its inherent, natural scale and the common scale µ.By solving their RGEs, they may be evaluated at their natural scales (where the logarithms are minimized) and then evolved towards a common scale µ, resumming all the large logarithms in the process.
The relevant degrees of freedom for the bulk regime for n = 2 angularities are shown in the middle panel of figure 2. The orange dot represents the (ultra)soft mode, the green dot the collinear mode, and the blue dot corresponds to a collinear-soft mode [40,42], which contributes to the measurement of both angularities.For our process of interest, e + e − → dijets, there are two distinct collinear directions corresponding to the two jets, and hence also two corresponding collinear and collinear-soft modes.The factorization formula for this regime was derived using SCET in refs.[40,45] and reads where we have defined convolutions between two functions f and g through

JHEP10(2019)130
Here the dots represent possible additional arguments.Furthermore, the short-hand notations The region of phase space represented by the left panel of figure 2 is reached through raising − α 1 /α 1 or lowering − α 2 /α 2 until the two are equal, joining the collinear-soft mode with the collinear mode in the process.The factorization formula of this region then no longer contains any collinear-soft functions, but instead involves jet functions depending on both angularities The consistency relation between the double-differential jet function and the convolution between the collinear-soft and single-differential jet function that this implies was verified explicitly at one-loop order in ref. [45].In this regime, soft radiation does not contribute to the measurement e α 2 and the factorization formula is simply a more differential version of the factorization formula for the sole measurement of e α 1 .Analogously, to obtain the factorization formula describing the region of phase space depicted in the right-most panel in figure 2, the soft and collinear-soft functions merge into a more differential soft function to yield This is a more differential version of the factorization theorem for the single-differential cross section in e α 2 and only resums large logarithms involving this angularity.The renormalization group equations of the perturbative functions occurring in the factorization formula in eq.(3.10) can be found in appendix B. Using the expressions for the anomalous dimensions given in eq.(B.2), the consistency relation of the factorization formula in eq.(3.10), given by is indeed found to be satisfied.
For the measurement of n angularities, the factorization formula for the cross section describing the bulk region R n (α 1 , . . ., α n ) follows from the modes appearing in the corresponding Lund plane.Specifically, there is a single soft mode, a single collinear mode (for JHEP10(2019)130 each of the two collinear directions) and there are n − 1 collinear-soft modes (per collinear direction), leading to the general factorization formula d n σ Rn(α 1 ,...,αn) dQ α 1 e α 1 . . .
When taking derivatives of this expression with respect to µ, the anomalous dimensions of all intermediate collinear-soft functions effectively combine into a single collinear-soft anomalous dimension involving the angularities e α 1 and e αn , leading to the conclusion that this factorization formula also obeys the corresponding consistency relation.Factorization formulas corresponding to regions that involve fewer angularities are again obtained by merging two degrees of freedom, i.e. merging two functions into a single, more differential function.These can involve more than two angularities, but are not required to obtain cross sections at NLL accuracy, since only tree-level expressions are needed for the functions and the anomalous dimensions smoothly merge.Specifically, the tree-level expression of each function is simply a product of delta functions of its arguments, and the more differential functions that arise due to the merging of modes thus do not give rise to any different results for the cross section.This then implies that the cross section for a specific region of interest does not depend on the total set of angularities that are measured, but instead only on the subset of angularities that occur in said region.

Resummation
The large logarithms in a factorization formula, such as eq.(3.16), can be resummed by evaluating each ingredient at its natural scale (where its logarithms are minimized) and then evolving them to a common scale.To solve the RGEs in eq.(B.1) that involve a convolution, it is convenient to switch to a conjugate space.Picking Laplace space, the transformation of a function f (t) is defined through where LT[. ..] denotes the Laplace transform and s is the variable conjugate to t.The plus distributions that appear are defined as The transformations of these distributions that are required up to NLL, are given by Solving the RGEs in Laplace space, inserting these results in the cross section in eq.(3.16) (at NLL) and transforming back to momentum space then yields the resummed cross JHEP10(2019)130 section.To ensure that our procedure allows the recovery of the inclusive cross section upon integration over all the angularities that the differential cross section depends on, we perform the resummation at the level of the cumulative cross section [67][68][69].To obtain the cumulative cross section, we integrate over each angularity e α up to some cut-off e cut α , which is the basis for our numerical implementation.Leaving this "cut" superscript implicit, the cumulative cross section is Qe Here we have defined in terms of the evolution kernels that can be found in appendix B. Furthermore, the notation has been employed to simplify the result.The natural scales of the functions at which the large logarithms are minimized depend on the (sub)set of angularities under consideration.Denoting this set of angularities by β 1 , . . ., β m , the various natural scales are given by where again we have suppressed the superscript "cut" on the angularities.To avoid the Landau pole in our numerical implementation, we freeze the value of α s below 2 GeV.

Power corrections
The power corrections to each factorization formula can be determined by considering the ratio of scales involved in the functions that are merged into a more differential function when transitioning towards a lower-dimensional region in phase space.The lowerdimensional region will be referred to as a 'daughter region' with respect to the higherdimensional 'parent region'.The measurement of three angularities will be used as an example to display this procedure.The various regions of phase space and their B-conditions, i.e. the equalities in eq.(3.9), can be found in table 1 for α 1 > α 2 > α 3 .In general, we denote the power corrections from an n-dimensional parent region R n (α 1 , . . ., α n ) towards an (n − 1)-dimensional daughter region R n−1 (α 1 , . . ., α i , α i+1 , . . ., α n ) by P n (α 1 , . . ., α n ; α i ), where the argument after the semicolon indicates the angularity that the daughter region lacks with respect to the parent region.Using this notation, the power corrections from the one-dimensional regions towards the fixed-order region are given by [45] P The power corrections of the factorization formula of a two-dimensional region towards the one-dimensional daughter regions are given by The powers denoted by # are (different) constants, which may be fixed by demanding that the power corrections should reduce to those of the one-dimensional region at the corresponding boundary, i.e.

JHEP10(2019)130
where the boundary conditions B 1 can be found in table 1. Plugging in the various scales then yields the complete set of power corrections of the two-dimensional region (3.27) The power corrections from the three-dimensional region towards any of the neighboring two-dimensional regions can be found in an analogous way.By again demanding that these power corrections should reduce to those of any other boundary theory, we obtain the set of power corrections of the three-dimensional region , e , . (3.28) This procedure is easily generalized to the case of n angularities.There are three different types of power corrections that need to be considered, all of them already present for n = 3.They are given by En(α 1 ,...,αn;αn) , P n (α 1 , . . ., α n ; α i ) = e En(α 1 ,...,αn;α 1 ) . (3.29) The powers E can then be found through demanding for 2 ≤ i ≤ n.

Matching phase space regions
The different regions of phase space that are found using Lund diagrams are each described by different cumulative cross sections.In order to obtain a combined prediction valid throughout phase space, these regions need to be matched to each another.For any given point in n-dimensional phase space, the combined cumulative cross section is defined as a linear combination of all possible regions that occur in phase space as σ(e α 1 , . . ., e αn ) = Rm a m (β 1 , . . ., β m ) σ Rm(β 1 ,...,βm) , (3.31) -14 -

JHEP10(2019)130
where again β 1 , . . ., β m is any subset of the full set of angularities α 1 , . . ., α n and the dependence of the transition variables and the cross sections on the full set of angularities has been suppressed.The sum runs over all possible regions with n ≥ m ≥ 0 and the set of coefficients is normalized as at every point in phase space spanned by the n angularities under consideration.In principle, this includes the matching to the fixed-order region, denoted by R 0 , but this matching is not performed here for simplicity, so that we simply set a 0 = 0. Following the approach in ref. [70], the specific admixture of transition variables a i is determined by the size of the power corrections to the factorization formula in each region.We start by defining a transition function that smoothly interpolates between 0 and 1 as where the constants c j are determined by demanding the continuity of f trans and its first and second derivative at both transition points x i and x f .The explicit expressions obtained in this way are given by The explicit values of the transition variables a m (β 1 , . . ., β m ) at a given point p in the n-dimensional space spanned by α 1 , . . ., αn are determined iteratively.All transition variables are initialized at 0. The region in which p lies is determined through the conditions given in eq.(3.9).If it lies outside of all regions, the transition variables are kept fixed at zero.If p lies inside a certain region R m (β 1 , . . ., β m ) involving m angularities, the following procedure is followed: • The set of daughter regions involving m − 1 angularities, obtained by removing any single angularity from R m (β 1 , . . ., β m ), is determined.
• The shortest Euclidean distance in the space spanned by α 1 , . . ., αn from the point p towards each daughter region is determined using the method of Lagrange multipliers [71].These distances are translated to a number between 0 and 1 through eq.(3.33),where the initial and final points x i and x f correspond to the distances where the power corrections are 10% and 50% respectively.The result of this procedure is denoted by ãm (β 1 , . . ., β m ; β j ), where the angularity after the semicolon again indicates the angularity that is involved in the parent region, but not in the daughter region.

JHEP10(2019)130
• The coefficient of the region R m (β 1 , . . ., β m ) is defined through . (3.36) • For each of the daughter regions R m−1 , the steps above are repeated in order to determine the transition variables a m−1 .The only notable difference is that the right-hand side of eq.(3.35) is to be multiplied by a factor 0 ≤ x ≤ 1, given by the sum of all the preliminary weights that the region under consideration might have inherited from all of its parent regions.For a region R m−1 (γ 1 , . . ., γ m−1 ), this factor is then given by This procedure is repeated until all regions from R m down to R 2 have been considered.The transition variables of the regions R 1 (β i ) are then given by the sum of the preliminary weights that they might have inherited from any of their parent regions.After all the coefficients have been determined, the cumulative distribution can be obtained through eq.(3.31).
In some cases in our numerical implementation, the cumulative distribution turns out to slightly decrease towards the fixed-order region due to the finite bin size.To ensure that this does not lead to negative spectra upon differentiation, any such bins are set equal to the average of their neighboring bins.

Results
This section contains results obtained through the reweighing procedure described in section 2. By default we show results from Herwig 7.1.4for leading order e + e − → dijets (excluding bottom and top quark jets) at center-of-mass energy Q = 1 TeV.The final-state parton shower is turned on, but the initial-state QED radiation and modeling of hadronization are switched off.The two jets are obtained via the exclusive k t algorithm [60] with the winner-take-all recombination scheme [61,62] using the FastJet package [72].We consider the set of angularities with exponent α i = 0.2 × s with s = 1, 2, . . ., 15, and use k = 4-body phase space for reweighing.We also show our analytic predictions, as well as those obtained from Pythia 8.240.We begin by showing a comparison between the Herwig, Pythia and analytic predictions for the single angularity distribution in figure 4. We find good agreement between the Herwig and Pythia results.The analytical result agrees very well with the numerical results for the angularity 2.6 , but shows some deviations for 1.2 .The reason for this is that the resummation region gets squeezed between the fixed-order region and the non-perturbative region. 6 In figure 5, we show results for two examples, obtained through reweighing with the best possible set of n = 1, 2, 3 angularities.In the left panel we show the results for α j = 1.2, where the performance for each n is comparable.In the right panel, where we show 6 While we do not include hadronization, this region is sensitive to the unphysical shower cut off.α j = 2.6, there is a dramatic improvement going from the best possible reweighing with n = 1 angularity to the best result for n = 2 reweighed angularities.We stress that the best set of angularities is obtained through a global minimization and is thus not optimized for any specific α j .The improvement from n = 2 to n = 3 is substantial, although not as dramatic.
Figure 6 shows the corresponding set of plots constructed by reweighing the NLL resummed results obtained from the calculation in section 3.In this case we have restricted the set of angularities that we examine to α i = 0.2 × s with s = 6, 7, . . ., 15 instead, and we have used the "best set" of one or two angularities obtained from the analogous procedure done with Herwig.The restriction on the angularity exponents that we consider is chosen such that it allows for a sufficient number of bins in the analytical resummation to be populated, which is otherwise not the case for lower values of α i .To focus solely on the differences that arise due to reweighing with a different number of angularities, all the distributions that enter in these plots are obtained from projecting the full three-dimensional distribution with angularity exponents {α j , 1.4, 2.8}.The reason for this is that the projection of an analytically resummed cross section involving a higher number of angularities down to a cross section involving a lower number of angularities does not exactly agree with the corresponding cross section obtained from a direct analytic calculation, i.e. without any projection.A more in-depth discussion is relegated to appendix A. With these comments in mind, we note that the reweighed results of figure 6 show a similar trend as those of figure 5, constructed using Herwig distributions.
To indicate the improvement obtained over the full domain of considered angularities, we show the goodness-of-fit, χ 2 α j ,min for the best set of n = 1, 2, 3 reweighed angularities as a function of α j in the left panel of figure 7.For n = 1, one would expect that χ 2 α j ,min goes to zero for the best single angularity.However, we take the median of 11 replicas and the optimal single angularity is not the same for each of these.On the other hand, α i = 0.2 is always part of the set of best angularities for n = 2.It is clear that the n = 2 case performs substantially better than n = 1, though there are a few angularities for which n = 2 performs worse.This happens because they are close to the best single angularity and therefore reproduced very well by n = 1, but not as well for n = 2 since the best two angularities are further away.For n = 3, there is a non-negligible, but less significant, improvement over n = 2. Finally, in the right panel of figure 7, we show the minimum goodness-of-fit χ 2 min for n = 1, . . ., 5. For n = 1, 2, 3 this is the global minimum, whereas for n = 4, 5 the results were obtained iteratively, as described in section 2. For comparison we include n = 0, which simply states how well pure phase-space predictions (without any reweighing) describe the angularity distributions in Herwig, thus providing a baseline.The black error bars correspond to roughly one standard deviation, having been constructed from the spread of the 7 most central replicas out of a produced total of 11 replicas.The dots represent the median of the replicas.One can again observe a substantial improvement from n = 1 to n = 2 and a smaller but still visible improvement from n = 2 to n = 3.The degree of improvement for going to n = 4 or n = 5 reweighed angularities is much smaller.
In appendix A.1 we provide plots that demonstrate the robustness of our results under different variations.These include the use of Pythia as the Monte Carlo for the reweighing, considering either 5-or 6-body phase space, restricting the values of angularity exponents, and considering lower or higher center-of-mass energies.We also show results that discuss the quality of the projections from higher-dimensional distributions and the impact of these on the reweighing procedure in appendix A.2.

JHEP10(2019)130
modes, and estimating the power corrections to connect them becomes more complicated for three (or more) angularities.This can be extended to processes with jets in hadronic collisions, in which case gluon jets also enter, and non-global logarithms [73] arise from soft radiation that simultaneously contributes to the angularities (measured on the jet) and the out-of-jet region.
Taking distributions obtained from this analytical resummation, as well as from the Herwig and Pythia Monte Carlo parton showers, we have studied whether employing them to reweigh a flat phase-space generator leads to improved predictions for other angularities (not used as input) via kinematic correlations.We have found an order of magnitude improvement when reweighing by distributions of two angularities over using only one, demonstrating the benefit of joint resummation.Reweighing with three or more angularities provides further improvement, albeit with a diminishing effect.The robustness of our conclusions is demonstrated by varying parts of our setup.
Our study shows that reweighing leads to improved predictions, particularly if the observable used in the reweighing procedure is similar to the observable of interest.Augmenting Monte Carlo parton showers by analytic resummation at NLL is probably not that useful, due to the sizable perturbative uncertainty at this order.However, this could be improved by matching the NLL to a fixed-order calculation.Furthermore, the factorization formulae presented here are not limited to a specific resummation order, and in principle all ingredients needed to NNLL are the same as for two angularities in ref. [45].The reason for this is that, apart from the anomalous dimensions, all ingredients in the factorization formulae are only needed at one-loop order.They can therefore depend on at most two independent variables, so any additional angularities can be expressed in those two.
We believe that this approach opens up a new route towards precise and detailed (i.e.differential) predictions for collisions at the LHC, supplementing current advances in Monte Carlo parton showers.Such predictions are particularly important in an era in which the Standard Model is subjected to ever more stringent tests, and Machine Learning techniques are developed in order to uncover faint signals through detailed features in the data.

A Additional plots
In this appendix we provide some additional plots that strengthen the universality of our conclusions, and briefly discuss the projections of the analytical resummed results.

A.1 Robustness of conclusions
First of all, we show the analogue of figure 7 using Pythia instead of Herwig.Specifically, the left panel of figure 8 shows the goodness-of-fit, χ 2 α i ,min , as a function of α i .In this case, there is a single best angularity at α = 0.8 for each of the replicas.Also, there is now a small region that does not improve from n = 2 to n = 3.The right panel shows the global minimum goodness-of-fit χ 2 min for n = 1, . . ., 5. The qualitative behavior is very similar to that in the right panel of figure 7 for Herwig, but the values of χ 2 α i ,min are somewhat larger.Figure 9 is the analogue of figure 7, but using 5-body flat phase space (i.e.k = 5) instead of 4-body phase space.While the qualitative behavior is largely the same, there are some numerical differences that are mostly driven by statistical fluctuations, most visible in the left panel.Specifically, the sampling of the collinear and soft regions of phase space relevant for angularities is worse for k = 5.The left panel of figure 10 is the analogue of the right panel of figure 7, showing the minimum goodness-of-fit χ 2 min for 6-body phase space.As with k = 5, the larger number of phase-space particles worsens the sampling of the angularity phase space and hence increases the statistical fluctuations, reflected by the larger uncertainties.The right panel shows χ 2 min where we repeated the analysis for k = 4, restricting to the angularity exponents in the interval α i ∈ [1.2, 3.0] (again in steps of 0.2).Because the total number of angularities is smaller, we only show the result of reweighing with n = 1, 2, 3 angularities.This is also the reason for the faster convergence as function of the number n of angularities used in the reweighing.
Finally, figure 11 is the analogue of the right panel of figure 7, but for different centerof-mass energies, Q = 200 GeV (left panel) and Q = 4 TeV (right panel).We do not show the reweighing with n = 4, 5 angularities.The qualitative behavior is similar as for Q = 1 TeV, but suggests that the reweighing procedure performs better for lower energies.This is in line with what one would expect from the increase in jet entropy with Q [74].gives the difference between the one-dimensional resummed angularity distribution and each of these projections.

A.2 Projections
In figure 6 we showed results obtained by applying the reweighing procedure to analytical resummed results, all derived from the same multi-dimensional calculation.Here we comment briefly on the issues leading to this approach, and show results for an alternative method.
Figure 12 shows the projections from higher-dimensional analytical results down to two specific angularity choices α j = 1.2 and α j = 2.6.The yellow dotted and green dashed curves represent the projection from the indicated three-dimensional and two-dimensional cross section respectively.The red curve shows the directly determined cross section differential in the single angularity.The figure clearly illustrates that the curves are not identical and that projections differ more in comparison to the one-dimensional distribution if one starts from higher-dimensional cross sections.This fact is quantified in the χ 2 between the one-dimensional resummed angularity distribution and each of the projections.Note that there is no corresponding issue for the Herwig and Pythia results, since they originate from fully exclusive events.
We suspect that this discrepancy is largely due to binning issues.As described in section 3.1, the number of distinct kinematic regions in phase space increases dramatically when cross sections differential in more angularities are considered.The result of the differential cross section in (the center of) each bin is obtained by determining the cumulative cross section on the edges of the bin7 and taking a numerical derivative.Due to the relatively small number of bins and the increasing number of kinematic regions, situations in which the edges and the center of a bin lie in different kinematic regions might occur.In these cases, the prediction of the spectrum (at the center of the bin) is obtained from input provided by cumulative distributions obtained from factorization formulas that are not valid at that point.As this is a binning issue, we expect the effect to diminish when a larger number of bins is considered. .Same as figure 6, but with the reweighing performed using the multi-dimensional analytical distributions calculated directly (rather than projections from a single multi-dimensional distribution).

JHEP10(2019)130
For completeness, we show in figure 13 the reweighing performed for the same restricted set of angularities as in figure 6, but using the n-dimensional analytical distributions calculated directly, and comparing to the one-dimensional analytic prediction.

JHEP10(2019)130
The MS cusp anomalous dimension to two loops is given by [76] Γ cusp (α s ) = α s 4π 4C F + α s 4π To NLL accuracy, the evolution kernels occurring in the resummed cumulative cross section in eq.(3.20) are given by where r = α s (µ)/α s (µ 0 ) has been defined.The required one-and two-loop coefficients of the beta function in the MS scheme are given by [77] Open Access.This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

Figure 1 .
Figure1.Illustration of the Lund plane with x = log 10 (1/θ) and y = log 10 (1/z).The crosses represent emissions and the red line the measurement of an angularity α , set by the dominant emission.There can be no emissions in the shaded area, while emissions above this only contribute at higher order.The green and orange dot denote the collinear and soft mode respectively.

Figure 2 .
Figure 2. The measurement of two angularities α1 and α2 represented in the Lund plane.Each panel describes a distinct region of phase space.The left and right panels involve only collinear (green) and soft (orange) modes, while the center panel contains an additional collinear-soft (blue) mode.Emissions in the shaded region are vetoed.

Figure 3 .
Figure 3.The Lund plane describing the region of phase space for the measurement of three angularities in which the most logarithms may be resummed.The various modes are denoted by the green (collinear), orange (soft) and blue (collinear-soft) dots.The dotted line serves to indicate the point C α1α2 α3 ) and a preliminary weight b m−1 (γ 1 , . . ., γ m−1 ; β i ) is assigned to each of the m daughter regions R m−1 (γ 1 , . . ., γ m−1 ).Here the β i after the semicolon in this case indicates the angularity that should be added to the set {γ 1 , . . ., γ m−1 } ⊂ {β 1 , . . ., β m } to obtain the full set of angularities {β 1 , . . ., β m } on which the parent region depends.The preliminary weights are given by b m−1 (γ 1 , . . ., γ m−1 ;

2 eFigure 4 . 2 =k = 4 HERWIG, αj = 2 . 6 PhSpFigure 5 .
Figure 4.The Herwig (red), Pythia (green dotted) and analytical NLL (yellow dashed) predictions for the 1.2 and 2.6 distributions.The fixed-order region has been grayed out and the analytical results have been normalised to the fraction of the area of the Herwig results that lies to the left of that region.

= 1 Figure 7 .
Figure 7. Left panel: we show the goodness-of-fit χ 2 αj from eq. (2.3) for the best set of n = 1 (red), n = 2 (green dotted) and n = 3 (yellow dashed) reweighed angularities as a function of the angularity exponent α j for Herwig.For each point we take the median value of the 11 replicas.Right panel: the global minimum χ 2 min from eq. (2.4) for n = 1, . . ., 5. The error bars represent the uncertainty as described in the text, with the central black dot representing the median over the replicas.The initial value of χ 2 (blue star) shows how well flat phase space reproduces the angularities prior to reweighing.

e= 1 Figure 10 .
Figure 10.Same as the right panel of figure 7, but obtained from k = 6-body phase space (left panel) or by restricting the domain of angularity exponents to the set α i ∈ [1.2, 3.0] (right panel).

e= 4 Figure 11 .
Figure 11.Same as the right panel of figure 7 but for Q = 200 GeV (left panel) or Q = 4 TeV (right panel).In this case, n = 4, 5 are not shown.

2 e 6 eFigure 13
Figure 13.Same as figure6, but with the reweighing performed using the multi-dimensional analytical distributions calculated directly (rather than projections from a single multi-dimensional distribution).