NLL′ resummation of jet mass

Abstract
 
 Starting from a factorization theorem in effective field theory, we present resummed results for two non-global observables: the invariant-mass distribution of jets and the energy distribution outside jets. Our results include the full next-to-leading-order corrections to the hard, jet and soft functions and are implemented in a parton-shower framework which generates the renormalization-group running in the effective theory. The inclusion of these matching corrections leads to an improved description of the data and reduced theoretical uncertainties. They will have to be combined with two-loop running in the future, but our results are an important first step towards the higher-logarithmic resummation of non-global observables.


Introduction
Up to now, higher-logarithmic resummations of collider observables have only been performed for the narrow class of global observables which constrain radiation uniformly over the entire phase space. This category includes very inclusive observables such as selected event shapes, but it excludes all observables with hard phase-space cuts or a fixed number of jets. In recent years, a lot of progress was made in the theoretical analysis of non-global observables [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. This includes work on the structure of higher logarithms as well as studies of leading logarithms beyond the large-N c limit.
In this paper we start the computation of higher-logarithmic terms for non-global observables by analyzing two simple observables, the jet mass and the interjet energy flow, and presenting resummed predictions which include the full one-loop corrections to the relevant hard scattering processes, as well as the associated jet and soft functions. In the effective-theory framework we use for resummation [6,8], these correspond to matching corrections and they will need to be supplemented by corrections to the renormalizationgroup (RG) running in the future to arrive at a complete higher-logarithmic treatment of the non-global part.
Our main goal in the present work is to develop the Monte Carlo methods to include these corrections as a step towards full higher-logarithmic resummation, but it is also interesting to study their numerical size, since they have never been computed for nonglobal observables and often dominate numerically in the global case. It is customary to add a prime to the logarithmic accuracy to indicate the presence of higher-order matching corrections. In this notation our next-to-leading-logarithmic results for the jet mass have NLL accuracy.
In Refs. [8,10] we have derived a factorization formula for interjet energy flow and lightjet mass. The key element is the presence of multi-Wilson-line operators which generate the intricate pattern of Non-Global Logarithms (NGLs). Explicitly, the result for interjet energy flow at a lepton collider has the form where Q is the center-of-mass energy, and Q 0 = βQ is the energy scale above which we veto energy in the gap outside the jet cones. For simplicity, we choose the jet axis along the thrust axis. The above factorization formula neglects power corrections from O(β) terms. The hard functions H m describe hard radiation inside the jet cone, and their characteristic scale is Q since radiation inside the cones is unrestricted. The index m represents the number of hard partons inside the jet, which propagate along the directions {n} = {n 1 , n 2 , . . . , n m }. Each of these sources soft radiation, which we describe by a Wilson line along the direction of the hard parton. The matrix elements of these Wilson lines define the soft functions S m ({n}, Q 0 , µ). To obtain the cross section, one integrates over the directions {n}, which is indicated by the symbol ⊗. The hard and soft functions are matrices in the color space of the m partons and one takes the color trace . . . after multiplying them. The operator definition for these functions and further explanations can be found in [8].
The second observable we consider is the jet mass distribution at a lepton collider. To define the jet mass, we use the thrust axis to split every event into two hemispheres. One can then (randomly) select one of the two jets and compute its invariant mass M , which is usually discussed in terms of the dimensionless variable ρ = M 2 /Q 2 . Alternatively, one computes the mass in both hemispheres and chooses the heavier mass ρ h or lighter one ρ . Obviously, there is a relation among the these observables: the jet mass distribution is simply the average of heavy-jet mass and light-jet mass one We will call the hemisphere we select to measure the mass the left one, which means that the radiation in the right hemisphere is unconstrained. 1 We introduce a light-like reference  4). The black lines represent hard radiation with typical scale Q which is constrained to be inside the cones, and the red lines depict soft radiation with a low energy scale Q 0 which is allowed to populate the full phase space. In the right figure, the blue lines in the left hemisphere represent collinear radiation which is described by the inclusive jet function in (1.4).
four-vector n µ = (1, 0, 0, 1) pointing to the right along the thrust axis and an opposite vectorn µ = (1, 0, 0, −1) pointing to the left. The hard partons in the right hemisphere then generate the complicated pattern of soft radiation and associated NGLs. The main difference to formula (1.1) is that one also needs the standard inclusive jet functions to describe collinear radiation in the left hemisphere. Resummation effects in the jet mass distribution have been discussed in Refs. [17][18][19][20][21], however only in [17] the leading NGLs were resummed. Our work is based on the factorization theorem for jet mass derived in [10]. The invariant mass of the left jet is obtained from the momentum pc of the energetic particles collinear ton and the soft partons in the left hemisphere, ρ Q 2 = M 2 = (pc + p s ) 2 = p 2 c + Qn · p s + O(p 2 s ) .
(1. 3) In the factorization theorem, the sum results in a convolution of the soft and jet functions.
To avoid this, one can work in Laplace space, where the factorization formula has the product formσ where τ is the Laplace conjugate variable of ρ, andj i is the inclusive jet function [22,23], which by now is known to three loops [24,25]. In (1.4) the index m indicates the number of partons in the inclusive (right) hemisphere, so that m = 1 at leading order (LO). As long as we consider large jet cone sizes of O(1), the leading-logarithms (LLs) in interjet energy flow at a lepton collider are of the form α n s ln n β. The interjet energy flow is a single logarithmic observable, because collinear logarithms cancel inside the large cone region and only soft logarithms remain. These logarithms arise from the multi-Wilson-line operators S m in (1.1) and one needs to use parton shower methods to resum the enhanced logarithms already at the LL level. In [15] we have written a dedicated parton-shower code to perform the resummation for such observables and have interfaced it with the MadGraph5 aMC@NLO event generator [26]. This provides an automated framework to perform the LL resummation for single-logarithmic observables. However, collider observables are typically double logarithmic. The leading logarithms in the jet mass distribution, for example, are α n s ln 2n ρ. Even for non-global observables, these double logarithmic terms have a simple structure, and they can be factored out and treated separately. In the parton shower framework, we therefore subtract these "global" contributions and exponentiate them manually, as Dasgupta and Salam did in their original paper on NGLs [27]. Given their different nature, it is interesting to analyze both the interjet energy flow and the jet mass as examples and we will present LL and NLL improved results for single logarithmic and double logarithmic observables, separately. A second motivation to also analyze the jet mass, is that there are LEP measurements to which we can compare to, in contrast to the interjet energy flow. Unfortunately, the typical jet mass at LEP jet is quite low M 10 GeV, which translates to a scale of the soft radiation of Q 0 ∼ M 2 /Q 1 GeV so that non-perturbative effects are very important in the peak region of the distribution.
Our paper is organized as follows. In the next section, we will discuss LL resummation for interjet energy flow and show how one implements the one-loop corrections to the hard and soft functions. We then move to the jet mass distribution in Section 3, focussing on the differences to the single-logarithmic case. We will in particular show how to subtract global logarithms in the parton shower and in the soft function. After presenting numerical results in Section 4 and comparing to LEP data and PYTHIA results, we conclude in Section 5.

Interjet energy flow at LL accuracy
The perturbative expansion of the interjet energy flow in (1.1) suffers from large logarithms of the ratio of the hard scale Q and the soft scale Q 0 . To resum these, one solves the RG equation of the hard function and evolves it from its characteristic scale µ h ∼ Q down to a soft scale µ s ∼ Q 0 . This yields the RG-improved expression [8] where the evolution factor is defined as a path-ordered exponential of the anomalous dimension The RG-evolution generates additional partons and maps the l-parton configuration along the directions {n } = {n 1 , . . . , n l } into an m-parton final state along the directions {n} = {n 1 , . . . , n l , n l+1 , . . . , n m }. The symbol⊗ in (2.1) indicates the integral over the directions of the additional m − l partons generated in the evolution. At the leading logarithmic level, we only need the one-loop anomalous dimension and can rewrite the exponent as In the last step, we have introduced the evolution time t ≡ t(µ h , µ s ). For a given µ h , there is a one-to-one correspondence of the evolution time to the low scale µ s . Obviously, for µ h = µ s , we have t = 0. During the evolution, t grows and goes to infinity as µ s hits the Landau pole. For µ h = M Z and two-loop running with a Landau pole at Λ = 0.230 GeV, the choice µ s = 1 GeV corresponds to t = 0.08. A plot connecting t and µ s for different values of µ h can be found in Figure 1 of our previous paper [15].
In [15] we implemented the RG evolution factor U ({n}, µ s , µ h ) in the large-N c limit using the parton shower method proposed by Dasgupta and Salam in [27]. We don't want to repeat the entire discussion here, but we give the algorithm in Appendix B, since we need to extend it to compute the soft functions, as discussed below. Let us also list the one-loop anomalous dimension, since its form will be relevant in the discussion of the jet mass below. It is given by [8] The entries R m and V m are angular functions associated with the emission of a real or virtual soft gluon and take the form where the color matrices T i,L act on the hard function from the left, i.e. on the amplitude, while T i,R acts on the conjugate amplitude. The sum runs over all unequal pairs (ij) of the m hard partons. The anomalous dimension involves the dipole radiator , (2.6) which is given by the product of the associated eikonal factors. In the virtual corrections, one integrates over the direction n k of the emission. We note that individually R m and V m suffer from collinear divergences, which cancel in the cross section. In the Monte Carlo implementation, one works with a collinear cutoff to regularize the divergences. As long as we choose the µ h and µ s properly, the hard and soft functions will be free of large logarithms and the large logarithmic terms are resummed in the evolution factor. Because they are free of large logarithms, the higher-multiplicity hard functions are suppressed by α s as H l ∼ α l−2 s H 2 . At LL level, we thus only need to include the hard function H 2 and the soft function is given as the unit matrix in the color space S m ∼ 1. At LL accuracy, the RG-improved result (2.1) simplifies to To extend these results to NLL, one needs two ingredients: the one-loop matching corrections and the corrections to the RG running due to the two-loop anomalous dimensions. The present paper focuses on the first set of corrections, i.e. LL accuracy. Specifically, we need one-loop corrections to H 2 , the tree-level result for H 3 and the one-loop soft functions S m . We write their perturbative expansions in the form (2.8) In this notation, the full LL resummed cross section takes the form (2.9) -6 - We used here that the leading-order soft function S (0) m is the unit matrix 1 in color space. The first line contains the LL result (2.7), and the remaining three lines show the different NLO corrections, which are depicted in Figure 2.
The hard functions H m include the momentum conservation and phase-space constraints on the hard partons. For two partons, these constraints render the integrals over the parton directions trivial. The momentum and jet direction constraints impose that the vectors n 1 and n 2 must point along the thrust axis and in opposite directions so that where we have used that also the color structure is trivial for two hard partons. The function H 2 (Q 2 , µ) is the standard dijet hard function which arises also for global observables such as the event shape thrust. In the large-N c limit, we should replace C F → N c /2. In [8] we have derived an expression for the hard function H 3 , which corresponds to the QCD process γ * → q(p 1 )q(p 2 )g(p 3 ). By definition H 3 only depends on angular information of the three partons, since their energies have already been integrated over. For convenience we split the phase space integration into different regions according to the direction of the thrust axis, which for three-parton final states points in the opposite direction of the most energetic parton. Due to momentum conservation, the three partons must be in a plane. Using invariance of the cross section under rotation around the thrust axis, in Region I only the angles θ 2 and θ 3 , between the partons and the thrust axis, are not fixed.
For convenience we parameterize these angles in terms of two variables u and v each going from 0 to 1 and defined aŝ where the variable v is directly related to the larger angle θ 3 , while u characterises the relative size of the angles. Please note that the variables u and v differ from the quantities of the same name used in [8], where we defined the variables such that v = 1 corresponded to the angle of the jet cone, rather than a 90 • angle as in (2.12). Because the same hard function H 3 also arises for the jet mass studied below, we prefer to not incorporate the specific phase-space constraint into its parameterization.
The bare hard function H 3 in terms of the anglesθ 2 andθ 3 was given in (4.4) of [10]. The corresponding representation includes a θ-function constraint imposed to prevent the thrust axis from flipping. For simplicity, we choose the jet opening half-angle α ≤ π 3 so that the axis constraint is automatically fulfilled. The hard function suffers from divergences when u and v go to zero. In dimensional regularization after performing MS subtraction, the contribution of Region I to the renormalized hard function H (1) 3 is given by The function Θ in (v) ensures that all hard emissions are inside the jet. For the interjet energy flow it is given by where α is the jet opening halfangle. In the large-N c limit, the color structure of the hard functions becomes trivial and we use non-bold symbols such as H (1) 3,I to indicate the scalar quantities which are relevant in this limit. The expression for the auxiliary function F (u, v) is given by (2.14) Similarly, in Region II we have Region III describes the situation, where the gluon is the most energetic particle and we parameterizeθ 1 = uv,θ 2 = v. The hard function reads . Left: Angular dependence of S 3 for fixed evolution time t = 0.08. Note that the angles θ q and θ g of the hard partons to the jet axis must be smaller than the cone angle α = π/3 ≈ 1.04. Right: Dependence on the evolution time t at fixed angles.
Next, we will discuss how to implement the above expressions into the parton shower code. We first rewrite the angular integral in the H which is the LL RG evolution or parton shower soft function. To implement this formula into a Monte Carlo framework, we will randomly generate u and v and then run the shower S 3 (u, v, µ h ) for the given configuration. There is, however, one complication, namely that the hard function is a distribution and can therefore not be integrated point by point. One way to solve this problem is to evaluate S 3 (u, v, µ h ) on a grid, interpolate and then perform the integrations over u and v. This works well because S 3 (u, v, µ h ) is a smooth function of the angles as can be seen from Figure 4. Note in particular that the limit v → 0, in which both angles go to zero and the two Wilson lines become collinear, is completely smooth. In this limit the quark and gluon Wilson lines combine and produce the same radiation as a single quark Wilson line, encoded in the function S 2 . The relation will lead to important simplifications below. In the right plot, we show the evolution time dependence of the soft function S 3 for fixed angles. One observes that the function falls off much faster when the hard partons approach the jet cone. In this configuration, more soft radiation exits the cone, explaining this suppression.
Interpolating the soft function S 3 gives accurate results, but is not efficient since the function depends on the phase-space constraints and thus needs to be recomputed when one changes the cone angle. It is much more natural to compute the convolution (2.19) directly in the Monte Carlo code. The simplest way to implement the plus distributions in the hard function into the Monte Carlo is to use a slicing method. To explain it in a simple setting, let us for the moment only consider the v dependence and forget about the variable u. Then the convolution (2.19) takes the form where B(v) represents a regular function. Thanks to relation (2.20) the A term can be combined with the LL parton shower result involving S 2 and the contribution from B(v) can be computed by randomly generating v-values and running the shower for each chosen configuration. The slicing method introduces a lower cutoff v 0 into the plus distribution integrals C i (v) to ensure that v can not go to zero. With the cutoff in place, we can integrate the subtraction term, e.g.
where one can use the same Monte Carlo method as for the B(v) terms to simulate the first term with the collinear cutoff v 0 , and then adds back the second term which is given by the LL parton shower result, multiplied by a logarithm of the cutoff parameter. The v 0 dependence will cancel out between the two terms up to power corrections. The power corrections in the artificial parameter v 0 can be neglected as long as one chooses it small enough. The slicing method involves large cancellations between the two terms on the right-hand side of (2.22), so for numerical stability reasons one should not choose v 0 too small. These two opposing requirements make slicing methods delicate, but we compared to the result using the interpolated soft function S 3 and found good consistence. The cutoff independence is demonstrated in Figure 12 in Appendix A.
Up to now we have disregarded the u dependence, but the Monte Carlo implementation of the full equations (2.13), (2.15) and (2.17) involves nothing beyond the above discussion, except that we have to consider both integrations. As (2.20) shows, the soft function becomes trivial for v → 0 and we can combine all δ(v) dependent terms with the parton shower for S 2 . We thus only need to apply the slicing method to the δ(u) ln i v/v + and (1/u) + (1/v) + terms. The corresponding cutoff dependent compensation terms are collected in Appendix A.
The final ingredient we need to implement is the one-loop soft function, which is defined as a sum over all dipoles where the sum runs over all unordered pairs (ij). In the large-N c limit only neighbouring legs give a contribution We evaluate the one-loop soft function numerically within our Monte Carlo code. It is well suited for this task since it generates emissions between neighbouring dipoles in an efficient way, by randomly choosing the rapidityŷ and azimuthal angleφ of the emission in the COM (center-of-mass) frame of the emitting dipole (n i , n j ). Here and in the following, we will use hats to indicate kinematic quantities in the COM frame. Our hard function shower keeps emitting additional hard partons until one of them enters the veto region at which point it terminates. In our implementation, we use this last parton in the veto region to obtain the NLO correction to the soft function. At NLO, the renormalized soft function can be expressed as with Θ lab out (ŷ,φ) constraining soft radiation to be outside of the jet cone in the lab frame. In the Monte Carlo implementation, the factor in square brackets is a weight factor for the corresponding emission. The auxiliary function f ij (φ,ŷ) connects the transverse momentumk T in the COM frame to the energy Q 0 in the lab frame,k T f ij (φ,ŷ) ≤ Q 0 , and is given by f ij (φ,ŷ) = 2 M −β cosφ + coshŷ , where M 2 = 2 n i · n j is the invariant "mass" of the dipole pair, and β = 1 − M 2 /4. The logarithm of | sinφ| arises from expanding the azimuthal angular integration in , which is related to the space-time dimension through d = 4 − 2 . A detailed derivation of expression (2.25) can be found in Appendix A.
While our slicing implementation of the hard function is simple but specific to the dijet processes and certainly not optimal, the above procedure to obtain the NLO soft function is simple, efficient and general. Compared to the LL parton shower code, including the one-loop soft function correction (2.25) yields where one evolves the hard function from hard scale to soft scale and multiplies it with the soft function S (1) m of the corresponding multiplicity. When running our Monte Carlo code we fill three histograms, one for the LL shower, one for the logarithmic part of (2.25) and one for the non-logarithmic part. Further details of the Monte Carlo algorithm, including the implementation of the one-loop soft function are given in Appendix B.
The computer time needed to run the shower including the one-loop corrections depends on the maximum evolution time needed in the computation. For the interjet energy flow, we run the shower until t = 0.08, corresponding to µ s ≈ 1 GeV. For a collinear cutoff at η cut = 4 (η cut = 5) in the parton shower we then end up with about 15 (30) hard partons per event on average. To resolve the peak region of the jet mass, discussed in the next section, we have to run to extremely low scales µ s = 0.275 GeV, corresponding to t = 0.3, near the Landau pole at Λ = 0.230 GeV. At this scale, hundreds of partons are generated in each event and we need a few days of computer time on a cluster to obtain our numerical results, which will be presented in Section 4 below.

NLL resummation for jet mass
Our second task is to perform the resummation for the jet mass distribution at electronpositron colliders. In contrast to the interjet energy flow, this observable suffers from soft-collinear double logarithms. These then constitute the LL results, while the nonglobal structure only arises at NLL. The resummation of jet mass including the leading non-global logarithms has been discussed in [10,17,27,28]. At NLL level, the non-global logarithms yield a simple overall factor which multiplies the cross section. Beyond NLL this simple factorization does not hold anymore, and one needs to include the corrections piece by piece. 2 The basic structure of the corrections is of course the same as for the interjet energy flow, see (2.9) and Figure 2, and we therefore mainly focus on the differences to this case. In addition to the double logarithms, the most important new element is that the factorization arises in Laplace space. We use the same notation as [10], where we presented NLL resummation results. For NLL accuracy we need to keep one-loop matching corrections in the factorization formula (1.4) and the theorem then reads In the first line we must include one-loop corrections for the quark jet functionj q , the hard function H 1 and soft functions S m . We do not include the O(α 2 s ) cross terms so that the first line turns into a sum of terms with the individual corrections. The hard function H i 2 in the second line includes two hard partons in the right jet. Since it involves a power of α s due to the hard emission, the remaining ingredients are only needed at LO. The second line also includes a gluon-jet contribution, for the case where the qq pair is in the right hemisphere. The one-loop hard functions are the same as for interjet energy flow, up to the different phase-space constraints. They are given in Appendix C.
In Laplace space, RG-evolution is multiplicative and we can factor out and exponentiate the double logarithms. Removing the double logarithmic part is important since our shower evolution, which also takes place in Laplace space, is purely soft. The subtraction of collinear contributions will also be needed for our numerical computation of the one-loop soft function. Using standard techniques introduced in [30], we can perform the inversion to momentum space analytically at the end and write a momentum space result directly in terms of Laplace-space ingredients.
The anomalous dimension Γ H in (2.4) which drives the resummation of the logarithms in interjet energy flow (1.1) can be viewed in two ways: As the hard anomalous dimension, used to evolve the hard functions to the soft scale, or as the soft anomalous dimension which evolves the soft functions to a higher scale. RG invariance of the cross section implies that the two evolutions must agree. The situation is more interesting for the light-jet mass (1.4) which involves three ingredients. In this case RG invariance translates into the statement where The Casimir C i for the quark-jet channel is C q = C F , while the gluon configuration has C g = C A . In our paper [10], we have analyzed the one-loop soft anomalous dimension and found that it has the form whereΓ lm is a regular non-logarithmic anomalous dimension, which takes the same form as (2.4), except for a subtraction to remove the collinear singularities, which give rise to the cusp piece in (3.4). The subtraction is achieved by replacing the diagonal elements in · n n · n k n k · n Θ L (n k ), (3.5) where Θ L (n k ) ensures that the emission is in the left hemisphere with the light jet. The trivial color structure arises from color conservation Note that V 0 is equal to the one-loop result (real plus virtual) for the case where there is only one hard parton on the right, which then, by momentum conservation, flies along n. The subtraction therefore removes the "global" one-loop part of the soft anomalous dimension. After this, the Monte Carlo result no longer involves collinear singularities. As before we regularize the collinear singularities in the individual entries ofΓ using a cutoff. The parton shower algorithm of Dasgupta and Salam [27] instead uses a veto algorithm to remove global logarithmic terms. Our subtraction of the global piece has the advantage that our Monte Carlo weights are always positive. Let us also note that the role of the subtraction is to separate out the collinear singularities, so that the same subtraction can be used for any process with the same double logarithmic structure, i.e. also in cases with more complicated geometry, where we cannot analytically compute the one-loop function.
To make use of the separation of the anomalous dimension into two pieces, we now factor the soft function as The splitting of the soft function into single and double logarithmic pieces is of course not unique. We have chosen the double-logarithmic "global" part S i G such that it includes the full one-loop result, so that the "non-global" remainder functionŜ i m starts at two loops for m = 1 partons in the right hemisphere. For the gluon case, we only need the tree-level result S g G = 1 since the hard function for this channel is suppressed by α s . The global piece fulfills a standard RG-evolution equation driven by the cusp piece of (3.4) which can be immediately solved in Laplace space. Using the technique introduced in [30], the associated momentum-space solution takes the form with η S = 2C i A γcusp (µ s , µ), where the logarithm L s has been replaced by a derivative operator with respect to η S . With the global function at hand, the Monte Carlo simulation only needs to provide the remainderŜ i m . Its single logarithmic RG-evolution is obtained by the subtracted parton shower described above and the one-loop correction for an m-parton configuration is given byŜ which, by construction, is free from collinear logarithms. We compute this difference in the large-N c limit by running the shower until it produces a parton in the left hemisphere, which is the veto region for the present case. The outside parton is the soft emission and we then compute the relevant one-loop weight factor precisely as in (2.25). The form of the Laplace space soft function can be found in the appendix in (C.7). When the emission arises from the first dipole, which involves the left parton along n 0 =n, we subtract the global part. For the quark-jet channel the subtraction is given by with a re-weighting factor The factor X is simply the ratio of the radiator (2.6) associated with the original (n, n) dipole and the one of the dipole (n, n j ) which emits the gluon and defines the frame in whichŷ andφ are generated. The subtraction removes the collinear divergence in the (n, n j ) dipole and yieldsŜ q m . The function g ij in (3.11) relates the momentum component n · k in the lab frame to the transverse momentumk T in the COM frame of the dipole (n i , n j ), analogously to the function f ij in (2.25). Its explicit form is given in the appendix in (C.8).
The final ingredients in (3.1) are the one-loop jet functions, which are well known. In Laplace space, the one-loop jet function is given bỹ with η J = 2C i A γcusp (µ j , µ). The relevant expressions for the ingredients are listed in Appendix E. Combining the global soft function with the jet function, we obtain where we define η = η J + η S . The full result is obtained after combining this with the subtracted shower evolution, the hard functions and the one-loop soft correction (3.10).
To implement this expression in practice, we run the shower, tabulate the results for the individual contributions to (3.1) and then replace the global function S i G (ω, µ h ) in (3.15) by the full result which includes the hard functions, evolution and one-loop corrections.
Up to NNLL, the integrated heavy-jet mass distribution is obtained as

Numerical results
In this section we will present numerical results, first for the interjet energy flow, then for the jet mass. For our plots, we work with Q = M Z and α s (M Z ) = 0.1181, and use two-loop α s (µ) running with n f = 5 quark flavors. To our knowledge, no measurements are available for the interjet energy flow, but we will compare our results for the jet mass to LEP measurements by ALEPH [31].

Interjet energy flow
For our numerical discussion we choose jet cone size parameter as α = π/3. This is equivalent to δ = tan α 2 = 1/ small cone angles, or equivalently large rapidity gaps, in order not to have to deal with large collinear logarithms. In our plots we show the gap fraction which is the fraction of events in which the soft radiation outside the jets has an energy E s below the cutoff Q 0 . By definition, the amount of energy in the gap must be below Q/2, otherwise the thrust axis, which defines our jet axis, would flip. The fixed order result is therefore R(Q 0 = Q/2) = 1 at any order in perturbation theory. The O(α 0 s ) result with just two back-to-back partons is of course R(Q 0 ) = 1, a nontrivial Q 0 dependence only arises at O(α s ) when the third parton is inside the gap. We will refer to the O(α s ) result as LO.
As a first step, let us check the size of the individual corrections and investigate whether the scale dependence is reduced after including them. In Figure 5 we show the hard and soft corrections separately and then plot the scale bands from varying the associated scales by a factor two around their default values µ h = Q and µ s = Q 0 . Compared to the LL scale bands shown in red, the scale dependence is reduced in both cases after including the corrections. We observe that the hard corrections are quite significant and positive, while the soft corrections are moderate and negative. The hard corrections have two sources, virtual corrections to H 2 and real emission contributions encoded in H 3 . The first of these is just a constant factor multiplying the LL result, while the second one comes together with the higher soft function S 3 . Both corrections are positive. At high values of Q 0 the three parton contribution from H 3 is about twice as large as the one from the one-loop correction to H 2 and it becomes more dominant at smaller values.
It is clear that the large hard function corrections at Q 0 Q/2 must be compensated by terms which are power suppressed in Q 0 /Q and are not captured by the resummation based on the factorization formula (1.1), which arises in the limit Q 0 → 0. One can obtain these power suppressed terms by matching to the fixed-order result. More precisely, one adds to the resummed result the fixed-order prediction minus its expansion around Q 0 . The subtraction removes the terms which are already included in the resummation. These power suppressed matching terms can be obtained as To evaluate this integral, one computes the cross section to find a parton inside the gap and subtracts from it its soft limit. The subtraction eliminates the virtual contributions and leads to a finite integral, which one can evaluate numerically. However, even after the matching to the fixed order result, the resummed result does not yet tend to R(Q 0 ) = 1 for Q 0 → Q/2 because we resum logarithms of µ s /µ h → 1/2 for µ s ≈ Q 0 and µ h = Q. To switch off the resummation, one can choose the soft scale in such a way that it approaches the hard scale µ h as Q 0 → Q max = Q/2. This can be achieved, for example with a profile function [32] of the form or even resum, the power corrections to resolve the difference. The first step would be to include the matching up to NNLO, which would in principle be possible since the fixedorder results are available [33][34][35]. In practice it would require some effort since we would need to compute the fixed-order expansion of our results (including the shower).
In Figure 7, we show an improved numerical result which includes the matching correction ∆R(Q 0 ), shown as a black dotted line, and uses the scale choice (4.3) to switch off the resummation at the end-point. The matching correction is negative and compensates the large hard corrections near the end-point. The LL corrections lead to a larger gap fraction R(Q 0 ). As mentioned earlier, there is unfortunately no experimental data to which we can compare our results, but we compare to PYTHIA [36]. While the two results are similar at very low Q 0 , PYTHIA is higher at intermediate values. We remind the reader, that the intermediate values heavily depend on the profile function used to switch off the resummation.

Jet mass
Let us now turn to the jet mass ρ. For interjet energy flow, we considered the integrated cross section, i.e. all events with energy in the gap below the veto, while we will look at the differential spectrum in the present case, since this is what was measured by the LEP experiments. We will however compute the spectrum by taking the derivative of the integrated cross section, which has the advantage that the spectrum is correctly normalized if the resummed prediction for the integrated cross section matches the fixed-order result at large ρ.
As a first step, we again separately plot the different ingredients and their scale dependence in Figure 8. In the first three plots we compare NLL to NLL with corrections from the jet, hard and soft functions. The red bands are the NLL result with scale variation, where we vary either the jet, hard or soft scale by a factor of two around the default values µ h ∼ Q, µ j ∼ √ ρ Q and µ s ∼ ρ Q. The blue curves show contributions at NLL accuracy from one of the three ingredients with its associated scale variation. Obviously, the scale  Figure 8. NLL corrections from the jet, hard and soft functions and their scale uncertainties. Each band comes from varying the scale associated with the correction by a factor of two around the default value. In the last plot we show LO power corrections from the fixed-order computation.
We have multiplied the distributions by ρ in order to make the results at larger ρ visible.
dependence is strongly reduced from NLL to NLL for jet and hard corrections. The soft scale dependence, on the other hand, is only modestly reduced after including one-loop soft function corrections. The scale bands mostly overlap with each other, which indicates that perturbative convergence is reasonably good in all the three cases.
In the last plot of Figure 8 we show the effect of adding the O(α s ) power corrections to the NLL results. The LO power corrections for the heavy-jet mass are known analytically and given in Appendix E. They are the same as for thrust, because the three-parton results for jet mass and thrust agree. Since the light-jet mass vanishes at O(α s ), we can immediately also obtain the LO power corrections for the jet mass distribution. From the plot, we observe that the difference between NLL and NLL + LO is very small, and that the contributions from power corrections will reduce the resummed result in the large jet mass region. In order to reproduce the full fixed order result, we use C F = 4/3 instead of the strict large-N c value C F = 3/2 for the hard, jet and soft one-loop corrections in the  Figure 9. Jet-mass distribution compared to PYTHIA results. On the left side we plot our default result, based on using the profile scale (4.3) and exponentiating the matching corrections. On the right-hand side, we do not perform these modifications such that we get a negative cross section at low ρ and hit the Landau pole at a nonzero ρ.
resummed results. We also use the exact color factors in the evolution factors of the global part (3.15).
The end-point of the jet mass distribution is at ρ max = 1/3 at O(α s ), corresponding to a symmetrical configuration of the three partons. We will work with the same profile function (4.3) to switch off the higher-order terms at the end point. To adapt it to the present case, we set Q 0 = ρ Q and Q max = Q/3. For simplicity, we will adopt the canonical value µ j = √ µ s µ h in the following and only indirectly vary the jet scale through the variations of µ s and µ h , which we vary independently by a factor of two around their default values. At very low values of ρ, the scale µ s (Q 0 ) hits the Landau pole at Λ = 0.23 GeV. Near the pole the soft corrections become large and negative, resulting in a negative cross section. To avoid this unphysical behaviour, we replace µ s (Q 0 ) → µ s (Q 0 )+Λ so that the pole occurs at ρ = 0. We also exponentiate the hard, jet and soft corrections to avoid the negative cross section. In the left plot of Figure 9 we show our result for the jet mass distribution after these modifications. In the right plot, we show the result with µ s (Q 0 ) = ρ Q and without exponentiation. We observe that the soft scale dependence changes sign at a point to the right of the peak. In this region the soft scale dependence becomes very small. With the modifications in µ s , we end up with quite small scale bands to the right of the peak, which are likely not an accurate characterization of the true uncertainties. The NLL peak in the right-hand plot is quite a bit higher because the cross section becomes negative below ρ = 0.004 and our distributions are by construction normalized. An important feature of our result is that peak occurs at a very low value ρ ≈ 0.006, which corresponds to µ s ≈ 0.5 GeV so that the peak region is strongly affected by nonperturbative effects. In Figure 9 we also show the PYTHIA [36] results, both on the parton level (dashed lines) and including hadronisation. The hadronisation effects shift the peak to the right by about ∆ρ ≈ 0.006, in accordance to what one expects from non-perturbative effects in the soft  Figure 10. Jet-mass distribution and comparison to ALEPH data [31] (green dots with error bars). The black curve represents the LO prediction for jet mass, where its analytical expression is given in (E.3). The red curve is the NLL resummation result and the band is from scale variation. The blue curve corresponds to NLL + LO results, in which we switched off resummation effects at large ρ using (4.3).
functions [37,38]. The parton-level PYTHIA result is quite close to the NLL result.
In Figure 10 we compare the NLL + LO jet mass distribution with ALEPH results [31], obtained by combining their measurements for the light-jet and the heavy-jet mass using (1.2) and adding the uncertainties on the individual measurements in quadrature. One immediately sees that the experimental peak shifted to the right from non-perturbative effects and the shift is compatible with the PYTHIA hadronization result. We also observe that the jet mass distribution falls off quite rapidly and to make the region of larger ρ visible, we include also a logarithmic plot in Figure 10. The plot also illustrates what motivated the profile function (4.3) with n = 4. The choices ensures that we start switching off the resummation fairly quickly about half-way to the endpoint and go over to the fixed-order result. The plots show that, compared the LO fixed-order result, resummation greatly improved the description of the experimental data. On the other hand there is -if at all -only a relatively narrow region in ρ in which both higher-order power corrections and non-perturbative corrections are small.
For completeness, we show in Figure 11 numerical results for the heavy-jet mass ρ h and the light-jet mass ρ . The heavy-jet mass is global and provides a reference variable at the same accuracy, but free from all the complications which arise for the jet mass. From the difference of the heavy-jet mass and the jet mass we obtain the light-jet mass. This is more sensitive to the non-global structure and also only has a nontrivial distribution at O(α 2 s ) so that there is no matching at the accuracy we work. The end-point for the NLO light-jet mass is at ρ max = 1/6, which is achieved when the four parton momenta form a tetrahedron, and we use this as the endpoint in our profile function (4.3). From the plot, one observes that also the heavy-jet distribution is affected by nonperturbative effects in the peak region, however, the peak is at a larger ρ value than for the jet mass itself. Not surprisingly, the worst description of the data arises for the light-jet mass distribution.  Figure 11. Light-jet and heavy-jet mass distribution in comparison to ALEPH data [31].
At larger ρ values the description is worse because the fixed-order result starts at O(α 2 s ) so that the matching corrections are beyond the accuracy of our computation. The peak region is not well described because it is in the nonperturbative regime and very narrow.

Conclusion and Outlook
In this paper we analyzed non-global observables and, for the first time, went beyond a resummation of only the leading non-global logarithms. Specifically, we analyzed the single-logarithmic interjet energy flow at LL and the double-logarithmic jet mass at NLL . The prime indicates that we included the full next-to-leading-order corrections to the hard and soft functions, as well as the jet function in the case of jet mass. The practical implementation of these corrections is the main result of the present paper. To achieve full NLL resummation for the interjet energy flow, and NNLL accuracy for the jet mass, we will need to also include the two-loop corrections the RG running, but we observe that the inclusion of the one-loop matching corrections already leads to an improved description of these observables. Since the jet mass peaks at a low value corresponding to a soft scale of M 2 J /Q ≈ 0.5 GeV for LEP energies, the peak region is strongly affected by non-perturbative effects, similar to what is observed for other event shapes.
Due to the intricate structure of the soft emissions, factorization theorems for nonglobal observables and the associated RG evolution are much more complicated than in the global case. Instead of analytical computations, one needs to resort to a numerical Monte Carlo framework to perform the resummation. While the global heavy-jet mass only involves a soft function with two Wilson lines, the shower evolution for jet mass produces additional legs, and for low jet masses we can end up with soft emissions from hundreds of hard partons. However, concerning the NLO soft function, this is a minor complication, since we only connect pairs of legs at this accuracy. Indeed, the inclusion of the NLO corrections to the soft function is a minor modification of the leading-logarithmic shower framework. Using the shower emissions which end up in the veto region, we are able to compute the next-to-leading-order correction to the soft function in a general way, with almost no additional computer time.
The more involved part is the implementation of the NLO hard functions. These are in essence the usual real and virtual fixed-order corrections to the Born-level process, but individually suffer from collinear divergences. Computing them in dimensional regularization and renormalizing, one ends up with distributions in the angles of the hard partons which must be implemented into the Monte Carlo framework. We do this with a simple slicing scheme, which works well for two-jet production in e + e − but is certainly not the most efficient method. The problem of combining a parton shower with fixed-order results arises of course also for general purpose showers and elegant solutions such as MC@NLO [39] and POHWEG [40] are available and have by now been fully automated. A complication in our case is that our shower systematically neglects small soft momenta and therefore does not conserve momentum. As a result, its kinematics is different from the one in the hard functions. While more work is needed on the NLO hard functions, let us note that we have achieved full automation for the leading-order hard functions in our previous paper [15] by working with Les Houches event files generated by the tree-level generator in MadGraph5 aMC@NLO. The same code also provides NLO shower matching and it would be very interesting to adapt it to our shower.
An important next step is of course the inclusion of second-order corrections into the RG-running to achieve the full resummation of subleading non-global logarithms. The corresponding anomalous dimension matrix involves three types of corrections: Double real emissions, real-virtual terms and fully virtual two-loop corrections. The relevant anomalous dimension matrix has been presented in a related framework by Caron-Huot [4]. We are working on determining the anomalous dimension also in our formalism. The implementation into a Monte Carlo framework will be nontrivial, because one needs to numerically handle the collinear singularities of the individual entries. There are a number of recent papers addressing the issue of double emissions in general parton showers [41][42][43][44].
A second interesting challenge is the inclusion of finite-N c effects, especially for nonglobal observables at hadron colliders. Our RG-evolution framework is in the general class of showers characterized in [45] and valid at finite N c , but implementing the interference effects and complex phases which arise beyond N c → ∞ is challenging. Interesting progress towards the computation of such corrections has been made in [14,46].
We have analyzed two simple non-global observables in the present paper. This is a first step, but our ultimate goal is of course to use the same methods to understand jet structure at the LHC. For narrow jets, the non-global structure actually factorizes into a structure for each separate jet [6,8,17]. Boosting our hemisphere jet mass result such that the left hemisphere transforms into a cone of radius R, one immediately obtains the non-global structure of the jet mass for an LHC jet of this radius. It will be interesting to analyze such observables in the future. where Ω d is the surface of the d-dimensional unit sphere and Ω 1 = 2. Introducing the auxiliary function f ij via we can perform the integral overk T . This integration yields a soft divergence, which is renormalized away in the MS scheme. After expanding in we then immediately arrive at expression (2.25) which only involves a finite angular integration which we perform with the parton shower, which generates its emissions using the variablesŷ andφ.

B Monte Carlo algorithm for the interjet energy flow
The inclusion of the NLO soft function is only a minor modification of the algorithm for LL resummation. In fact, the first three steps are identical to what was shown in Appendix B of [15]. The only difference arises in the last step, where we also compute the soft function.
To record the results of the shower, we fill three histograms: h U contains the LL evolution, h L the coefficient of the logarithm of the soft function (2.25) and h c its non-logarithmic part.
The shower algorithm for the evolution of the function H 2 ({n 1 , n 2 }, Q, µ h ) to lower scales involves the following steps: 1. Start at evolution time t = 0 from an initial event E with vectors {n 1 , n 2 } and weight w = 1.
2. Generate a random time step ∆t according to the probability distribution P E (t) = V E exp(−V E ∆t), and insert the event weight w into the histogram h U at time t + ∆t.
3. Choose a dipole associated with a pair of neighbouring vectors n i and n j in E with probability V ij /V E . Generate a new random vector n k and multiply the weight by the factor R k ij /V ij , expressed in the random variables chosen to generate the direction of the new vector n k .
4. If n k is outside the veto region, add this new vector to the event which then becomes E = {n 1 , · · · , n i , n k , n j , · · · , n 2 }, multiply the weight by a factor V E /V E and return to Step 2. Otherwise, add the weight factors w and ln 2 | sinφ| to h L and h c at time t, go to Step 1 and start a new event.
In terms of these histograms, the soft function correction reads α s (µ s ) 4π where the factor 1 2 has the same source as in (C.1) and the interjet functions were given in (2.13), (2.15) and (2.17). The anti-quark hard function Hq is equal to the quark function. For the gluon function, there is also a region θq > θ q which is parameterized analogously and gives an identical contribution.
As explained in [10], the soft function for the light-jet mass is directly related to the coft function in Sterman-Weinberg dijet cross section defined in [8]. In Laplace space, we have T i,L · T j,R d d k (2π) d−1 δ(k 2 )θ(k 0 )e −n·k/(τ e γ E ) n i · n j n i · k n j · k θ(n · k −n · k) . (C.6) The evaluation of this expression proceeds along the same lines as for the interjet energy flow case derived in detail in Appendix A. If both emitting partons are in the right hemisphere, the renormalized one-loop result is given by with the measurement function Θ lab L (ŷ,φ) constraining the soft radiation to the left hemisphere, and a function g ij (φ,ŷ) = 1 βM 2β coshŷ + βeŷ tanh y i + βe −ŷ tanh y j − cosφ 2β 2 + tanh y i + tanh y j + sechy i sechy j sinφ sin(φ i − φ j ) . (C.8) If one of the two partons is on the left, the function has a collinear divergence, which can be subtracted, as detailed in Section 3. The subtraction was given in (3.11).

D Monte Carlo algorithm for the jet mass distribution
In this appendix we provide the Monte Carlo algorithm used for jet mass resummation, which is also applicable for other non-global observables with soft-collinear double logarithms. Compared to interjet energy flow, we need to subtract the global anomalous dimension and the one-loop global soft function. As for the interjet energy case, we fill three histograms: h U contains the LL evolution, h L the coefficient of the logarithm of the soft function (C.7) and h c its non-logarithmic part. The algorithm for evolving H q 1 to lower scales involves the following steps: 1. Start at evolution time t = 0 from an initial event E with vectors {n, n 1 } and weight w = 1.
2. Generate a random time step ∆t according to the probability distribution P E (t) = V E exp(−V E ∆t), and insert the event weight w into the histogram h U at time t + ∆t.
3. Choose a dipole associated with a pair of neighbouring vectors n i and n j in E with probability V ij /V E . Generate a new random vector n k and multiply the weight by the factor R k ij /V ij , expressed in the random variables chosen to generate the direction of the new vector n k .
4. If n k is in the right hemisphere, add this new vector to the event so that E = {n, · · · , n i , n k , n j , · · · , n 1 }, multiply the weight by a factor V E /V E and return to Step 2. If n k is in the left hemisphere and was emitted from dipole (n, n j ), we need to subtract the global one-loop soft function S q (1) G in equation (3.11). This is achieved with the weight factors 1 − X(φ,ŷ) w and ln 2 | sinφ| g 0j (φ,ŷ) 1 − X(φ,ŷ) w , (D. 1) which are added to the histograms h L and h C at time t. After filling the histograms go to Step 1 and start a new event. Otherwise, add the unsubtracted weight factor w and ln 2 | sinφ| to the respective histograms, go to Step 1 and start a new event.
The quantity V E denotes the subtracted global anomalous dimension V E = V E − V 0 , where V 0 is the large-N c result for the subtraction (3.5) obtained by replacing the Casmir operator C i in this equation by N c /2 for a quark jet, or N c for a gluon jet, respectively.

E Ingredients for jet mass resummation
For convenience, we collect here the perturbative results for ingredients used in the resummation formula for jet mass distribution. The evolution factors at NLL accuracy are given by The LO integrated jet mass distribution is written as The integrated light-jet mass distribution is trivial at this order because the light jet has zero mass for three partons.