Phenomenology of jet angularities at the LHC

We compute resummed and matched predictions for jet angularities in hadronic dijet and Z+jet events with and without grooming the candidate jets using the SoftDrop technique. Our theoretical predictions also account for non-perturbative corrections from the underlying event and hadronisation through parton-to-hadron level transfer matrices extracted from dedicated Monte Carlo simulations with SHERPA. Thanks to this approach we can account for non-perturbative migration effects in both the angularities and the jet transverse momentum. We compare our predictions against recent measurements from the CMS experiment. This allows us to test the description of quark- and gluon-jet enriched phase-space regions separately. We supplement our study with SHERPA results based on the matching of NLO QCD matrix elements with the parton shower. Both theoretical predictions offer a good description of the data, within the experimental and theoretical uncertainties. The latter are however sizeable, motivating higher-accuracy calculations.


Introduction
The LHC data-taking campaigns during run I and II recorded huge amounts of high-quality data which can be explored to perform high-precision tests of the Standard Model. Recent studies have demonstrated the particular potential of jet cross-section measurements. For example, jet data can be used in determinations of the strong coupling constant, see e.g. [1][2][3][4][5], and form an important input for constraining parton distribution functions (PDFs) [6][7][8][9][10][11][12][13][14][15]. Information on the physical processes underlying jet formation can be gained by studying so-called jet substructure observables (for recent reviews see e.g. [16][17][18]). In this context, a popular jet substructure technique used to remove radiation in phase-space regions dominated by non-perturbative physics is SoftDrop grooming [19], including the modified Mass Drop Tagger (mMDT) as a special case [20,21]. Physics observables measured on SoftDrop jets have been shown to be less affected by non-perturbative corrections and to be amenable to precision calculations in QCD [22][23][24][25][26][27]. This has triggered the attention of the jet substructure community, both theorists and experimentalists [28][29][30], over the past few years.
One of the central goals of jet substructure analyses is to determine differences between jets originating from quarks versus gluons. This information can be exploited, for instance, to enhance sensitivity in searches for new physics. Consequently, the consistent modelling of quark and gluon jets has been an active area of research over the past years, see for example [31][32][33][34][35][36][37][38][39]. An often-studied benchmark, widely used in many of the above references, is the so-called generalised jet angularities [40,41]. A sound theoretical understanding of these angularities is thus critical to the task of quark-gluon discrimination. These observables, and closely related ones, have hence been studied in a variety of theoretical frameworks [40,[42][43][44][45][46][47][48][49][50][51][52]. Recent measurements of jet angularities at the LHC have for example been presented in [53][54][55]30].
In this paper we focus on the wealth of data presented by the CMS experiment in [30], which measured angularities, with and without SoftDrop grooming, on jets produced either in association with a Z boson or in dijet events. In an earlier article, [51], we already presented NLO + NLL accurate predictions for the Z+jet case, derived using the resummation plugin to the SHERPA event generator framework [56,57].
Here, we aim to extend this study to the case of dijet production. In this context, we introduce an improved prescription to apply non-perturbative corrections, extracted from Monte Carlo (MC) simulations, to our resummed analytic predictions. Since Z+jet and dijet events receive different contributions from quark-and gluon-initiated jets, a combined study of (groomed) jet shapes such as the angularities in both channels is expected to shed light on our fundamental understanding of jet formation from hard quarks and gluons. To that aim, besides providing results for the angularity differential distributions, we also consider the angularities mean values in quark-and gluon-dominated phase-space regions of Z+jet and dijet final states.
The article is organised as follows: in Sec. 2 we introduce the considered angularity observables and define the fiducial phase space used in our study of jets in Z+jet and dijet production. In Sec. 3 we introduce the theoretical methods used for making predictions at NLO + NLL accuracy, based on the implementation of the CAESAR formalism in SHERPA. Here, we also present our new approach to account for non-perturbative corrections based on transfer matrices extracted from MC simulations, see Sec. 3.3. We present our final NLO + NLL + NP results in Sec. 4, alongside with MC predictions from SHERPA, and compare those against data from the CMS experiment. We compile our conclusions and give an outlook in Sec. 5.

Phase space and observable definition
To be able to directly compare to the CMS measurement presented in [30], based on proton-proton collisions at a centre-of-mass energy of √ s = 13 TeV, the exact definition of angularities and the eventselection criteria have to be carefully matched to the experimental setup. For completeness, this section summarises the final definitions and cuts and we refer the reader to Ref. [30] for additional details and motivations.
The events considered are conceptually triggered either by an approximately on-shell Z boson, decaying into a pair of muons, accompanied by a hard jet (Z+jet case), or by two high-p T jets (dijet case). The precise phase space for the two cases is defined by the selection criteria summarised in table 1. Substructure observables are then calculated for the jet with the highest p T in the Z+jet case. For dijet events, the two jets with the largest transverse momenta are considered and the measurement is performed separately for the more central (smaller rapidity |y|) and more forward (larger |y|) of those jets. For the Z+jet process p T,jet bins in the interval [50,1500] GeV were considered, while in the dijet case this range was extended to [50,4000] GeV.
The substructure observables we are interested in belong to the family of jet angularities. These probe both the angular and the transverse momentum distribution of particles within a given jet. They are defined from the momenta of jet constituents as follows: (2.1) Here is the azimuth-rapidity distance of particle i from the jet axis, which is determined using the Winner-Take-All (WTA) recombination scheme [46] for α ≤ 1, and using the default E-scheme axis otherwise. R denotes the radius parameter used for the anti-k t [58] jet clustering. We restrict our analysis to the infrared-and-collinear-safe case κ = 1 with α = 0.5, 1, 2, while Ref. [30] in addition considered the infrared-unsafe combinations κ = 0, α = 0 and κ = 2, α = 0. For the different choices of α we follow the common naming convention [31,59,30], and refer to λ 1 0.5 as Les Houches Angularity (LHA), λ 1 1 as jet width, and λ 1 2 as jet thrust. The CMS collaboration also measured the same jet angularities based on the jet constituents obtained after SoftDrop grooming, as implemented in the FASTJET [60] contrib package. The SoftDrop parameters are chosen as β = 0 and z cut = 0.1. All results are provided either using all the particles in the jet (both Process Selection cuts

Dijet
Require at least two R = 0.4 (or R = 0.8) anti-k t jets with p T > 30 GeV. The two leading (highest p T ) jets have to satisfy |y jet | < 1.7, p T,jet > 30 GeV, |φ j1 − φ j2 | > 2 and |p T,j1 − p T,j2 | / |p T,j1 + p T,j2 | < 0.3 . Table 1: Event-selection cuts used by the CMS collaboration in [30]. The muons used to reconstruct the Z boson are considered as "bare", i.e. without QED emissions, and are excluded from the input to the jet clustering. charged and neutral), or using only the charged jet constituents. Note that the p T,jet bin is always determined from the full jet, without any grooming or restriction on the charge of particles.

Theoretical framework for jet-angularity predictions
We here consider two sets of theoretical predictions for jet angularities in Z+jet and dijet production to confront with the CMS data. The first are simulations with the SHERPA multi-purpose event generator. The jet-angularity distributions presented in [30] were compared against several event-generator predictions. The highest accuracy was reached in two simulations based on LO merging [61,62] of tree level matrix elements with additional jets dressed with parton showers, one using MADGRAPH5 aMC@NLO [63] in conjunction with PYTHIA [64], the other one using SHERPA [65,57]. In addition, results obtained from the parton showers of PYTHIA, HERWIG [66,67] and SHERPA, based on Born matrix elements, were presented. We here provide SHERPA predictions at NLO accuracy for the Z+jet and jj processes, which were not covered in [30].
The second set of predictions is obtained from NLO + NLL semi-analytic calculations, based on the implementation of the CAESAR resummation formalism [43] in the SHERPA framework. In [30] the CMS results for Z+jet production were already compared against our NLO + NLL predictions from [51]. We here extend this comparison to the case of dijet production. In addition, we introduce a more physical way to apply non-perturbative corrections to our calculations, that we refer to as transfer-matrix approach which we also apply to both the Z+jet and dijet cases.
In the following, we describe our setup for the MC simulations with SHERPA, then introduce the necessary additions to our resummation framework to carry out the NLO + NLL calculation for the dijet case, and finally introduce the transfer-matrix approach.

Monte Carlo simulations with SHERPA
We compile hadron-level predictions for jet angularities with the SHERPA [65,57] event generator, version 2.2.11, using the NNPDF-3.0 NNLO PDF set [68]. We consider inclusive dijet and µ + µ − +jet production based on the SHERPA implementation of the MC@NLO formalism [69], dubbed SH-MC@NLO in what follows, matching the NLO QCD matrix elements for two-jet and µ + µ − +jet production with the SHERPA Catani-Seymour dipole shower [70]. The involved QCD one-loop amplitudes we obtain from OPENLOOPS [71], using the COLLIER library [72] for the evaluation of tensor and scalar integrals. The central values for the perturbative scales entering the calculation are set to To estimate the perturbative uncertainties of our predictions, we perform on-the-fly [73] 7-point variations [74] of the factorisation (µ F ) and renormalisation (µ R ) scales in the matrix elements and the parton shower, while we keep the parton-shower starting scale (µ Q ) fixed as it is related to the dipole kinematics of the parton shower. As error estimate we then consider the envelope of the results for The SHERPA Underlying Event (UE) simulation is based on its implementation of the Sjöstrand-Zijl multiple-parton interaction model [75]. To account for the parton-to-hadron transition we use the SHERPA cluster fragmentation [76]. Later, when we will use SHERPA to account for non-perturbative corrections to the analytic calculation, we include an estimate of the uncertainty related to the hadronisation model. This uses the interface to PYTHIA 6.4 [77], provided with SHERPA, with the Lund string fragmentation model [78,79] as an alternative. In both models the default set of tuning parameters is used, see [57] for details. However, we additionally vary the default for the rescale-exponent parameter α MPI min used to model the dependence of the cutoff regulator p T,min for the 2 → 2 cross-section integration in the UE model on the hadronic centre-of-mass energy, i.e. 244. This is motivated by a preliminary study comparing SHERPA predictions for the jet angularities in particular for the low-p T,jet region against the data from [30]. As systematic up-and down-variations of the MPI activity we consider α MPI min,up = 0.08 and α MPI min,down = 0.24.

NLO + NLL resummation for dijet final states
Throughout this paper, we target a description of the angularity distribution at the NLO+NLL accuracy, whereby the exact distribution at order α 2 s (NLO) * is matched to the all-order resummation, which is relevant for small values of the observable at hand, carried out at the NLL accuracy.
The NLL resummation is performed using the SHERPA implementation of the CAESAR formalism [43]. This plugin to SHERPA was first presented in [56] and further developed in subsequent work, e.g. to introduce new matching schemes [80], to include the SoftDrop phase-space constraints [81], or to support some non-global configurations in Ref. [51]. The master formula for the all-order resummation can be written as the following sum over partonic channels δ: In this expression, dσ δ / dB is the Born-level cross section, obtained through the COMIX matrix element generator [82]; the collinear radiators R l are directly obtained from the CAESAR paper [43] or extended to include the effect of SoftDrop (see e.g. Refs. [19,81]); S is the soft function including global and nonglobal colour evolution; P includes corrections associated with PDF evolution for initial-state radiation (which are absent in our case, i.e. P = 1); F is the multiple-emission function, which for our application * Note, in this counting the powers of αs from the underlying Born process, i.e. α 2 s for dijet and α 1 s for Z+jet, are not included. is simply given by e −γ E R /Γ(1 + R ) with R = ∂ L R since angularities are additive observables. Finally, Θ hard implements the cuts on the Born-level phase space, cf. table 1. The phase-space integration as well as the separation over different partonic channels is greatly facilitated by the integration of the plugin into the SHERPA framework. With this approach, one reaches NLL accuracy including an estimate of the resummation scale uncertainty (introducing a rescaling parameter x L ), as well as a treatment of the endpoint corrections and associated uncertainties. Modulo an enhanced treatment of non-global logarithms to account for the more complicated colour structure of dijet events that we discuss below, our usage of the resummation plugin to SHERPA is the same as in our study of angularities in Z+jet events and we therefore refer the reader to Ref. [51] for details. In particular, through the iterative flavour-k t clustering procedure [83] (referred to as BSZ algorithm in what follows), our automated resummation plugin allows us to keep track of the jet flavour (see again [51] for a detailed description of the algorithm), yielding the so-called NLL resummation accuracy.
Let us now discuss our treatment of soft logarithms and, in particular, of non-global logarithms. In practice, the soft function S can be factorised in a global contribution S global and a non-global [84,85] part S non-global . Although we perform the resummation of non-global logarithms in the large-N C limit, we can account, to all orders, for the full set of colour correlations induced by the first emission. In order to achieve this we note that, for a given Born-level partonic channel B δ , the overall contribution from soft logarithms can be written as where α s = α s (µ Q ). Note that these contributions arise only at NLL accuracy and so the actual scale choice for α s is subleading. In practice, we evaluate α s at µ R . Here, H and c are the hard function and colour metric for this channel. The I ij t coefficients in tΓ B δ correspond to global contributions, stemming from soft-wide-angle radiation. These depend only on the jet radius R and are computed in the small-R limit. They are taken from [86] and we note that, compared to the Z+jet case, the dijet case also has to include dipoles involving the recoiling jet. The non-global contributions to S are encoded in the functions f ij (t) in tΓ B δ . They are computed independently for each dipole using a Monte Carlo approach as described in [84]. † The main difference compared to our previous work on Z+jet events is the fact that one also has to include the f ij (t) corresponding to dipole configurations involving the recoiling jet. On top of the standard dependence on t(L) and on the jet radius, these dipoles involving the recoiling jet also depend on the rapidity separation between the measured jet and the recoiling one (see e.g. [86]). This approach allows us to compute the nonglobal soft function for ungroomed angularities. For SoftDrop groomed distributions, the non-global factor remains the same for v ≥ z cut and saturates at that value, i.e. S non-global (max(v, z cut )) . Finally, we note that while for Z+jet processes the matrix structure of Eq. (3.8) disappears, the situation for dijet events is more complex and the full matrix structure needs to be taken into account.
The above-described resummed observable calculation is matched to exact distributions at NLO. These are obtained using SHERPA with RECOLA [87,88] and OPENLOOPS [71] as one-loop generators in conjunction with the COLLIER library [72], and using Catani-Seymour dipole subtraction [89,90]. All tree-level contributions are obtained from the built-in COMIX [82] generator. For a generic angularity λ, the matching is done at the level of cumulative distributions Σ δ (λ) defined for each flavour channel δ. For this, we expand the (fixed-order) jet cross section σ fo (p T ) as well as the fixed-order and resummed cumulative distributions, Σ fo (λ) and Σ res (λ) in series of α s : fo + . . . , (3.11) Σ δ res (λ) = σ δ(0) + Σ δ(1) res + Σ δ(2) res + . . . , (3.12) where the p T dependence is left implicit, the number in parentheses indicates the respective order in α s (without including the Born-level factors of α s and, possibly, α EW ), and σ δ(0) is the Born-level jet cross section in the specified flavour channel. Our final matched expression for the cumulative distribution reads: (2) can be calculated without the need for the two-loop virtual corrections to the jet cross section.
To estimate the perturbative uncertainties of our predictions, we again consider 7-point variations of the factorisation and renormalisation scales in the matrix elements and the resummation, cf. Eq. (3.3).
The central values in the Z+jet case are taken from [51], i.e. the transverse momentum of the Z boson, multiplied by the jet radius in the case of µ Q . Similarly, in the dijet case we use µ R = µ F = µ Q /R = H T /2, closely matching the scale choice in the SH-MC@NLO simulations see Eq. (3.2). When implementing the matching between our fixed-order and resummed results, we treat the endpoint of the resummed expression and of its expansion by setting with v max the kinematic endpoint of the fixed-order distribution. Varying x L between 1 2 and 2 then allows us to quantify the resummation uncertainty (we keep p = 1 fixed). As total perturbative uncertainty estimate we consider the envelope of all (µ R , µ F ) and x L variations.

Non-perturbative corrections in the transfer-matrix approach
So far we have explained our analytic description of jet angularities in perturbative QCD. Clearly, these observables probe QCD dynamics in the infrared regime and are therefore also sensitive to nonperturbative (NP) corrections. At a hadron collider such as the LHC, this includes the Underlying Event i.e. remnant-remnant interactions -on top of the parton-to-hadron transition. Because jet angularities are IRC-safe observables, NP effects can be considered as power corrections and are therefore expected to become less important as one increases the transverse momentum of the jet they are measured on. However, at moderate p T,jet their effect can be sizeable and it is crucial to include them in any realistic phenomenological prediction. Additionally, the size of the jet-radius parameter impacts the susceptibility to NP corrections.
In previous works [22,23,91,51], we have adopted a rather simple procedure to account for NP effects. We have computed predictions for a given jet-angularity distribution in the considered fiducial phase space, i.e. applying the phase-space constraints on jets and restricting to a given p T,jet bin, both at parton and hadron level, using a general-purpose MC generator. The ratio between these two distributions was then used to correct our resummed and matched predictions for that particular observable and phase-space region. While this model captures the dominant NP effects of redistributing the partonic jet-constituent momenta onto hadrons, it is still rather crude, as it does not fully account for alterations of the underlying parton-level event kinematics due to NP effects. In particular, both hadronisation and UE affect the transverse-momentum distribution of jets, potentially leading to migration between the considered p T,jet -bins. By neglecting this important feature of NP corrections, this approach turns out to be rather sensitive to cutoff-effects in the employed parton shower. In particular, NP corrections can become pathological for observable values corresponding to scales below the shower cutoff. This generator dependence can be partially overcome by considering several MC programs and using the spread of the extracted NP ratios as an uncertainty on the NP corrections, as was done in Ref. [51], and consequently used in the comparison to CMS data in [30].

The transfer-matrix approach
In this new study, we attempt to overcome the limitations of the method for extracting NP corrections we previously used. To this end we develop and implement a more realistic and detailed model to include NP corrections, in which we account for the change of parton-level event kinematics due to hadronisation and the UE. In particular, considering double-differential measurements in jet-transverse momentum and the angularity observables, we derive non-perturbative transfer matrices which account for the alteration of both p T,jet and the angularity. Through the migration in transverse momentum we significantly reduce the sensitivity to phase-space restrictions and non-perturbative parameters in the underlying MC simulations, in particular to the parton-shower cutoff parameter. We here present our approach in full generality, applicable to an arbitrary set of observables, measured in a multi-differential way. We thereby aim to keep track both of the migration in the underlying event-kinematical variables, used to define the fiducial phase space, that get partially integrated over, as well as changes in the actual observable of interest, e.g. a specific angularity variable.
Let us consider a scattering process which results in a partonic configuration P. Through NP effects the set of parton momenta is then mapped onto a hadron-level configuration H (P). The map H, which does not need to be fully specified at this point, accounts for hadronisation and UE corrections. It could be derived from field-theoretical considerations (see for example Refs. [92,93] for recent work on SoftDrop observables) or it can be extracted from any given parton-shower simulation interfaced to a model of NP phenomena. For a given configuration P or H (P), we then measure a set of m observables, V (P) or V (H (P)). ‡ We define the transfer operator as the conditional probability to measure a hadron-level set of observables v h , evaluated on H (P), given that the parton-level observables were v p : .

(3.15)
This way, the multi-differential distribution for the set of hadron-level observables v h can be written as When performing numerical studies, we often work with binned distributions, i.e. we consider binned cross sections obtained by integrating the multi-differential distribution over hypercubes in the observables' space. If we consider, for instance, the parton-level case, the cross section in any given hyper-bin p is written as (3.18) ‡ Here, we have chosen the same set of observables V on the parton-and hadron-level configurations for simplicity. It is of course trivial to extend this to differing sets of observables for parton and hadron level, for example using additional auxiliary observables to parameterise the parton-level phase space, or adding selection criteria like, for instance, particle charge, at hadron level.
If we now consider a binned distribution at hadron level, the transfer operator from parton-level bin p to a given hadron-level bin h becomes a matrix of the form Consequently, the final hadron-level distribution in the hyper-bin h is obtained by the weighted sum of all parton-level contributions We note that the Θ-functions introduced in Eqs. (3.18) and (3.20) can also be used to represent eventselection cuts, such as the ones reported in Tab. 1, that define the fiducial region of a measurement, at parton-and hadron level, respectively. We further note that the parton-and hadron-level bins do not necessarily have to be the same. For example, one would typically define underflow and overflow bins at parton level so as to include their contribution to the final hadron-level predictions.
The elements of the transfer matrices appearing in Eq. (3.21) can easily be extracted from a multipurpose generator in a single run, given that events are accessible at different stages of the simulation process. Note however that, while the parton shower and hadronisation are treated in a factorised form in all multi-purpose event generators [94], this is not necessarily the case for the UE. In particular PYTHIA [77,64] makes use of an interleaved evolution of the initial-state shower and the secondary interactions [95,96]. Accordingly, in a full event simulation within such model there is no notion of an intermediate parton-level final state that is directly comparable to a resummed calculation. However, in SHERPA the parton showers off the hard process and the simulation of multiple-parton interactions are fully separated, i.e. the UE is simulated only after the shower evolution of the hard interaction is completed. The secondary scatterings then get showered and ultimately the partonic final state consisting of the showered hard process and multiple-parton interactions gets hadronised. In what follows we therefore determine the transfer matrices using the SHERPA generator.
Of course, the matrix T hp depends on the model or the specific generator used to derive them. However, we argue that a milder dependence than with the naive parton-to-hadron level ratios we previously adopted can be expected. As the matrices are based on the conditional probabilities for transitions between bins, e.g. in p T,jet , they only depend on the parton shower insofar as it determines what exact phase space those probabilities are averaged over.
It is interesting to note that, the task of correcting PL predictions for NP effects closely resembles the situation of accounting for detector effects in an experimental analysis. For this purpose proposals have recently been made (see for example Refs. [97][98][99][100][101]) to use machine-learning techniques to accomplish such a mapping for either binned or unbinned distributions. One can envisage that similar methods can be applied for the problem of non-perturbative corrections, possibly even in an invertible way, relating the observed hadronic final state to an underlying parton-level prediction. §

Transfer matrices for jet-angularity observables
For our study of jet angularities we follow the above strategy to determine T hp from MC simulations with SHERPA, and use them to correct our binned perturbative NLO + NLL predictions. To this end we employ the RIVET analysis tool [102,103], which allows us to directly access the HepMC [104, 105] event record and extract the intermediate truth-level information on the event evolution after parton showering, prior to UE and hadronisation. This enables us to consistently derive the desired migration matrices for all angularity observables (with and without grooming, both jet radii, and based on all/charged hadrons) for a particular UE parameter set and hadronisation model in a single generator run for Z+jet and dijet production, respectively. To assess the systematic uncertainties related to the UE activity and hadronisation, we derive alternative transfer matrices for the up and down variation of the α MPI min parameter and using the string fragmentation model, cf. Sec. 3.1. In total, we consider the envelope of (model, α MPI min ) = {(cluster, 0.16), (string, 0.16), (cluster, 0.08), (cluster, 0.24)} , (3.22) using cluster hadronisation with α MPI min = 0.16 as our default. Ideally, the uncertainties on the parton-to-hadron transfer matrix should include a contribution associated with the choice of the parton-shower cutoff scale, especially if one then wants to use the transfer matrix with analytic parton-level results which do not have an explicit cutoff. Varying the parton-shower cutoff would however require a full re-tuning of the non-perturbative parameters of the Monte Carlo, a task which is clearly beyond the scope of the present study. In practice, our use of two hadronisation models, cluster and string fragmentation, should at least partially include the effect of varying the parton-shower stopping scale, at a much lower cost.
In order to illustrate our transfer-matrix method, we discuss some concrete examples. First, we consider the change in transverse momentum for the leading jet in Z+jet production. We apply the transfer matrix extracted from the SH-MC@NLO simulation for the Z+jet event selection for ungroomed jets, based on all hadrons, to our parton level (PL) NLO+NLL resummed calculation. Thus, we consider v p and v h in Eq. }. We first focus on the jet-p T distribution. Integrating out the angularity at hadron level, our transfermatrix approach yields the following result: dσ HL dp HL where we consider binned distributions, and multiply our result for λ 1 2 at NLO + NLL accuracy by the matrix form of T ({p HL T,jet }|{p PL T,jet , λ 1,PL 2 }). We note that Eq. (3.23) still non-trivially depends on the assumed double-differential distribution of σ PL , even though in itself it is only differential in one observable. This happens because the above expression uses different parton-level distributions for the determination of the transfer matrix T , here SH-MC@NLO simulations, and for the parton-level cross section σ PL , given by the NLO + NLL prediction.
It should be noted that there is an additional contribution from events where the parton level does not pass some of the event-selection cuts, like the lowest transverse momentum bin edge but also any other cut in principle, whereas the corresponding hadron-level event actually does pass the cuts. We have found that this contribution is basically negligible starting from the p T,jet ∈ [88 − 120] GeV bin. For the Z+jet case, we added an underflow bin, dependent on a technically necessary cut p T,jet > 30 GeV on matrix-element level, and confirmed that this contribution is indeed dominated by lower p T,jet regions. In the following, we will consider distributions starting from p T,jet > 88 GeV and neglect this additional contribution.
The result is illustrated in Fig. 1, where we show the set of transverse-momentum bins measured by CMS and show the effect of bin migration from parton-to hadron level on the respective NLO+NLL cross section for anti-k t jets with R = 0.4. For example, the hadron-level cross section for p T,jet ∈ [65, 88] GeV receives non-negligible contributions from parton-level events with smaller transverse momentum, e.g. through an additive contribution from the UE, as well as higher slices, e.g. through the widening of jets  due to hadronisation. It is evident that these effects, which were completely ignored in our previous study, can be rather sizeable, especially at low and moderate p T,jet , and are therefore important to consider in the derivation of NP corrections. Note that we do not illustrate the flow in or out of the underflow bin at p T,jet < 50 GeV for simplicity. The net flow over this border appears to be negligible. In principle, a similar study could be performed for all the other variables that are used to characterise the event kinematics and to define the fiducial region, for example the jet rapidity. However, we expect that their impact is much less pronounced and, therefore, in order to keep the dimensionality of the transfer matrices as low as possible, we have ignored them.
As second step in our assessment of NP corrections, we study their effect on angularity distributions. While at parton level the angularities are measured on all the partons in the jet, at hadron level they can either be measured on all the particles or using charged particles only. To illustrate the effect of NP corrections, we consider in Figs. 2 and 3 the LHA, λ 1 0.5 , measured on ungroomed jets, using different approximations for the p T,jet migration from parton to hadron level: including only the contribution of the transfer matrix for the same p T,jet bin (dashed orange curve), including as well migration from the two neighbouring bins (dashed red curve), and including the full transfer matrix (solid black curve). Each time, the same approximation is used for both the numerator and the denominator of the normalised distribution. For reference, the parton-level result is shown by the solid blue line. Fig. 2 concentrates on a medium-p T,jet bin, p T,jet ∈ [120, 150] GeV, for both the Z+jet and the central dijet selection and Fig. 3 shows the corresponding results for the highest p T,jet interval, p T,jet ∈ [1000, 4000] GeV.
In Fig. 2, we see that including only the same p T,jet bin already accounts for the bulk of the nonperturbative correction, and that NP effects cause sizeable migrations between the different angularity bins. However, we see a non-negligible contribution of events from the neighbouring p T,jet bins. The effect of p T,jet migration beyond the two closest bins differs from the full correction only at large values of the angularity. For the final prediction, using the full transfer matrix, we also derive an uncertainty estimate, shown as a grey band, corresponding to the variation of the UE α MPI min parameter as described above. We observe that, despite the large overall effect of the NP corrections, they are fairly stable under the considered variations. The observed uncertainty bands do not fully cover up the case where we discard the effect of p T,jet migration from the neighbouring bins, in particular for large λ 1 0.5 values. However, as one might expect, NP effects decrease when considering the higher p T,jet bin in Fig. 3. Despite being also much smaller at high p T,jet , the migration between different observable bins is still by no means negligible. The effect of migration of events from lower p T,jet regions does however appear to be very small in this case, and is mostly covered by the UE variations.
To better illustrate the effect of hadronisation in our transfer-matrix approach, we directly plot the transfer matrix itself in Fig. 4. We select jets with R = 0.8 and p T,jet ∈ [120, 150] GeV in central dijet events and measure the jet width, λ 1 1 , either on ungroomed jets (left plot) or on groomed jets (right plot). At hadron level, the jet width is measured on charged particles only. The transfer matrix for ungroomed jets shows clear signs of bin migration, especially at low-to-mid angularities where the parton-level values predominantly get pushed to larger values of the jet width at hadron level. If we instead apply grooming, the transfer matrix appears to be much more diagonal, albeit with longer tails away from the diagonal. ¶ Finally, it is interesting to compare our new method with the simpler approach we have used in earlier studies, which consisted in accounting for NP contributions via a bin-by-bin HL/PL ratio. In Fig. 5 we perform such comparison for an observable known to be rather sensitive to NP contributions, namely the groomed LHA λ 1 0.5 measured on charged tracks on the hardest jet in Z+jet events [51]. We again consider the moderate p T,jet region, here both for R = 0.4 (left), and R = 0.8 (right). The old approach, used in [51], is shown in blue, ‖ while the results obtained with the new transfer-matrix approach are shown in red. For both predictions we estimate theoretical uncertainties, illustrated by the bands, corresponding to the envelope of the 7-point scale variations, the alternative x L -parameter settings, and variations of the α MPI min parameter of the UE model. We see that the difference in the nominal predictions is rather substantial, and that the new treatment of NP corrections yields a significantly better agreement with the CMS data, shown in black, although visible differences remain for the R = 0.8 case. One notices that the uncertainty estimates for the results based on the transfer-matrix approach are somewhat larger than for the old ratio method. A source of this increase is the larger range of kinematics probed when allowing for migration from lower and higher p T,jet slices.
We conclude this section by noting that a similar procedure to account for NP effects was employed by the ALICE collaboration in Ref. [54], where a measurement of various angularity observables in inclusive jet production in pp collisions at √ s = 5.02 TeV was reported. These data, based on SoftDrop- ¶ At larger p T,jet , not shown here, the transfer matrix for groomed jets shows clear signs of migration from relatively large partonic values of the angularity to very small hadron-level values. These are likely related to subjets passing the SoftDrop condition at parton level being pushed below the SoftDrop cut after hadronisation, an effect which is typically included in analytic treatments of non-perturbative corrections of SoftDrop observables [21,23,92]. ‖ Note that, for a meaningful comparison between the two methods, both the approach based on the HL/PL ratio and the one based on transfer matrices are derived using the SHERPA generator with variations of the UE parameter α MPI min .  groomed charged-tracks jets, was compared to analytical predictions at NLL accuracy obtained from SCET [48,49]. Two techniques for correcting the analytic predictions for NP effects and the selection of charged-particle jets were considered. The first being based on "folding" the NLL result with a response matrix extracted from MC simulations that maps parton-level jets with (p PL T , λ PL ) to hadron-level jets with (p HL T , λ HL ), thereby accounting for hadronisation corrections only. To incorporate UE effects an additional bin-wise correction has been applied. As an alternative a NP shape-function approach [49] to simultaneously correct for hadronisation and the UE has been employed. We refer to [54] for additional details.

Results for jet angularities in dijet and Z+jet production
In this section, we present the results obtained from the calculation detailed in section 3, i.e NLO + NLL accurate predictions accounting for NP corrections through the transfer-matrix approach. We start in section 4.1 with a few considerations at the purely perturbative level. In section 4.2, we then discuss our results with NP corrections included, and present full hadron-level predictions at NLO QCD accuracy from SHERPA.

Selected parton-level results
As discussed in section 3.2, the implementation of the CAESAR resummation in the SHERPA plugin [56] relies on a separation into different flavour channels δ. Since we obtain fully exclusive parton-level events from COMIX, we can easily apply the (BSZ) flavour-clustering algorithm from Ref. [83] so as to identify the flavour assignment of a given phase-space point. We use this algorithm iteratively, following the procedure introduced in [51] for the leading jet in the Z+jet case. This can trivially be extended to the central/forward jet in dijet production. It allows us to talk about the flavour of a particular anti-k t jet of radius R in an IR-safe way. Note that in practice we run the bland variant of the algorithm whereby each jet gets identified as either quark-or gluon-like. While flavour identification is a necessity in the context of our matching scheme (at least at LO in the angularity distribution), these results also provide a well defined way to analyse the flavour decomposition of jet cross sections. Such truth-level flavour Ratio to Born (c) Z+jet Figure 6: Fraction of gluon jets as a function of the jet transverse momentum p T,jet in dijet (left and middle) and Z+jet (right) processes. For the dijet production we present results for the leading central and forward jet, respectively. We show predictions at Born (red bands) and at one-loop level (blue bands). Flavour tagging for the considered anti-k t jets relies on the algorithm described in [51], that employs the BSZ flavour-k t algorithm, and in turn defines the partonic channels δ in our implementation of NLL CAESAR resummation in SHERPA.
assignment could for example be used as well-defined input, e.g. for machine-learning based methods, or as benchmark for flavour-tagging algorithms [106][107][108][109][110][111][112]. We start our discussion by considering the fractions of gluon and quark anti-k t jets in the Z+jet and dijet final states. These are compiled in Fig. 6 for the central-and forward jet in dijet events, as well as in Z+jet production, as a function of the jet transverse momentum, p T,jet . On each plot, two curves are shown together with their respective scale uncertainties from the 7-point scale variation: gluon fractions at Born level (red bands), i.e. O(α 2 s ) for dijets and O(α 2 EW α s ) for Z+jet * * , and for the NLO QCD matrix element (blue bands), at O(α 3 s ) and O(α 2 EW α 2 s ), respectively. Note that the (NLO) one-loop result corresponds to the LO accuracy for jet-angularity distributions. One sees that for low-p T dijet events, both the central and forward jets have a large fraction of gluon jets (∼ 70%), and that high-p T dijet events and Z+jet events at any p T,jet , with gluon fractions of about 20−30%, are dominated by quarks.
A striking feature of Fig. 6 is the size of the NLO QCD corrections. For dijet production, we see an (absolute) decrease of about 15% of the gluon fraction, regardless of the transverse momentum and rapidity of the jet, marginally larger than the estimated one-loop perturbative uncertainty. For Z+jet events, NLO corrections cause a 2−10% (absolute) increase with a clearly-visible dependence on the jet p T . This is in line with earlier observations that the NLO corrections to the gluon fraction in Z+jet events are both larger than in dijet events and show a stronger jet-p T dependence (see e.g. [113,24]). It is also worth commenting on the relative size of the theoretical uncertainty, which appears to increase when going from Born-level to NLO QCD. The reason for this counterintuitive behaviour can be traced back to the fact that the considered gluon fractions are determined by cross-section ratios. At Born-level the dependence on α s , and hence on the renormalisation scale, exactly cancels between numerator and denominator, leaving only a rather weak dependence on the factorisation scale.
Our results can be directly compared to the corresponding gluon-jet fractions reported by CMS in [30] (see Fig. 2 therein and the corresponding discussion), although that study was done using a different definition of gluon jets based on MADGRAPH5 aMC@NLO +PYTHIA simulations as well as a series of other generators. Qualitatively, the gluon fractions presented in [30] show the same pattern as the ones obtained here. However, with the BSZ-based approach we find gluon fractions that are generally smaller than the ones reported by CMS. This is especially the case in Fig. 6c for the Z+jet process where CMS found gluon fractions reaching almost 40% at large p T,jet . This is in agreement with their comment that other generators predict up to 25% smaller gluon fractions in the Z+jet sample, which is also where the largest NLO corrections are observed. Next, Figs. 7 and 8 show examples of our matched NLO + NLL results for the LHA λ 1 0.5 distribution, for different p T,jet and event selections. In each case, we have separated the total cross section into contributions from quark and gluon jets. The total result is given by the sum of both components, indicated by the solid (black) line. The fraction of gluon jets in each bin is shown in the lower panels. The selected p T,jet ranges are chosen such that they coincide with the gluon-and quark-enriched samples studied by CMS in [30] (see also table 2 in the following section). The results confirm the findings from Fig. 6: we indeed see that the high-p T dijets and the Z+jet sample for p T,jet ∈ [120, 150] GeV are dominated by quark jets, and the low-p T dijet events are instead dominated by gluon jets. It is also clearly visible in these figures that the gluon distributions contribute at larger values of the angularity variable than the quark distributions, indicative of their potential as a quark-gluon discriminator. Groomed LHA λ 1 0.5 Theory/Data Figure 9: NLO + NLL + NP and hadron-level SH-MC@NLO predictions for the differential cross section in the Les Houches Angularity λ 1 0.5 for ungroomed (top row) and groomed (bottom row) R = 0.4 anti-k t jets with p T,jet ∈ [120, 150] GeV, compared to data from CMS [30]. The left and middle panel correspond to the central and forward jet in dijet events, respectively, the right one to the leading jet in Z+jet production. The NP corrections to the perturbative NLO + NLL prediction have been obtained with the transfer-matrix approach.

Hadron-level results
In this section we present our final NLO + NLL + NP predictions, with NP corrections implemented through the transfer-matrix approach, for jet angularities in both dijet and Z+jet production and compare them to the results obtained by the CMS collaboration. To avoid the proliferation of plots, we concentrate mostly on the LHA λ 1 0.5 measured on (groomed or ungroomed) R = 0.4 anti-k t jets at moderate transverse momentum, namely p T,jet ∈ [120, 150] GeV. As before, we consider separately the central and forward jet in dijet events, as well as the leading jet in Z+jet events, we thereby restrict to the data for all hadrons. Corresponding results for the jet width λ 1 1 and jet thrust λ 1 2 we compile in App. A. The full set of predictions, i.e. for all the p T,jet slices considered by CMS in [30], for jet radius R = 0.8, as well as based on charged tracks, and their comparison to the data from CMS are also publicly available [114] † † . To estimate the theoretical uncertainty of our NLO + NLL + NP predictions we consider 7-point variations of the factorisation and renormalisation scales, as well as variations of the x L parameter, see section 3.2 for details. Besides the NP corrections corresponding to the default parameter set in SHERPA, we also consider up-and down variations of the UE activity and use Lund string fragmentation as alternative hadronisation model, cf. section 3.1. This way we derive four alternative hadron-level distributions for each combination of scale-and x L -variations. The final uncertainty is then obtained by taking the envelope of all those variations.
Alongside the NLO+NLL +NP results we present hadron-level simulations with SHERPA at NLO QCD accuracy, as described in section 3.1. The perturbative uncertainty of the SH-MC@NLO predictions are taken as the envelope of 7-point variations of µ F and µ R in the matrix elements and the parton shower.
In Fig. 9 we present NLO + NLL + NP and SH-MC@NLO predictions for normalised differential cross sections in λ 1 0.5 and compare them with the CMS experimental data. The top row of plots thereby corresponds to ungroomed jets, while the bottom row ones are obtained with SoftDrop (β = 0, z cut = 0.1) applied to the jets prior to the observable evaluation. The left-hand plots correspond to LHA measurements on the central jet, the middle ones on the forward jet in dijet events, respectively, and the right-hand ones on the leading jet in Z+jet production.
Overall, our resummed and matched predictions when corrected of NP effects (shown in red) provide a good description of the hadron-level data, with the exception of the last (and in some cases the second to last) bin at large values of the angularity. The corresponding region of phase space is outside the jurisdiction of the all-orders calculation and one might have hoped that it would be well-described by the NLO contribution. However, the last bin contains the kinematic endpoint of the fixed-order calculation, accordingly, this part of the distribution is very sensitive to the effect of multiple emissions. Indeed the SH-MC@NLO predictions (in blue) are able to populate this region of phase space through additional parton-shower real radiation, resulting in a better description of the data. In the future, it would be interesting to see if higher-order, e.g. NNLO, corrections will yield an improved description in this largeangularity region. For the groomed distribution, we also see that the MC simulation provides a better description of the peak region than the analytic calculation. One should however bear in mind that all the λ 1 0.5 bins, except the lowest one, are in a region with λ 1 0.5 ≥ z cut , i.e. not directly affected by SoftDrop at NLL accuracy. In this region, non-trivial subleading effects can have a sizeable impact (see also Refs. [91,115] for discussions and potential improvements).
A crucial difference between the resummed and the SH-MC@NLO results is the size of the theoretical uncertainty, which is much larger for the NLO + NLL + NP calculation. We have investigated this feature by decomposing the total uncertainty into its various perturbative and non-perturbative contributions. This clearly led us to the conclusion that this effect is dominated by the resummation-scale variation, i.e. the variation of x L in Eq. (3.14). A systematic reduction of this uncertainty would require to improve the accuracy of the resummation, i.e. to include NNLL contributions. We will comment on the feasibility of this calculation in the conclusions. We furthermore note that in the here shown MC simulations no corresponding variation is performed (see section 3.1). It would be interesting to study systematic variations when using a different shower model, e.g. the DIRE cascade [116] available in the SHERPA framework.
Following the analysis performed by the CMS collaboration in Ref. [30] we also consider the mean value of the angularity distributions, as a function of p T,jet . Corresponding results for the case of the λ 1 0.5 angularity are shown in Fig. 10. The layout is organised as in Fig. 9 with ungroomed jets in the top row and SoftDrop-groomed jets in the bottom row and with, from left to right, the central and forward jet in dijet events and the leading jet in Z+jet events. Analogous results for the jet-width and jet-thrust angularities are collected in App. A.
The comparison between the resummed and the MC predictions highlights the feature previously discussed, namely that the estimated uncertainties of the NLO+NLL +NP results are significantly larger than the ones obtained for the SH-MC@NLO simulations, with the former being dominated by the x L variation. We observe that both theoretical predictions consistently underestimate the experimentally measured mean values. From the ratio plots in the bottom panels we can read off that the theoryto-data ratios for the central values are almost constant, exhibiting only a very mild dependence on the transverse momentum. This holds in particular for the SH-MC@NLO predictions that undershoot the data by about 5 − 10%. Larger deviations, reaching up to 18% at the highest p T,jet values, are seen for the NLO + NLL + NP results in the ungroomed case. In the groomed case, although our NLO + NLL + NP predictions systematically undershoot the data, we observe agreement, within the theoretical uncertainties. In Ref. [30], similar results have been found for the LO MC generators considered there. Part of this effect can be explained by the observed underestimation of large angularity values in the theoretical predictions when comparing to data. However, in this region the cross section is quite low, thus it contributes rather little to the average value λ 1 0.5 . As we have seen in section 3.3, the angularity distribution is significantly affected by NP corrections. Accordingly, it would clearly be interesting to include the new data on the jet angularities in tunes of the NP models, thereby further investigating which parameters affect their description. In turn, improved NP transfer matrices could be derived and applied to the NLO + NLL predictions. Detailed studies of jet angularities allow us to quantitatively assess how well MC simulations, and more generally theoretical calculations, describe the particle distribution inside jets. In turn, they offer potential to analyse the description of QCD radiation off quarks and gluons, separately. This issue was, for instance, studied in Refs. [32,33]. The measurement of CMS presented in [30] now allows us to test theoretical predictions on jets from gluon-enriched and quark-enriched samples. For this purpose, the analysis identified five interesting phase-space regions and parameter choices, detailed in table 2, that can be classified as having an enhanced gluon-or quark-jet contribution in the dijet or Z+jet process, respectively. We refer to section 4.1 and Fig. 6 in particular for a discussion on flavour fractions in these specific channels. SoftDrop (β = 0, z cut = 0.1) R = 0.4 [120,150] dijet central Z+jet For these five configurations we consider the three IRC safe angularities λ 1 α with α ∈ {0.5, 1, 2}, i.e. the LHA, jet width, and jet thrust. We present in Fig. 11  Both theoretical approaches yield results that are in fair agreement with the data. The deviations of the central values from the measured ones are in fact rather similar for all three observables, and for all the configurations under consideration. Also, we do not observe a qualitative difference in the description of the g-or q-enriched samples. Note that, both in the data and in the predictions, we observe smaller jet angularity mean values for the q-enriched samples compared to the g-enriched case. This is theoretically anticipated given that gluons carry more colour charge and accordingly radiate more. For the highp T configuration (2) this effect is, as expected, rather marginal, since the fractions of gluon jets are relatively similar in the central and forward cases. Our theoretical calculations nicely capture this effect even quantitatively. As seen for the LHA already, the theoretical uncertainties on the NLO + NLL + NP predictions appear to be much larger than the SH-MC@NLO ones, motivating to consider the evaluation of the NNLL, and eventually NNLO, corrections. From these considerations we can conclude that for the calculations considered here, i.e. NLO + NLL + NP and SH-MC@NLO, QCD radiation off both hard quarks and hard gluons is well-modelled. This is interesting as it was noted before [32,33] that generalpurpose MC event generators do not always agree in their description of QCD radiation off gluons, while they largely do for radiation off quarks, heavily constrained by tunes on LEP data. However, it has to be noted that, although the considered samples are certainly q-and g-enriched, they still contain significant contributions from the respective other flavour channel, cf. Fig. 6.
To further study the quality in the modelling of gluon and quark jets, we next consider ratios of the angularity mean values in the g-and q-enriched samples for the five phase-space selections. In Fig. 12 we present our corresponding theoretical predictions, normalised to the result reported by CMS. While the left-hand plot contains the results for the NLO + NLL + NP calculation, the right-hand plot shows corresponding predictions obtained from SH-MC@NLO. values agree remarkably well with the measured ones i.e. have a ratio to data centred close to unity. For the SH-MC@NLO simulation we observe the largest deviation, about 10%, for the jet-thrust angularity for the case of SoftDrop-groomed jets of moderate-p T , see also Figs. 15 and 16 in App. A.
This good description reflects the fact that, despite of the notorious underestimation of the angularity mean values seen in Fig. 10, both our calculational approaches seem to treat quark and gluon jets similarly well, and do not introduce an artificially-enhanced discriminative power between both hypotheses. This is to some degree in contrast to the generator predictions studied by CMS in [30]. There, in particular for the low-p T,jet regions, predictions based on the Born matrix elements for Z+jet and dijet production overestimated the difference between gluon and quark jets by up to 20 − 30%. This might to some extent originate from the different gluon-and quark-jet decomposition they predict for the various phase-space regions (see Fig. 6 and the discussion of flavour fractions at Born and NLO QCD accuracy). We finally point out that, as noted before, the uncertainty estimate for the NLO + NLL + NP predictions is quite large, in particular for jet thrust, λ 1 2 .

Conclusions
Following a recent measurement accomplished by the CMS collaboration [30], we have performed a detailed phenomenological analysis of IRC safe jet angularities, considering both ungroomed and SoftDrop jets. Our theoretical predictions account for all-order effects at next-to-leading logarithmic accuracy, including non-global logarithms, and matching to next-to-leading order QCD calculations. Thanks to a flavour-sensitive matching procedure, based on the BSZ flavour-k t algorithm, we are able to reach NLO + NLL accuracy. Furthermore, our predictions are supplemented by non-perturbative corrections, that account for hadronisation effects and the Underlying Event. With respect to our previous work, see e.g. [51], our theoretical predictions feature two improvements. First of all, on the perturbative side, we now include resummed and matched calculations for dijet events, in addition to the ones for Z+jet production recently presented in [51] and confronted with experimental data in [30] already. The resummation for dijet final states features more complicated colour structures, for both the global and the non-global part, which is dealt with by our largely-automated implementation of colour matrices in the SHERPA resummation framework [56,80].
Secondly, we employ a more sophisticated approach for dealing with non-perturbative corrections that we extract from fully exclusive hadron-level simulations at MC@NLO accuracy with the general-purpose event generator SHERPA. Rather than relying on a simple bin-by-bin rescaling, we have developed and implemented a rather general transfer-matrix approach that can account for migrations across different kinematical variable and observable bins. In the present case, we account for migrations both in the jet transverse momentum and in the angularity observable when going from parton-to hadron level. We thereby largely reduce the sensitivity to differences in the perturbative predictions from the analytic resummation and the employed MC simulations. This cures deficiencies of our previous approach, that showed pathological features for observable values corresponding to scales below the parton-shower cutoff.
We compare our final NLO + NLL + NP results to a plethora of distributions measured by CMS, both in dijet and Z+jet production. We thereby consider differential cross sections for the LHA, jet width and thrust, as well as their respective mean values measured as a function of the jet transverse momentum. For comparison, we also include results from SHERPA at MC@NLO accuracy. Overall, we find a fair agreement between the data and the theoretical predictions. Our results however tend to systematically underestimate the mean values measured by CMS. With respect to the comparisons presented in Ref. [30] for the Z+jet selection, we clearly achieve an improved description of the data, thanks to the new transfer-matrix method used to include non-perturbative corrections. Finally, still following the CMS analysis, we compare our theoretical predictions to the measurements performed on quark-and gluon-enriched samples. We find that both the resummed calculation and the SH-MC@NLO prediction are able to describe the data well. This is interesting because jet angularities have been proposed as IRC-safe quark/gluon taggers [38,39] that can be applied not only in the context of searches for new physics but also to Standard Model measurements that aim for extractions of fundamental parameters, such as the strong coupling, or parton distribution functions. The better description of the data that we find originates to some degree from a more accurate modelling of the flavour fractions, i.e. the relative contributions of partonic channels, that receive sizeable NLO QCD corrections, included in both our theoretical predictions.
Throughout this paper, we have noted that our resummed calculations suffer from rather large theoretical uncertainties, dominated by variations of the resummation scale. An obvious way to systematically improve on this is to promote the resummed calculation to NNLL accuracy. Many of the contributions that are relevant for NNLL resummation of jet shapes have already been computed [44], predominately using SCET techniques, and have reached a considerable degree of automation, see e.g. [117][118][119]. Furthermore, first calculations for groomed observables exist even at N 3 LL [25] and with an improved description of the transition between groomed and ungroomed regimes [115]. The resummation framework we are employing in our plugin also has already been extended, for global observables, to NNLL [120]. The bottleneck in this enterprise is the inclusion of subleading non-global logarithms. The first resummation of these effects has been achieved very recently [121,122] and we look forward to investigating whether this method can be easily interfaced with our framework, in order to perform higher-precision phenomenology of jet angularities.
Another obvious improvement of our calculation would be to upgrade the fixed-order component of the angularity distributions to NNLO accuracy. This requires computing the hadronic production of three QCD partons, respectively Z+2 partons, at two-loop accuracy. In this context, we note that the first NNLO evaluation of three-jet observables has recently been presented [123]. Finally, the two-loop corrections to the dijet and Z+jet processes (see [124,125] and references therein) are needed to achieve NNLL accuracy.
In the context of assessing the impact of NP corrections, it would be interesting to compare the transfer-matrix approach developed here to results obtained using first-principle field-theoretical arguments (see for instance Refs. [92,93] for recent work on hadronisation corrections for SoftDrop jets). In a more generic context, it would be interesting to see the impact that the measurement of jet angularities has on tuning NP parameters of general-purpose Monte Carlo generators. One could then also provide independent tunes for different parton-shower cutoffs. This could be, in turn, added to our transfer-matrix approach of NP corrections to parton-level analytic calculations. We close by noting that it would clearly be interesting to apply our calculational methods to jets arising from different production modes, e.g. involving top-quarks [55,53], or at different centre-of-mass energies [54]. We leave this for future work.

Acknowledgements
We would like to thank Andreas Hinzmann for fruitful discussions and comments on the manuscript. We are grateful to Vincent Theeuwes for collaboration at an early stage of the project. This work is supported by funding from the European Union's Horizon 2020 research and innovation programme as part of the Marie Sk lodowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104). SS and DR acknowledge support from BMBF (contracts 05H18MGCA1 and 05H21MGCAB). SS acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -project number 456104544. The work of SC, OF and SM is supported by Università di Genova under the curiosity-driven grant "Using jets to challenge the Standard Model of particle physics" and by the Italian Ministry of Research (MUR) under grant PRIN 20172LNEEZ. We also acknowledge support from the Royal Society through the project 580986 "Resum(e) The Path To Discovery". SC would like to thank the Institute for Particle Physics Phenomenology (Durham University) for hospitality during the course of this work.

A Results for jet width and jet thrust
In this appendix, we collect additional results obtained with our NLO + NLL + NP calculation and SH-MC@NLO for the jet-width (λ 1 1 ) and jet-thrust (λ 1 2 ) angularities, for anti-k t , R = 0.4 jets. In Fig. 13 we show the jet-width distributions for p T,jet ∈ [120, 150] GeV, for the central and forward jet in dijet events, as well as the leading jet in Z+jet production, for ungroomed and SoftDrop jets. In