Jet Angularities in Z+jet production at the LHC

We present a phenomenological study of angularities measured on the highest transverse-momentum jet in LHC events that feature the associate production of a $Z$ boson and one or more jets. In particular, we study angularity distributions that are measured on jets with and without the SoftDrop grooming procedure. We begin our analysis exploiting state-of-the-art Monte Carlo parton shower simulations and we quantitatively assess the impact of next-to-leading order (NLO) matching and merging procedures. We then move to analytic resummation and arrive at an all-order expression that features the resummation of large logarithms at next-to-leading logarithmic accuracy (NLL) and is matched to the exact NLO result. Our predictions include the effect of soft emissions at large angles, treated as a power expansion in the jet radius, and non-global logarithms. Furthermore, matching to fixed-order is performed in such a way to ensure what is usually referred to as NLL$'$ accuracy. Our results account for realistic experimental cuts and can be easily compared to upcoming measurements of jet angularities from the LHC collaborations.


Introduction
During the first two runs of the CERN Large Hadron Collider (LHC) the experimental collaborations have accumulated a vast amount of high-quality data of proton-proton collisions at energies as high as 13 TeV. In the current shutdown period, the collaborations refine their search and measurement strategies, in order to fully exploit the physics potential of the acquired data, and to prepare for the upcoming third run of the LHC. However, no significant increase in collision energy is planned for the foreseeable future. Accordingly, the focus of the theory and experimental communities should be on devising new effective and robust analysis techniques to interrogate the data, so that no stone is left unturned. In this context, the use of deep-learning algorithms to augment performance is now becoming standard practice (see e.g. [1][2][3] for recent reviews). Complementary to this effort, measurement campaigns that deliver unfolded data that can be compared to state-of-the-art theoretical calculations have been, and will further be, pursued in order to stress-test our understanding of the Standard Model.
-1 -non-global logarithms, in the limit of large number of colours (N c ). 1 Furthermore, our initial detailed study of MC predictions allows us to supplement our perturbative NLO + NLL predictions with non-perturbative corrections extracted from hadron-level MC simulations. Our theoretical predictions are fully differential in the kinematics of the jets and of the leptonic decay products of the Z boson, so that we can impose realistic fiducial cuts and compare them to the distributions obtained from MC event generators. Our predictions are obtained with the resummation plugin of the SHERPA generator framework [39] and they can be directly compared to unfolded measurements from the LHC collaborations, once these become available. 2 The paper is organised as follows: in Section 2 we introduce the class of angularity observables to consider and detail our event selection cuts. In Section 3 we present and compare hadron-level predictions for groomed and ungroomed jet angularities from the HERWIG and PYTHIA event generators, based on the leading-order matrix elements, and SHERPA, using multijet merging at leading and nextto-leading order. Section 4 is devoted to the all-orders evaluation of jet angularities at NLO + NLL accuracy using the SHERPA resummation framework. In Section 5 we compare our resummed predictions with the parton-level results of our MC simulations. We extract non-perturbative corrections from the full particle-level simulations and apply them to our NLO+NLL predictions. Our conclusions are presented in Section 6. Additional results are collected in Appendix A.

Observable definition
In this study we concentrate on the production of a leptonically decaying Z boson associated with one or more jets at the LHC. Jets are defined using the anti-k t clustering algorithm [44] with radius parameter R 0 and standard E-scheme recombination, i.e. the momenta of objects that are paired together are simply added so that the resulting jet momentum is given by the sum of its constituents' momenta.
On the hardest jet, i.e. the jet with the largest transverse momentum p T,jet , we measure the angularities [34,45,46] λ κ α = i∈jet p T,i where is the Euclidean azimuth-rapidity distance of particle i from the jet axis. Our study is limited to a subset of the above defined general angularities, namely the ones that exhibit infrared and collinear (IRC) safety. This poses the restrictions κ = 1 and α > 0. Furthermore, it is well-known [29,47] that IRC safe angularities with α ≤ 1 are sensitive to recoil against soft emissions, leading to a rather complicated resummation structure. To circumvent these additional complications, we compute for these angularities the distance measure in Eq. (2.2) with respect to the jet axis obtained by reclustering the jet constitutents with the anti-k t algorithm but using the Winner-Take-All (WTA) recombination scheme [48]. We are also interested in computing angularities on groomed jets. In this case, we recluster the jet with the Cambridge-Aachen (C/A) algorithm and consider the SoftDrop grooming algorithm with parameters z cut and β [21]. The angularity is then computed on the resulting SoftDrop jet, i.e. both sums in Eq. (2.1) are restricted to the particles that survived the grooming. For groomed jets we also adopt the WTA prescription for angularities with α ≤ 1.
In order to highlight the behaviour of the angularity distributions in different energy regimes, standard and groomed angularities are considered in different jet transverse momentum bins. This will allow us, for instance, to better study the impact of non-perturbative corrections, e.g. from the parton-to-hadron transition. We note that, in order to avoid issues related to bin-migration, which greatly complicates the structure of the resummation of the SoftDrop β = 0 (modified MassDrop) angularities [49], we shall always choose the reference transverse momentum p T,jet to be the one before grooming. Our theoretical predictions, both the ones obtained with matrix-element improved parton showers and with all-order resummation matched to fixed-order calculations, account for fiducial cuts and therefore can directly be compared to unfolded measurements.

Event selection cuts
We close this section by detailing the fiducial volume for the jets and leptons in our subsequent studies. Our choices follow the selection of a recent (preliminary) CMS measurement [50]. We consider the inclusive production of a pair of oppositely charged muons in proton-proton collisions at 13 TeV centre-of-mass energy. We require all final state particles to have pseudo-rapidity |η| < 5. For both muon candidates we require p T,µ > 26 GeV , and |η µ | < 2.4 . (2. 3) The lepton pair has to pass the additional conditions 70 GeV < m µ + µ − < 110 GeV , and p T,µ + µ − > 30 GeV . (2.4) In what follows we will refer to the lepton pair as Z boson, where nevertheless we imply the inclusion of off-shell effects. We then select events that exhibit at least one anti-k t jet [44] with |y jet | < 1.7 , and R 0 = 0.8 3 , (2.5) and consider several bins in p T,jet starting at 50 GeV following [50]. We compiled results for all bins used there 4 , but will in the following limit the discussion to p T,jet ∈ [120, 150] GeV and p T,jet ∈ [408, 1500] GeV that we find representative for the lower and higher jet scale selections. In order to gain better control over NLO QCD corrections for the Zj production process, that become large when the Z-boson transverse momentum is significantly smaller than p T,jet [51], we require the transverse momenta of the lepton pair and the leading jet to be largely balanced. To this end we impose the constraint Finally, we require the Z-boson and the leading jet to be well separated in azimuthal angle, i.e.
3 State-of-the-art Monte Carlo study of jet angularities The Born-level contributions to Zj production are given by the two channels pp → µ + µ − + q/g, where q represents a massless quark or anti-quark. Jet angularities become non-trivial upon inclusion of additional radiation, which in general purpose MC generators is accounted for by parton shower simulations [52]. However, to properly include non-logarithmic corrections, e.g. from hard real emissions, these should be matched to exact higher-order matrix elements, as obtained, for example using event generators such as MCFM [53][54][55], POWHEG [56][57][58] or MADGRAPH5_aMC@NLO [59,60]. Here we consider simulations including the complete set of NLO QCD corrections to the Zj and Zjj production processes matched to the SHERPA parton shower in the MEPS@NLO formalism [61]. To address the sensitivity of the angularities to non-perturbative corrections we also account for UE contributions to the collision final states and model the parton-to-hadron transition. In this section we describe the calculational setups for our Monte Carlo simulations and then present corresponding predictions for the groomed and ungroomed LHA, Width, and Thrust angularities.

Monte Carlo generator setup
We compile MC predictions for the jet angularities with the SHERPA [62,63] event generator, version 2.2.10, using the NNPDF-3.0 NNLO PDF set [64]. We consider the inclusive Zj production process in the MEPS@NLO multijet merging formalism [61], thereby combining the NLO QCD matrix elements for µ + µ − j and µ + µ − jj production, matched with the SHERPA Catani-Seymour dipole shower [65]. We set the merging-scale parameter to Q cut = 30 GeV (which is of the order of the jet p T cut used in our event selection). We obtain the QCD one-loop amplitudes for the one-and two-jet processes from RECOLA [66,67], using the COLLIER library [68] for the evaluation of tensor and scalar integrals. To assess the impact of the QCD one-loop corrections, we furthermore compile MEPS@LO [69] predictions based on merging the one-and two-jet leading-order matrix elements, using Q cut = 30 GeV as well.
The perturbative scales entering the calculation are defined according to the CKKW-style scale setting prescription [69,70]. In this procedure the hard-process partons are clustered into a Born-like 2 → 2 configuration that defines the core process with an associated scale µ core . For the production channels considered here the event-wise determined core process, of order α ncore s , can correspond to jj → µ + µ − (n core = 0), jj → Zj (n core = 1), or jj → jj (n core = 2). For the three possible cluster configurations the corresponding core-process scale is given by The core-process scale is then used to define the factorisation scale and the parton shower starting scale of the core process, i.e. µ F = µ core , µ Q = R 0 µ core . The effective renormalisation scale, µ CKKW , of the n-parton hard matrix elements corresponds to with t i the reconstructed shower-branching scales. To estimate the perturbative uncertainties of our MC predictions, we do on-the-fly [71] 7-point variations [72] of the factorisation and renormalisation scales in the matrix elements and the parton shower. The uncertainty bands given later on correspond to the envelope of the settings To model the parton-to-hadron transition we use the SHERPA cluster fragmentation model [73]. The UE simulation relies on the SHERPA implementation of the Sjöstrand-Zijl multiple-parton interaction model [74]. In both models the default set of tuning parameters is used, see [63] for details.
To obtain further independent predictions, in particular for the modelling of non-perturbative effects, we compile additional results from HERWIG [75,76], version 7.2.1, and PYTHIA [77], version 8.303. For both generators we consider the µ + µ − j-production process at leading order, 5 invoking the respective default models for the QCD parton shower, hadronisation and UE simulation.
For event selection and analysis we employ the RIVET analysis package [80,81]. For jet reconstruction we use the FASTJET package [82]. For the SoftDrop grooming we rely on the implementation in the RecursiveTools class which is a part of the FASTJET contrib package.

Hadron-level Monte Carlo predictions
To set the stage we begin by directly comparing the hadron-level predictions obtained from the various event generators, i.e. calculational setups, for the three considered jet angularities, both, with and without SoftDrop grooming. These will later serve as means to extract non-perturbative corrections and estimates for related uncertainties for our resummed predictions, see Section 5.
In Fig. 1 and Fig. 2 we collate the generators' predictions for the lower and higher p T,jet slice, respectively. For each angularity we provide results without (left column) and with (right column) SoftDrop grooming. The broken x-axes indicate the fact that the first bins start at zero and, hence, they appear unbounded on a logarithmic scale. For each individual panel we include a ratio plot, with the SHERPA MEPS@NLO result taken as reference. For the SHERPA MEPS@LO and MEPS@NLO predictions we include uncertainty bands, green solid and red hatched, respectively, that reflect the envelope of the 7-point scale variations in both the matrix-element and parton-shower component. The inclusion of the exact NLO QCD one-loop corrections in the MEPS@NLO method results in a significant reduction of scale uncertainties, for all angularities both in their groomed and ungroomed variants.
Considering the ungroomed angularities for the lower p T,jet window first, we observe that with increasing α the distributions peak at lower observable values. While the LHA angularity distribution exhibits a Sudakov peak around λ 1 1/2 ≈ 0.25, for the Thrust variable the maximum is at λ 1 2 ≈ 0.05. All generator predictions largely agree on the peak position. However, for lower observable values quite sizeable deviations can be observed, reaching and partially exceeding 50%. Notably, the two SHERPA predictions agree quite nicely in this observable range, that in fact is significantly affected by non-perturbative corrections, described by the same models in the MEPS@LO and MEPS@NLO simulations. However, HERWIG and PYTHIA use alternative models and parameter tunes for hadronisation and UE. Both generators predict somewhat narrower distributions in comparison to SHERPA. For large values of the angularities, i.e. towards the kinematical endpoints, the MEPS@NLO calculation predicts somewhat larger event fractions than what is obtained using LO matrix elements, with the largest deviation appearing in comparison to the PYTHIA LO result. This systematic effect which barely exceeds 15% is, however, within the scale uncertainty band. It is interesting to note that a similar pattern between LO and NLO matrix elements was already observed by the CMS collaboration [83] when comparing their measurement of the groomed jet mass to analytic QCD predictions.
Upon invoking SoftDrop grooming of the jet constituents, the spread in the generator predictions is sizeably reduced. In particular, in the region of very small observable values, dominated by nonperturbative effects in the ungroomed case, the predictions now agree to within 10%. It is clearly visible that grooming decreases the jet-angularity values, and in fact significantly sculpts the observable distributions. With increasing α, a secondary peak emerges below the grooming transition point that moves towards smaller values of λ. For large observable values grooming has much smaller impact and the spread in the generator predictions observed in the ungroomed case basically remains unaltered.
For higher jet transverse momenta, i.e. p T,jet ∈ [408, 1500] GeV presented in Fig. 2, the picture somewhat changes. Both for the groomed and ungroomed angularities we observe a larger variation of the predictions for λ 0.1. The differences between the SHERPA MEPS@NLO predictions and the PYTHIA LO results reach up to 30%. This increase with respect to the lower-p T,jet slice can be expected as the inclusive K-factors, with the fiducial cuts adopted here, are also found to be increasing functions of jet transverse momentum. It is interesting to note that the SHERPA MEPS@LO results agree very well with the HERWIG predictions, for all three considered angularities, with and without grooming. The inclusion of the full set of one-loop virtual corrections for the one-and two-jet processes in the MEPS@NLO calculation translates into corrections of O(10 − 15%) for larger values of λ. Despite the larger sensitivity to higher-order perturbative corrections, the relative impact of non-perturbative effects, most relevant for small angularity values, is reduced for jets with high transverse momentum.
We note that measurements of (groomed) jet angularities differential in p T,jet provide means to probe the modelling of perturbative as well as non-perturbative corrections in MC generators. In Sec. 5 we will use the presented simulations to extract non-perturbative corrections for the resummed calculations. There we will also present a comparison of the generators' parton-level predictions, i.e. without the inclusion of UE and hadronisation, with the NLO + NLL results.   -9 -To perform the NLL resummation of the jet angularities we use the SHERPA implementation of the CAESAR resummation formalism [47], first presented in [39]. This framework has recently been employed to obtain resummed predictions for SoftDrop thrust [84] and multijet resolution scales [85] in electron-positron collisions, as well as NLO + NLL predictions for SoftDrop groomed hadronic event shapes, in particular groomed transverse thrust [86]. For a generic observable, the master formula for the all-order cumulative cross section for observable values up to v, with L = − ln(v), can be written as a sum over partonic channels δ: where dσ δ dB δ is the fully differential Born cross section for the partonic channel δ and H implements the kinematic cuts applied to the Born phase space B. F denotes the multiple emission function which, for additive observables such as the angularities considered in this paper, is simply given by The ratio of partondistribution-functions (PDFs) P takes into account the true initial-state collinear scale. The soft function S implements the non-trivial aspects of colour evolution. The collinear radiators R l for the hard legs l were computed in [47] for a general observable V scaling for the emission of a soft-gluon of relative transverse momentum k (l) t and relative rapidity η (l) with respect to leg l as The CAESAR resummation plugin to SHERPA hooks into the event generation framework, facilitating the process management, and providing access to the COMIX matrix-element generator [87], as well as phase-space integration and event-analysis functionalities. The SHERPA framework is also used to compile all the required higher-order tree-level and one-loop calculations. For the NLO QCD computations we use the SHERPA implementation of the Catani-Seymour dipole subtraction [88] and the interfaces to the RECOLA [66,67] and OPENLOOPS [89] one-loop amplitude generators. The plugin implements the building blocks of the CAESAR master formula Eq. (4.1), along with the necessary expansion in α s used in the matching with fixed-order calculations. The building blocks are evaluated fully differentially for each Born-level configuration B δ of a given flavour and momentum configuration. 6 In Ref. [86] the CAESAR formalism and the corresponding implementation in the SHERPA plugin were extended to include the phase-space constraints given by SoftDrop grooming with general parameters z cut and β. Note that, as already pointed out in [21], differential distributions computed from Eq. (4.1) at NLL are discontinuous at L = L z = − ln z cut . A way to cure this behaviour was proposed in Ref. [84]. However, as this discontinuity is a subleading (NNLL) effect, we decided to leave it untouched and checked that the two approaches were in agreement within our theoretical uncertainties. Furthermore, we note that our resummation is strictly valid in the limit of small z cut , i.e. λ α z cut 1. However, for β = 0 SoftDrop , finite-z cut corrections are already present at the (leading) single-logarithmic accuracy [19]. For the specific choice of z cut = 0.1 adopted in this paper, these corrections have been studied in [49] in the context of the groomed jet mass in dijet processes, and found to have a negligible effect around one percent. The treatment of the kinematic endpoint is implemented in the same way as in Ref. [86] by shifting the relevant logarithms and adding power-suppressed terms to achieve a cumulative distribution that approaches one at the kinematic endpoint and has a smooth derivative approaching zero. In particular, we introduce the additional parameters p, x L , v max and modify all logarithms by power-suppressed terms according to Here we set v max to the numerically determined endpoint of the NLO distribution, and by default use p = 1.
We fix the renormalisation and factorisation scale to the transverse momentum of the Z boson, p T,µ + µ − , and the resummation scale to µ Q = p T,µ + µ − R 0 . Note, for the Born configuration, on which the resummation is performed, it holds p T,µ + µ − ≡ p T,jet . To evaluate the perturbative uncertainties of our results, we vary µ R and µ F according to a 7-point variation, (1,2), (2, 1), (2, 2)}, simultaneously in the fixed-order calculation and the resummation. The argument of α s in the resummation is always taken to be µ R , with a compensating term for the LL dependence to remain NLL accurate and ignoring the NNLL ambiguity this introduces in pure NLL terms. Additionally, we calculate the resummed distribution with The ingredients of the master formula Eq. (4.1) are readily available in the case of global event shapes, with or without SoftDrop grooming. However, some adaptations are needed if we want to apply this approach to the resummation of non-global jet angularity distributions measured on the leading jet in Z+jet events. These are described in Section 4.1 below. We then discuss aspects related to the matching to NLO fixed-order calculations in Section 4.2 and present numerical results in Section 4.3.

Jet angularities resummation within the SHERPA framework
Compared to the formalism described above and used in [39,85,86], some adjustments are necessary to the SHERPA resummation framework. In particular we have to take into account the fact that the angularities are sensitive to radiation within the hardest jet in the event only. Let us label, for definiteness, the initial state legs as l = 1, 2 and the measured final state leg as l = 3. As our observable is not sensitive to radiation collinear to the initial state legs, we are allowed to set P = 1 and R 1 = R 2 = 0 in Eq. (4.1). Keeping in mind that the Born phase space is given by final states consisting of a single parton plus the Z boson, only one radiator is left. In terms of Eq. (4.2), the angularity observables λ 1 α can be parametrised by where y 3 ≡ y jet denotes the rapidity of the hard jet. Note that for our choice of µ Q = p T,jet R 0 the last factor in g 3 d 3 equals unity.
-11 -A number of complications arises when considering the presence of a jet boundary. First, when originally computing the radiators, one conventionally divides the phase space for each dipole at η = 0. An arbitrary rapidity condition with respect to leg l, η l > η min,l , leads to an additional contribution where we have introduced Note that the last equality holds with α s = α s (µ R ) at NLL accuracy. We can use this to reflect the rapidity requirement given by the jet condition by setting η min,3 = ln(2 cosh(y 3 )/R 0 ). Here Q is an arbitrary scale, which we choose to identify with the partonic centreof-mass energy, that ultimately cancels in the final expressions. For global observables this cancellation happens with the soft function S, however, the restriction to radiation inside the jet additionally affects the structure of soft emissions. The details were worked out in [90] for the jet-mass observable. Those results apply to the general class of angularities as well, as they depend on the scaling of the observable with k t only, i.e. on the parameter a in Eq. (4.2). This implies (using a = 1) In a general way, the global contribution can be written as where c and H are the colour metric and hard function, respectively, see [39,85] for details of the notation. We can write the matrix Γ in colour space as a sum over dipoles In our case only three dipoles, namely 12, 13 and 23, contribute to the soft function. As a consequence, by exploiting colour conservation, all colour operators T i · T j can be written as linear combinations of the SU (3) Casimir invariants C F and C A and, consequently, the matrix structure of Eq. (4.10) disappears. The coefficients of the single-logarithmic function t, can be obtained by integrating the matrix elements for each dipole ij over the appropriate phase space. Because the collinear contributions have already been included in the radiators R l in Eq. (4.1), the I ij coefficients are purely due to soft wide-angle radiation. Consequently, they are functions of the jet radius R 0 only, and are independent of the angularity exponent α. In particular, they coincide with those computed in [90] for the jet mass: Note that by including the above dipoles we account for initial-state radiation at NLL accuracy, as an expansion in powers of the jet radius. Higher order terms in the jet radius expansion can be computed, but in addition to the suppression in R 2 0 the corresponding coefficient drops rapidly, and even for R 0 = 1 the first terms are an appropriate approximation [90].
non-global (t(L)). The two curves correspond to the case where the final-state parton is either a quark (solid red) or a gluon (dashed blue). The Sudakov is plotted as a function of the scaled variable t(L) (cf. Eq. (4.8)). For illustrative purpose, the top axes show the corresponding values of λ 1 α (irrespective of α) for two representative jet transverse momenta.
The non-global part S B δ non-global (t) is computed numerically, in the large-N c limit, following the algorithm highlighted in [17], as also done in [90]. 7 This can straightforwardly be applied to 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 )). We also note that the non-global contribution only depends on v (and the jet transverse momentum through α s (µ R )) but not on the angularity parameter α. Similarly, for the groomed case, the contribution from non-global logarithms only depends on z cut and not on β. The main reason behind this is that, at singlelogarithmic accuracy, non-global logarithms originate from configurations that are strongly ordered in energy with angles commensurate with the jet radius R 0 .
In practice, S non-global (t) is computed separately for each possible colour dipole, either incomingincoming or incoming-final. These contributions are independent of the jet rapidity meaning, in particular, that the two possible incoming-final configurations are equal. One can then reconstruct the two Born channels δ from the dipole contributions as done in [90]. This procedure allows one to recover the full-N c result at least for the first non-trivial term of the expansion in α s . We finally note that the numerical procedure described in [17] introduces an angular cut-off θ min to regulate the collinear divergence. We have performed 4 independent runs with θ min = {0.008, 0.004, 0.002, 0.001} and extrapolated the result to θ min → 0 using the same approach as in [91]. Fig. 3 shows the resulting non-global contributions S B δ non-global (t) for both (Born-level) quark and gluon jets. The upper set of axes indicate what value of t corresponds to the observable value λ 1 α for fixed p T,jet . We have checked that they agree within ∼ 1% with the results obtained in [90] at t 0.4.
However, the results presented here extend to larger values of t, as appropriate for the kinematic range investigated in this paper.
To gauge the numerical impact of the global and non-global soft function on our NLL predictions for the jet-angularity variables, we present in Fig. 4 results for the ungroomed case. As both S-function contributions contain terms expected to vanish in the R 0 → 0 limit, 8 we here show results for the two considered jet radii, R 0 = 0.4 and R 0 = 0.8. Besides the complete NLL distribution (in red), we depict results for S non-global = 1 (in blue), allowing us to quantify the numerical effect of non-global logarithms, and S global = S non-global = 1 (in green), effectively neglecting soft wide-angle emissions entirely. We observe that these contributions have a bigger impact for larger jet radii as theoretically anticipated. Furthermore, we observe a somewhat milder (relative) impact of S in the α = 0.5 case, when comparing to α = 2. For the latter a quite significant shift of the distribution is found. The case α = 1 lies in between the other two. This is understood from the fact that, for a given value of the angularity, the soft function is independent of α while the effect of the collinear contribution increases when decreasing α. When invoking grooming, the soft function remains unchanged for observable values larger than the transition point. However, below the transition point the softfunction contributions become constant, given by S non-global (z cut ) and S non-global (z cut )S global (z cut ), respectively. This is illustrated and confirmed in Fig. 15 in App. A.1. We have checked explicitly that our observations on the impact of the soft function both for the groomed and ungroomed case also hold after matching to NLO, as described in what follows.
We close this discussion of the NLL resummation, by noting, had we defined jets with algorithms that differ from anti-k t , for instance, C/A or k t , a new class of so-called clustering logarithms would have appeared at NLL, see e.g. [92][93][94][95][96].

Matching to NLO and achieving NLL accuracy
In order to achieve a faithful description of the angularity distributions across their whole spectrum, we need to match our all-order results with fixed-order predictions, computed here at NLO accuracy. Our matching procedure will keep track of the jet flavour, so that we can obtain what is usually referred to as NLO + NLL accuracy.
The matching procedure is defined for the cumulative distribution (for simplicity, we use λ ≡ λ 1 α ) We also introduce the notation Σ δ res,fo,match (see also [85,86]) in order to denote the results for the cumulative distribution for the (family of) flavour channel(s) δ computed, with resummation, at fixed order, or matched between the two, respectively. In the following, the channel label is omitted in general expressions and we use the shorthand σ = Σ(1). We denote the expansion of any cumulative distribution Σ to order α 2 s relative to the pp → µ + µ − j Born process as In practice, at least for Σ (2) fo , we only calculate We employ a multiplicative matching scheme, in which, for every partonic channel δ, we have . (4.17) By this procedure we automatically include the correct coefficients We use the anti-k t algorithm with our choice of parameters for defining jets, and employ the flavour-k t algorithm from [97] (BSZ) to assign the flavour channels. However, the BSZ algorithm is not immediately applicable here, due to the non-global structure of the angularities. In particular, at any perturbative order there exist contributions with only a single jet constituent, such that the angularity variable equals zero. In principle we need to perform the resummation around multi-jet configurations. However, to achieve our target accuracy, it is sufficient to guarantee that configurations occuring at LO are dressed with the proper LL Sudakov factor. An explicit example is given in Fig. 5(b). Our matching scheme achieves this by multiplying these Z +jj events with the full NLL Σ δ res corresponding to the channel δ obtained after dropping the second jet, which will contain the correct LL factor. It Figure 5: Examples for partonic input configurations to the flavour assignment algorithm. Configuration (a) needs to be identified as a gluon jet in the collinear limit, whereas for (b) we need to perform the LL resummation for a quark jet. IR safety requires that configurations (c) and (d) in the limit where p g → 0 are identified as gluon and quark jet, respectively.
is hence crucial that the flavour channel of those contributions is defined by the single object inside the leading jet alone. At the same time, in the example shown in Fig. 5(a), the quarks need to be clustered in the collinear limit and the jet identified as a gluon jet. While at this stage we could just use the combined flavour of all anti-k t jet constituents, NLO contributions would spoil the IR safety of this procedure. Examples of relevant configurations are given in Fig. 5(c) and (d), where an IR-safe algorithm needs to ensure that the g → qq splitting is undone first in the soft gluon limit. Here we propose and employ the following algorithm to assign flavour to the leading anti-k t jet: 0. Start with the list O of all coloured final-state objects, containing particle four-momenta and flavour labels, and the beams B,B with their respective flavours.
1. Run the standard anti-k t algorithm with radius parameter R 0 on O, and obtain the objects in the leading, i.e. highest p T , jet J ⊂ O.  (a) If d BSZ = d BSZ ik , update O by removing i and k and adding a new object with momentum p i + p k and the combined flavour of objects i and k. We finally can define bland versions of this algorithm, vetoing clusterings that would lead to jets with multiple flavours by setting the corresponding distance measures to infinity. We use this bland variant and identify every event as having a leading jet that is either quark-or gluon-like.

If
Note that the only contribution that matters are the cross terms between C δ,(1) and the leading logarithms, which only depend on this flavour assignment and not on details of partons outside the jet. Therefore, we can sort any configuration with the same flavour assignment of the leading jet into the same channel. With this procedure we can ensure a correct assignment of the LO real corrections, while maintaining an infrared-safe definition of jet flavour also at NLO.

NLO + NLL validation and predictions
Splitting up the full fixed-order calculation into separate flavour channels in an infrared safe way, enables us to validate our results, in particular the effective inclusion of the C δ,(1) coefficients, for each partonic channel δ individually. We hence compare the leading-order contributions obtained respectively for the exact matrix element and for the expansion of the resummed calculation.
Similarly, at NLO, we compare where Σ δ exp NLO (λ) is the NLO expansion of the all-order result, while Σ δ exp NLO+C (λ) also includes the contribution from the C δ, (1) In Figs. 6 and 7 we present the differential distributions dσ/d log λ ≡ dΣ(λ)/d log λ in λ, separately for the gluon (Fig. 6) and quark (Fig. 7) channel, identified by the BSZ flavour-k t algorithm, for the representative jet-p T slice p T,jet ∈ [408, 1500] GeV. We have checked that analogous results are obtained for the other transverse momentum slices. In all these validations we assume a jet radius of R 0 = 0.8. As before, we consider the angularities λ 1 α with α ∈ {1/2, 1, 2} either without grooming or with SoftDrop grooming using z cut = 0.1, β = 0. For each observable we include an additional panel that contains the differences between the fixed-order result, i.e. Eqs. (4.20) and (4.22), and the expansions of the resummation up to order α s and α 2 s , i.e. Eqs. (4.21) and (4.23), (4.24), respectively. We first note that for all angularities -groomed or ungroomed -the difference between the derivatives of Σ LO and Σ exp LO vanishes for λ 1, as expected in the ungroomed case. For the groomed angularities there are in principle, for β = 0, deviations of O(z cut ). Those appear to be too small to observe numerically here, as was also observed for example for event shapes in [86]. In comparing dΣ NLO /dL and dΣ exp NLO /dL (with L ≡ log(1/λ)) we observe a linearly rising difference, indicating missing terms of order α 2 s L 2 in Σ exp NLO . By including the C δ,(1) coefficient in Σ exp NLO+C , the derivatives only differ by a constant, confirming that missing terms in the cumulative distribution are reduced to order α 2 s L. Having validated our resummation calculations by comparing their expansions to the corresponding fixed-order predictions, we finally present our matched resummed NLO + NLL results and compare them to their NLO counterparts. Figs. 8 and 9 contain our predictions for the two considered transverse-momentum slices, using the same set of observables and grooming parameters as above. These figures also illustrate how the full result is obtained as the sum of the two flavour channels -quarks and gluons -identified by our BSZ flavour-assignment procedure. We note that, for all observables, the contribution from gluon jets is increased for the higher p T,jet slice compared to the case p T,jet ∈ [120, 150] GeV. However, for both slices we observe a more significant contribution from gluon jets for larger values of the angularities, while the low-λ 1 α tails are entirely dominated by quark jets. This is a confirmation that a cut on the jet angularity can serve as a theoretically well-defined and IRC safe quark-gluon discriminant, as pointed out, for instance, in Refs. [29,[34][35][36]. We leave further investigation on this topic to future work. In addition to the matched resummation, we show the full fixed-order predictions at NLO accuracy. All-order effects turn out to be important essentially over the entire observable range. In particular, only for the highest values of the observables do NLO and matched predictions start to look similar. This appears to be the case for both ungroomed and groomed distributions. After a detailed analysis, we have concluded that this effect, the size of which is rather surprising, is driven by the large constant contribution of Eq. (4.18), which in turn originate from the large perturbative corrections that characterise the Z+jets process. Similar observations have been reported for the groomed jet mass in Zj production in Ref. [98].
-19 -     In this section we provide our final theoretical predictions for the leading-jet angularities. We start the discussion with a comparison of our NLO + NLL results with parton-level MC predictions from the PYTHIA, HERWIG and SHERPA event generators. We then discuss the size of non-perturbative corrections due to hadronisation and UE, and demonstrate how to account for these effects in our final NLO + NLL predictions.

NLO + NLL at parton level
Let us start by comparing our NLO+NLL results for the jet angularities against different parton-level MC predictions. As in Section 3 we consider the LHA (λ 1 1/2 ), Width (λ 1 1 ), and Thrust (λ 1 2 ) observables. In Figs. 10 and 11 we compare the NLO + NLL distributions against the parton-level MEPS@NLO predictions of SHERPA and the LO results from PYTHIA and HERWIG. As before we consider two p T,jet slices, namely p T,jet ∈ [120, 150] GeV in Fig. 10 and p T,jet ∈ [408, 1500] GeV in Fig. 11. The theoretical uncertainty for the resummed and matched predictions is obtained by varying renormalisation (µ R ) and factorisation (µ F ) scales, as well as the parameter x L that determines the resummation scale, as detailed in Sec. 4, while the uncertainty band of the SHERPA MEPS@NLO predictions corresponds to variations of µ R and µ F only, however, both in the matrix element and parton-shower component. For this reason, the former yields a more conservative assessment of the uncertainty than the latter.
A useful viewpoint to understand our results is to consider the different types of contributions that affect the distributions. Roughly speaking, even at parton level, the behaviour of the distributions at low values of the observables strongly depends on the details of the treatment of radiation in the infrared region, e.g. the parton-shower cutoff or the Landau pole in the resummation. This region is characterised by large non-perturbative corrections, mostly due to hadronisation. The size of the region depends on the angularity exponent α, the presence of grooming and the considered transverse momentum region. Following the analysis of, e.g., Refs. [19,21] we can estimate the region of large non-perturbative corrections to be ungroomed jets or SoftDrop jets with α ≤ 1: λ 1 α µ NP p T,jet R 0 min(α,1) , where µ NP is a non-perturbative scale that we take to be 1 GeV. For α ≤ 1, the non-perturbative corrections first come from the region of hard-collinear radiation which is mostly unaffected by the grooming procedure. Conversely, when applying SoftDrop for angularities with α > 1, the factor in the second bracket of Eq. (5.2) is smaller than 1 and the boundary of the non-perturbative region is pushed towards smaller observable values. 9 Above the thresholds identified by Eqs. (5.1), (5.2), we expect (resummed) perturbation theory to provide a good description of the underlying physical processes. In this second region we can expect our NLO + NLL predictions and the MC ones to agree best. Finally, we can identify a third region, which is characterised by very large values of the observables, close to the kinematic endpoint of the distribution. In our resummed and matched prediction this region is under the jurisdiction of the NLO calculation, however the actual value of the kinematic endpoint is strongly affected by the presence of additional emissions, as captured by the parton-shower simulations. For the same reason, we also expect this very last bin to be rather sensitive to non-perturbative contributions, especially the UE, as we shall discuss below. We can now move to analyse in more detail the various cases. First of all we consider the ungroomed λ 1 1/2 distribution for the lower transverse momentum bin in Fig. 10. From the qualitative argument of Eq. (5.1), we expect the boundary between the first two regions to be located around λ 1 1/2 0.1. Indeed, in the region to the left of this value the agreement between the NLO + NLL and the three parton-level MC predictions is rather poor, with the NLO + NLL and SHERPA MEPS@NLO uncertainty bands not overlapping. Above this value, the agreement improves, as expected, until we end up in the third region, close to the endpoint, where the predictions from the two approaches are in strong disagreement. The situation does not improve with SoftDrop, as we have anticipated earlier.
If anything the agreement between the two predictions deteriorates at the boundary between the first and second region. We attribute this to the presence of transition-point effects (remember that we have chosen z cut = 0.1). However, as these are formally NNLL effects, they are not very-well modelled by our calculation. In the case of the higher transverse momentum bin, the boundary between the first two regions is pushed to smaller values of the observable (λ 1 1/2 0.05) and, consequently, the region where we find agreement between resummation and MC simulations should widen. This appears clearly in the top-left plot of Fig. 11, which represents the ungroomed case. In fact, the NLO + NLL result is in very good agreement with the PYTHIA prediction. Changes in behaviour of the groomed distribution with respect to the lower-p T case are instead minimal. In both cases, we still see large deviations in the endpoint region. Therefore, we must conclude that the distribution of the LHA observable is not well theoretically controlled in our NLO + NLL calculation. Since the different parton-shower simulations agree among themselves much better, a more systematic study of the differences with respect to the analytic approach seems well motivated. In this context we note that Ref. [99] presented such comparison for the closely related BKS observables [45,100] and fractional-energy correlations [47] in e + e − collisions, in an approach that could be extended to groomed observables, also in hadronic collision. In any case, as we shall see shortly, the other angularities considered in this study do not suffer from this problem to the same extend. Now let us examine the remaining two cases, i.e. the λ 1 1 and λ 1 2 distributions. 10 Similarly to the λ 1 1/2 case it is convenient to consider the three different λ regions. Let us concentrate on the p T,jet ∈ [120, 150] GeV slice in Fig. 10 first. We observe differences between our NLO + NLL results and MC predictions in the first region for the ungroomed distributions, i.e. λ 1 α 0.01. However, we notice that grooming now significantly improves the overall agreement between the different predictions, especially for the λ 1 2 case, as expected from Eq. (5.2). Differences between NLO + NLL and the parton-shower simulations in the endpoint region are also seen in this case, but this time all predictions remain consistent within the estimated uncertainties. We also note the abrupt change of behaviour in the SoftDrop distributions around λ 1 α z cut = 0.1, which marks the boundary between the groomed and ungroomed regions. As already explained, this discontinuity is beyond the accuracy of our resummation and it is therefore reassuring that the theoretical uncertainty in this region is rather large. Finally, we note that the overall agreement between our NLO + NLL calculation and the parton showers further improves for high transverse momentum jets, evidenced by Fig. 11. Thus, unlike the case of the LHA observable, the NLO+NLL predictions for jet width and jet 10 Note that, unlike in the λ 1 1/2 case, the binning for ungroomed and groomed distributions is chosen in a different way. More precisely, the groomed distributions have finer bin-spacing than the ungroomed variants. Also note that we use different binning for the λ 1 1 and the λ 1 2 distributions.
-25 -thrust in general agree with the corresponding results from parton-shower simulations. This allows us to conclude that the λ 1 1 and λ 1 2 distributions are well under theoretical control. It is interesting to point out that, in the context of quark-gluon discrimination, the ATLAS collaboration has also observed [101] a worsening of the agreement between the experimental data and Monte Carlo simulations, for smaller values of α although with large uncertainties.

Impact of non-perturbative corrections
Both the NLO + NLL and the MC predictions in Section 5.1 were given without including nonperturbative corrections due to the UE and hadronisation. In this section we provide an estimate of these contributions and demonstrate how one can introduce corresponding corrections to our NLO + NLL results.
Let us consider the parton-to-hadron transition first. One way to add hadronisation effects to our NLO + NLL results is to perform a field-theoretical analysis, based on non-perturbative matrix elements. However, the jet structure may not only be significantly affected by hadronisation but also by the UE. The underlying event, unlike hadronisation, should not be seen as a non-perturbative correction to our Zj matrix elements, but rather as independent semi-hard processes, contributing with additional final-state partons that ultimately hadronise. Therefore, the UE cannot be straightforwardly included into the resummation framework beyond the simple model of uniform radiation [102,103]. Moreover, hadronisation may lead to non-trivial clustering of partons from different processes, i.e. the hard scattering and the UE, into final-state hadrons. Therefore, we prefer to rely on dynamical nonperturbative models as implemented in the PYTHIA, HERWIG and SHERPA MC event generators [52]. Each of these programs is using its own default combination of non-perturbative models for hadronisation and UE. More precisely, in both the PYTHIA and SHERPA frameworks the UE is described by the Poisson-based Sjöstrand-Zijl model presented in [74,104,105] whereas in HERWIG the UE is simulated based on the eikonal-model for soft emissions [106,107]. The hadronisation effects are simulated according to the Lund string model [108,109] in PYTHIA and according to cluster-fragmentation models both in HERWIG [110,111] and SHERPA [73].
In order to estimate the overall impact of non-perturbative effects we consider the ratios of hadronto-parton-level predictions. To this end, each MC simulation with PYTHIA, HERWIG and SHERPA is performed twice: first at parton level (without hadronisation and UE) and second at hadron level (with both UE and hadronisation switched on). This approach allows us to consider the ratio dσ HL /dλ / dσ PL /dλ evaluated per λ-bin for each MC generator. The dσ HL /dλ / dσ PL /dλ ratios obtained with PYTHIA, HERWIG and SHERPA are then combined into an envelope, with extreme points given by the maximal and minimal predictions per given angularity bin. The central value of the envelope is determined by the respective arithmetic mean value. The resulting ratios, together with their uncertainties are shown in Figs. 12 (a) and (b), differential in λ 1 α , for the two representative p T,jet bins, respectively. Experimental studies of jet structure may be based on all final-state particles, i.e. jet constituents, or purely on the tracks left by charged hadrons only. Therefore, it is convenient to consider two types of dσ HL /dλ / dσ PL /dλ ratios: the first one where all final-state hadrons in the jet are considered (red solid areas in Fig. 12) and the second one where only the charged jet constituents are used in the evaluation of the angularities (blue hatched areas). We note that the shape and the size of non-perturbative corrections in the two cases are rather similar. By comparing Figs. 12 (a) and (b), we observe that the behaviour of the non-perturbative correction factors obtained here feature all the expected properties. The overall size of the non-perturbative corrections decreases as the transverse momentum increases, in line with our expectation for IRC safe observables. This remains true for the charged-hadrons only case, despite IRC unsafety. Applying SoftDrop typically reduces the size and the onset of non-perturbative corrections, although this feature is less prominent in the LHA case, for which non-perturbative corrections remain rather large also for groomed jets.
Finally, we use the results in Fig. 12 to correct our NLO + NLL predictions for non-perturbative effects. In order to do that we take our NLO + NLL distributions as in Figs. 10 and 11 and multiply them on a per-bin basis by the corresponding central value of the dσ HL /dλ / dσ PL /dλ ratios shown in Fig. 12. The final uncertainty bands are obtained by summing in quadrature the perturbative and non-perturbative uncertainties.
The NLO + NLL + NP distributions, represented by black solid lines, in Figs. 13 and 14 present the main result of this paper. They feature our NLO + NLL perturbative predictions and include an MC-based estimate of non-perturbative corrections, here for the case where all hadrons in the jet are considered, for the three angularities λ 1 α , for the two representative transverse momentum bins. Corresponding results for the case where only charged tracks are considered in the observable calculation are presented in Appendix A.2 in Figs. 16 and 17. The gray uncertainty bands represent the perturbative, (µ R , µ F , x L ) variations, and non-perturbative, δ NP , systematics added in quadrature. For comparison, we also report the SHERPA MEPS@NLO hadron-level predictions, red dashed-dotted lines, including systematic variations of the factorisation and renormalisation scale, red hatched band. We see a good agreement between the two predictions for the jet Width and Thrust. We expect that a comparison to upcoming experimental data would find agreement with these distributions. However, the size of the experimental uncertainties will tell us whether the theoretical calculation needs to be improved, in order to perform precision phenomenology. In this context, these observables could be exploited to extract jet properties, such as the flavour, or, even more ambitiously, parameters of the Standard Model such as couplings and masses. The situation for the LHA distribution is instead in stark contrast. Here the resummed calculation and event-generator predictions are not in agreement, especially in the groomed case. This is true both at parton-and hadron-level. In this case upcoming LHC data may help us to shed light on the discrepancy. The data may favour the MC prediction, indicating that modelling used in the resummed calculation, which for instance neglects recoil, is not appropriate. On the other hand, the data may favour the resummation, indicating the need for a better understanding of the logarithmic structure that is achieved by the parton-shower models.
We conclude this section by reporting that most of the conclusions reached above, also hold for smaller jet radii. Explicit resummed and matched results for R 0 = 0.4 , as well as their comparison to SHERPA MEPS@NLO predictions, are collected in Appendix A.3.
-27 -         -32 - In this paper we have performed a detailed phenomenological analysis of variables that describe the internal structure of jets produced in association with a Z boson decaying into muons. Our study has focussed on IRC safe jet angularities that are characterised by an angular exponent α, the variation of which allows us to probe QCD radiation in different ways. In particular, inspired by an upcoming measurement by the CMS collaboration, we have derived results for α = 1/2, 1, 2. Furthermore, we have considered standard (ungroomed) and groomed jets. Our grooming algorithm of choice has been SoftDrop, which is characterised by the momentum fraction parameter z cut and the angular exponent β. We have shown explicit results for the combination z cut = 0.1, β = 0, corresponding to the default choice of the CMS collaboration in their ongoing study.
Comparisons between unfolded experimental measurements and accurate first-principle theoretical predictions can allow us to test, and improve, our description of jet angularities. These observables are interesting for many reasons. For instance, they can be employed to distinguish quark-like from gluonlike jets, in a theoretically well-defined way [35][36][37]. They can also be used to extract Standard Model parameters, such as the strong coupling constant [38]. Furthermore, jet-angularity measurements can form an important input to constrain Monte Carlo event generators in general and the tuning of non-perturbative model parameters in particular [37].
The first part of the paper has been devoted to a detailed comparison of different MC event generators that provide an increasing sophistication in the way they include fixed-order matrix elements. In the region of small angularities, where soft radiation, both of perturbative and non-perturbative origin, dominates, we have confirmed that the spread in the predictions provided by widely used general-purpose event generators is largely reduced when grooming is considered. In our study we have also investigated the impact of fixed-order corrections to the parton shower. In particular, we have compared LO predictions from PYTHIA and HERWIG, to the ones obtained with merged samples from SHERPA, both at LO and NLO. The inclusion of full NLO QCD corrections in the SHERPA MEPS@NLO simulations leads to a significant reduction of systematic uncertainties, estimated by consistent variations of the factorisation and renormalisation scales in both the matrix element and parton-shower components. When comparing to the LO simulations, we have observed noticeable effects at large values of the angularities, ranging from ∼ 10% at moderate transverse momentum to ∼ 30% in the high-p T,jet region. This applies to both, standard and groomed jets. This behaviour is perhaps related to an analogous trend observed in the inclusive K-factors, which are increasing functions of the jet transverse momentum. Our findings underline the importance of using state-of-the-art generators for jet phenomenology, especially when considering processes characterised by fairly large radiative corrections, as the one considered here, i.e. hadronic Zj production.
The second part of this study has instead focussed on all-order predictions obtained with resummed perturbation theory at single-logarithmic accuracy. Thanks to a flavour-dependent matching to fixedorder predictions, we have been able to achieve NLO+NLL accuracy. This includes, for the ungroomed case, the resummation of non-global logarithms, which we perform in the large-N c limit. The resummed and matched calculations are implemented in the resummation plugin to SHERPA. This allows us to exploit the main features of the event generation framework, including phase-space integration, matrix-element evaluations and matching, obtaining NLO + NLL predictions that include fiducial cuts and can be faithfully compared to unfolded measurements. To this end, we complemented our perturbative predictions with non-perturbative corrections, derived from MC simulations, both for angularity measurements based on all or charged-only jet constituents. As expected, groomed jets benefit from smaller corrections, which however remain rather large for the λ 1 1/2 angularity. Although -33 -we have focussed on one particular combination of SoftDrop parameters, our resummed calculation is more general. It applies for all β ≥ 0 values and can be extended to negative β if SoftDrop is used in tagging rather than grooming mode. If one wanted to consider much smaller values of z cut , in addition the systematic resummation of related logarithmic corrections should be included. Similarly, for large values of z cut , power corrections can become important and should be included.
The main deliverable of this paper are NLO + NLL predictions for jet angularities in Z+jets events in proton-proton collisions that can be immediately exploited for data-theory comparisons. However, this work has already shown us the path towards higher precision that we need to pursue. In this context, we can identify four complementary lines of research. First of all, although we do not expect large corrections, one should improve the accuracy of the resummation, in order to reduce the theoretical uncertainty. The framework for NNLL calculations has been worked out using SCET, although explicit results only exist, to the best of our knowledge, for the case α = 2 [98,112]. Second, it would be desirable to match the resummation to NNLO distributions. However, we recognise that this endeavour is much more formidable as it requires calculating Z plus two partons at NNLO accuracy. Some first steps have however already been taken in this direction [113][114][115][116]. Next, given that the size of non-perturbative corrections can remain large even for groomed angularities, e.g. for α < 1, it would be interesting to study additional grooming strategies, such as complementing the SoftDrop procedure with a filtering step [117] or using the Recursive SoftDrop algorithm [118], although these tools would most likely further complicate the structure of the analytic resummation. Finally, a SCET framework to deal with non-perturbative effects has been recently developed [119,120] and it would be interesting to quantitatively compare its predictions to the phenomenological model employed here.

A Additional results
In this Appendix, we collect supplementary results.

A.1 Numerical effects of the global and non-global soft function for groomed observables
In Fig. 15 we show the same results as in Fig. 4 in the main text, but with grooming.    The comparison yields very similar findings to the R 0 = 0.8 case reported in the main text. However, for the groomed width in the lower transverse momentum slice we here observe very large non-perturbative corrections (with corresponding large uncertainties). This signals the breakdown of our approximate treatment for including non-perturbative corrections. Here we probe a region of phase space close to the parton-shower cutoffs, below which the Monte Carlo simulations are purely determined by non-perturbative effects and, consequently, the hadron-to-parton-level ratios get large. This effect becomes more important as we lower p T,jet and/or the jet radius R 0 . Its onset can be estimated using Eqs. (