Leading jets and energy loss

The formation and evolution of leading jets can be described by jet functions which satisfy non-linear DGLAP-type evolution equations. Different than for inclusive jets, the leading jet functions constitute normalized probability densities for the leading jet to carry a longitudinal momentum fraction relative to the initial fragmenting parton. We present a parton shower algorithm which allows for the calculation of leading-jet cross sections where logarithms of the jet radius and threshold logarithms are resummed to next-to-leading logarithmic (NLL$'$) accuracy. By calculating the mean of the leading jet distribution, we are able to quantify the average out-of-jet radiation, the so-called jet energy loss. When an additional reference scale is measured, we are able to determine the energy loss of leading jets at the cross section level which is identical to parton energy loss at leading-logarithmic accuracy. We identify several suitable cross sections for an extraction of the jet energy loss and we present numerical results for leading subjets at the LHC. In addition, we consider hemisphere and event-wide leading jets in electron-positron annihilation similar to measurements performed at LEP. Besides the average energy loss, we also consider its variance and other statistical quantities such as the KL divergence which quantifies the difference between quark and gluon jet energy loss. We expect that our results will be particularly relevant for quantifying the energy loss of quark and gluon jets that propagate through hot or cold nuclear matter.


Introduction
Highly energetic jets play a major role at high-energy collider experiments such as the Large Hadron Collider (LHC) and the Relativistic Heavy Ion Collider (RHIC), as well as the the future Electron-Ion Collider (EIC) [1]. In the past years significant progress has been made in performing high-precision calculations for exclusive and inclusive jet production as well as jet substructure observables. Aside from being a means to constrain parton distribution functions (PDFs) of the proton [2][3][4][5][6], an integral part of searches for new physics [7][8][9], and a sensitive probe of the strong coupling constant [10,11], another fundamental concern in these studies is how exactly energy is distributed into the states registered in the detector. These states can be considered at different levels of resolution, from the irreducible individual hadrons, to large radius jets which may or may not have multi-prong substructure [12][13][14]. The ability to resolve the final state of a collision at multiple scales is critical in being able to test our understanding of the dynamics that lead to these states. For example, the inclusive jet cross section pp → jet + X has been calculated to nextto-next-to leading order (NNLO) [15,16], and is an important observable to constrain the gluon PDF. An inclusive jet sample is obtained by measuring the transverse momentum p T of all the jets in a given rapidity range. The factorization in QCD can be formulated in terms of hard-scattering functions and (semi-)inclusive jet functions [17][18][19][20][21]. The formation and evolution of jets described by the inclusive jet function is illustrated in the left panel of Fig. 1. Here all jets are taken into account that are produced by the QCD fragmentation process and which are identified with a given jet algorithm. The jet functions satisfy DGLAP evolution equations which allow for the resummation of logarithms of the jet radius R. One can tune R to capture various stages of the shower, and eventually as R → 0, jet production would merge into the traditional observable of inclusive hadron production [22].
However, inclusive jet production forms only one part of the set of observables one can probe in a fragmentation process, where one wishes to know the dynamical means by which the object of concern (the initial quark or gluon) is randomly broken up. Asking more differential questions about the fragmentation process, or probing more exclusive observables, can reveal the underlying mechanism of fragmentation based on general considerations of probability theory alone [23]. In QCD scattering, the object we are concerned with is the total momentum in the underlying hard process, and how the resulting fragments are possibly labeled according to polarization or flavor composition. While in vacuum QCD, the underlying dynamical process of fragmentation can be claimed to be qualitatively and even quantitatively understood, the propagation of partons through a strongly interacting medium has required a more careful theoretical treatment.
Thus it is critical to move beyond the consideration of inclusive jets, and in this work we focus on leading jet production. That is, we consider the cross section when only the leading jet is measured in a given rapidity interval per event (analogous to the largest fragment considered in [23]). The corresponding leading jet function only takes into account the formation and evolution of most energetic jet resulting from an active parton which is illustrated in the right panel of Fig. 1. The renormalization group (RG) equation turns out to be a non-linear DGLAP-type evolution equation which was first introduced in [18] at leading-order (LO) and leading-logarithmic (LL) accuracy using a generating functional approach. In [24], these results were extended by including the full fixed order jet function at next-to-leading order (NLO). In addition, the jet functions were incorporated using a complete factorization formula at NLO which was obtained within Soft Collinear Effective Theory (SCET) [25][26][27][28][29]. Here we further extend the work of Ref. [24] by evolving the entire NLO jet function using a parton shower Monte Carlo approach and we include the resummation of threshold logarithms [30,31] which dominate the cross section when the momentum fraction carried by the leading jet relative to the initial parton approaches unity. In addition, we focus specifically on cross sections where an additional reference scale Q is measured such that we can directly measure the momentum fraction of the leading jet relative to Q. Vital to our approach is that the leading jet functions constitute normalized probability densities, even outside the Sudakov region. Thus the leading jet is a (theoretically) well-defined object of the event, whose evolution we can track. The probability distribution allows us to calculate the mean and variance of this distribution. The mean corresponds to the average energy contained inside the leading jet relative to the fragmenting parton i which we denote by z 1i . Correspondingly, z i,loss = 1 − z 1i is the average out-of-jet radiation or the leading jet energy loss.
Given that leading jets form a well-defined object distributed probabilistically by the fragmentation process, they allow for a well-defined notion of jet energy loss at the jet function and cross section level. We identify the following three criteria that allow for a meaningful definition of the energy loss distribution and its average: • Different than inclusive jets, the leading jet constitutes a well-defined object which has lost energy relative to the initial parton due to out-of-jet emissions. The corresponding jet functions are normalized probability densities which allow for a perturbative evaluation of the (average) energy loss. We will also discuss a possible extension of the present work to leading hadrons which also allows for a well-defined, but nonperturbative definition of energy loss. We stress that it is not possible to construct the corresponding probability density for inclusive jets since the number of inclusive jets is not fixed but generated dynamically event-by-event through the QCD fragmentation process.
• To quantify the lost energy, we not only need to know the energy of the leading jet but also a reference scale Q with respect to which we define the energy loss. We consider different observables where the reference scale is given for example by the center-of-mass (CM) Q = √ s in e + e − collisions. Other examples include jet substructure measurements, Semi-Inclusive Deep-Inelastic Scattering (SIDIS) or photon/Z-jet correlations.
• Lastly, we require that the measured jet energy loss at the cross section level agrees with the (average) parton energy loss at LL accuracy. Higher order effects give corrections to this direct relation which, however, are calculable order-by-order in perturbation theory. The analogy between parton and jet energy loss here is similar to the identification of the variable Bjorken x B in Deep Inelastic Scattering (DIS) and the parton momentum fraction x at LL accuracy.
The concept of jet or parton energy loss has played an important role in theoretical calculations of jet quenching in heavy-ion collisions [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49]. Typically, the notion of energy loss is defined in the soft gluon approximation where flavor-changing processes are suppressed, and calculations are performed in the lowest non-trivial fixed order or resummed expansions. Ideally, one would like to measure the energy of the parton before and after the interaction with the quark gluon plasma (QGP). The difference is the energy loss due to vacuum and medium-induced emissions and allows for the extraction of properties of the QGP. However, the concept of an energetic parton that exits the hard interaction, losing radiation only due to soft emissions, that then emerges from the scattering is tenuous even in the pure vacuum evolution case. The work presented here using leading jets provides the closest connection to this idealized scenario of energy loss measurements. In the vacuum, the parton/jet energy loss is calculable perturbatively, and nonperturbative effects may be modeled via shape functions. In the medium the average energy loss can be determined experimentally and compared to theoretical model calculations, and then tied to the underlying physics of the QGP.
Beyond simply calculating the average energy loss (which corresponds to the mean of the leading jet distribution), since we have the full probability distribution at hand, we also consider for the first time the variance of the jet energy loss that characterizes eventby-event fluctuations. We analyze in particular the change with the jet radius parameter R for these moments of the leading jet distributions, and we explore differences between the quark and gluon energy loss. Moreover, we introduce additional statistical quantities to help understand the differences between leading quark and gluon jets and their fragmentation processes. For example, we compute the Shannon entropy and KL divergence. The latter quantifies the difference between the quark and gluon leading jet distributions. In addition, we present receiver operator characteristic (ROC) curves for different values of the leading jet radius. This represents a first step toward assessing the potential impact of this observable for quark/gluon jet tagging.
A reliable quantitative understanding of leading jets in the high-energy collider experiments is also necessary for the calibration of the jet energy using Z/γ-tagged jets [50,51]. See also Refs. [52][53][54]. In addition, it can be advantageous from an experimental point of view to measure additional quantities such as jet substructure observables on the leading and first subleading jet in the event instead of an inclusive jet sample. See for example Refs. [55][56][57]. A reliable evaluation of the quark/gluon fractions requires corresponding factorization formulas analogous to those developed in this work.
One of the main novelties of our work is the development of a Monte Carlo parton shower framework which solves the non-linear DGLAP-type evolution equations of leading jets, while including the complete (to NLL ) threshold resummed hard-scattering and jet functions. We find full agreement with analytical results for inclusive jets, validating that the parton shower framework which allows for a systematic extension to leading jet cross sections. In addition, we present numerical results for leading partons where we run the shower down to the nonperturbative scale ∼ 1 GeV which is a first step toward understanding the fragmentation spectrum of leading hadrons. We leave more detailed studies of leading hadrons for future work. While the shower introduced here should be considered as a single or few purpose Monte Carlo event generator [58][59][60], we expect that it allows for systematic extensions to other observables consistent with analytical results obtained within QCD factorization. See Refs. [61][62][63][64][65][66] for recent developments of parton shower algorithms.
The remainder of this work is organized as follows. In section 2, we introduce the main theoretical concepts of leading (and subleading) jets and compare them to inclusive jet production. We discuss the evolution equations of leading jets and factorization formulas at fixed order and in the threshold limit. In addition, we discuss the connection between leading and subleading jets and inclusive single-, di-and tri-jet functions. In section 3, we discuss the setup of the parton shower Monte Carlo framework and we present first numerical results. We discuss the resummation at LL accuracy and the extension beyond LL by including (threshold resummed) hard and jet functions in the shower algorithm. In section 4, we derive the threshold resummed hard and jet functions for e + e − hemisphere leading jets and leading subjets. We discuss how nonperturbative effects can be included in the threshold limit which is phenomenologically relevant for leading jets and present numerical results for both processes at the cross section level. In section 5 we calculate the average leading jet energy loss and the variance at NLO by taking moments of the leading jet function. We present numerical results for the mean and variance of the leading jet/energy loss distribution and focus in particular on quark/gluon differences. In addition, we present numerical results for the Shannon entropy and the KL divergence. In section 6, we study the discrimination power of leading (sub)jets for quark/gluon jet tagging. In section 7, we discuss further applications of our framework such as event-wide leading jets in e + e − collisions similar to existing data from LEP. In addition, we consider leading jets in SIDIS and photon-jet correlations in proton-proton collisions. Both processes also allow us to perform jet energy loss measurements. Lastly, we present results for leading partons at the nonperturbative scale. We conclude in section 8.

Fixed order, evolution and factorization
We start by reviewing the NLO jet function, its evolution and the factorization formula for inclusive jets pp → jet + X in section 2.1. Correspondingly, we discuss the leading jet cross sections pp → jet 1 + X in terms of their NLO jet function, and evolution equations in section 2.2. In addition, we introduce the relevant jet functions for subleading jets. Subsequently, we discuss the structure of relevant factorization formulas in section 2.3. In section 2.4, we extend the inclusive jet function for single jets to di-and tri-jet functions and discuss their relation to leading and subleading jet functions. Here we refer to Ref. [23] where these relations were proposed in the context of statistical properties of randomly broken objects and spin glasses.

Review of inclusive jet production
The inclusive jet function was calculated in Refs. [19][20][21] for k T -type algorithms and in Ref. [67] for cone algorithms. See also [68,69] for the extension to massive quarks. We can write the inclusive jet function J i (z, QR, µ) for i = q, g in terms of the momentum fraction of the inclusive jets z, the jet radius R, the large reference scale Q and the renormalization scale µ. For inclusive jet production in proton-proton collisions the large reference scale is given by the transverse momentum of the initial fragmenting parton Q =p T = p T /z. Here p T denotes the transverse momentum of the final observed jet. We note thatp T is not an observable quantity in proton-proton collisions which is why we convolve the jet function with a hard-scattering function in the factorization formula of the cross section. The range of the convolution integral is determined by the allowed range ofp T which is determined by the jet's transverse momentum and rapidity p T , η and the CM energy √ s, see Eqs. (2.5) and (2.6) below. In subsequent sections we will focus on cross sections where we have access to the initial scale Q which is necessary to define the lost energy of a leading jet as mentioned in the Introduction. Therefore, we keep the general notation here and denote the hard reference scale by Q. For quarks and gluons we find the following results at NLO The inclusive jet function measures the momentum fraction z of all jets that are produced from the initial parton. The evolution equations of the inclusive jet functions are the standard timelike DGLAP evolution equations [17,18,20,21] µ d dµ 3) similar to fragmentation functions. Here P ji (z) denote the time-like Altarelli-Parisi splitting functions which allow for a perturbative power expansion in terms of the QCD strong coupling constant The inclusive jet functions are then evolved from their characteristic scale µ ∼ QR, which removes logarithms of the jet radius in Eq. (2.1), to the hard scale µ ∼ Q. This DGLAP evolution resums large logarithmic corrections of the jet radius to all orders. The factorization formula for the inclusive jet cross section pp → jet + X differential in the transverse momentum p T and rapidity η of the jet is given by Here f i,j denote the parton distribution functions (PDFs) for the two incoming protons and H ijk are the hard-scattering functions ij → k similar to inclusive hadron production. They are known analytically at NLO [70,71]. For hadron production, the same factorization formula as in Eq. (2.5) can be used except that the perturbative jet functions are replaced by fragmentation functions. The factorization in Eq. (2.5) holds up to power corrections O(R 2 ) which are usually found to be small even for large values of R [24,72]. The first two symbols ⊗ in Eq. (2.5) denote appropriate integrals over the momentum fractions x i,j . The third integral denoted by ⊗ is a convolution integral in terms of the momentum fraction z. To make this structure more explicit, we rewrite Eq. (2.5) as where we absorbed the PDFs and corresponding integrations over the variables x i,j in the new hard functions ffH i . Here the lower integration limit is given by z 0 = 2p T / √ s cosh η. Note that in Eqs. (2.5) or (2.6) we do not have access to the partonic transverse momentum p T which appeared in the jet functions in Eq. (2.1). Therefore, we write both the hard and jet function in terms of the observed jet transverse momentum p T instead of the initial partonic transverse momentum using the relationp T = p T /z as in Ref. [20]. This is valid as long as we are not in the regime where z 1 which would require an additional resummation of small-z logarithms. See Refs. [73,74] as well as earlier work in Refs. [75][76][77][78][79] We end this section by summarizing some key features of the inclusive jet functions. Similar to PDFs and fragmentation functions, the inclusive jet functions constitute number densities in the sense that an integral over the jet function (first Mellin moment) yields the event averaged number of jets N jet which originate from the fragmenting parton. This number is generated dynamically through the QCD fragmentation process and depends on the jet radius, the jet algorithm and the scale Q. Analogously, the same integral over fragmentation functions gives the average number of hadrons or the hadron multiplicity. The first moment (second Mellin moment) is related to momentum conservation in the sense that the initial scale Q has to be recovered in the inclusive jet sample that is produced resulting from the initial parton. Of course, in practice not the entire inclusive jet sample resulting from a highly energetic quark or gluon is reconstructed by the experiment. The limited range of detectors is taken into account in the factorization formula in Eq. (2.5) which depends on the jet rapidity η. At the level of the jet function, we thus have the following two sum rules which hold for quarks and gluons separately. Note that the quantity N i,jets introduced here is related to the entropy of a jet [80,81]. In order to evaluate the first Mellin moment analytically, the resummation of small-z logarithms is required as mentioned above. Related experimental results can be found in Ref. [82] where the subjet multiplicity was measured. The momentum sum rule for an inclusive jet sample in Eq. (2.8) is illustrated on the left side of Fig. 1 where three jets are reconstructed resulting from an initial fragmenting quark. In this case, the momentum fractions of the three jets have to add up to unity z 1 + z 2 + z 3 = 1. An important aspect is that the momentum sum rule in Eq.

Leading and subleading jet functions and their evolution
We are now going to introduce the jet functions for leading and subleading jets analogous to the inclusive case discussed in the previous section. The LO and LL resummation for leading jets was first introduced in Ref. [18]. The O(α s ) correction of the leading and subleading jet functions was first discussed in Ref. [24]. Here we denote the leading jet function by J i (z 1 , QR, µ) which depends on the momentum fraction z 1 of the leading jet relative to the initial scale Q of the fragmenting quark or gluon, see the right panel of Fig. 1. Analogously, the jet function that describes the formation of the leading and the first subleading jet is given by J i (z 1 , z 2 , QR, µ). It depends on the momentum fractions z 1 and z 2 of the leading and first subleading jet, respectively. It is a more differential version of the leading jet function. Note that throughout this work we write jet and hard-scattering functions associated with leading and subleading jets in script font to distinguish them from the corresponding inclusive jet quantities. We start with the fixed order calculation of the leading and subleading jet function. At LO there is just one parton and we trivially find At NLO there are at two partons which can be clustered into a single jet or two separate jets depending on their distance and the jet algorithm. If both partons are clustered into the same jet, the contribution is the same as for inclusive jets and it is proportional ∼ δ(1−z 1 ). If the two partons are clustered into separate jets, we only take into account the jet which contains the more energetic parton. Instead, for inclusive jets we always take into account both jets independent of how energetic they are. It turns out that at NLO we can write the leading jet functions in terms of the inclusive ones by including a theta function which requires z 1 > 1/2. We thus find for the leading and subleading jet functions up to NLO Note that the corresponding jet functions J i (z 1 , . . . , z n , QR, µ) which take into account the dynamics down to the n-th leading jet can be constructed in a similar way. We also note that the leading jet function is obtained from the subleading jet function upon integration 1 0 dz 2 J i (z 1 , z 2 , QR, µ) = J i (z 1 , QR, µ) . (2.14) The leading jet function at NLO is only non-zero for z > 1/2 where it agrees with the inclusive jet function. At next-to-next-to-leading order (NNLO), we need to consider three particles in the final state which gives a lower bound for the leading jet function of z > 1/3. In general, at N n LO, the minimal non-zero value of the leading jet function is thus given by 1/(n + 1). Next, we consider the evolution equation of the leading jet functions. Different than for inclusive jets, the evolution equations for leading and subleading jets are non-linear. At NLO there are only two partons and it is sufficient to follow the leading parton at the 1 → 2 splitting. If there are subsequent splitting processes as well, one needs to know the value of the leading jet function at each branching point. This is illustrated in Fig. 2 which shows an exemplary branching tree, where the green line leads to the leading jet which carries a momentum fraction z 1 . In order to obtain the correct path, we need to known at each branching point the value of the branching fraction z of the splitting i → jk as well as the leading jet functions J j (z j1 ) and J k (z k1 ) of each branch. This feature of leading jets makes Figure 2. Illustration of an exemplary branching tree: Starting from parton i at the scale Q, four partons/jets are produced which carry a momentum fraction z n relative to the initial scale. In order to determine the correct path (green) that leads to the leading jet, we need to know the momentum fractions and leading jet functions of each branch at a given splitting i → jk. This results in a non-linear evolution equation for observables involving leading jets. Instead, for inclusive jets we only need to know the final momentum fractions z of each jet irrespective of the they compare to the rest of the event since we sum over all contributions. the evolution equations non-linear. Instead, for inclusive jets we only need to sum over all possible paths of the branching tree and we obtain the usual linear DGLAP evolution equations, see section 2.1 above. The measurement of leading jets requires knowledge of all the other jets when the evolution terminates. Only if we know all other jets we can determine which jet is the leading jet. Instead, for an inclusive measurement we simply sum over all possible jets at the end which does not require simultaneous knowledge about the other jets that are produced.
We can write the non-linear evolution equation for the leading jet functions as 1 Here we follow Ref. [83] and use the notation P i→jk (z) = P (0) ji (z) for the LO Altarelli-Parisi splitting functions and we regulate both endpoints 1/(1 − z) + and 1/z + . We thus have P q→qg (z) = P q→gq (1 − z) and P g→gg (z) = P g→gg (1 − z). We can then check explicitly that the leading jet functions in Eq. (2.12) satisfy the non-linear evolution equations. Similar non-linear evolution equation were obtained in the context of fractal jet substructure observables [83], the jet charge [84,85], track functions [86,87] and di-hadron fragmentation functions [88][89][90].
We also note that the evolution equations here are different than for the central subjets or the jet shape with the winner-take-all axis (WTA) [91] which was considered in Refs. [67,92,93]. Even though the fixed order expressions of the jet functions contain a similar term ∼ Θ(z > 1/2)P ji (z), the use of the WTA axis leads to a linear DGLAP-type evolution equation with modified kernels.
Similar to the evolution equation for the leading jet functions in to Eq. (2.15), we can write down an evolution equation for the subleading jet functions. We have Different than for inclusive jet functions in Eq. (2.7), we find that the leading jet functions are normalized to unity (first Mellin moment). The first moment (second Mellin moment) corresponds to the average energy or transverse momentum fraction z 1 which is contained in the leading jet. Recall that for inclusive jets, the second Mellin moment is unity due to momentum conservation, see Eq. (2.8). We thus find the following expressions for the leading jet functions In order to interpret the leading jet functions as probability densities for the leading jet to carry the longitudinal momentum fraction z 1 , they need to be i) normalized to unity and ii) positive for all values of z 1 . The first requirement is satisfied as can be seen from Eq. (2.17). For example, we can check this result using the NLO expression. In addition, this normalization is conserved under the non-linear evolution, i.e., we have which follows from Eq. (2.15). The conservation of the first Mellin moment for leading jets is analogous to the conservation of the second Mellin moment for inclusive jets, see Eq. (2.9) and only holds when the jet function is written in terms of Q =p T instead of the transverse momentum of the observed jet. The second requirement that the leading jet functions need to be positive for all values of z 1 can be violated at finite perturbative accuracy. For example, at NLO the leading jet function can become negative especially for large values of z 1 and small values of the jet radius R. However, we observe that this problem can be solved by including threshold resummation which we discuss in more detail in section 4 below. Therefore, we can indeed treat the leading jet functions as probability densities which allow us to calculate the expectation value or the average energy contained in the leading jet z i1 as in Eq. (2.18) above. We note that probability densities are rather unusual in this context since we typically obtain number densities as for inclusive jets or inclusive hadrons. Similar sum rules as in Eqs. ( Having defined the average energy contained in the leading jet z i1 , we can now define the average energy loss of a leading jet which is given by all the energy which is not contained in the leading jet We discuss z i,loss and other statistical quantities that quantify the probability distribution of jet energy loss in section 5 for both the jet functions and at the level of cross sections. We also would like to stress that the notion of a leading or subleading jet does not directly correspond to the leading jet in the entire event in proton-proton collisions which is not possible due to the finite detector coverage and the incoming particle beams. Instead, we calculate the leading jet in a given rapidity interval. This dependence is taken into account in factorization formulas which we discuss in the next section similar to inclusive jets. We can now extend Eq. (2.18) to subleading jets. For example, we can use the subleading jet function to calculate the the average energy which is contained in the first subleading jet as Here, J i2 is the jet function of the second leading jet which is indicated by the additional subscript 2. It constitutes the probability density for finding the first subleading jet with momentum fraction z i2 and we can thus calculate the average quantity z i2 . Note that the full subleading jet function is needed for the evolution in Eq. (2.16). However, the information of the leading jet momentum fraction z i1 can be integrated out after the evolution and we can calculate the probability density for the first subleading jet. Analogous equations hold for the n-th leading jet The average values of the energy fractions contained in the subleading jets is also relevant for studies of the jet energy loss. They contain the information of how the lost energy of the leading jet z i,loss , as defined in Eq. (2.20), is distributed amongst the subleading jets. We have We note that the non-linear evolution equations for leading and subleading jets in Eqs. (2.15) and (2.16), respectively, are difficult to solve analytically. In Ref. [18,24] an iterative approach was used starting with the LO jet function as the initial condition of the evolution. However, this approach is impractical when the whole NLO jet function is evolved and when additional threshold logarithms are resummed to all orders which is discussed in more detail in section 4 below. In addition, in order to calculate the (sub)leading jet spectrum at small-z many iterations would need to be computed analytically. Ref. [18] also explored a strictly leading order parton shower method. However, we introduce a new Monte Carlo approach to solve the non-linear evolution equations for leading and subleading jets which we discuss in detail in section 3, that allows for the incorporation of threshold resummation effects and jet function contributions, moving beyond the LO process. Different than general purpose Monte Carlo event generators our approach is a single purpose (or few purpose) parton shower which is designed to specifically solve the above evolution equations. It allows for a well defined perturbative accuracy similar to analytical calculations and it allows for systematic improvements to yet higher perturbative accuracy. Before introducing the parton shower algorithm, we first discuss relevant factorization formulas for leading and subleading jets in the next section.

Factorization
Different than for inclusive jets, see Eq. (2.5), the form of the factorization formula for leading jets pp → jet 1 + X depends on the perturbative accuracy we are working at. The structure of the hard functions changes and additional jet functions need to be included as the perturbative accuracy is increased. Here we consider the cross section differential only in transverse momentum of the leading jet p T 1 for η = 0. The extension to η = 0 is straightforward since only the hard functions depend on the rapidity and the overall theta functions which enforce momentum conservation. Following Ref. [24], we can write the leading jet cross section at LO and LL accuracy as where we sum over all contributing channels ij. Note that we denote the appropriate LO hard-scattering functions by ff H ij describe the production of two partons with momentâ p T i,j in a hard-scattering event which are back-to-back at LO in the transverse plane. In principle we could also write the LO hard function only as a function ofp T i =p T j . The hard functions also include the initial state parton distribution functions and appropriate integrals over the momentum fractions of the incoming partons. The jet functions J i,j take into account the formation and evolution of the two leading jets which originate from partons i and j, respectively. The delta function in the last line of Eq. (2.26) then picks one of the two leading jets from partons i, j. The one with the larger transverse momentum is measured and denoted by p T 1 . The two theta functions in the last line ensure momentum conservation. We can also rewrite Eq. (2.26) in terms of the measured jet transverse momenta instead of the partonic quantities p T i = z i1pT i which is valid as long as we are sufficiently far away from the region where z i1 1. We find which is closer to the factorization formula for inclusive jets in Eq. (2.5). In either case the evolution equations only resum logarithms of the jet radius. Next, we consider the leading jet cross section with NLO hard-scattering functions. At NLO a third hard parton can be radiated which requires us to introduce an additional jet function. The final jet with the highest transverse momentum can result from the fragmentation process of either of the three hard partons produced in the hard-scattering process. We thus have ijk denotes the NLO hard function with three hard partons ijk which subsequently fragment into jets. This result can be generalized to higher orders.
Next, we consider the cross section where we not only measure the transverse momentum of the leading jet but also the first subleading jet. We denote the transverse momentum of the fist subleading jet by p T 2 . The result obtained here can also be extended to the measurement of further subleading jets. At LO and LL, we find the following result Here the leading jet p T 1 is given by the transverse momentum of the leading jet originating either from parton i or j as above. The transverse momentum of the first subleading jet p T 2 is given by the smaller one of the leading jets from partons i and j or one of the subleading jets. The theta functions in the last line do not change as they are written in terms of the initial partonic transverse momentum. The factorization for leading and subleading jets in Eq. (2.29) can be generalized to higher perturbative accuracy similar to Eq. (2.28).
It is instructive to consider how the factorization for inclusive jets is recovered by summing over the leading jet and all subleading jets. After carrying out the sum over all jets, the factorization structure simplifies significantly and it has the same structure to all orders. To make this connection more explicit, we work with a factorization formula as in Eq. (2.26). The delta function in the last line of that equation, which specifies the measurement, needs to be replaced by a sum over delta functions that measure the transverse momentum of all jets. Schematically, for two initial partons we have where the first two terms which are written out explicitly correspond to the two leading jets originating from partons i and j. As an example, we consider only the first delta function. We find that we can rewrite the corresponding contribution to the inclusive cross section as where the lower integration limit here is given by z 0 i1 = 2p T / √ s as expected from Eq. (2.5) for η = 0. The jet function J i in the second and third line is written in terms of the observed jet p T instead of the parton transverse momentum similar to Eq. (2.27). After performing the integral over z j andp T j and implicitly defining ff H i , the last equation has the typical convolution structure as it is found for inclusive jets. Similar steps hold for the other jets that are produced and eventually we find The right hand side of the equation is written in terms of the inclusive jet and hard functions of Eq. (2.5). Note that if we consider the inclusive cross section where only the two leading jets are measured for every event, it would approximate the inclusive jet cross section at high p T but differ at low p T where subleading jets are important. The factorization of the corresponding cross section would involve the subleading jet functions as in Eq. (2.29). These kind of differences may also be relevant in the context of assessing QCD scale uncertainties of jet cross sections, see for example [94,95]. We also note that a factorization formula similar to Eq. (2.29) is relevant for the precise calculation of jet substructure observables in order to compare to recent measurements from ATLAS [55,56] and CMS [57]. In this case the jet substructure measurements were performed only on the leading and the first subleading jet. In addition, a requirement was imposed on the ratio of the transverse momenta of the two leading jets. These constraints affect the quark/gluon fractions which are not be properly taken into account when a factorization formula for inclusive jets is used instead.

Leading and subleading jets from inclusive single-, di-and tri-jet functions
The (single-)inclusive jet function J i (z) was already introduced in section 2.1. Analogously, we can introduce the inclusive di-and tri-jet functions which we denote by (2.33) The inclusive di-jet function is similar to di-hadron fragmentation functions considered for example in Refs. [88][89][90]. Analogously, we denote the inclusive n-jet function by J i (z 1 , . . . , z n ). Different than for the leading and subleading jet functions no ordering of the momentum fractions z i is implied. For example, at NLO we can write the di-and tri-jet functions in terms of the (single-)inclusive jet function Note that the NLO di-jet function in Eq. (2.34) is similar to the NLO subleading jet function in Eq. (2.13) but without the theta function Θ(z 1 > 1/2). Higher order results for the diand tri-jet functions can be obtained with the parton shower which will be introduced in the next section. In the following, we leave the arguments QR, µ of the jet functions implicit. Following Ref. [23], we should be able to compute the leading and subleading jet functions from the inclusive n-jet functions. If the observed jet has a momentum fraction z > 1/2 it is the largest jet. We thus have In the range of 1/3 < z < 1/2 the inclusive jet function is given by the sum of the leading and subleading jet functions.
Analogous to the notation introduced in Eq. (2.23) we include a subscript n to denote the jet function which depends only on the momentum fraction of the n-th leading jet. For example, we have More generally, we have Here we limit ourselves to confirm the observation made in Ref. [23] in the range of 1/3 < z < 1/2. For lower values of z a similar analysis is possible as outlined in Ref. [23]. We have And the leading jet function is given by These relations were derived in [23] from general considerations. In a general fragmentation process, one would need to know the whole set of leading, sub-leading, or sub-sub-leading fragment functions, etc., or the whole suite of inclusive multi-fragment distribution functions to completely characterize the fragmentation process, unless a more dynamical rule for generating these functions can be found. In the jet case, the parton shower constitutes such a dynamical rule, at least to leading logarithmic accuracy.

Monte Carlo setup
In this section, we introduce the new Monte Carlo setup which solves the non-linear evolution equations of leading and subleading jets including the threshold resummed hard and jet functions. We start by presenting the Monte Carlo setup at LO/LL accuracy in section 3.1 following Ref. [18]. For inclusive jets we compare the Monte Carlo result to an analytical solution of the DGLAP evolution equations in Mellin space where the contour integral of the inverse transformation is performed numerically. In section 3.2, we then discuss in detail how the Monte Carlo code can be extended beyond LL accuracy by including the (threshold resummed) hard and jet functions which brings the perturbative accuracy of the shower to the same level as analytical calculations of inclusive jets.

The parton shower at leading log
At LO/LL order we follow the parton shower setup of Ref. [18]. See also Ref. [96]. We have at any time a set of partons S that represents the state of the shower: Here, z i is the energy fraction carried by the i-th parton, and f i = q,q, or g is its flavor. We choose a uniformly distributed random number rnd ∈ [0, 1] and solve for a time step ∆t according to The P f i denote the final state summed splitting functions for a quark to split to a quark and gluon (f i = q orq) or a gluon to split to anything (f i = g).
Here is a small parameter which cuts off the integral at both endpoints. The sum in the exponent in Eq. (3.2) runs over all the generated partons in the event. We then advance the Monte Carlo time t → t + ∆t and solve for R in That way the MC time and R = tan(R/2) are ordered variables and the shower terminates when t > t max = t(Q, R). For strict LO/LL comparisons, we take K = 0, otherwise, K is a factor which makes the threshold contribution consistent with the two-loop cusp, using the so-called CMW trick [97]. The minimum angle R is set by the jet radius R, which is where the shower algorithm terminates. We then pick a parton to split with a probability given by its relative contribution to the decay rate of Eq. (3.2). The momentum fraction z of the splitting of the chosen parton is sampled according to the appropriate DGLAP splitting functions, and in the case of a splitting gluon, the flavor is assigned based on the relative contribution of the final state particles to the gluon's phase space. We now histogram the inclusive momentum fraction z of the partons produced by the shower relative to the initial quark or gluon. The result is the LO evolved jet function. To establish the consistency of the parton shower result and the analytical Mellin space evolution of Ref. [20], we compare the two results in Fig. 3. As an example, we show the results for Q = 91.2 GeV and five exemplary jet radii R = 0.05 − 0.8. Here we use the LO expression for the running coupling constant and α s (M Z ) = 0.1187. We observe that the two results agree exactly for the five jet radii and over the entire range of z.

Extension beyond leading-logarithmic accuracy
In order to improve the accuracy of the shower, it is helpful to examine the factorization theorem for inclusive jet production. We write Here C i represents a hard-scattering/coefficient function of a particular process. We can then describe how the parton shower generates U . First, we introduce: • Notation: Let f (S) i R denote the shower average on partonic states S of some function f on those states at the angular scale R, when the shower was initiated by the parton of flavor i. S = {{z 1 , f 1 }, {z 2 , f 2 }, ...., {z n , f n }} is the final state of the shower given above, which is the stochastic variable being averaged over.
Then the DGLAP evolution kernel U (x; R; Q) to leading logarithmic accuracy, evolved from a hard scale Q down to QR, is given by (3.6) Where we have introduced the energy fraction z k of a parton in the event. In practical algorithmic terms, what Eq. (3.6) actually means is that we start with a histogram H ij x , indexed by the momentum fraction x ∈ [0, 1], and labeled by the initial and final partons i, j. Each x falls into some bin of size ∆(x) 1 (which can also depend on where we are in the distribution). Initialize H x = 0. Then: 1. Run shower to get the set of momenta S at angular scale R.

For each
After a sufficient number of events, divide all bins in H x by the number of events. Then H x = U (x; R) in the limit ∆(x) → 0. Note that this effectively approximates the delta function in Eq. (3.6) by the step function: Substituting in Eq. (3.6), we then get (3.10) We wish to implement this last integral completely stochastically in the shower. We proceed as follows. First we introduce the cumulant function and its functional inverse Then we change variables If we can normalize CJ(x) on some interval, so 0 < CJ (x) < 1, then we can take σ to be a random number between 0 and 1, and thus perform this final integral via a Monte Carlo Here [·] is the final averaging process over the random variable rnd ∈ [0, 1] uniformly distributed on the unit interval. The ∝ arises from normalization issues. Given a coefficient and jet function, we can compute the stochastic averaging as follows: we start with a histogram H ij x , indexed by the momentum fraction x ∈ [0, 1], and labeled by the initial and final parton flavors i, j. Each x falls into some bin of size ∆(x) 1. Initialize H ij x = 0. Then: 1. Run shower starting with initial parton i to get the set of momenta S at angular scale R.
2. For each k ∈ S, draw a random number uniformly between 0 and 1, rnd, and calculate CJ −1 if k (rnd) (recall that f k is the flavor of the k-th parton in the shower), then add After a sufficient number of events, divide all bins in H ij x by the number of events. We can then improve the accuracy of the shower by using instead of the leading order coefficient and jet functions their threshold resummed expressions, as described in the next section. The algorithm for leading jet production with threshold corrections is then as follows: 1. Run shower starting with initial parton i to get the set of momenta S at angular scale R.
2. For each k ∈ S, draw a random number uniformly between 0 and 1, rnd, and set x k = z k × CJ thr−1 if k (rnd) (recall that f k is the flavor of the k-th parton in the shower). This gives a set of jet momenta {x 1 , ..., x n }. Then add 1 to the bin in After a sufficient number of events, divide all bins in H ij x by the number of events. We note that the inverse function CJ thr is calculated from the threshold resummation expressions for C and J.
Finally, we note that in principle the NLO DGLAP evolution kernels can be included in a parton shower following Ref. [98]. We leave a combination of the NLO DGLAP evolution and the NLO/threshold resummed hard and jet functions for future work.

Threshold resummation for leading jet observables
For leading jets the threshold region z → 1 [30,31] is phenomenologically important. Therefore, we need to include the corresponding double logarithmic corrections to all orders at NLL accuracy. We note that the jet at threshold is the leading jet -up to power corrections, and so we do not need to consider the difference between the inclusive and leading jets when performing the resummation. We consider first the threshold resummation for e + e − hemisphere jets which will be the standard reference process for the subsequent section. Second, we consider the threshold resummation for subjets in proton-proton collisions. The two processes are illustrated in Fig. 4 and the corresponding threshold resummation is discussed in sections 4.1 and 4.2. Different than event-wide leading jets in proton-proton or e + e − collisions, for both processes considered here we only have one initial parton at LO/LL accuracy. Different than for fixed order calculations of leading jets, as discussed in the previous section, the structure of the factorization does not change order-by-order in the threshold region. By including the threshold resummed hard and jet functions in the Monte Carlo approach discussed above, we can obtain numerical results for leading jets. We discuss nonperturbative effects in section 4.3 which can be included in the threshold region by convolving the perturbative spectrum with a shape function. We then present numerical results for both processes in section 4.4. Parton showers should naturally, at least to leading logarithmic accuracy, resum final state threshold logarithms, such as treated in this paper. We also note that in Refs. [99,100] initial state threshold resummation was also included in the Deductor parton shower.

e + e − hemisphere jets
For e + e − hemisphere jets 3 there is only one initial parton at LO which is different than for e + e − event-wide leading jets which will be discussed in section 7.1 below. Note that for inclusive jets this distinction is not relevant at the level of the factorization theorem since we sum over all possible jets in the final state. Per-event and per-hemisphere inclusive jets only differ in terms of their normalization. We consider the cross section dσ/dz where additional angular dependencies are integrated out and z = 2E/Q is the hemisphere jet energy relative to half of the CM energy Q = √ s. The collinear factorization using ln R resummed jet functions is analogous to the results discussed above and it's form depends on the perturbative accuracy. At LO/LL accuracy we have Therefore, this process is ideally suited to study jet energy loss since all three criteria listed in the Introduction are satisfied. At LL accuracy, we find a direct connection between the leading jet cross section and the average parton energy loss see Eq. (2.18). Here the average energy loss is a function of the energy scale Q and the jet radius R. At higher perturbative accuracy, additional jet functions need to be introduced in Eq. (4.1). However, in the threshold limit z → 1, where additional emissions are soft, we can obtain a closed form of the factorization formula. We derive the threshold resummation in Mellin transform space and then perform the inverse transformation analytically. Throughout this work we adopt the following convention for the Mellin transform and it's inverse In this section we derive the threshold resummed cross section for inclusive e + e − hemisphere jets. The leading jet cross section is then obtained by including the threshold resummed hard and jet functions in the parton shower as discussed in the previous section. In the threshold region, the cross section can be refactorized as Here H i is the hard function [103], J i accounts for the recoiling soft radiation [104,105], S i is a soft-collinear function [106,107] taking into account soft radiation in the direction of the observed jet [108,109] and J i takes into account unmeasured collinear radiation inside the observed jet [110,111]. In addition, S NGL i accounts for non-global logarithms (NGLs) due to correlations between the in-and out-of-jet region [112]. For our numerical results presented below we include the fit to the Monte Carlo algorithm of Ref. [112]. To the accuracy we are working at, the NGL factor can be included multiplicatively. For completeness, we summarize the NLO order expressions of the relevant functions in the Appendix A. For phenomenological applications at LEP, the sum in Eq. (4.5) runs over quark and anti-quark channels. For our numerical results presented below we also consider leading jets initiated by gluons. In e + e − collisions they can be obtained from a hardscattering process with an intermediate scalar φ → gg, where φ could be an intermediate Higgs boson. The natural scales which eliminate the large logarithms of the different functions at fixed order and set the initial scale of their RG evolution are given by (4.6) We note that both within collinear factorization, Eq. (4.1), and when threshold resummation is included as in Eq. (4.5), the parton shower resums logarithms between the scales QR → Q which is illustrated in Fig. 5. The additional resummation of the threshold logarithms is carried out analytically as described in this section and included in the parton shower algorithm. As an example, we consider the RG equation of the soft-collinear function S i . In Mellin space the evolution equation is multiplicative and can be written as The anomalous dimension is given by  The cusp Γ i and non-cusp γ i anomalous dimensions allow for the following perturbative expansion (4.10) For the soft-collinear function S i the non-cusp term vanishes to NLL, γ S,0 i = 0. The relevant terms of the cusp anomalous dimension at NLL are given by T F N f . Instead of performing the inverse transformation numerically using for example the minimal prescription of Ref. [113], we follow here the method introduced in Ref. [114,115] and perform the inverse transformation analytically. Following the notation of Ref. [116], we introduce the functions K i , η i and η i,γ as 14) Evaluating these expressions at NLL, we find with r = α s (µ)/α s (µ 0 ) and similarly for η i,γ . The function η i,γ vanishes for the soft function S i at NLL. We can then write the evolved soft function S i in Eq. (4.12) in Mellin space as Next, we consider the Mellin inverse transformation back to z-space. Note that here we did not make a scale choice for µ S . Therefore, the N -dependence appears both in the factor N −2η i and also the initial condition of the evolution S i (QR/N, µ S ) which we need to take into account at NLL accuracy. The initial condition of the evolved soft function S i is a polynomial of L = ln(µ 2 we can then write the initial condition of the soft function in Eq. (4.18) as With this replacement, the only remaining dependence on the Mellin variable N is the factor N −2η . We calculate the following Mellin transformation (4.21) In the threshold limit we can drop terms which are power suppressed as O(1/N ). With this result we can now write the NLL resummed soft function in z-space as .

(4.22)
In order to implement the threshold resummed jet function in the Monte Carlo code discussed in section 3 above, we need the cumulant instead of the differential result. We take the cumulant to be the integral of Eq. (4.22) from z to 1. We find .

(4.23)
Here the superscript c indicates that S c i is the cumulative result of the soft function. Next, we extend our result to NLL accuracy where we need to include the initial condition S i (−∂ η , µ S ) of the evolved soft function at O(α s ). Taking the first and second order derivatives of Eq. (4.23) with respect to η, we find Substituting the derivatives (−1) m ∂ m η for the operators O m in S i (−∂ η , µ S ), we find the following expression for the resummed cumulant of the soft function at NLL accuracy .

(4.25)
Lastly, we need to calculate the convolution of the resummed soft function and the jet function J i . At NLL, we find . (4.26) The relevant anomalous dimensions can be found in the Appendix A. The extension to NLL can be obtained in analogy to Eq. (4.25) above. Besides the NGLs, the remaining functions in the refactorized expression in Eq. (4.5) are independent of z and satisfy multiplicative evolution equations. Their anomalous dimensions are also listed in the Appendix A.

Subjets
In proton-proton and heavy-ion collisions we do not have access to anx event-wide reference scale like in e + e − collisions to define the energy loss of the leading jet. However, we can construct a reference scale by first identifying an inclusive jet sample with jet radius R. The transverse momentum of the identified jets can be used a reference scale. We then recluster the particles inside a given jet with a jet radius r < R and measure the momentum of the identified leading subjet p r T relative to the initial jet z r = p r T /p T . A calculation of the inclusive subjet distribution within collinear factorization was discussed in [21,67]. See also Refs. [117,118]. We briefly review the factorization for inclusive subjets within collinear factorization and we then extend it to include the resummation of threshold logarithms which can also be implemented in the parton shower framework introduced above. We consider the cross section 27) where p T and η is the transverse momentum rapidity of the inclusive jet with radius R. See also Eq. (2.6). Up to power corrections O(R 2 ), the cross section can be factorized in terms of parton distribution functions f i,j , hard-scattering functions H ijk and a jet function G k (typically denoted by G k e.g. in Ref. [67]) which depends on the jet substructure observable z r . The symbols ⊗ denote appropriate integrals over the involved longitudinal momentum fractions x i,j and the fraction z contained in the observed jet with radius R. We change the index of the jet function to G i from here on. At NLO, the jet function for a quark is given by and similarly for a gluon, see Ref. [67]. Note that here z denotes the momentum fraction contained in the inclusive jet with radius R relative to the initial parton and z r the fraction of the observed jet contained in the subjet with radius r. The logarithm ln(µ 2 /p 2 T R 2 ) in the first line is resummed through DGLAP evolution which is satisfied by the entire jet function G i . In order to resum the logarithm ln(R 2 /r 2 ) in the second line of Eq. (4.28), the jet function G i can be further factorized in terms of a hard matching coefficient and a jet function for the subjet ∼H ij ⊗ J j which was carried out in [67]. Here we extend this calculation and include also the resummation of threshold logarithms which is phenomenologically important similar to the e + e − hemisphere jets discussed above. At threshold z r → 1, the jet function G i section can be refactorized as which allows for the joint resummation of threshold logarithms and ln(R 2 /r 2 ) similar to Eq. (4.5) above. The hard functions H ij can be combined with the remaining functions in Eq. (4.27) which altogether can be considered effective quark/gluon fractions. Here S j is a collinear-soft function -the same as for hadon-in-jet production at threshold which was discussed at NLL accuracy in Ref. [119]. See also Ref. [120]. The soft-collinear function S j and the jet function J j are the same as for e + e − hemisphere jets in Eq. (4.5) except that they are evaluated at the subjet radius r. We note that the NLO expressions of the soft-collinear and collinear-soft funtion are the same up to an overall minus sign (and the different jet radii). For the subjet refactorization at threshold, we find two types of NGLs. They arise independently at the boundary of the initial jet due to a correlation of the functions H ij , S j and at the boundary of the subjet due to the correlation of the functions S j , J j . The fixed order expressions of all relevant functions and their anomalous dimensions can be found in the Appendix A. The natural scales of the different functions in the factorization at threshold in Eq. (4.29) are given by The solution of the evolution equations can be obtained in analogy to the calculation outlined for e + e − hemisphere jets above. Similar to leading e + e − hemisphere jets, Eq. (4.2), we can relate the first moment (second Mellin moment) of the leading subjet cross section to the average partonic energy loss Here the result is weighted by appropriate quark/gluon fractions f q,g and σ tot is the inclusive jet cross section (without the substructure measurement).

Nonperturbative effects
Finally, we must deal with possible nonperturbative effects in the shower. To regulate the Landau pole, we shift the argument of the running coupling as: We take m reg ∼ 500 MeV. To model the nonperturbative contribution to the cross section, we focus on the case of leading jets in e + e − . Then we replace the soft function in Eq. (4.5) as The function S c i (z, QR, µ) pert explicitly refers to the perturbative resummed result of Eq. (4.25). Here S is a shape function that will shift the spectrum near z = 1, similar to what happens in the case of event-shapes [121][122][123]. The 1/R scaling of the nonperturbative correction was first identified in Ref. [124].

Numerical results
In this section we present numerical results for e + e − hemisphere leading jets and leading subjets in proton-proton collisions. We implement the threshold resummed hard and jet functions in the Monte Carlo parton shower discussed above which allows us to calculate both the inclusive and leading jet cross section at NLL accuracy. In Fig. 6, we show the results for e + e − hemisphere jets for quarks and gluons separately. 4 As an example, we choose the jet radius of R = 0.5 and the hard scale Q = √ s = 91.2 GeV. The inclusive and leading jet spectra agree for z > 1/2. For e + e − hemisphere jets, a jet with momentum fraction z > 1/2 is automatically the leading jet. Note that this does apply to event-wide leading jets in e + e − collisions as discussed in section 7.1 below. We observe that both spectra peak at large values of z which indicates that it is very likely to find a jet that carries a large momentum fraction of the initial quark or gluon. See also Ref. [125]. The peak is less pronounced for an initial gluon than for quarks which is expected due to the different color factors. The peak structure at large values of z confirms that the identified leading jet is a good proxy of the underlying parton level degrees of freedom. We note that the peak arises due to the threshold resummation. At LO/LL accuracy the numerical result diverges near the endpoint, see Fig. 3. Therefore, it is phenomenologically important to include threshold resummation for leading jet measurements. Note that the suppression of the cross section for z → 1 is unusual since threshold resummation is typically associated with an enhancement of the cross section [30,31]. For z < 1/2 the inclusive and leading jet spectrum differ due to the subleading jets which contribute only to the inclusive jet spectrum. The leading jet cross section drops significantly below z = 1/2 indicating that it is very unlikely to find a leading jet that carries only a small momentum fraction z.
Another intriguing feature of the results in Fig. 6 is the shape of the leading jet spectrum around z = 1/2. At LO/LL accuracy there is a sharp kink or cusp which is now smeared out. In addition, we observe that the leading jet spectrum now extends to relatively small values of z which is numerically quite different compared to the result at LO/LL accuracy. We would like to emphasize that the resummation here is critical to obtain the full leading jet spectrum. For all values of R, the tail of the distribution extends to small values of z. Instead, any fixed order computation at N n LO would not give finite values below z = 1/(n + 1). Next, we study the dependence of the e + e − hemisphere leading jet cross section on the jet radius R. The results for R values in the range of 0.05 − 0.8 are shown in Fig. 7. For small values of R, the leading jet can only capture a small fraction of the initial energy. Indeed, we observe that the curves shift toward smaller values of z as we decrease the jet radius R. This shift differs significantly for quark and gluon jets. In the next section we consider the average z value of these distributions which is related to the jet energy loss.
In Fig. 8, we show numerical results for inclusive and leading subjets in proton-proton collisions at √ s = 13 TeV. We choose exemplary jet kinematics as indicated in the figure and include appropriate quark/gluon fractions. The leading subjet cross section is the most straightforward possibility to directly quantify jet energy loss in proton-proton and heavy-ion collisions. We note that the size of the scale uncertainty band could be reduced by including higher order corrections.

Quantifying jet energy loss
As discussed in section 2, the leading jet functions constitute probability densities that allow us to quantify the energy loss of leading jets. In this section we discuss different statistical quantities both at the level of the (threshold resummed) jet functions and and cross sections. In particular, we focus on e + e − hemisphere leading jets. However, due  to universality of the leading jet functions, the same arguments apply to subjets or other suitable observables discussed in section 7 below. In section 5.1 we start with the mean z i1 and variance σ i of the leading jet energy distribution. The mean is directly related to the average energy loss which we define as all the energy that is not contained in the leading jet z i,loss = 1 − z i1 . In section 5.2 we focus on differences between quark and gluon leading jets and in section 5.3, we study the Shannon entropy and the KL divergence to further quantify differences between quarks and gluons.

Mean and variance
We start by studying the mean and variance which are two fundamental quantities that quantify parton/jet energy loss. The mean or average energy of the initial parton which is contained inside the leading-jet z i1 is given by the first moment (second Mellin moment) of the leading jet function J i as introduced in section 2.2. We repeat the relevant equation here for convenience. For quarks and gluons, we find The average energy fraction of the leading jet depends on the scale Q and the jet radius R which we omit on the right hand side. In Ref. [18], an expansion of the average out-of-jet radiation in α s ln R was performed. Here we perform the complete expansion in powers of the strong coupling constant α s which requires knowledge of the entire jet function. The mean or average energy loss of the leading jet is given by z i,loss = 1 − z i1 . This relation holds to all orders. At NLO, z i,loss coincides with the average energy fraction contained in the first subleading jet. At higher orders, the average lost energy z i,loss is shared amongst the different subleading jets. We consider both cone [126,127] and k T -type [128][129][130][131] jets.
We start with an NLO computation of the average momentum fraction which is contained in the leading jet at the jet function level. From Eq. (5.1), we find for quarks and gluons We note that the terms ∼ ln(1/R 2 ) are the same for both k T and cone-type jets and they agree with the result in [18]. The remaining finite O(α s ) corrections are reported here for the first time. As expected those terms depend on the jet algorithm. An important aspect of the full O(α s ) result is that it leads to a finite energy loss even for R → 1, where the logarithms of the jet radius vanish. At small jet radii the logarithmically enhanced terms in Eq. (5.2) dominate and need to be resummed to all orders. In principle analytical results at higher orders in α s could be obtained from the non-linear evolution equations in Eq. (2.15). Instead, here we are going present numerical results using the parton shower framework which was introduced in the previous sections which also includes the resummation of threshold logarithms. We present numerical results not at the jet function level but for a cross section which can be measured in experiments, i.e. including also the (threshold resummed) hard function. We consider e + e − hemisphere leading jets as an example. The average energy loss for quark and gluon jets is shown in Fig. 9 using the parton shower algorithm described above including threshold resummation at NLL . We plot the numerical result for z i,loss as a function of the jet radius on a logarithmic scale. We choose two exemplary hard scales of Q = 91.2 GeV and 500 GeV which set the initial parton energy. The rightmost R values correspond to QR = 1 GeV, where the average energy loss for both Q values turns out to be very similar. We observe that the energy loss of gluons is larger than for quarks which is generally expected due to the different color factors and DGLAP splitting functions. In addition, the average energy loss of leading jets is larger at smaller scales Q. For large jet radii, the average energy loss of a quark is around 10% and 15-20% for an initial gluon. This observation quantifies why jets are generally considered to be excellent proxies of parton level dynamics compared to hadrons. For the small jet radii, the average energy loss increases up to around 70% for quarks and 80% for gluons. As expected, a relatively small jet radius can only capture a small fraction of the initial momentum of the quark or gluon ln (1/R) Q = 500 GeV NLL 1/2 < ζ i < 2 Figure 9. Average energy loss of e + e − hemisphere leading jets z k T i,loss for k T -type jets. We separately show the quark and gluon results for an initial energy of Q = 91.2 GeV and 500 GeV as a function of the jet radius. and the energy loss grows significantly. We note that hadrons correspond to jets with vanishing radius where an additional nonperturbative parton-to-hadron transition needs to be taken into account which makes the average energy loss of hadrons even larger than what is shown in Fig. 9. We would like to stress again that such a statement is not possible for inclusive jets since any emission outside the leading jet constitutes another jet which is also taken into account when an inclusive jet sample is measured as discussed in section 2 above. The result for the parton/jet energy loss presented here depends significantly on the perturbative higher order corrections which we include in the parton shower 3. In this sense, the results presented here constitute the first quantitative calculation of vacuum parton/jet energy loss which allows for a meaningful connection to experimental measurements.
Besides the mean of the leading jet probability distribution, we can also calculate the variance or the fluctuations of jet energy loss. The variance for quarks and gluons is defined as Similar to the mean, the variance σ i depends on the jet algorithm, the jet radius R and initial scale Q. We omit the explicit dependence of σ i on those quantities here for notational convenience. However, we study the dependence of σ i on those variables numerically below.
Here z 2 i is the second moment (third Mellin moment) of the leading jet function. At the jet function level it is given by We start again with an NLO calculation of the variance at the jet function level. For quarks and gluons, we find Similar to the mean in Eq. (5.2), the ∼ ln(1/R 2 ) term is independent of the jet algorithm. Next, we study the variance using the full parton shower for e + e − hemisphere leading jets.
The results for quarks and gluons are shown in Fig. 10 as a function of the jet radius. We observe that the variance for both quark and gluon jets is in the range of ∼ 0.1 − 0.2 for Q = 91.2 GeV and 500 GeV. For large R, the variance is almost independent of the scale Q but the scale dependence becomes more visible toward smaller R. For gluons the variance peaks at intermediate values of the jet radius R. The leading jet distribution for gluons peaks either at large-z (large R) or large-z (small R) which leads to relatively small values of σ k T g . Only in the intermediate R region, the gluon z-distribution is broad which leads to the maximum of σ k T g that we observe in Fig. 10. For quarks the z-distribution is even more peaked at large R than for gluons which leads to a smaller variance. It evolves more slowly toward small-z than the gluon distribution which is why the variance continues to increase toward small R and eventually becomes even larger than for gluons. At small R the variance for quarks has a maximum and (close to the nonperturbative region) becomes smaller again. It will be interesting to study how the mean and variance of leading jets/the energy loss is modified in heavy-ion or electron-nucleus collisions where the notion of (medium induced) parton/jet energy loss plays an important role.
We end this section by noting that there is no unique definition of energy loss. For example, we could adopt the definition that all the energy which is not contained in the first two leading jets is "lost energy". We can also calculate the average energy loss z i,loss for this alternative definition from the subleading jet function as ln (1/R) Q = 500 GeV NLL 1/2 < ζ i < 2 Figure 10. Variance σ k T i of the energy loss of e + e − hemisphere leading jets for quark and gluons. As in Fig. 9, we choose a k T -type jet algorithm, two initial reference scale of Q = 91.2 GeV and 500 GeV and we plot the result as a function of the jet radius.
Also other definitions are possible as long as we consider a fixed number of jets. For inclusive jets it is not possible to construct a probability density since the number of jets is not fixed but it is generated dynamically event-by-event.

Quark/gluon energy loss
The dependence of the energy loss mechanism on the parton flavor has been discussed extensively in the literature -in particular in the context of jet quenching studies in heavy-ion collisions. While we only focus on vacuum energy loss in this work, our results set the baseline for studies of energy loss in the nuclear medium. In this section we compare the LO estimate of the quark/gluon energy loss to the full result from the parton shower. First, we consider the LO emission spectrum in the soft gluon approximation. For an initial quark/gluon, we find dI q,g dz ∼ C F,A 1 − z . (5.14) This relation implies that in this limit the ratio of the average jet energy loss of a quark and gluon is z q,loss z g,loss which is sometimes referred to as Casimir scaling in the literature. The ratio of the quark/gluon energy loss z q,loss / z g,loss is shown in Fig. 11 as a function of the jet radius R for Q = 91.2 GeV and 500 GeV. Here we choose again e + e − leading jet production as a representative example. We observe that the ratio of the average energy loss of quarks and gluons is almost identical for the two different Q values at large R but differs at intermediate and smaller R. Interestingly, the ratio of the average energy loss of quarks and gluons is significantly closer to 1 compared to the LO estimate in We explore this feature in more detail in the context of quark/gluon tagging below. It turns out that the value of the ratio of the quark/gluon average energy loss agrees for QR = 1 GeV, which corresponds to the rightmost points in Fig. 11 and is given by ≈ 0.9.

Shannon entropy, KL divergence
Since the leading jet cross section constitutes a probability distribution, we can compute various other statistical quantities besides the mean and variance. In this section we consider the Shannon entropy and the KL divergence. In order to quantify the average uncertainty of the leading jet/energy loss probability distribution, we consider the Shannon/information entropy. We can write the continuous version of the Shannon entropy h i at the jet function level as In addition, we consider the KL divergence D KL to quantify the difference between quarks and gluons . The KL divergence is not symmetric under J q ↔ J g . Nevertheless, it is a very useful measure to quantify the similarity of two probability distributions. An alternative measure would be the Jensen-Shannon divergence which is symmetric. Note that here we introduced both quantities in Eqs. (5.16) and (5.17) for a continues variable z. However, experimental measurements are binned which is why we replace the integral versions of these metrics in Eqs. (5.16) and (5.17) by their corresponding discrete versions. For our numerical results presented below we choose a binning of N = 1000 steps in z. The measurement of both of these quantities will also shed new light on the energy loss mechanism in the hot or cold nuclear matter environment. Here we only consider the Shannon entropy and the KL divergence for the probability densities of the leading jet but they can be extended to subleading jets as well.
The Shannon entropy at the cross section level for e + e − hemisphere leading jets is shown in Fig. 12. It peaks at intermediate values of the leading jet radius. At large values of R, the uncertainty of gluon jets is larger compared to quark jets and a larger value of Q leads to a lower value. At intermediate to small values of R the ordering changes. The KL divergence is shown in Fig. 13. We observe that the KL divergences peaks at intermediate values of R which indicates that there is an optimal value of R to distinguish quark and gluon leading jets. The maximum value depends on the scale Q and is shifted toward smaller R for larger values of Q. Interestingly, this maximum is in the perturbative regime indicating that leading (sub)jets may be a good observable for quark/gluon discrimination which is under perturbative control. We further explore the potential of leading (sub)jets as a quark/gluon discriminant in the next section. We also notice that the rightmost points of the two curves in Fig. 13 agree, which correspond to QR = 1 GeV in both cases. This is likely due to the relative simplicity of the non-perturbative model that we are using, without tuning the it to account for differences between quark or gluon initiated jets. 6 Quark/gluon discrimination with leading (sub)jets A typical task of jet substructure observables is the discrimination between quark and gluon jets. In this section we study the tagging performance of leading (sub)jets. It is instructive to compare leading (sub)jets to established quark/gluon jet tagging techniques in the literature [132]. Using the largest momentum fragment in a jet as a classifier is in the same class of observables known as "fractal jet observables" [83], see also Refs. [84,133]. These observables are tuned to the energy flow patterns generated in DGLAP evolution of the final state. Another set of useful classifiers is given by the so-called generalized angularities [134] (see also [110,[135][136][137]) Here the sum runs over all particles inside the jet with radius parameter R, and z i , θ i are their momentum fractions and angles with respect to the jet axis, respectively. For example, the average momentum fraction of the leading jet can be obtained from Eq. (6.1) as follows. First, instead of summing over individual particles in the jet, we sum over all subjets which are obtained by reclustering the initial jet with a jet radius r < R, see also section 4.2 above. Second, we choose β = 0 which is similar to jet multiplicities and the jet p D T observable [138,139]. Third, take the expression to the power of 1/κ and we take the limit κ → ∞ which singles out the leading jet momentum fraction.
Typically, the performance of a classifier is quantified by studying the ROC curve. The ROC curve shows the quark/gluon true positives (cumulative distribution function (CDF)) vs. the false positives rates for a given decision threshold. The result for (sub)leading jets and two values of the energy scale Q and various values of the jet radius R are shown in Fig. 14. We observe that the discrimination power changes significantly as a function of the jet radius.
In order to obtain a single value to quantify the performance of leading (sub)jets as a quark/gluon discriminant, the area under the ROC curve (AUC) is commonly used which is shown in the right panel of Fig. 15. We note that the location of the peak differs slightly between the AUC and the KL divergence in Fig. 13 above. As mentioned above, in both cases, the peak of the distributions is in the perturbative region. Instead, other observables that perform very well are often nonperturbative such as the particle multiplicity. Since leading (sub)jets correspond to the most energetic part of a jet, it may capture relevant information with little overlap with e.g. the particle multiplicity. We plan to explore the information content of leading (sub)jets relative to other observables in future work. See for example Ref. [140].

Further applications
In this section we discuss several cross sections involving leading jets and hadrons. We start by considering event-wide leading jets in e + e − collisions in section 7.1 in contrast to hemisphere leading jets. We present numerical results and compare to Pythia 8 [58] simulations. Similar measurements were performed at LEP. We then discuss leading jets in SIDIS (section 7.2), photon-jet correlations (section 7.3) and eventually extend our parton shower framework toward leading hadrons instead of jets (section 7.4).

e + e − event-wide leading jets
Instead of the e + e − hemisphere leading jets, which we discussed in the previous sections, we are now going to consider the leading jet in the entire event. The necessary factorization structure here is similar to leading jets in proton-proton collisions which was discussed in section 2.3, except that here we have direct access to the initial reference scale Q = √ s. We consider the cross section differential in the energy fraction of the leading jet z 1 = 2E 1 /Q. At LL accuracy, we have qq (Q, µ) dz q dzq J q (z q , QR/2, µ) Jq(zq, QR/2, µ) × δ(z 1 − max{z q , zq}) . (7.1) The quark and anti-quark each fragment into a leading jet with energy fraction which we denote by z q,q . The event-wide leading jet momentum fraction is then obtained by picking the larger value of z q,q which is taken into account by the delta function Eq. (7.1). The topology of event-wide leading jets in e + e − collisions is illustrated in the left panel of Fig. 16. As discussed in section 2.3, at higher perturbative accuracy we need to take into account additional jet functions as well as integrals over the hard function. We stress again that the factorization structure here is very different compared to inclusive jets. In this case any jet is taken into account independent of their energy fraction which allows for a factorization structure in terms of a convolution integral and which is independent of the perturbative order. We include again the threshold resummation for e + e − event-wide leading jets. Working at NLL accuracy, the threshold resummed jet and hard functions can be obtained analogous to the results in section 4. We calculate the cross sections in two different ways. First, we evolve the threshold resummed hard and jet functions with the parton shower, i.e. for one initial parton, and we then compute the integral in Eq. (7.1). Second, we initialize the parton shower with two (back-to-back) partons and determine the event-wide leading jet directly from the output of the parton shower. We find that both methods give the same result as expected. We note that we cannot directly connect event-wide leading jets in e + e − collisions to parton energy loss due to the structure of the factorization in Eq. (7.1). The main difference compared to e + e − hemisphere leading jets is that now we have two instead of one parton at LO, see criterion 3 in the Introduction. We show our numerical results for e + e − event-wide leading jets in Fig. 17 for two exemplary values of the jet radius. We observe that the spectrum peaks close to z ≈ 1 and falls off steeply toward smaller values of z. We note that the spectrum looks significantly different compared to e + e − hemisphere jets and is more peaked at large values of z since there are now at least two jets produced in each hemisphere and we pick the more energetic one. For comparison, we also show the inclusive jet spectrum in Fig. 17. The inclusive spectrum starts to deviate from the leading jet result around z ≈ 0.94 which differs from the hemisphere jet case where the two spectra only start to differ for z < 0.5. In addition, we show Pythia 8 [58] results for the anti-k T [131] and C/A [129,130] algorithm. Pythia is tuned very well to LEP data and can serve as a benchmark for e + e − jet observables. We find good agreement with our numerical results where the nonperturbative parameter is chosen as Λ QCD = 0.2 GeV. The two jet algorithms used for the Pythia simulations give identical analytical results in the threshold region at NLL . Indeed, we observe that the Pythia results for the two algorithms are very similar.
Interestingly, the OPAL Collaboration at LEP provided data for the leading and the first two subleading jet cross sections in Ref. [141]. OPAL used an algorithm that finds, by construction, three jets in every recorded event. This procedure does not correspond to our definition of jets making a one-to-one comparison impossible. We show the OPAL data in Fig. 18. The OPAL measurement imposed an intra-jet angle which roughly corresponds to the jet radius chosen in the right panel of Fig. 17. While the shape of the spectrum is quantitatively different compared to our results in Fig. 17, we do find qualitative agreement. We note that for event-wide jets, the second leading jet is not required to carry an energy fraction of z 2 < 0.5 (as for e + e − hemisphere jets) since we have at least two jets in the event. In fact, from Fig. 18 we find that the first subleading jet energy extends up to z 2 ∼ 0.94. The OPAL results demonstrate that it is possible to experimentally measure leading and subleading jet spectra at LEP. We expect that experiments like the LHC and EIC can measure leading jet observables equally well.

Semi-Inclusive Deep-Inelastic Scattering
SIDIS measurement allow for a clean measurement of jet/parton energy loss since we only have one quark which fragments into the observed leading jet at LL accuracy. We consider the process e(k)+p(P ) → e (k )+jet 1 (P jet 1 )+X where both the final state electron and the leading jet are observed. The reference scale with respect to which the leading jet energy is measured is set by the virtuality Q 2 = −q 2 = −(k − k ) 2 of the exchanged photon. We consider the cross section dσ ep→e +jet 1 +X dx B dy dz 1 .
where we introduced the usual variables for SIDIS with hadrons in the final state: Bjorken x B and the inelasticity y which are given by We have Q 2 = x B ys, where √ s is the electron-proton CM energy. The momentum fraction of the leading jet z 1 is defined as Here the last equality holds in the Breit frame and P + jet 1 denotes the large lightcone momentum component of the leading jet. In the target rest frame we have z 1 = E jet 1 /Q which OPAL √ s = 91.2 GeV Figure 18. OPAL results for the energy spectrum of the leading jet (Jet 1) and the first two subleading jets (Jet 2,3) using 3-jet events. Data taken from [141]. makes this asymmetric process very similar to e + e − hemisphere leading jets. Suitable jet clustering algorithms for this type of observable in the Breit frame were discussed in Ref. [125]. Within collinear factorization, the cross section for inclusive jets can be written in terms of F T,L , the transverse and longitudinal structure functions [103,142] dσ ep→e +jet+X dx dy dz = 4πα 2 Q 2 1 + (1 − y) 2 2y F T (x, z, Q) + 1 − y y F L (x, z, Q) , (7.5) where α is the fine structure constant and we dropped the subscript B of Bjorken x B for simplicity. In the current fragmentation region, the structure functions a = T, L can be written as Here, ⊗ denote convolution integrals in x and z, f i are the PDFs and J j are the inclusive jet functions. Therefore, also for SIDIS at LL accuracy we can consider the first moment as a direct measure of the average quark energy loss Here σ tot is the inclusive DIS cross section dσ/dx/dy. The threshold resummation for jets in SIDIS follows from Refs. [125,[143][144][145][146] which can also be implemented in the parton shower algorithm discussed in section 4. We leave more detailed numerical studies especially for kinematics at the EIC for future work.

Photon-jet correlations
In this section we consider the production of a direct photon and measure the momentum of the leading jet in the opposite hemisphere. See Refs. [53,54,147,148] for related recent experimental results. We consider the cross section differential in the photon's transverse momentum and rapidity p γ T , η γ . In addition, we measure the recoiling leading jet's transverse momentum p T relative to the photon x J1γ = p T /p γ T . At LO/LL accuracy the factorization of the cross section can be written as [149] dσ Here we consider the transverse momentum of the direct photon p γ T as the reference scale with respect to which we measure the energy loss of the recoiling leading jet. This observable can also give direct access to the weighted average quark/gluon energy loss However, different than the processes considered above, x Jγ can generally be > 1 making a clear interpretation in terms of energy loss more difficult. We may nevertheless get a handle on the energy loss of the leading jet recoiling against the photon by considering the difference between inclusive and leading jets. At the jet function level, we find that we can rewrite the average momentum fraction of the leading jet as z = We can therefore also get access to the average jet energy loss by measuring the difference between the inclusive and leading jet spectra recoiling the direct photon. The inclusive and leading spectra from Pythia 8 [58] are shown in Fig. 19. In practice, jets can only be reconstructed down to a certain transverse momentum. Therefore, we need to introduce a small cutoff z cut in Eq. (7.10), which gives z (z cut ) = 1 + 1/2 zcut dz z J i (z, p γ T R, µ) − J i (z, p γ T R, µ) . (7.11) We leave more detailed numerical studies using our parton shower framework for future work.

Toward leading hadrons
In order to extend our calculation from leading jets to leading hadrons, we need to evolve the shower down to the nonperturbative scale µ NP = O(1GeV) and include hard-scattering functions and fragmentation functions in the shower algorithm. The leading and subleading hadron fragmentation functions are currently not known and we are not aware of existing data sets that could constrain them. Similar to jet functions, the leading and inclusive parton-to-hadron fragmentation functions agree for z > 1/2 but differ for smaller values of z. We expect that the necessary experimental measurements are feasible which can provide important new insights into the QCD fragmentation mechanism. We leave more detailed phenomenological studies for future work. Here we only focus on the parton shower evolved to the nonperturbative scale without including fragmentation functions and hard-scattering functions. Analogous to leading jets discussed above, we can use the Monte Carlo setup introduced in section 3 to calculate the spectrum of leading partons since the corresponding leading hadron fragmentation functions satisfy the same non-linear evolution equations. Now we evolve the shower down to a scale of order µ NP = O(1 GeV) which we consider as the onset of nonperturbative physics. To be specific, the Monte Carlo algorithms terminates when t > t max = t(Q, 1/Q) or equivalently QR = O(1 GeV) in Eq. (3.3).
In the upper two panels of Fig. 20, we show the result for the z-spectrum for leading D i and inclusive D i partons evolved down to µ NP starting from Q = 91.2 GeV. For the numerical results shown here we choose R = 0.01. Even without sampling from the nonperturbative fragmentation function, the gluon distribution (left) already has the typical shape of a hadron fragmentation spectrum which falls off steeply toward z → 1. Instead, the quark spectrum (right) still has a peak at very large values of z similar to the evolved LL jet functions discussed above. See Fig. 3. The leading and inclusive parton spectra start to deviate around z ≈ 0.4. We note that all the leading jet cross sections discussed in previous sections can also be calculated for hadrons. In particular, we can also define a nonperturbative version of energy loss in terms of leading hadrons instead of jets. We leave more detailed studies of the leading hadron energy loss for future work.

Conclusions
In this work we discussed leading and subleading jet cross sections which probe fundamental aspects of the QCD fragmentation process. Different than inclusive jets, the formation and evolution of leading jets is described by jet functions with non-linear DGLAP-type evolution equations. These leading and subleading jet functions constitute normalized probability densities to find a (sub)leading jet with a given longitudinal momentum fraction from an initial quark or gluon. Instead, the jet functions relevant for inclusive jet production are number densities where the total number of jets per event is not fixed and produced dynamically event-by-event. Motivated by results in probability theory [23], we established relations between leading/subleading jets and inclusive single-, di-and tri-jet functions which we plan to explore in more detail in future work.
We focused in particular on cross sections where we have access to an additional reference scale Q with respect to which we can define the longitudinal momentum fraction z 1 which is contained in the leading jet. We are then able to define the (average) out-of-jet radiation or the energy loss of leading jets as z loss = 1 − z 1 which can be computed orderby-order in perturbation theory and which can be accessed directly through experimental measurements. For observables where we only have one parton at leading order, the experimentally accessible leading jet energy loss can be identified with parton energy loss at leading-logarithmic accuracy. We identified criteria of suitable observables which allow for a direct measurement of the leading jet energy loss which include e + e − hemisphere jets, subjets in proton-proton collisions, jets in Semi-Inclusive Deep Inelastic Scattering and photon-jet correlations. From these cross sections we can obtain the average energy loss of leading jets z loss which can be compared to our theoretical results.
One of the main new developments of our work is a parton shower framework that allows us to compute threshold resummed leading jet cross sections at next-to-leading logarithmic accuracy (NLL ). Hard-scattering functions, jet functions and also fragmentation functions can be included directly in the parton shower framework. The results of the parton shower agree exactly with analytical results for inclusive jets and allow for a well-defined extension to leading jet cross sections. The threshold resummation which we include for both the hard and jet functions is phenomenologically important for leading jet observables. We derived the threshold resummation for leading jets in e + e − collisions and leading subjets in proton-proton collisions. While the developed framework is a "few purpose" parton shower, we expect that it can be extended systematically to other observables.
We presented numerical results for e + e − hemisphere and event-wide leading jets as well as leading subjets in proton-proton collisions at NLL accuracy using the parton shower framework. For e + e − event-wide leading jets we compared to Pythia 8 results and found good agreement. Interestingly, the OPAL Collaboration at LEP measured similar leading jet cross sections in Ref. [141]. While a one-to-one comparison to the existing data is not possible since their definition of leading jets differs from ours, these measurements demonstrate that leading jet and hadron measurements are generally feasible.
We investigated the differences of the average energy loss between leading quark and gluons jets. We observed that the differences are surprisingly small compared to a leadingorder estimate and that they are rather independent of the jet radius R. Besides the average energy loss, the mean of the energy loss probability distribution, we also considered for the first time the variance which quantifies event-by-event fluctuations of energy loss. In addition, we computed the Shannon entropy and the KL divergence. The latter quantifies the difference between quark and gluon jet energy loss. We further explored the potential of leading (sub)jets to discriminate between quark and gluon jets. We presented ROC curves and the AUC for different values of the jet radius R. Interestingly, the best discrimination power is achieved for a perturbative value of R. In the future, we plan to explore the tagging performance of leading (sub)jets and study their complementarity to other typical observables such as particle multiplicity.
In addition, we outlined how our work can be extended to leading hadrons. Similar to leading jets, leading hadron fragmentation functions are normalized probability densities which allow us to establish a well-defined but nonperturbative notion of energy loss. Our results constitute the first quantitative calculation of leading jet energy loss which can compared directly to experimental data. We expect that our results will be particularly useful to study leading jets which traverse hot or cold nuclear matter in heavy-ion collisions or electron-nucleus collisions at the future EIC. Medium induced emissions can generally increase the energy loss of leading jets and the corresponding energy loss spectrum introduced here can provide important information about the Quark-Gluon Plasma/cold nuclear matter and its interaction with hard probes.
where P ji (N ) denote the Altarelli-Parisi splitting functions in Mellin space. In addition, we have γ H ji = α s /πP ji (N ) for i = j.