Jet mass distribution in Higgs/vector boson + jet events at hadron colliders with $k_t$ clustering

We address the issues of clustering and non-global logarithms for jet shapes in the process of production of a Higgs/vector boson associated with a single hard jet at hadron colliders. We perform an analytical fixed-order calculation up to second order in the coupling as well as an all-orders estimation for the specific invariant mass distribution of the highest-$p_t$ jet, for various jet algorithms. Our results are derived in the eikonal (soft) limit and are valid up to next-to-leading logarithmic accuracy. We perform a matching of the resummed distribution to next-to-leading order results from MCFM and compare our findings with the outputs of the Monte Carlo event generators Pythia 8 and Herwig 7. After accounting for non-perturbative effects we compare our results with available experimental data from the CMS collaboration for the Z + jet production. We find good agreement over a wide range of the observable.


Introduction
The invariant mass of a jet is a typical example of a jet shape that plays an important role in the study of the substructure of jets, testing QCD, and identifying newphysics signals. Being sensitive to soft and/or collinear emissions from the parton initiating the jet and from the other incoming and outgoing partons, this observable provides an indispensable mean for probing various aspects that are relevant to achieving better accuracy in QCD calculations. Examples of such aspects include, on the non-perturbative side, hadronisation corrections, underlying event, pile-up interactions, and on the perturbative side, initial and final-state radiation, colour flow, resummation of large logarithms, etc. Analytical calculations for these aspects pave the way for a deeper insight into QCD processes, a better control of theoretical uncertainties, and a precise quantification of missing higher-order contributions and their significance, all of which are issues not very clear in Monte Carlo event generators.
In this paper we shed light on the resummation of large logarithms that arise due to a miscancellation of soft a e-mail: naima.ziani@univ-batna.dz b e-mail: kamel.kkhelifa@iu.edu.sa c e-mail: yazid.delenda@univ-batna.dz (corresponding author) and collinear singularities between real emissions and their corresponding virtual corrections. The convergence of the perturbative series, in the invariant jet mass (m j ) distribution, is spoiled by the presence of large logarithms in the ratio of the jet mass and its transverse momentum p t , L = ln(m j /p t ), at each order in perturbation theory. In the exponent of the integrated distribution, these logarithms take the form α n s L m , with α s being the strong coupling constant and m ≤ (n + 1), and thus they require an all-orders resummation. A next-to-leading logarithmic (NLL) resummation ensures that all single logarithms of the form α n s L n are resummed, in addition to the leading (double) logarithms (LL) α n s L n+1 .
The jet mass is a non-global observable, i.e., an exclusive observable that is sensitive only to gluon emissions which end up inside the jet. To ensure a proper NLL resummation then its distribution must carefully be treated for a class of large single logarithms known as non-global logarithms (NGLs), which are related to secondary non-Abelian emissions of soft gluons [1,2]. Furthermore, another type of large single logarithms known as clustering logarithms (CLs) [3,4], related to primary gluon emissions off the hard Born configuration, needs to be resummed when the jets are reconstructed using jet algorithms such as k t [5,6] and Cambridge-Aachen (C-A) [7,8]. The anti-k t clustering algorithm [9] is known to cause no CLs (see for instance refs. [10,11]). The full resummation of both NGLs and CLs has thus far proven to be a formidable challenge. The resummation of NGLs is usually achieved numerically via a Monte Carlo approach [1,2] in the large-N c limit (N c being the number of quark colours), though full-colour numerical resummation has been provided in refs. [12,13] based on an analogy between small-x BFKL resummation in Regge scattering and the Weigert equation [14]. Additionally, NGLs may also be resummed via an evolution equation known as the Banfi-Marchesini-Smye (BMS) equation [15] valid at large N c .
In this work we study the jet mass distribution in the process of production of a single jet associated with a vector boson (γ, Z or W ) or a Higgs boson H at the Large Hadron Collider (LHC). In ref. [11], the jet mass distribution was calculated at NLL accuracy combined with nextto-leading order (NLO) results in Z + jet and di-jet processes at hadron colliders, 1 for jets defined with the antik t clustering algorithm. The NGLs were computed therein analytically at fixed order (at O(α 2 s )) and numerically to all orders in the large-N c approximation. In the context of soft-collinear effective theory the jet mass distribution was also studied in ref. [16] for di-jet events, in ref. [17] for γ +jet events, and in ref. [18] for H + jet events. We elaborate herein on the work of ref. [11] by considering the jet mass distribution when jets are reconstructed using k t or C-A clustering algorithms. We additionally consider other vector bosons, namely γ and W , as well as Higgs boson + jet production processes. On the experimental side, the jet mass distribution in W/Z + jet events at the LHC was studied by the CMS collaboration [19], where the jets were reconstructed using various jet algorithms. Additional jet substructure techniques such as trimming, filtering, and pruning, were also addressed in the same work [19]. We do not address these techniques in the present paper.
We compute NGLs and CLs at fixed order, specifically at O(α 2 s ) where they first appear, for the invariant mass distribution of the highest-p t jet. We provide results for the following three jet algorithms: k t , C-A and anti-k t , where we note that for the latter algorithm NGLs were first computed in ref. [11] and that CLs are absent. Moreover, we approximate the all-orders resummed CLs and NGLs by an exponential of the O(α 2 s ) result in the case of k t and C-A algorithms. This is justified by the fact that for the anti-k t algorithm the said exponential approximates the all-orders numerical result very well as we shall demonstrate. 2 We then compare the NLL-resummed and NLO-matched result for the jet mass distribution, which includes the resummed global and non-global (NGLs and CLs) form factors convoluted with the Born cross-section 1 Note that our convention for LO, NLO, etc, is different from that in ref. [11]. In our convention, the LO differential distribution is proportional a delta function.
2 While the all-orders numerical resummation of NGLs for the anti-kt algorithm may be computed using the Monte Carlo code of ref. [1], as was done in ref. [11], we found that this code produces unreliable results for some dipoles in the case of kt clustering. Note that the C-A algorithm is not implemented in the code of ref. [1]. and corrected for NLO effects for each of the four V /H + jet processes, with results from Pythia 8 [20] and Herwig 7 [21,22] parton showers. Finally, we estimate the nonperturbative corrections to this distribution and compare our predictions with experimental data from the CMS collaboration [19] for the jet mass distribution in Z + jet events at the LHC. This paper is organised as follows. In section 2 we discuss kinematics of the processes under consideration and define our observable. We calculate, in section 3, the distribution of the jet mass at leading order and construct the resummed global form factor up to NLL accuracy in the exponent. In section 4 we compute the leading CLs at O(α 2 s ) for both k t and C-A clustering algorithms, which happen to give identical results at this particular order. We also calculate NGLs at O(α 2 s ) for the aforementioned jet algorithms in addition to the anti-k t . We are then able to assess the impact of the various clustering algorithms on NGLs. In section 5 we discuss the all-orders resummation of NGLs and CLs. In section 6 we compare our NLL-resummed result including NLO corrections for the jet mass distribution with the outputs of Pythia 8 and Herwig 7 parton showers. In section 7 we estimate the non-perturbative corrections, which include hadronisation corrections and the underlying event, on the distribution, and compare our results with the experimental data. Finally, in section 8, we draw our conclusions.

Processes and kinematics
In this paper we are interested in the calculation of both CLs and NGLs at single logarithmic accuracy, for the jet mass distribution in the process of production of a single jet associated with a vector (W/Z/γ) or Higgs boson at hadron colliders. For this purpose, it suffices to consider the eikonal (soft) approximation in the squared matrix elements for the emission of gluons. The emitted gluons are assumed to be strongly ordered in transverse momenta, i.e., k tn ≪ · · · ≪ k t2 ≪ k t1 ≪ p t , where k ti is the transverse momentum of the i th emission and p t is that of the outgoing hard jet. The latter ordering simplifies the calculations of the emission amplitudes while being sufficient for capturing the single logarithmic CLs and NGLs.
For a vector boson + one jet production in hadron collisions, there are three partonic channels that contribute to the Born process, namely: qq → g, qg → q, andqg →q. For W ± production, flavour changing needs to be taken into account at the Born level, but this does not affect the QCD structure of initial and final-state radiation. As for the Higgs + one jet process there are four partonic channels to be considered. These are: qg → qH,qg → qH, qq → Hg, and gg → Hg. As far as QCD calculations are concerned all mentioned channels, whether for Higgs or vector boson production, are in fact identical as they all involve three hard coloured (QCD) partons and a colour-neutral boson. This means that the resummation of the jet mass distribution is essentially identical in all of the said channels, with differences pertaining to just the Born cross-section and the associated colour factors for the various channels. We note that the relevant total cross-sections have been calculated up to next-to-next-toleading order (NNLO): Higgs + jet in refs. [23][24][25][26], Z + jet in refs. [27,28], W + jet in ref. [29], and γ + jet in ref. [30].
In the current work, we henceforth consider the three partonic channels shown in figure 1: (δ 1 ) : qq → g + X, (δ 2 ) : qg → q + X, and (δ 3 ) : gg → g + X, 3 where X refers to the colour-neutral boson (γ, Z, W ± or H). We label the incoming partons with (a) and (b) and the outgoing parton initiating the hard jet with (j). The four-momenta of the three hard Born partons and the emitted soft gluons are given by where η i and φ i are the rapidity and azimuth of the i th emission and y and ϕ are those of the outgoing hard jet, measured with respect to the beam axis. The incoming partons a and b carry momentum fractions x a and x b of the incoming protons, and √ s is the collision centreof-mass energy. We shall be ignoring recoil against soft emissions throughout, as it is beyond single logarithmic accuracy.

Jet mass observable and jet algorithms
We study the normalised (squared) invariant mass of the outgoing hard jet j defined by where the sum is over all emitted soft gluons which end up inside the hard jet after the application of a jet algorithm on the final state partons. Notice that we are considering massless quarks and that the soft approximation has been assumed in the above equation whereby p j · k i ≫ k ℓ · k m . The k t , C-A and anti-k t jet algorithms work as follows. For each pair (im) of hadrons in the final state one defines a distance and for each single hadron a beam distance for some fixed jet radius parameter R. Here the parameter p = 1, 0, −1 for k t , C-A, and anti-k t clustering, respectively. If the smallest of all of these distances is d im , then particles i and m are combined into a single particle with four-momentum p i + p m , whereas if the smallest is d i then particle i is considered as a jet and is removed from the list of particles. This procedure is iterated until one is left only with jets in the final state. For the k t algorithm, and in the regime of stronglyordered emissions, the clustering of particles starts with the softest real gluon. Then, in a given event this softest gluon is dragged towards the next-to-softest real parton within a circle of radius R in the (η, φ) plane. If no such harder parton exists then this softest gluon is considered as a jet and is removed from the list of partons. The process is then repeated until no particles are left. When clustering two partons together, the resulting pseudo-jet is essentially aligned along the direction of the harder, and its four-momentum is just that of the harder parton.
For the anti-k t algorithm, on the other hand, clustering starts with the hardest particle, and hence it works in an apposite way to k t clustering. For the C-A algorithm, only geometric distances between partons in the (η, φ) plane decide how clustering happens. Particles which are closest to each other get clustered first.

Jet mass distribution
In what follows we calculate at NLL accuracy the jet mass distribution for a given channel δ, defined by (following the notation of refs. [11,31]) where d 2 σ δ /dB δ d̺ is the differential cross-section with respect to both the Born configuration B δ and the jet mass observable ̺. Details of the differential Born configuration dB δ are discussed further in appendix A. The integrated jet mass distribution is obtained by integrating dΣ δ (ρ)/dB δ over B δ with some chosen kinematical cuts (which we denote by Ξ B ), and summing over all Born channels. That is Following ref. [11], we write eq. (4) in the region ρ ≪ 1 in the factorised form where dσ 0,δ /dB δ is the differential partonic Born crosssection for channel δ (see appendix A) and the factor C B,δ (ρ) depends on the Born kinematics and has the perturbative expansion where C (n) B,δ (ρ) are channel-dependent terms that correct the resummation for non-logarithmically-enhanced terms. The ρ-dependent function f B,δ (ρ) resums all the large logarithms. It has the form [31] where the function Lg 1 resums the leading (double) logarithms (LL) of the form α n s L n+1 ,g 2 resums next-toleading (single) logarithms (NLL) of the form α n s L n , and α sg3 resums next-to-next-to-leading logarithms (NNLL) of the form α n s L n−1 , and so on, with L = ln(R 2 /ρ). The LL functiong 1 receives contributions from soft-collinear emissions from the parton initiating the jet and depends on its colour Casimir scalar. The NLL functiong 2 receives contributions from various sources: (a) hard-collinear emissions from the outgoing hard parton, (b) soft wide-angle emissions from all hard partons, (c) starting at O(α 2 s ), NGLs from soft wide-angle correlated secondary emissions, and (d) CLs, when jet algorithms other than anti-k t are implemented for jet reconstruction, from soft wide-angle primary emissions off the hard partons. These again appear starting from O(α 2 s ). The whole functiong 1 and parts ofg 2 , namely contributions (a) and (b) stated above, have been determined in ref. [11] for the anti-k t algorithm. The exact same result also applies for the case of k t and C-A clustering as the effect of jet algorithms first appears at O(α 2 s ). Our task is to determine the other two contributions tog 2 , namely (c) NGLs and (d) CLs, for k t and C-A algorithms. Before doing so, we review in the next section the basic calculations that lead to the determination ofg 1 and contributions (a) and (b) ofg 2 .

.1 Fixed-order calculation
In this section we compute the jet mass distribution at leading order in QCD and present the all-orders resummed result. Our calculations are valid in the eikonal approximation and accurate up to NLL accuracy. First, we define the following antenna functions relevant for the squared matrix elements for the emission of soft gluons It is worth noting that these antennae are purely angular functions, i.e., they involve no energy or momentum dependence.
Consider the process of emission of a single soft gluon off the three-hard-legs Born configuration (abj), i.e., the process a+b → j+k 1 shown in figure 2. The corresponding factorised eikonal amplitude squared is given by with ∆ δ = {(ab), (aj), (bj)} denoting the three dipoles formed from the partons in channel δ. The colour factor C iℓ is defined as where T i are the generators of the SU(N c ) group with Casimir scalar given by for quarks (and anti-quarks) and T 2 i = C A = N c for gluons. Conservation of colour implies that for our leading-order process a + b → j + k 1 we have [32]: T a + T b + T j = 0, where the generators are taken as if all partons were incoming. Explicitly written, the colour factors relevant to our dipoles are: The term W R 1,δ is the eikonal amplitude squared for the emission of a real soft gluon in the partonic sub-process δ. The corresponding virtual correction in the eikonal limit is simply W V 1,δ = −W R 1,δ . Notice that we are adopting the notation used in our previous work on eikonal amplitudes for e + e − → di-jet process [32]. In our recent paper [33] we have generalised the latter to the case of hadron collisions, specifically considering three-hard-legs Born processes. The corresponding phase-space factor is given by whereᾱ s = α s /π = g 2 s /4π 2 , g s is the strong coupling, and ξ 1 = k t1 /p t . The running of the coupling is irrelevant at one loop and only becomes important at higher orders. The full resummation that we present later will include running-coupling effects.
Following the procedure of measurement operators (see for instance ref. [34]), we write the jet mass distribution at one loop as where the function Ξ in (k 1 ) ensures that the angular integration region for gluon k 1 is such that it gets clustered to the hard jet when the jet algorithm is applied. At this order all jet algorithms essentially work in the same manner, and Ξ in (k 1 ) is then a simple Heaviside step function; At higher loops, as we shall see, this is not as simple. Substituting the expression of the eikonal amplitude squared (10) into eq. (13) we obtain with the antenna function w 1 iℓ for each emitting dipole in ∆ δ given by . (15c) Note that the upper limit of the k t1 integral is the renormalisation scale µ R = p t , which translates into an upper limit 1 on ξ 1 . In order to perform the angular integrations we introduce the polar variables (r 1 , θ 1 ) such that and make a change of variables in the integration such that dη 1 dφ 1 = 1 2 R 2 dr 2 1 dθ 1 . One may expand the jet mass ̺ 1 defined in eq. (2) as a series in R as follows In fact, at single logarithmic accuracy it suffices to keep just the first term in this expansion, and thus we write the step function in eq. (14) as Θ ξ 1 R 2 r 2 1 − ρ . We now perform the integrations for each dipole.
-The dipole (ab): The contribution of the in-in dipole (ab) to eq. (14) at single logarithmic accuracy may be written as follows with L = ln(R 2 /ρ) being the large logarithm that we aim to resum. This contribution corresponds to soft wide-angle radiation from the in-in dipole into the interior of the measured outgoing jet, and is thus free from collinear logarithms.
For the in-jet dipole (aj) eq. (14) reads Note here that the step function Θ ξ 1 R 2 r 2 1 − ρ which restricts ξ 1 > ρ/(R 2 r 2 1 ) also implies that R 2 r 2 1 > ρ since ξ 1 < 1. This serves as a collinear regulator for the integral over r 1 , which would otherwise diverge, resulting in an overall double logarithm as well as a single logarithm. Evaluating the ξ 1 integration yields We perform the integration over θ 1 by expanding the integrand as a series in R and neglecting higher-order terms that have small coefficients. Thus we find The first term in this expansion corresponds to soft and collinear emissions from the outgoing hard leg (j) into its own jet. It contributes at the double logarithmic level, giving the result which is independent of the jet radius (other than in the argument of the logarithm). The other terms in the expansion (21) are purely soft wide-angle contributions, hence we can set ρ → 0 in the lower limit of integration over r 2 1 , and throw away the sub-leading ln r 2 1 term in the integrand. Performing the integration we obtain (23) We note that the coefficient of R 8 in this expression is vanishingly small (O(10 −7 )).
-The dipole (bj): For the other in-jet dipole (bj) the only differences relative to the dipole (aj) are the colour factor C aj → C bj and a minus sign to be inserted in the exponent of the exponential in the integrand of eq. (19), i.e., exp(R r 1 cos θ 1 ) → exp(−R r 1 cos θ 1 ). This is equivalent to a change R → −R (the rest of the integral is invariant under this change). This actually does not produce any differences in the integration since only even powers of R appear in the results (22) and (23).
We can therefore write the assembled soft-collinear double-logarithmic result as and the soft wide-angle single-logarithmic contribution as with This result was first derived in ref. [11], and it actually exponentiates to all orders. However, the running coupling, whose argument is the invariant transverse momentum κ 2 t1,(iℓ) = k 2 t1 /w 1 iℓ of the emission k 1 with respect to the emitting dipole (iℓ) [35], contributes at higher orders and modifies the single logarithmic contribution f where β 0 is the one-loop coefficient of the QCD β function. Accounting for the running coupling for the double logarithmic contribution f is more subtle. In fact the running coupling introduces additional single logarithmic components which depend on the renormalisation scheme. We discuss the all-orders resummed result in the following subsection.

Resummed global result
The full NLL-resummed global form factor f global B,δ (ρ) has been computed in ref. [11] (eqs. (3.3), (3.11) and appendix C therein). The interested reader is referred to the latter reference, together with ref. [31], for details. Here we only state its form, which is given by [11] f global with γ E the Euler-Mascheroni constant (γ E ≈ 0.577) and Γ denotes the Gamma function. The radiator R and its derivative with respect to L, R ′ , are presented in appendix B. We note that the global form factor is identical for all jet algorithms. We also note that the expression of f global B,δ (ρ) may be deduced from the general form presented in ref. [31] as we show in appendix B.
In the next section we treat the case of two-gluon emission where clustering and non-global logarithms first popup.

Two-gluon emission
In the eikonal approximation, the factorised squared amplitude for the emission of two real gluons k 1 and k 2 off the three-hard-legs Born configuration is given by [33] where the one-loop amplitude squared W R i,δ , which builds up the reducible part of the above two-gluon squared amplitude (first term on the right-hand side), is given in eq. (10), and the irreducible contribution W RR 12,δ reads The virtual corrections at this order are Following ref. [34], and implementing the measurementoperator method, we write the jet mass distribution at this order as with phase-space factor dΠ 12 Here, the first integral produces the primary-emission contribution, which contains CLs, and the second integral gives NGLs. The functions Ξ p , where p stands for primary, and Ξ NG , where NG stands for non-global, result from the application of the jet algorithm and restrict the angular integration regions for gluons k 1 and k 2 .

Clustering logarithms
In this subsection we focus on the primary-emission integral in eq. (32) and leave the treatment of the correlatedemission NGLs term to the next subsection. The primaryemission contribution may be split into two parts. The first is the global component which results from integrating both gluons within the measured jet region. This has, however, been accounted for by the all-orders resummed formula (28) discussed in the previous section, and will thus be skipped here. The second part is related to the way jet algorithms cluster gluons and results in large single logarithms that are referred to as clustering logarithms [3,4,36]. These logarithms are a result of miscancellation between real emissions and virtual corrections. The key point is that while real gluons may be dragged into/out of the jet by other real gluons and thus get clustered together, virtual gluons can neither drag nor get dragged. We note that CLs are totally absent when jets are reconstructed using the anti-k t algorithm. At two loops, the C-A and k t algorithms produce identical CLs, but they start to differ at higher orders as was shown in ref. [36].
To perform the first integral in eq. (32) we begin by simplifying the clustering function Ξ p (k 1 , k 2 ). To this end, we introduce the same change of variables as in eq. (16) such that (η 1 , φ 1 ) → (r 1 , θ 1 ) and (η 2 , φ 2 ) → (r 2 , θ 2 ). Note that the upper limit of r i is π/ (R | sin θ i |) since we have −π < φ i − ϕ < π. We then have for the k t clustering algorithm [4,36] where the algorithm distances d im have been defined in eq. (3). The first term exactly reproduces half the square of the one-loop result (14), i.e. 1/2! [f (1) B,δ ] 2 , and persists at higher orders as 1/n! [f (1) B,δ ] n . This signifies that the oneloop result simply exponentiates into the global form factor discussed before. It is the second term, the CLs term, that we shall focus on in the remainder of this subsection.
We write the CLs contribution at this order as follows The expressions inside the square brackets are the oneloop eikonal amplitudes squared (10) for gluons k 1 and k 2 , respectively. To single logarithmic accuracy the ξ integrations factor out from the rest of the integrals yielding the result L 2 /2!, and we are left with where where the first term represents contributions from independent dipoles, that is, each dipole consecutively emits softer gluons at each order independently of the other dipoles. This situation is analogous to that in e + e − annihilation to di-jet process (see for instance ref. [36]). The second term in eq. (36) represents contributions arising from the interference of dipoles in channel δ.
To carry out the integrations we expand the integrand as a power series in R and use the change of variable θ 1 − θ 2 → θ 1 for the angular integrations. We obtain the following result for the dipole-interference part. Notice that the interference term F 2,int is not symmetric under the interchange of the dipoles (aj) and (ab), or (bj) and (ab), as apposed to the dipoles (aj) and (bj). This stems from the fact that integrands such as w 1 aj w 2 ab and w 2 aj w 1 ab are not identical, though symmetric under (r 1 , θ 1 ) ↔ (r 2 , θ 2 ). Since the angular restrictions on k 1 and k 2 are not identical then the results one obtains for the two mentioned terms are different. This boils down to the effect of the k t algorithm which does not treat the two gluons symmetrically. Furthermore, independent and interference terms involving the in-in (ab) dipole vanish in the limit R → 0.
Substituting the results (37) back into eq. (36) we obtain the corresponding CLs coefficients for each channel.

They read
for channel qq → g + X, for channel qg → q + X, and for gg → g + H. We show in figure 3 a plot of the CLs coefficient 1 2! F δ 2 as a function of R for the various channels δ. We notice that gluon-initiated jets have larger CLs co- Fig. 3. Two-loops CLs coefficient as a function of jet radius R in the kt and C-A algorithms for the three channels. efficient than quark-initiated jets, mainly due to the corresponding colour factors (C A = 3 and C F = 4/3, respectively). These series expansions in R converge, and at small values of R it suffices to keep only the leading terms. At very small values of R we observe that This result for CLs obtained here for H/W/Z/γ + jet events at hadron colliders coincides with that found in refs. [10,36,37] for jet mass distribution in e + e − → di-jet events. It does, however, deviate from it as R increases due to initial-state radiation from the incoming partons. Inline with the findings of refs. [4,36] we expect the term (35) to simply exponentiate to all orders. Nonetheless, there will be new CLs terms at each order that are not captured by the latter exponential and that are highly non-trivial to deduce (see ref. [36]). Moreover, we expect that at higher orders the small-R limit of the CLs coefficient in H/V + jet events at the LHC will coincide with that in e + e − → di-jet events found in ref. [36].
In the next subsection we compute NGLs at two loops.

k t and C-A clustering algorithms
We now turn to the evaluation of the correlated secondaryemission contribution in eq. (32) for the k t and C-A clustering algorithms. To this end we write where the clustering function reads As before, the integration over ξ 1 and ξ 2 yields L 2 /2!, and we may write with NGLs coefficient Performing the integration, as in the previous subsection, we obtain the results for G (iℓ) 2 for each dipole as a series in R In terms of channels we have for channel qq → g + X, for channel qg → q + X, and for gg → g + H. Moreover, in the small-R limit we observe that The result for channel δ 2 is exactly the same small-R limit found in the case of jet shapes in e + e − → di-jet events (see for instance ref. [37]). Results for channels δ 1 and δ 3 in the limit R → 0 are also the same and differ from those for channel δ 2 only in the colour factor. In figure 4 we plot the NGLs coefficient 1 2! G δ 2 at this order as a function of jet radius R. Once again we notice that gluon-initiated jets have larger NGLs coefficient due to their large gluonemission colour factor (C A ). We observe from the plots in figures 3 and 4 that the CLs coefficient for the gg → g channel grows larger with R while that for NGLs does not change much.
Moreover, in order to assess the overall impact of CLs and NGLs at this order, we plot in figure 5 the combined coefficient of the single logarithmᾱ 2 s L 2 resulting from the non-global nature of our observable. We note that at large jet radii (R 1.0) and for all partonic channels the CLs and NGLs tend to balance each other out, but not entirely though. For small jet radii the said single logarithmic CLs + NGLs coefficient is quite large in magnitude especially for gluon-initiated jets.

Anti-k t clustering algorithm
For the sake of assessing the effect of clustering on NGLs, we report the results for the NGLs coefficient in the anti- k t algorithm. Note that there are no CLs in this case. The corresponding integral is identical to that in eq. (40) except for the clustering function. It reads The results we obtain for each dipole are In terms of channels we have G δ1,akt These results are in agreement with those reported in ref. [11]. Notice again that the R → 0 limit of the above expressions produces a result (which is proportional to 1.645 = ζ 2 ) that is identical to that reported in ref. [37] for e + e − → di-jet process. We plot in figure 6 the NGLs coefficient 1 2! G δ,akt 2 with anti-k t -clustered jets as a function of the jet radius R for the various partonic channels. As is clearly evident from the plots, NGLs in the anti-k t algorithm are much larger compared to those in the C-A or k t clustering case. This is made clearer in figure 7 where NGLs coefficients for each dipole are plotted for both k t and anti-k t algorithms. This observation was also made in previous studies of NGLs with k t clustering [4,36,38]. While k t clustering induces another tower of large single logarithms, namely CLs, it actually diminishes the impact of NGLs. Additionally, as we observed in the previous subsection 4.2.1, the induced CLs play a role of further reducing NGLs since their coefficients have opposite signs. This may hint at a (phenomenological) favour for the k t (or C-A) clustering algorithm over the anti-k t algorithm.

All-orders treatment of CLs and NGLs
Including the resummation of NGLs and CLs together with the global form factor (28) then the all-orders NLLresummed jet mass distribution may be cast into where S δ (ρ) and C δ (ρ) account for the resummation of NGLs and CLs, respectively. We note that, unlike the global form factor, the factors S δ and C δ are algorithmdependent.
In the anti-k t algorithm, the NGLs form factor results from multiple correlated gluons outside the jet that coherently emit the softest gluon into the jet. For the k t and C-A clustering algorithms, gluons can be moved into and out of the jet by the clustering, thus NGLs can be induced when more than one gluon is emitted within the jet region from an ensemble of harder gluons. The NGLs factor S δ can be computed numerically and in general only in the large-N c limit [1,15]. For the e + e − → di-jet process, finite-N c results do exist though [12,13]. Moreover, the CLs form factor results from multiple independent (primary) emissions that are clustered by the k t or C-A algorithm. Just like NGLs, the latter CLs can also be resummed numerically.
For the anti-k t algorithm, the all-orders numerical resummation of NGLs may be obtained from the dipoleevolution Monte Carlo code of ref. [1] as reported in ref. [11] for the various dipoles. We see from figure 8 that the Fig. 8. The full resummed differential jet mass distribution in the anti-kt jet algorithm with NGLs factor as an exponential of the two-loops result and as an all-orders numerical result obtained from ref. [11]. We explain the in the next section how these plots are obtained. exponential of the two-loops result (46) approximates very well the all-orders numerical result for the NGLs factor in the Z + jet process. The same is observed for the other processes. Hence we shall confine ourselves to simply using the exponential of the two-loops result for the k t and C-A algorithms. To this end we write, for a given channel δ, where G δ 2 , for the k t and C-A algorithms, is given in eq. (44), and the evolution parameter t is defined by Note that at fixed order t reduces to justᾱ s L.
As for CLs, it was shown in refs. [4,36] that the perturbative CLs series exhibits a pattern of exponentiation, and that the exponential of the two-loops result is a very good approximation to the numerically-resummed CLs factor obtained from the code of ref. [1]. Therefore, and just as we did with NGLs, we shall be using the exponential of the two-loops result for the CLs resummed factor C δ (ρ). Thence where F δ 2 is given for k t and C-A algorithms in eq. (38).

Comparison to Pythia 8 and Herwig 7 parton showers
In this section we present comparisons of our results for the jet mass distribution with those obtained from Pythia 8 [20] and Herwig 7 [21,22] parton showers (PS), where the jets are clustered with FastJet [39]. The resummed result is obtained by convoluting dΣ δ /dB δ given in eq. (6) with parton distribution functions (we use MSTW 2008 (NLO) PDFs [40] and µ F = µ R = 200 GeV). For doublechecking we perform the convolution using two different methods. In one method we simply use a Monte Carlo code to integrate over the momentum fractions of the partons x a and x b and over the transverse momentum p t and rapidity y of the jet, as explained in detail in appendix A.
In the other approach we generate a set of unweighted parton-level Born events using MadEvent from MadGraph [41,42] in the "Les Houches Event File" format [43], with the cuts Ξ B being applied. We then weigh each event by the resummed form factor S δ (ρ) C δ (ρ) f global B,δ (ρ), sum over all events, and divide by the effective luminosity L = N tot /σ 0 , with N tot the total number of events and σ 0 the Born cross-section calculated with MadGraph. This results in the integrated distribution given in eq. (5) from which the differential distribution can straightforwardly be obtained. To avoid low-p t resummation we impose a cut on p t of the final-state jet, e.g., p t > 200 GeV, i.e., we only consider high-p t jets, at a centre-of-mass energy √ s = 7 TeV.
In our resummed result we also include an approximation to the NLO effects on the distribution through the NLO factor C (1) B,δ (ρ). The full NLO distribution may ideally be analytically calculated using the full squared amplitude with two partons in the final state as well as virtual corrections to the Born cross-section. Though possible this is a delicate task. The alternative numerical approach would be to exploit fixed-order programs and obtain the factor C (1) B,δ (ρ) as a fully differential distribution in the Born configuration, and then perform the integration including the resummed form factor over the Born kinematics. Practically this is not feasible. Instead, one could obtain an NLO factor C (1) B,δ (ρ) that is averaged over the Born configuration [44] and insert it in eq. (6) as if it were unintegrated over dB δ . In this paper we employ this method and estimate the Born-configurationaveraged factor C (1) δ (ρ) as was done in refs. [11,44], using the NLO jet mass distribution obtained from the fixedorder program MCFM [45,46].
In refs. [11,44], the NLO factor C (1) δ (ρ) was calculated in the small-ρ limit as a constant, and then the ρ-dependence of the NLO contribution to the jet mass distribution was included at the stage of matching. This is equivalent to using the full ρ-dependence of C At NLO there are new channels that open up, specifically processes with incoming qq ′ or two gluons, that are not present at the Born level. These channels are not logarithmically enhanced and only contribute a small correction to the distribution. 4 In figure 9 we show plots for the differential jet mass distribution 1/σ dΣ/d √ ρ, where Σ(ρ) is defined in eq. (5), in Z + jet events at the LHC with k t clustering. We choose two values for the jet radius, one for which the size of NGLs + CLs is expected to be small, R = 1.0, and another where NGLs + CLs are expected to be important, e.g., R = 0.6. The global and pure-resummed distributions are normalised to the Born cross-section, while the resummed + C (1) , Pythia 8, and Herwig 7 distributions are normalised to the total cross-section.
We observe from the R = 1.0 plot in figure 9 that the global and full-resummed distributions are quite close to the Pythia 8 PS result, indicating the smallness of the effect of NGLs and CLs factors in this case. For the R = 0.6 plot, there is a clear difference between the global and Pythia 8 PS curves, and our full resummation, which is based on the exponential of the two-loops NGLs and CLs result, seems to do better. We also note that the NLO term C (1) slightly modifies the peak and tail of the distribution especially for R = 0.6, bringing it even closer to the Pythia 8 PS result.
As is clear from the plots, the Pythia 8 PS result seems to be in better agreement with our resummed distribution near the peak than Herwig 7. This observation was also made in ref. [11]. It should be noted, however, that a more comprehensive comparison is feasible only when one includes non-perturbative effects, where different event generators are then expected to be in agreement. We do this in the next section.
In figure 10 we plot the same distribution employing the C-A algorithm. We recall that up to two loops both k t and C-A algorithms produce identical results. This means that the resummed formula that includes the exponential of the two-loops NGLs and CLs as well as the C (1) δ terms for all channels are the same in both algorithms. We expect, however, differences between the two cases when one performs an all-orders NGLs and CLs resummation, and also when one includes the higher-order C (n) terms. We compare, in figure 10, the resummed + C (1) result with the Pythia 8 PS result employing both k t and C-A algo-rithms, where we notice that the peak of the distribution is slightly higher in the latter algorithm. We additionally show in figure 11 the differential jet mass distribution in the process gg → Hg for jet radii R = 1.0 and R = 0.6. Notice from figures 5 and 6 that, for this channel, the combined effect of NGLs and CLs at two loops in the k t algorithm is small compared to that of NGLs with anti-k t clustering, and so we expect our resummed distribution to fit well with the PS result in the case of k t clustering. This is indeed the case as is clear from figure 11, particularly for Pythia 8 PS.
Finally, in figure 12, we plot the resummed differential jet mass distribution in the processes W + jet and γ + jet at the LHC with k t clustering and R = 1.0. Our results are in general in good agreement with Pythia 8 results particularly near the peak of the differential distribution. The discrepancy between the results of Pythia 8 and Herwig 7 may be lifted when nonperturbative effects are included, as we do in the next section.

Matching to fixed order
Before we end this section we discuss the matching of the resummed result with the NLO fixed-order distribution. In fact, including the constant term C (1) δ (ρ) in eq. (6), the expansion of the resummed distribution now agrees with the fixed-order result over the entire range of ρ, except for the small correction due to the missing channels at the Born level (specifically the channel with incoming gg).
Additionally, as was shown in ref. [11], the NLO distribution has a kinematical end point of ρ max = tan 2 (R/2), which the resummed distribution does not have. In order to match the resummed distribution to the NLO result, specifically at the end point, we introduce the following change of the large logarithm [10] such that the large logarithm L ′ vanishes when ρ → ρ max , and L ′ → L when ρ → 0. We then use the simple matching formula where now both Σ NLL and Σ NLL,αs include the C (1) term. The subtracted term Σ NLL,αs cancels both the large logarithms and the C (1) terms in Σ NLO (ρ), leaving only corrections due to the channels missing at the Born level. We show in figure 13 a plot of the matched differential jet mass distribution compared to the fixed-order result from MCFM for the Z + jet process at the LHC. In this plot the resummed curve is plotted with the standard definition of the large logarithm (L = ln(R 2 /ρ)), and thus does not posses the end-point character, while the matched curve does have an end point exactly as in the MCFM curve. We note from this figure that the matched curve coincides with the resummed curve at small ρ, indicating a perfect cancellation of the large logarithms between the expanded result Σ NLL,αs and the fixed-order MCFM result Σ NLO .

Comparison to CMS data
In order to compare our results with the experimental data we first need to account for non-perturbative effects from hadronisation corrections and the underlying event. One commonly used numerical approach to extract these corrections is to compute the ratio of the results obtained from Monte Carlo event generators with nonperturbative effects switched on and off. In this paper we include these corrections analytically by considering the mean value of the change in the jet mass δm 2 j due to these non-perturbative effects. This change was computed in ref. [47] to be with and Here µ I is an arbitrary matching scale (chosen to be of order of a few GeV) and α 0 is the averaged coupling over the non-perturbative low-k t region, α 0 = 1/µ I µI 0 α s (k t ) dk t . The A(µ I ) is rescaled by the so-called Milan factor (M = 1.49 for anti-k t clustering and M = 1.01 for k t clustering [48]) to account for gluon decay. The constant K is defined in the appendix.
Non-perturbative effects are dominated by the contributions of the dipoles involving the outgoing jet, which scale like O(R), and which account for hadronisation corrections, while the smaller O(R 4 ) contributions from the incoming legs account for the underlying event. Since the mean value of δm 2 j depends both on the Born channel and kinematics, then we perform the shift on the mass of the jet on an event-by-event basis, that is we make the change m 2 j → m 2 j − δm 2 j in the resummed form factor and then perform the convolution. Furthermore, we shift the terms C (1) (ρ) accordingly.
We compare, in figure 14, the NLL+NLO resummed result (with the C (1) term), including non-perturbative corrections, with experimental data from the CMS collaboration [19,49] (obtained with integrated luminosity L = 5 fb −5 ), in the Z + jet process at the LHC with anti-k t clustering and R = 0.7. We also include in this figure the Monte Carlo results obtained from interfacing MadGraph with Pythia 8 [50,51] and Herwig 7 including hadronisation corrections and the underlying event. The plots in this figure are for the un-normalised jet mass variable m j rather than the normalised one √ ρ = m j /p t . In this figure we have 300 GeV < p t < 450 GeV. CTEQ6L parton distribution functions [52] have been used both in the convolution and MadGraph/Pythia 8/Herwig 7/MCFM results. For best fit we choose µ I = 3.5 GeV. This plot shows a good agreement between the data and the resummed prediction over the entire range of the jet mass, as well as with the Monte Carlo simulation. Resummed differential mj distribution with antikt clustering and R = 0.7 in Z(→ ℓ + ℓ − )+ jet events, with ℓ = e, µ, compared to experimental data from CMS [19] and MadGraph+Pythia 8 and Herwig 7 results. The experimental data are taken from ref. [49].
We note that the NLL+NLO+NP curve is cut off at around 40 GeV, which is just a manifestation of the shift of the resummed distribution to the right, as explained above. The value of the NLL+NLO+NP distribution at say m j = 40 GeV is related to the value of the resummed distribution at m 2 j − δm 2 j , with δm 2 j (see eq. (56)) varying from 20 GeV to around 40 GeV depending on the p t of the jet and partonic channel. Hence, we have no result for the non-perturbative distribution below this value (∼ 40 GeV) of the jet mass. Additionally, due to the Landau-pole singularity at small values of the jet mass, the distribution is unreliable in the region to the left of the Sudakov peak.

Conclusions
In this paper we have presented state-of-the-art detailed fixed-order calculations as well as all-orders estimates of distributions of important observables at the LHC. Specifically we have considered a typical jet-shape observable that has been studied quite substantially in the literature, namely the invariant jet mass. It is a member of a large class of observables known as non-global observables, that have so far proven to be quite delicate to treat. The subtleties in the analysis of such observables stem from the fact that they are defined for a restricted phase-space region. This is unlike global observables which are defined over the whole phase space. The former non-global observables receive contributions that are totally absent for their global counterparts. These contributions appear at each higher order in perturbation theory and have so far shown no pattern of iteration.
We have extended the work of ref. [11] from various angles: (a) we have implemented two clustering algorithms, k t and C-A, instead of just the anti-k t considered in the latter reference. Generally, computations in k t and C-A algorithms are much more difficult to handle than in antik t case; (b) we have computed CLs, which are completely absent for anti-k t and thus not treated in [11]; (c) we have investigated the jet mass distribution in various different processes, namely W/Z/γ/H + one jet, while only Z + jet was considered in the said reference; and (d) we have provided analytical expressions for our results in the form of power-series expansions in the jet radius R. As the experimental data [19] for the jet mass distribution in Z + jet events at the LHC with anti-k t clustering were available only after the publication of ref. [11], we made the comparison of the resummed result with these experimental data herein.
We have confirmed previous results that were arrived at in studies of e + e − annihilation processes. These include, for instance, the observation that NGLs are decreased by the application of jet clusterings other than anti-k t . In other words, NGLs are more significant when anti-k t is used. This may hint at the advantage of using other jet clustering algorithms in order to bypass the difficulties posed by NGLs. Additionally, we showed that in the limit of very small jet-radius parameter the NGLs and CLs at hadron colliders coincide with those at e + e − colliders. Moreover, we have been able to identify new features that are not present in the simple e + e − annihilation case such as the significance of initial-state radiation and its impact on the jet mass distribution. The jet mass provides a tool to discriminate gluon and quark-initiated jets as their corresponding jet mass distributions were shown to be quite different.
It is worth, as a continuation to this project, investigating other crucial hadronic processes at the LHC such as di-jet production. The latter represents an important background for numerous potential new physics signals. Another issue that is also worth tackling is performing calculations beyond two-gluon emission. This will provide a deeper insight into the nature of QCD hadronic processes that have not been fully understood so far.
where τ = 7.65 and F = 0.09. The total Born cross-section σ 0 is simply the integral of dσ 0,δ /dB δ (including the kinematical-cuts function Ξ B ) over B δ , summed over possible δ, where we have where f i denotes the parton density function for the corresponding incoming parton (i) evaluated at a factorisation scale µ F . Substituting the Mandelstam variables into the delta function in the integrand we obtain This delta function can be used to perform the integration over one of the x's, say x b , and thus we set and multiply the integrand by Since x b > 0 then x a > e y p t / √ s and y < ln( √ s/p t ). Additionally, since x b < 1 then y > − ln( √ s/p t ) and x a > e y p t / √ s + M 2 /s 1 − e −y p t / √ s .
The latter inequality overrules x a > e y p t / √ s, and furthermore, since x a < 1 we deduce that which also overrules the condition |y| < ln √ s/p t . We perform the integration over p t , y and x a , either in the Born cross-section or in the jet mass distribution, numerically via Monte Carlo method.

B Resummed Global form factor
The Sudakov global form factor that resums global logarithms is given in eq. (28). The radiator R is composed of contributions from double-logarithmic soft-collinear and single-logarithmic hard-collinear emissions from the outgoing leg (j), and from single-logarithmic soft wide-angle emissions from all legs. The leading-order soft wide-angle contribution that we calculated in section 3.1 simply exponentiates to all orders. Additionally, the non-global and clustering corrections appear as a factorised part that multiplies the resummed global form factor. The remaining soft-collinear and hard-collinear contributions may be obtained using the general formalism of resummation introduced in ref. [31] as we show below.
First we write the definition of the observable ̺ in terms of the transverse momentum k (ℓ) t , rapidity η (ℓ) , and azimuth φ (ℓ) of a single soft-collinear emission with respect to the direction of the hard leg (ℓ). Emissions that are collinear to the incoming legs (a) and (b) do not end up inside the jet, so they do not contribute to its mass, hence ̺ (a) = ̺ (b) = 0. For emissions that are collinear to the outgoing leg (j), we introduce a coordinate rotation that takes the momentum of leg (j) to the z axis, where the momenta of the jet and the soft emission become p (j) j = p t cosh y (1, 0, 0, 1) , (72a) The jet mass observable (being invariant under rotations) is then given by Comparing the definition of the normalised invariant jet mass (73) to the general parametrisation of observables from ref. [31] where Q = p t is the hard scale, we see that a j = b j = g j = 1 and d j = 2 cosh y.
Employing the master formula for resummation (eq. (3.6) from ref. [31]) we obtain the expression of the radiator R δ (ρ), for a given Born channel δ, in the MS renormalisation scheme R δ (ρ) = C j [L g 1 (α s L) + g 2 (α s L) + g 2,coll (α s L)] + + g 2,wide (α s L) C ab R 2 2 + (C aj + C bj ) h(R) , (75) with h(R) being given in eq. (26). Here C j is the colour factor associated with leg j, C j = C F for outgoing (anti-) quark jet and C j = C A for outgoing gluon jet, and C iℓ is the colour factor for dipole (iℓ) introduced in the main text. We have with λ = α s (R p t ) β 0 ln(R 2 /ρ). The factor B j accounts for corrections due to hard-collinear emissions to the outgoing jet j and is given by for quark jets , with T R the Dynkin index (normalisation constant) for the SU(N c ) generators, T R = 1/2, and n f = 5 the number of active quark flavours. Additionally we have K = C A 67 18 − π 2 6 − 5 9 n f , β 0 = 11 C A − 2 n f 12 π , (78) In the master formula we excluded the single-logarithmic soft wide-angle term referred to in ref. [31] as ln S(T ), and we calculated it manually in section 3.1. It appears as the last term in the radiator (75). The derivative of the radiator R with respect to L, relevant in the expression (28), is given by

C Fixed-order expansion
For the sake of matching we need the fixed-order expansion of the resummed form factor (49). We can cast the latter in the form dΣ δ (ρ) dB δ = dσ 0,δ dB δ C δ (ρ) exp where the expansion coefficients in the exponent, G nm , up to O(ᾱ 2 s ), are