Energy correlation functions for jet substructure

We show how generalized energy correlation functions can be used as a powerful probe of jet substructure. These correlation functions are based on the energies and pair-wise angles of particles within a jet, with (N + 1)-point correlators sensitive to N -prong substructure. Unlike many previous jet substructure methods, these correlation functions do not require the explicit identification of subjet regions. In addition, the correlation functions are better probes of certain soft and collinear features that are masked by other methods. We present three Monte Carlo case studies to illustrate the utility of these observables: 2-point correlators for quark/gluon discrimination, 3-point correlators for boosted W/Z/Higgs boson identification, and 4-point correlators for boosted top quark identification. For quark/gluon discrimination, the 2-point correlator is particularly powerful, as can be understood via a next-to-leading logarithmic calculation. For boosted 2-prong resonances the benefit depends on the mass of the resonance.

In this paper, we introduce generalized energy correlation functions that can identify N -prong jet substructure without requiring a subjet finding procedure. These correlators only use information about the energies and pair-wise angles of particles within a jet, but yield discrimination power comparable to methods based on subjets. As we will see, (N + 1)point correlation functions are sensitive to N -prong substructure, with an angular exponent β that can be adjusted to optimize the discrimination power. To our knowledge, the 2-point correlators-schematically i,j E i E j θ β ij where the sum runs over all particles i and j in a jet or event-first appeared in Ref. [54] and independently in Ref. [55], with no previous studies of higher-point correlators. 1 Besides the novelty of not requiring subjet finding, a key feature of the generalized energy correlation functions is that the angular exponent β can be set to any value consistent with infrared and collinear safety, namely β > 0. In contrast, observables like angularities [9,51] are required to have β > 1 to avoid being dominated by recoil effects. 2 By choosing values of β 0.2, the correlators are able to more effectively probe small-scale collinear splittings, which will turn out to be useful for quark/gluon discrimination.
To put our work in perspective, it is worth remembering that the basic idea for using energy correlation functions to determine the number of jets in an event is actually quite old. As we will review, the C-parameter for e + e − collisions [61,62] is essentially a 3-point energy correlation function that can be used to identify events that have two jets. However, the C-parameter is defined as a function of the eigenvalues of the sphericity tensor and therefore 1 Our definition of the energy correlation function should not be confused with Refs. [56][57][58][59], which refer to an ensemble average of products of energies measured at fixed angles. Here, energy correlation functions are measured on an event-by-event basis. 2 As will be discussed in Sec. 2.2 and a forthcoming publication [60], N -subjettiness may or may not have recoil sensitivity depending on how the axes are chosen.
only gives sensible values for systems that have zero total momentum and for events that are nearly dijet-like. In contrast, our generalized energy correlation functions give sensible results in any Lorentz frame and can be used to identify any number of jets in an event (or subjets within a jet). In addition, they can be defined in any number of spacetime dimensions. The remainder of the paper is organized as follows. In Sec. 2, we introduce arbitrarypoint energy correlation functions and define appropriate energy correlation double ratios C (β) N (built from the (N + 1)-point correlator), which can be used to identify a system with N (sub)jets. We also contrast the behavior of C (β) N with N -subjettiness ratios. We then present three case studies to show how these generalized energy correlation functions work for different types of jet discrimination.
• Quark/gluon discrimination. Using C (β) 1 (built from the 2-point correlator) in Sec. 3, we perform both an analytic study and a Monte Carlo study of quark/gluon separation. Through a next-to-leading logarithmic study, we explain why quark/gluon discrimination greatly improves as the angular exponent approaches zero (at least down to β 0.2), highlighting the importance of working with recoil-free observables.
• Boosted W /Z/Higgs identification. Using C (β) 2 (built from the 3-point correlator) in Sec. 4, we will see that the discrimination power between QCD jets and jets with two intrinsic subjets from a colour-singlet decay depends strongly on the ratio of the jet mass to its transverse momentum. This occurs because a QCD jet obtains mass in different ways depending on this ratio. In particular, we will see that the energy correlation function performs better than N -subjettiness in situations where the jet mass is dominated by soft wide-angle emissions.
• Boosted top quark identification. Using C (β) 3 (built from the 4-point correlator) in Sec. 5, we find comparable discrimination power to other top-tagging methods. While one might worry that the 4-point correlators would face a high computational cost, we find that a boosted top event can be analyzed for a single value of β in a few milliseconds.
We conclude in Sec. 6 with an experimental and theoretical outlook. The energy correlation functions are available as an add-on to FastJet 3 [63] as part of the FastJet contrib project (http://fastjet.hepforge.org/contrib/).

Generalized Energy Correlation Functions
The basis for our analysis is the N -point energy correlation function (ECF) (2.1) Here, the sum runs over all particles within the system J (either a jet or the whole event). Each term consists of N energies multiplied together with N 2 pairwise angles raised to the angular exponent β. This function is well-defined in any number of space-time dimensions as well as for systems that do not have zero total momentum. Note that it is infrared and collinear (IRC) safe for all β > 0. Moreover, ECF(N, β) goes to zero in all possible soft and collinear limits of N partons. As written, Eq. (2.1) is most appropriate for e + e − colliders where energies and angles are the usual experimental observables. For hadron colliders, it is more natural to define ECF(N, β) as a transverse momentum correlation function: 3 where R ij is the Euclidean distance between i and j in the rapidity-azimuth angle plane, In this paper, we will only consider up to 4-point correlation functions: If a jet has fewer than N constituents then ECF(N, β) = 0. Note that the computational cost for ECF(N, β) with k particles scales like k N /N !. From the ECF(N, β), we would like to define a dimensionless observable that can be used to determine if a system has N subjets. The key observation is that the (N + 1)point correlators go to zero if there are only N particles. More generally, if a system has N subjets, then ECF(N + 1, β) should be significantly smaller than ECF(N, β). One potentially interesting ratio is which behaves much like N -subjettiness τ N in that for a system of N partons plus soft radiation, the observable is linear in the energy of the soft radiation. 4 Of course, this is but one choice for an interesting combination of the energy correlation functions, and one can imagine using the whole set of energy correlation functions in a multivariate analysis. 3 We will continue to use the notation ECF, though we will mainly use the transverse momentum version in this paper. 4 Unlike N -subjettiness, this ratio scales like γ 1−N β under transverse Lorentz boosts γ, which is somewhat undesirable when considering systems with several subjets.
In this paper, we will work exclusively with the energy correlation double ratio which is dimensionless. 5 One way to motivate this observable is that we already know that N -subjettiness ratios τ N /τ N −1 are good probes of N -prong substructure [49,50]. As we will see, the notation "C" is motivated by the fact that this variable generalizes the C-parameter [61,62]. One should keep in mind that C (β) N involves (N + 1)-point correlators, and when clear from context, we will drop the (β) superscript.
The energy correlation double ratio C N effectively measures higher-order radiation from leading order (LO) substructure. For a system with N subjets, the LO substructure consists of N hard prongs, so if C N is small, then the higher-order radiation must be soft or collinear with respect to the LO structure. If C N is large, then the higher-order radiation is not strongly-ordered with respect to the LO structure, so the system has more than N subjets. Thus, if C N is small and C N −1 is large, then we can say that a system has N subjets. In this way, the energy correlation double ratio C N behaves like N -subjettiness ratios τ N /τ N −1 , with key advantages to be discussed in Sec. 2.2.

Relationship to Previous Observables
While the definition of the energy correlation double ratio C N is new, it is related to previous observables for e + e − and hadron colliders that have been studied in great detail.
An energy-energy correlation (EEC) function for e + e − events was introduced in Ref. [54] for its particularly nice factorization and resummation properties. It is defined as where the sum runs over all particles in the event and n T is the direction of the thrust axis. This variable is IRC safe for all a < 2. The Θ-function is only non-zero if the pair of particles is in the same hemisphere. This removes the large correlation of the two initial hard partons which would otherwise dominate the sum, and means that EEC a behaves much like the jet angularities [9,51] with the same angular exponent a. The EEC was introduced because it is insensitive to recoil effects and has smooth behavior for all allowed values of a. In particular, EEC a has a smooth transition through a = 1, whereas angularities exhibit non-smooth behavior and also are increasingly sensitive to recoil effects as the angular power a increases. If one considers only one hemisphere of a dijet event, then EEC a is approximately the same as C (β) 1 in our notation with β = 2 − a. Both observables are sensitive to 1-prong (sub)structure, and we will discuss the issue of recoil further in Sec. 2.2. A related two-particle angular correlation function was introduced in Refs. [21,55,64] for discrimination of jets initiated by QCD from jets from boosted heavy particle decays. The angular correlation function is defined as where the Θ-function only allows pairs of particles separated by an angular scale of R or less to contribute to the observable. The behavior of the observable can be studied as a function of R, and jets that are approximately scale invariant should have an angular correlation function that scales as a power of R. For a fixed value of R, the properties of the angular correlation function are very similar to that of EEC a and C (β) 1 . As mentioned above, the notation C (β) N was chosen because of its relation to the Cparameter from e + e − collisions [61,62]. The C-parameter is used to identify two-jet configurations without recourse to a jet algorithm or explicit jet axes choice. It is defined as which can also be expressed in terms of the eigenvalues of the sphericity tensor. At first glance, this looks very much like C (2) 1 in the sense that the numerator looks like a 2-point correlation function with β = 2. There is a crucial difference between the behavior of sin 2 θ ij and θ 2 ij , however, such that the C-parameter vanishes for dijet configurations when the jets are back-to-back (i.e. θ ij = π). If we expand around the dijet limit, then the C-parameter really behaves like a 3-point correlation function (i.e. like C 2 ). To see this, note that for e + e − → qqg, the C parameter has the simple form where x i = 2p i · Q/Q 2 and Q is the total four-vector of the system. 6 The angle between final-state particles i and j in the e + e − → qqg system is (2.14) Thus, if we change θ → 2 sin θ 2 in the definition of ECF(N, β), then C 2 can be expressed as which, up to normalization, is the traditional C-parameter. Of course, at higher orders in perturbation theory the definitions of the C-parameter and C (2) 2 diverge. Both observables are sensitive to 2-prong (sub)structure, though C (β) 2 gives sensible answers even for systems with non-zero total momentum and has an adjustable angular exponent β.
Higher-point energy correlation functions have been studied very little in the literature. Two early studies for e + e − collisions are in Refs. [62,65]. However, both define observables that only make sense for systems with total momentum equal to zero and explicitly use operations only defined in three-dimensional space, such as cross-products and the properties of momentum tensors with rank greater than 2. Thus, these observables cannot be easily generalized to determine if a (boosted) system has N (sub)jets. Historically, observables like the D-parameter [61,62,66] have been used to identify peculiar phase space configurations such as a planar configuration of particles. However, this is not directly related to the number of jets in the event. Recent substructure variables like planar flow [7,9], Zernike coefficients [52], and Fox-Wolfram moments [53] are similarly sensitive to peculiar phase space configurations rather than prong-like substructure. Planar flow, for example, vanishes if the constituents of the jet lie along a line, which is a good probe for some (but not all) 3-prong configurations. The energy correlation double ratio C N is designed to directly probe N -prong configurations, though the high computational cost of ECF(N + 1, β) likely limits the practical range to N ≤ 3 (i.e. up to three-prongs).

Advantages Compared to N -subjettiness
The variable N -subjettiness [49,50] (based on N -jettiness [67]) is a jet observable that can be used to test whether a jet has N subjets, and it has been used in a number of theoretical [14,20,[68][69][70][71][72] and experimental [27,31] substructure studies. Since both N -subjettiness and the energy correlation double ratio C N share the same motivation, it is worth highlighting some of the advantages of the energy correlation double ratio.
First, a quick review of N -subjettiness. It is defined in terms of N subjet axesn A as 7 τ (β) where the sum runs over all particles in the jet and R A,i is the distance from axis A to particle i. There are a variety of methods to determine the subjet directions, with arguably the most elegant way being to minimize τ N over all possible subjet directionsn A [50]. If a jet has N subjets, then τ N −1 should be much larger than τ N , so the observable that is typically used for jet discrimination studies is the ratio (2.17) As discussed above, this ratio is directly analogous to the energy correlation double ratio C (β) N −1 . One immediate point of contrast between N -subjettiness and the energy correlation double ratio is that C N does not require a separate procedure (such as minimization) to determine 7 In Refs. [49,50], N -subjettiness was defined with an overall normalization factor to make it dimensionless.
Here, we remove the normalization factor so it has the same dimensions as Eq. (2.8). Figure 1: Example kinematics with soft wide-angle radiation. Left: recoil of the jet axis (dashed) away from the hard jet core (E 1 ) due to soft wide-angle radiation (E 2 ), which is relevant for small values of β. Right: a three-particle configuration that highlights the difference between C 2 and τ 2,1 .
the subjet directions. While novel, this by itself does not necessarily imply that C N will have better discrimination power than τ N,N −1 , though it does mean that C N is a simpler variable to study. 8 We now explain two test cases where C N can perform better than τ N,N −1 : insensitivity to recoil for C 1 and sensitivity to soft wide-angle emissions for C 2 .

Insensitivity to Recoil
A recoil-sensitive observable is one for which soft emissions have an indirect effect on the observable. In addition to the direct contribution to the observable, soft radiation in a recoilsensitive observable changes the collinear contribution by an O(1) amount. An example of a recoil-sensitive observable is angularities for the angular exponent a ≥ 1 (β ≤ 1), which was studied in Ref. [54]. Because C N is insensitive to recoils, it is better able to resolve the collinear singularity of QCD.
For 1-prong jets, the effect of recoil on an observable is illustrated in Fig. 1a. Because of conservation of momentum, soft wide-angle radiation displaces the hard jet core from the jet axis. Angularities (i.e. 1-subjettiness) are sensitive to this displacement since they are measured with respect to the jet center. For a jet with two constituents separated by an angle θ 12 (using the notation in Eq. (2.1) for simplicity), In particular, β serves two different roles for N -subjettiness. As in C (β) N , β controls the weight given to collinear or wide-angle emissions. In addition, when the minimization procedure is used, β controls the location of the axes which minimize N -subjettiness. When trying to determine the optimal value for β for subjet discrimination, it is difficult to disentangle these two effects.
Taking the E 2 E 1 limit one can view the first term as the contribution directly from the emission E 2 , while the second term comes about because particle 1 recoils when it emits particle 2. The dependence of τ (β) 1 on the energies and emission angle is different according to the value of β. For β > 1, the second term is negligible, and the angularities become such that τ 1 is linear in the soft radiation E 2 . However, for smaller values of β, the expression for angularities changes because recoil effects become important. For β = 1, both terms are identical in the E 2 E 1 limit and angularities become For β < 1, the first term is negligible in the E 2 E 1 limit, and the angularities are dominated by the effect of recoil of the hard radiation By contrast, the observable C has the same behavior for all values of β: which is dominated by the splitting angle and energy of the softer particle in the jet for all values of β > 0. Because N -subjettiness is essentially a sum over N subjet angularities, τ N can suffer from the same recoil-sensitivity as angularities for β ≤ 1, depending on how the subjet axes are determined. For example, if N -subjettiness is defined using k T subjet axes, then τ N is recoil sensitive. N -subjettiness is also recoil sensitive if the subjet axes are aligned with the subjet momenta. A related issue is that if the subjet axes are determined using the minimization procedure, then the minimization is only well-behaved for β ≥ 1. 9 In all of these cases, the useful range of β is limited to β ≥ 1. In contrast, the energy correlation double ratio is recoil-free and well-behaved for the whole IRC-safe range β > 0. As we will see in Sec. 3 (and demonstrated recently in Ref. [73]), this is relevant for quark/gluon discrimination, where β 0.2 for C 1 is the optimal choice.
It should be noted that one can construct a recoil-free version of N -subjettiness where the subjet axes are always chosen to minimize the β = 1 measure (see forthcoming work in Ref. [60]), regardless of which β is used in Eq. (2.16). We refer to axes defined in this way as "broadening axes", since β = 1 is closely related to the jet broadening measure [74]. We will make use of this fact later when comparing C 1 to 1-subjettiness in Sec. 3.3.

Sensitivity to Soft Wide-Angle Emissions
Another point of contrast between C N and τ N,N −1 is in how the two variables behave in the presence of emissions at multiple angular scales. The way N -subjettiness is defined, every jet is partitioned into N subjets, even if there are fewer than N "real" subjets. For example, when a jet has a soft subjet separated at large angle (as one might expect from the radiation off a quark or gluon), N -subjettiness will still identify that soft subjet region, yielding a relatively low value of τ N,N −1 (and therefore making the jet look more N -prong-like than it really is). In contrast, because the energy correlation function is sensitive to all possible soft and collinear singularities, C N takes on a relatively high value in the presence of a soft wide-angle subjet, making the jet look less N -prong like (as desired).
We can show this concretely for C 2 using the configuration in Fig. 1b where there is the following hierarchy of the energies and angles: 10 (2.23) Again using the notation in Eq. (2.1), the energy correlation functions are (2.25) For N -subjettiness with three jet constituents, it is consistent to choose axes that lie along the hardest particle in a subjet. For 1-subjettiness, the axis lies along particle 1. For 2subjettiness, one axis lies along particle 1 and the other axis lies along particle 2 or particle 3, depending on the relationship between E 3 θ 13 and E 2 θ 12 . This gives Regardless of the ordering of E 3 θ 13 and E 2 θ 12 we see that: so in the presence of a soft subjet at large angle θ 12 , C 2 yields a larger value than τ 2,1 (i.e. more background-like as desired). As we will see in Sec. 4, this allows C 2 to perform better than τ 2,1 for background rejection in regions of phase space where soft wide-angle radiation plays an important role.
One way to understand the improved performance of C 2 with respect to τ 2,1 is to consider the concrete example of β = 2 at fixed jet mass m. 11 Using the kinematic limit above, the jet mass-squared is given approximately by and it is convenient to define z as the energy fraction of the emission that dominates the mass (e.g. . For fixed jet mass, QCD backgrounds tend to peak at small values of z, but we see from Eq. (2.26) that τ 2,1 does not have any z-dependence for fixed jet mass. For C 2 , if particle 2 dominates the mass (i.e. if a soft wide-angle emission dominates the mass), then so C 2 penalizes small values of z. In this way, C 2 acts similarly to taggers that reject jets if the kinematics of the dominant splitting of the jet is consistent with background [3-6, 11, 12]. In contrast, τ 2,1 only exploits the degree to which radiation is collimated with respect to the two subjet directions, and does not take into account the z-dependence at fixed jet mass. If particle 3 dominates the mass (i.e. if the mass is dominated by a hard core of energy), then C 2 is constant in the energy fraction z, and so is no longer affected by the kinematics of the emission that generated the mass. However, there is still the potential for improved performance in identifying boosted color singlet resonances like Z bosons. For a boosted Z boson, emissions at wide angle with respect to the angle between decay products are suppressed by color coherence. As one goes to higher boosts where the ratio of jet mass to jet p T decreases for fixed jet radius, the volume of phase space for allowed emissions decreases, which can also be seen as a consequence of angular ordering. It is therefore less likely for a Z boson signal to generate final state radiation at large θ 12 , while background QCD jets will emit at large angle independently of the p T . Because radiation at large angles has an enhanced effect on C 2 as compared to τ 2,1 , cf. Eq. (2.27), we expect C 2 to be more effective at discriminating color-singlet signals from background QCD jets.

Quark vs. Gluon Discrimination with C 1
Our first case study is to use the energy correlation functions to discriminate between quark jets and gluon jets. The observable C 1 contains the 2-point energy correlation function ECF(2, β) and so is sensitive to radiation in a jet about a single hard core. 12 This case study is simple enough that we can predict the quark/gluon discrimination power through an analytic calculation, which we will subsequently validate with Monte Carlo simulations. In our later case studies involving higher-point correlators, we will rely on Monte Carlo alone. 11 We thank Gregory Soyez for helpful discussions on these points. 12 The CMS experiment uses an observable they call pT D = i p 2 ti /( i pti) 2 for quark versus gluon discrimination [39,75]. It is related to the β = 0 limit of C In any discussion of quark-gluon discrimination, one should start with a reminder that defining what is meant by a quark or a gluon jet is a subtle task, since the one existing infrared-safe way of defining quark and gluon jets [76] works only at parton level. Existing work on practical aspects of quark-gluon discrimination in Refs. [39,73,75,77,78] has not entered into these issues. Instead the discussion has relied on Monte Carlo simulations, defining a quark (gluon) jet to be whatever results from the showering of a quark (gluon) parton. We will adopt a variant of this methodology in our Monte Carlo studies. Our analytic approach will instead define a quark or gluon jet in terms of the sum of the flavors of the partons contained inside it. It is based on resummation and therefore contains similar physics to the Monte Carlo parton shower.

Leading Logarithmic Analysis
We begin our analysis by considering the leading logarithmic (LL) structure of the cross section for the observable C 1 . With L equal to the logarithm of C 1 , we define LL order as including all terms in the cross section that scale like α n s L 2n , for n ≥ 1. At LL order, quark versus gluon jet discrimination can be understood as a consequence of quarks and gluons having different color charges. To LL order, the strong coupling constant α s can be taken fixed and only the most singular term in the splitting function need be retained. With only one soft-collinear gluon emission, the normalized differential cross section for any infrared and collinear safe observable e has the same form for both quark and gluon jets: where C is the color factor, R 0 is the jet radius, 13 z is the energy fraction of the emitted gluon, θ is its splitting angle, andê is a function of z and θ. Recall that C F = 4/3 for quarks and C A = 3 for gluons. At this order, the observable C which takes a maximum value of We identify the logarithm L as which we use in the following expressions for compactness. This distribution can be resummed to LL order by exponentiating the cumulative C (β) 1 distribution. The resummed distribution 13 We use this somewhat non-standard notation because R will later be used with a different meaning.
that follows is then Because the quark color factor is smaller than the gluon color factor, the Sudakov suppression is less for quarks. Thus, the C (β) 1 distribution for quark jets is peaked at smaller values than for gluon jets.
To figure out the quark/gluon discrimination power from this C (β) 1 resummed distribution, we will make a sliding cut on C (β) 1 and count the number of events that lie to the left of the cut. Adjusting this cut then defines a ROC curve relating the signal (quark) jet efficiency to the background (gluon) jet rejection. To LL accuracy, the (normalized) cumulative distributions for quarks and gluons are: (3.6) Note that at LL order, there is a simple relationship between these cumulative distributions: Thus, if a sliding cut on C (β) 1 retains a fraction x of the quarks, it will retain a fraction x C A /C F of the gluons. The quark/gluon discrimination curve is then which (perhaps surprisingly) is independent of β. This LL discrimination result holds for a wide class of IRC safe observables sensitive to the overall jet color factor, including the jet mass. Only beyond LL order does the discrimination curve depend on β.

Next-to-Leading Logarithmic Analysis
We continue our analysis to next-to-leading logarithmic (NLL) order, which we define as including all terms that scale as α n s L n+1 and α n s L n in ln Σ. In addition, we will also include the non-logarithmically enhanced term arising at O(α s ). At NLL order, there are several new effects that must be included, which together turn out to improve the quark/gluon discrimination power of C (β) 1 compared to the LL estimate. The dominant effects are subleading terms in the splitting functions and phase space restrictions due to multiple emissions. In addition, one must account for the running of α s , fixed-order corrections, and non-global logarithms [79] arising from the phase space cut of the jet algorithm. We will consider how these affect the discrimination power of C (β) 1 , ultimately showing that small values of β improve quark/gluon discrimination. We will work in an approximation of small jet radius, R 0 1, which will allow us to consider only the effects of radiation from the jet, while neglecting modifications associated with the full antenna structure of initial and final-state partons.
The resummation to NLL for generic (global) observables was carried out in Ref. [54]. The central result of that analysis was an expression for the NLL cumulative distribution for an arbitrary observable (satisfying certain basic conditions, e.g. recursive infrared safety). From Ref. [54], the probability that the value of an observable is less than e −L takes a form such as where N is a matching factor to fixed order, N = 1 + O (α s ), and γ E 0.5772 is the Euler-Mascheroni constant. In a fixed-coupling approximation, the "radiator" function R for the observable C where C is the color factor of the jet and B encodes subleading terms in the splitting functions. 14 For quark jets B q = − 3 4 and for gluon jets where n f is the number of light quark flavors. The specific NLL resummed formula in Eq. (3.9) holds for observables that are global, recursively infrared and collinear safe (rIRC), and additive. The last two conditions are satisfied by C (β) 1 . The general expression for R with running α s appears in Ref. [54]. The scale at which α s is evaluated is p T R 0 , and we will use the shorthand unless an explicit scale is used as the argument of α s . Because C for a jet is non-global, it is necessary to include an extra factor in the resummation, discussed in detail in Sec. 3.2.3. We will also include information obtained by matching to the O(α s ) fixed-order cross section, where the matching procedure is described in App. A.
Armed with the matched NLL cumulative distribution, including the non-global and O (α s ) corrections, we can now determine the quark versus gluon discrimination curve by numerically inverting Σ q and plugging it into the expression for Σ g . This is shown for various values of β in Fig. 2a. In Fig. 2b, we fix 50% quark efficiency and show the gluon rejection rate (i.e. one minus the gluon efficiency) as a function of β for R 0 = 0.6. Also on this plot is an approximate analytic expression for the rejection rate as a function of β that we derive below in Eq. (3.22). We see that the discrimination power improves as β decreases. It is, however, not sensible to take β too small: for β = 0 our observable is collinear unsafe, and large non-perturbative effects can be expected as β approaches zero. Furthermore for β α s the convergence of our calculation breaks down (cf. App. B).
To understand the behavior of Fig. 2b semi-analytically, we will study the impact of different physical effects on the discrimination. To do so, we will again express Σ g in terms of 14 To obtain Eq. (3.10), we used the fact that, for a general jet observable that takes the form where Ri is the angle of the emission, Eq.   1 , calculated at NLL order matched to fixed order for various values of β. The β-independent LL prediction is shown for comparison. Right: Gluon rejection rates at 50% quark efficiency, as a function of β, demonstrating that β 0.2 is optimal at NLL order (for smaller values of β, non-perturbative effects become important). Also shown is an analytic approximation from Eq. (3.22) (C 1 Approx.) that includes the most important physics that enters at NLL. Σ q so as to determine the discrimination power of a cut on C 1 . In fact, we are most interested in the exponent relating Σ g to Σ q (as in Eq. (3.7)), so we will actually relate the logarithms of the two cumulative distributions to one another. We are interested in the regime where ln 1/Σ ∼ 1, which, from Eq. (3.6), implies that α s L 2 ∼ 1. The logarithm of the cumulative distribution has the schematic expansion

Quark vs. Gluon NLL LO
With the power counting of α s L 2 ∼ 1, we will consider all terms from Eq. (3.12) that scale as α 0 s , α 1/2 s , or α 1 s . This corresponds to all terms at order α s from Eq. (3.12), as well as the terms at α 2 s L 3 , α 2 s L 2 , and α 3 s L 4 . To illustrate this power counting, consider, for example, the term α s L, which scales as α 1/2 s as one varies α s while keeping α s L 2 fixed and of order 1. In what follows we will pay special attention to the terms at order α s L and α 2 s L 2 , which turn out to be the most relevant ones when establishing deviations from our LL analysis and whose dominant contributions have clearly identifiable physical origins. The terms at order α 2 s L 3 and α 3 s L 4 are simply proportional to the LL color factor, multiplied by powers of the β-function, and so do not significantly modify the LL analysis.

Subleading Terms in Splitting Functions
We first consider the effect on the discrimination from the subleading terms in the splitting functions. In the observable C (β) 1 , β controls the weight given to collinear and wide-angle emissions in the jet. At large values of β, wide-angle emissions are given greater weight, and at small values of β, collinear emissions are given greater weight. Wide-angle soft radiation is controlled by the term in the splitting function that diverges as the energy fraction goes to zero; i.e., the term dz/z. Both quarks and gluons have the same functional form for the soft limit of the splitting function, with the only difference being the overall color factor. By contrast, collinear emissions are controlled by the subleading terms in the splitting function, which differ for quarks and gluons (i.e. different values of the B coefficient). Therefore, as β goes to zero and the collinear emissions become more important in C (β) 1 , the differences between the quark and gluon splitting functions are accentuated.
To see this behavior directly from Eq. (3.9), we can ignore the R -dependent prefactor and focus on the e −R factor. We can write which is 1 9 for n f = 5. We then have (3.14) This last form allows us to relate the cumulative distribution for gluons to that of quarks, in the same spirit as Eq. (3.7): This implies that the separation between the quark and gluon distributions increases as β decreases and so smaller values of β result in better discrimination. Because this effect first arises at O √ α s , there will be corrections at O (α s ) due to the running coupling. Note also that the coefficient δB is quite small in QCD, and so the total effect from the subleading terms in the splitting functions on the discrimination power is minimal.

Multiple Emissions
Next, consider the effect of multiple emissions. The Sudakov logarithm corresponds to the integral of the area (in ln k t , ln θ space) over which emissions are forbidden. At LL, any number of emissions can lie arbitrarily close to the lower boundary of the phase space region without changing the value of the observable. At NLL, one must consider the cumulative effect of the emissions that lie near the phase space boundary. Multiple emissions tend to increase the value of the observable C 1 , and so, for a fixed value of C 1 , they must be suppressed. This introduces an extra degree of discrimination between quarks and gluons; there are likely to be more such emissions for gluons than quarks and so it costs more to "accept" a gluon jet.
For a given LL Sudakov factor, the extent of the boundary region is effectively increased as β is decreased, leading to better quark versus gluon discrimination at small β.
In Eq. (3.9), the effect of muliple emissions is seen in the R -dependent prefactor. For small values of R , the prefactor has the expansion (3.16) We will drop terms at O α 2 s LR/β 2 and higher, which constrains us to consider β α s L. The cumulative distribution can then be written approximately as which allows us to relate Σ g in terms of Σ q as This again suggests an increase in discrimination power for relatively small β. While this effect appears at order α s rather than √ α s , it has a substantially larger coefficient.

Non-Global Logarithms
Because jets are defined in a restricted phase space, non-global logarithms may contribute to the quark versus gluon discrimination power. The effect of non-global logarithms on the cumulative distribution can, for our purposes, be approximated in the large-N C limit as [79][80][81] Σ with NG = e −CC A π 2 12 α 2 s π 2 L 2 Σ = e −C A π 2 12 αs π βR Σ . This neglects some contributions starting at order α 3 s L 3 in the exponent, but these would not affect the quark-gluon discrimination at our accuracy. Recently a first numerical calculation has been performed including the full-N C structure [82] and it suggests that finite-N C corrections are small.
If we temporarily ignore the R -dependent prefactor in Eq. (3.9), the inclusion of nonglobal logarithms leads to (3.20) All quark/gluon dependence resides in the color factor inside R, so we still have the property from the LL calculation (again, ignoring the prefactor and setting δB = 0) Hence non-global logarithms do not modify the above arguments in any significant way. This analysis holds for the anti-k T jet algorithm, whose boundary is unaffected by soft radiation at angles ∼ R 0 . For other algorithms of the generalized-k T family, which have irregular, soft-emission-dependent boundaries, there are additional terms, clustering logarithms [83,84], which also appear starting from order α 2 s L 2 . Some of the O α 2 s clustering logarithms involve color factor combinations such as C 2 F and C 2 A for quarks and gluons respectively, and so presumably would have an impact on quark-gluon discrimination at our accuracy. We leave the study of these terms for future work.

Summary of NLL Result
Using the results of Ref. [54] and App. A to include all effects up through O(α s ) in the logarithm of the cumulative distribution we find This expression includes two terms beyond those discussed in the subsections above. The one proportional to b 0 , where b 0 = 11 3 C A − 2 3 n f is the one-loop β-function coefficient, has two origins: it comes from the running coupling corrections to the contribution from the subleading terms in the splitting functions and from the running-coupling corrections to the relation between the logarithm L and ln Σ q . The last term in parentheses comes from O (α s ) matrix element corrections, discussed in detail in App. A. It depends on the choice of the jet definition, including the procedure by which one defines quark versus gluon jets at parton level. Specifically, we assume any algorithm is equivalent to the generalized k T family of jet algorithms at order α s , and at this order define a quark jet to be one that contains a quark and a gluon, while jets containing gg or qq are considered to be gluon jets. 15 Beyond O (α s ), the calculation assumes that the algorithm maintains a rigid circular boundary in the presence of multiple soft particles at angles of order R 0 , i.e. that it behaves like the anti-k T algorithm.
Note that every subleading term in Eq. (3.22) is proportional to a difference of color or quark number factors and so the discrimination power depends sensitively on these differences. The overall quark versus gluon discrimination power increases as β is decreased (even though the last term favors larger values of β for n f > C A ). Numerically, this behavior is dominated by the subleading terms in the splitting functions and the multiple emissions effect. The effect of the subleading terms in the splitting functions goes like √ α s and so is formally more important than the multiple emissions effect which is O(α s ). However, the effect of 15 In contrast to the situation with LO studies, at O (αs) it does not makes sense to discuss jet flavor based on the flavor of the parton that "initiates" the jet, since interference effects between diagrams mean that the initiating parton cannot be uniquely identified. The question of quark-gluon jet definition at fixed order is discussed further in App. A.
the subleading terms is multiplied by the small number δB and so is numerically smaller than the contribution from multiple emissions. Running-coupling and fixed-order effects are significantly smaller. Robustly, then, smaller values of β lead to better discrimination between quark and gluon jets. One explicitly sees that we have an expansion in powers of α s /β, and so it can only be trusted for β substantially larger than α s ; in practice, perhaps β (2 ∼ 4) × α s (see App. B). It is interesting to comment also on traditional angularities: for β > 1 most of Eq. (3.22) still holds, and only the last term in parentheses would be modified. However, for β ≤ 1 angularities are dominated by recoil effects, with a structure that is independent of β, and so we expect that the discrimination should saturate. Because the energy correlation double ratio C (β) 1 is recoil-free for all values of β, it is better able to probe the collinear singularity and multiple emission effects that distinguish quarks from gluons.

Monte Carlo Study
We now use a showering Monte Carlo simulation to validate the above NLL analysis of C (β) 1 . A similar study of the EEC function appears in Ref. [73], where it was called the two-point moment. 16 Through this paper, jets are identified with the anti-k T algorithm [85] using FastJet 3.0.3 [63]. No detector simulations are used other than to remove muons and neutrinos from the event samples before jet finding, as was done in analyses for the BOOST 2010 report [1].
We generate pure quark and gluon dijet samples from the processes qq → qq and gg → gg in Pythia 8.165 [86,87] at the 8 TeV LHC using tune 4C [88]. While Pythia is not fully accurate to NLL, it does include subleading terms in the splitting functions and multiple emissions, so not surprisingly we find improved discrimination at smaller values of β, in agreement with Sec. 3.2. We scan over various jet radii and p T cuts to study the dependence of the quark/gluon discrimination on these parameters. For this study, we only use the hardest reconstructed hadron-level jet in the event with a transverse momentum in the ranges of p T ∈ [200, 300] GeV, [400, 500] GeV, or [800, 900] GeV. 17 If the hardest jet in the event lies outside the p T range of interest, the event is ignored. In addition, we scan over jet radii values 16 Ref. [73] examined the C1 quark-gluon discrimination for a range of β values and reached a conclusion that is consistent with ours. While their initial analysis naïvely suggests that C1, Fig. 18, performs worse than jet broadening ("girth", or equivalently τ1 with β = 1), Fig. 13, that comparison involves different Monte Carlo event samples. Table 1 of Ref. [73] compares the observables on equal footing, which shows that C1 indeed has better discrimination power than jet broadening, consistent with our discussion here. 17 The reason for focusing only on the leading jet is that we want to minimize ambiguities related to defining quark and gluon jets. The subleading jet is the one more likely to have undergone radiation, and with radiation, quark jets may change into gluon jets, and vice-versa. Additionally, the local emission environment is changed (e.g. non-global logs may become more important). The probability that an event has radiation in the vicinity of the subleading jet is O (αs), while it is O α 2 s near the leading jet. As a cross-check on the flavour composition of our events, we have clustered the parton-level showered events with the flavor-kt algorithm [76]. We find that the flavor of the leading jet is consistent with expectations except in a small fraction of events, between a few percent and ten percent depending on the generator. Quark vs. Gluon Pythia8 p T 400,500 GeV, R 0 0.6 for several values of β in Pythia. Also plotted is the leading log approximation for the discrimination curve, Eq. (3.8). of R 0 = 0.4, 0.6, and 0.8. Because our broad conclusions hold for all samples generated, we only show representative plots to illustrate the quark/gluon performance of C 1 .
In Fig. 3a, we plot the distribution of C . The jet sample is the same as used in Fig. 3b. Right: Gluon rejection rate for 50% quark efficiency as a function of β, for angularities, 1-subjettiness measured with respect to the broadening axis, and C (β) 1 . The broadening axis is defined as the axis which minimizes the β = 1 measure in N -subjettiness. The latter two observables are recoil-free, and therefore give better discrimination power for small values of β.
transverse momentum in the range [400, 500] GeV and jet radius R = 0.6 in Pythia. As expected, the gluon curve lies at larger values than the quark curve because of the greater Sudakov suppression in gluon jets. The quark/gluon discrimination curves for different values of β are shown in Fig. 3b, which are directly comparable to the NLL results in Fig. 2, up to jet contamination effects included in Pythia such as underlying event and initial-state radiation. Again, we see that β 0.2 is the optimal value. In Fig. 4, we show the gluon rejection rate for 50% quark efficiency as a function of β, comparing different p T ranges and R 0 values, all of which favor small values of β. Note that the gluon rejection power degrades as the jet radius is increased, exhibited in Fig. 4a. This may be associated with the increase in the amount of underlying event and initial-state radiation captured in the jet as the jet radius increases. This radiation is uncorrelated with the dynamics of the quark or gluon which initiates the jet. The degradation is most prominent at large values of β, where wide angles in the jet are emphasized (which is where most of the uncorrelated radiation resides).
To compare the discrimination power of C to other IRC safe observables, we consider 1-subjettiness τ  in Herwig++ (directly comparable to Fig. 5b). We also tested Pythia 6.425 and Herwig 6.520, whose results lie in between Pythia 8 and Herwig++. and β = 1 is jet broadening or girth. Among the angularities, Ref. [77] found that jet broadening (β = 1) was the most powerful angularity for quark/gluon discrimination, and so is a natural benchmark to compare to C is a recoil-free observable and is therefore expected to behave similarly to C (β) 1 . In Fig. 5a we plot the discrimination curves for angularities (i.e. 1-subjettiness measured with respect to the jet axis) for several values of β, as well as the discrimination curve for C (0.2) 1 in Pythia. Indeed, for most of the range, the most discriminating angularity is β = 1, but the performance of all angularities is roughly comparable to and only somewhat better than the LL expectation. By contrast, C (0.2) 1 yields a quark to gluon efficiency ratio that is about twice as large as any of the angularities over much of the range. In Fig. 5b, we highlight the importance of working with recoil-free variables, by plotting the gluon rejection rate at a fixed 50% quark efficiency. For β ≥ 1, C (β) 1 and 1-subjettiness have essentially the same performance. As β approaches 0, however, the discrimination power for the angularities degrades, while the two recoil-free observables (C (β) 1 and 1-subjettiness with respect to the broadening axis) have improved performance, as expected from the NLL analysis. 18 To verify the claims made about the performance of C (β) 1 as a quark/gluon discriminator, we also simulate quark and gluon dijet samples in Herwig++ 2.6.3 [89,90]. We use the same kinematic cuts and jet algorithm parameters as in the Pythia samples. As the same quali- 18 The reason for the mismatch between C1 and τ1 with respect to the broadening axis at very small values of β has not yet been determined. tative conclusions hold in the Herwig++ samples as in Pythia, we only reproduce Figs. 3b and 5b for the Herwig++ sample. In Fig. 6a, we plot the quark versus gluon discrimination curve with C (β) 1 . While the discrimination power of C (β) 1 in the Herwig++ sample is not as great as in the Pythia sample, the behavior that the discrimination increases as β decreases is robust. In Fig. 6b, we plot the gluon rejection rate at a fixed 50% quark efficiency for the three observables considered earlier. As in Fig. 5b, the discrimination power of the recoilfree observables increases as β decreases and degrades for the recoil-sensitive angularities (though the β-dependence is once again weaker than with Pythia 8). We also tested C (β) 1 using Pythia 6.425 [86] with tunes DW and Perugia 2011 [91] and Herwig 6.420 [92] with JIMMY [93], which exhibit discrimination power that is intermediate between Pythia 8 and Herwig++. We also checked that the behavior is robust as underlying event, initial-state radiation, and hadronization are sequentially turned off.
Of course, there is a substantial numerical difference between Pythia 8 and Herwig++ for quark versus gluon discrimination. Some of distinction between Pythia 8 and Herwig++ could be due to the fact that different evolution variables are used: Pythia 8 is a p T -ordered shower whereas Herwig++ is angular ordered. This could in turn affect the flavor content of the quark and gluon jets, thus leading to variations in their ability to discriminate quark and gluon jets. The energy correlation function observables seem to be particularly sensitive to these differences, especially at relatively small values of the angular exponent β. This suggests that any detailed study of the properties of quark and gluon jets should measure C 1 at multiple values of β. Beyond discriminating quark and gluon jets, measurements of energy correlation functions at both e + e − and hadron colliders could be useful for Monte Carlo tuning, especially given the current differences between generators.
In these studies, C 1 has been measured on jet samples which include both charged and neutral hadrons (and we have not applied realistic smearing of energies and angles). In order to exploit the discrimination power of C 1 with β 0.2, one needs excellent angular resolution, so in an experimental context, it may be advantageous to measure C 1 using only charged hadrons. A track-only calculation of C 1 , using e.g. the methods of Ref. [94], is beyond the scope of this work, but we did verify in Monte Carlo that the quark/gluon discrimination power only degrades by a few percent when using a track-only version of C 1 . We also note that CMS has successfully made use of p t D, related to the β → 0 limit of C 1 , using both charged and neutral hadrons [39,75].
Because we observe significant differences in the absolute scale of the quark versus gluon discrimination between different Monte Carlo generators, the performance of C 1 in an experiment may not be as optimistic as computed to NLL. However, the increase in the discrimination power as β → 0 seems robust and would be important to verify in data. As discussed in App. B, perturbative control over C 1 ceases to exist for β 0.2 − 0.4. While hadronization will then become significant, separately on the distributions of C 1 for both quark and gluon jets, it is not entirely trivial to relate this to the question of its expected impact on the quark/gluon discrimination performance. In any case, this kind of questions deserves further investigation, both theoretically and experimentally.

Boosted Electroweak Bosons with C 2
Our next case study is using C 2 to discriminate between QCD jets and jets with two intrinsic subjets, such as boosted W , Z, or Higgs bosons. 19 Recall that C 2 involves the 3-point correlator. To identify a boosted resonance, one first looks for jets whose mass is compatible with the resonance of interest. Then C 2 can be used to determine if the jet has two hard subjets, in which case the jet is tagged as coming from that boosted resonance. While we have not carried out analytic calculations to guide us to understand the performance of C 2 (and higher point correlation functions), it is still instructive to study its discrimination power in Monte Carlo. We use Pythia 8 to demonstrate the qualitative behavior and performance of C 2 , though one must of course be mindful of the quantitative differences in Monte Carlo programs seen already in Sec. 3.3. A calculation of C 2 will be left to future work.
The key finding from this section is that this tagging procedure is very sensitive to the ratio of the jet mass to the the jet transverse momentum. This arises because the structure of the QCD background depends strongly on the jet mass requirement, and the behavior of C (β) 2 differs depending on what type of radiation contributes dominantly to the jet mass. For a fixed jet p T , we will find that small values of β 0.5 are better for high mass resonance discrimination, whereas large values of β 2 give the optimal separation at lower masses. In both regimes, C 2 offers better discrimination power than τ 2,1 , with the difference being more pronounced for small m/p T . After describing this physics for a generic boosted 2prong resonance, we will specialize to the case of the Higgs boson, where additional b-tagging information is available.

Dependence on the Mass Criterion
Consider a quark or gluon jet with invariant mass close to the boosted resonance of interest (which we will call a Z boson for concreteness). For jets with mass comparable to their transverse momentum, the mass is dominated by a single, relatively hard, perturbative splitting. Thus, one expects that the QCD jets that can fake a Z boson are those with two relatively hard cores of energy surrounded by soft radiation. These jets are straightforward to analyze in fixed-order perturbation theory (to generate the jet mass) matched to resummed perturbation theory (to generate the radiation pattern for C 2 ), since there is a clear ordering of emissions in the jet. In particular, QCD jets with large mass should appear similar to jets initiated from heavy resonance decay, with differences controlled mainly by the color of the decay products and the phase space of the hard splitting.
For many systems of interest, however, the above analysis is not appropriate. Once the jet mass is less than around a fifth of the jet transverse momentum times R, the mass no longer arises dominantly from a hard perturbative splitting. For jets in the low to intermediate mass ranges, a significant mass can be generated by a single soft emission from a single hard core. At lower masses, the mass of a jet is generated by multiple soft emissions. Jets in the low and intermediate mass regimes require resummation of these soft emissions to accurately model their dynamics as fixed-order perturbation theory is no longer accurate. For this reason, we expect QCD jets in this mass regime to look very different from boosted resonances with two hard cores, and the discrimination power of C 2 should improve as the ratio of the jet mass to the transverse momentum decreases. In addition, as discussed in Sec. 2.2.2, C 2 is better able to exploit the color singlet nature of the Z boson when m/p T is small. To illustrate this, we generate a mixed sample of quark and gluon jets from pp → Zj with the Z decaying to leptons. These are simulated at the 8 TeV LHC in MadGraph5 1.5.0 [109], showered in Pythia 8.165 [86,87], and we identify the hardest anti-k T R 0 = 1.0 jet. In Fig. 7 we plot the invariant mass spectrum of QCD jets in three different transverse momentum bins, p T ∈ [200, 300] GeV, [400, 500] GeV, and [600, 700] GeV. We see that the mass distributions in each p T bin have steeply falling tails extending to masses of about p T /2. In the tail region, we expect fixed-order perturbation theory to accurately describe the origin of mass of the jet. At lower masses, below about p T /5, Sudakov suppression becomes important as the distributions peak and decrease toward zero mass. In this mass regime, fixed-order perturbation theory is no longer adequate to describe the distribution.
This differing origin of the jet mass is reflected in the C (β) 2 distributions. Because small values of C (β) 2 correspond to 2-subjet-like jets, the C (β) 2 distribution moves to lower values as the mass of a QCD jet increases, as shown in Fig. 8a for β = 2 in the p T range [400, 500] GeV. 20 In contrast, for a boosted heavy particle that decays to two partons, the C (β) 2 distribution is relatively insensitive to the resonance mass, since the mass of such a jet comes mostly from two partons from the decay regardless of the boost factor. Shown in Fig. 8b is the signal C  [70,90] GeV, [80,100] GeV, [100,120] GeV, [110,140]  distributions are remarkably similar. 21 In Fig. 9a, we show the QCD jet versus Z boson discrimination curve for m Z = 91 GeV with p T ∈ [400, 500] GeV for several values of β. To see how the physics changes as the resonance mass changes, we plot the QCD rejection rate for 50% boosted Z efficiency in Fig. 9b as a function of β, for m Z = {80, 91, 110, 125, 150, 200} GeV. At low masses, the most powerful discriminant is β 1.5 − 2. This is expected, since large values of β emphasize soft wide-angle emissions where there is more of a penalty for QCD jets in the Sudakov peak. However, we do not have a quantitative way to understand why the discrimination power saturates at β 2, as opposed to even higher values. At intermediate masses, a wide range of β values yield very similar results. At the high masses where QCD jets are in the tail region, the discrimination dependence on β inverts, with the most powerful discrimination for β 0.5. This is likely to do with the same quark/gluon color factor discrimination as in Sec. 3. In particular, high mass QCD jets are formed by a hard perturbative splitting, which is most likely to be a gluon, whereas the Z jet has two hard quark subjets. That said, we have not yet performed a NLL calculation to understand why β 0.5 is preferred, as opposed to even smaller values.
Finally, it is instructive to compare the discrimination power of C 2 to N -subjettiness. The ratio of 2-subjettiness to 1-subjettiness τ identify Z bosons decaying to two jets. To eliminate ambiguities in minimum axes finding at small values of β, we choose to define the subjet axes by those that minimize the β = 1 measure (i.e. the broadening axes). The discrimination curves of τ (β) 2,1 for m Z = 91 GeV is plotted in Fig. 10a, with the C (β) 2 curve with the most discriminating value from Fig. 9a shown for comparison. We also show the QCD rejection rate for 90% boosted Z efficiency in Fig. 10b. At low masses, C 2 performs as well as or better than τ (β) 2,1 over the entire range of β, except at very small values of β. At high masses, the discrimination power of τ (β) 2,1 becomes comparable to C (β) 2 , since both observables lock onto the hard subjets in the Z decay of the massive QCD jet. The increase in the relative discrimination power of C 2 with respect to τ 2,1 as the ratio m/p T decreases is expected from the discussion of Sec. 2.2.2. As m/p T decreases, soft wide-angle subjets become more important for determining the structure of the jet and C 2 emphasizes these emissions more than τ 2,1 .

Boosted Higgs Identification
One key application for 2-prong jet substructure observables is for identifying boosted Higgs bosons in the decay H → bb. Compared to the case of Z bosons, there is additional information from the presence of b quarks (and the resulting B hadrons) in the final state, which can be used to mitigate QCD backgrounds. Thus, to identify boosted Higgs bosons decaying to bottom quarks, we employ three criteria. First, we require the jet to have a mass comparable curve with the best discrimination (β = 1.7). The subjet axes for N -subjettiness are defined as those that minimize the β = 1 measure (broadening axes). Right: QCD rejection rate for 90% boosted Z efficiency as a function of β, sweeping the value of the Z boson mass to m Z = {91, 125, 200} GeV. Because these curves are with 90% Z efficiency, they are not directly comparable to Fig. 9b. Note that as m/p T decreases, the performance of C 2 improves relative to τ 2,1 .
to the Higgs boson. Second, we demand that two B hadrons are tagged in the jet. Third, we use a sliding cut on C (β) 2 to test for two hard subjets in the jet.
Because we demand that the jet have two B hadrons, the largest QCD background to Higgs decays to bottoms is gluon splitting to bottoms. The splitting function g → bb does not have a soft singularity, so the bottoms from this splitting will have comparable energies. This is also the case for Higgs decay, so we do not expect the same discrimination power for Higgs bosons compared to Z bosons studied above. That said, because of the difference in the total color of the jets, there is an additional handle on Higgs versus gluon discrimination; the bottom quarks from the gluon splitting will be in a color octet state, so there will be significantly more radiation at wide angles compared to Higgs jets.
This color octet versus color singlet distinction can be exploited in two ways. First, more wide-angle radiation can be included in the jet if the jet radius is increased. Larger jet radii improve the contrast for C (β) 2 , since more wide-angle radiation is included in the (octet) gluon jets compared to the (singlet) Higgs jets. Second, the value of β can be set to accentuate the importance of wide-angle emissions in the jet. As β increases, more weight is given to wide-angle emissions, further penalizing gluon jets compared to Higgs jets when using C is beyond the scope of this work, but we can get a sense for the discrimination power of C signal pp → ZH to the leading QCD background of pp → Zbb where both bottom quarks happen to be clustered into the same jet. We generate both the signal and background distributions for the 8 TeV LHC in MadGraph5 1.5.0 [109] plus Pythia 8.165 [86,87], with all ground state B hadrons stable to allow for naïve b-tagging of the jets (with 100% efficiency and no mistags). The mass of the Higgs is set to 125 GeV, and the Z is decayed to leptons and the Higgs is decayed to bb. We consider anti-k T jets with various values of the jet radius R 0 = {0.6, 0.8, 1.0, 1.2}. The leading jet is required to have transverse momentum in the range [400, 500] GeV with exactly two B-hadrons as constituents. To approximate realistic b-tagging within jets, we recluster the jet with the k T algorithm to find two exclusive subjets, and we require that each subjet contain exactly one identified B-hadron. Finally, the mass of the jet is required to be in the window of m J ∈ [110, 140] GeV (i.e. within 15 GeV of the Higgs mass). From the leading jet, we compute C (β) 2 for various values of β and determine the discrimination power.
In Fig. 11, we plot the distributions of C 2 for Higgs jets and the QCD background for various jet radii. As expected, the C 2 distributions dramatically increase as the jet radius increases for QCD jets, while they only increase slightly for Higgs jets. In Fig. 12a, we plot the discrimination curves of Higgs jets versus QCD using C (β) 2 for several values of β for the jet radius R 0 = 1.0. As expected, the discrimination power increases both as the angular exponent increases, again, a consequence of the greater amount of wide-angle radiation in the QCD jets. Fig. 12b shows the β dependence on the QCD rejection rate for 50% boosted Higgs efficiency for jet radii of R 0 = {0.6, 0.8, 1.0, 1.2}. The rejection rate increases dramatically as the jet radius increases. At small jet radii, large values of β lead to the best discrimination, as large β emphasizes wide-angle emissions which differ for QCD and boosted Higgs jets. As This may be because a large jet radius will tend to include more initial state radiation or underlying event, which is independent of the dynamics of the jet.
In Fig. 13a, we compare the discrimination performance of the N -subjettiness ratio τ (β) 2,1 to the most discriminating C (β) 2 (β = 2.0) with jet radius equal to R 0 = 1.0. Over the entire range of signal efficiencies, C (2) 2 performs better than τ (β) 2,1 for any value of β. In Fig. 13b, we plot the QCD rejection rate for 50% boosted Higgs efficiency at several jet radii. For a jet with a small jet radius, C (β) 2 performs significantly better than any N -subjettiness, with the distinction decreasing as the jet radius increases.

Boosted Top Quarks with C 3
Our final case study tests the discrimination power of even higher-point correlation functions, namely using C 3 to distinguish boosted top quarks from QCD jets. 22 Unlike the previous two case studies, this observable is significantly more challenging than lower point correlation functions. From a computational point of view, C 3 involves a 4-point correlator, so its computational cost is expensive since it scales like k 4 , where k is the number of particles in the system. That said, our FastJet add-on only requires a few milliseconds to analyze a boosted top event at one value of β. From an analytical point of view, each term in ECF(4, β) involves a product of 4 energies and 4 2 = 6 angles, complicating an understanding of how C 3 behaves in various limits.
We will find that C 3 performs significantly worse than one might expect from the strong performance seen in C 1 and C 2 . While it is possible that this is an artifact of choosing the particular double ratio combination in the definition of C 3 , we suspect that the proliferation of energy and angular factors in ECF(4, β) is reducing the sensitivity of C 3 to any individual soft emission. In particular, for a soft-collinear emission, C 1 and C 2 are independent of the kinematics of the hard structure of the jet. By contrast, even for a soft-collinear emission, C 3 retains dependence on the hard kinematics of the jet. This is because the correlation functions in the ratio defining C 3 are dominated by possibly different subsets of the hard emissions in the jet. Nevertheless, it is illustrative to see that even with these limitations, there is still discrimination power in C 3 .
To study the performance of C 3 as a top tagger, we use the boosted top and QCD background event samples created for the BOOST 2010 workshop [1]. 23 These events come from 7 TeV LHC collisions simulated with Herwig 6.510 [118] with underlying event simulated with JIMMY [93] with an ATLAS tune [119]. The event samples consist of 2 → 2 QCD processes, either all hadronic tt production or dijet production. For direct comparison to other top tagging procedures, we follow the analysis procedures used in Ref. [1]. We identify anti-k T jets with radius R 0 = 1.0 and demand that the jets have p T ∈ [500, 600] GeV. No detector simulation is performed at this stage, other than removing muons and neutrinos before jet finding. We impose three cuts to discriminate top jets from QCD. First, we demand that the 3 to other methods studied in the BOOST 2010 report [1]. The efficiency curves for N -subjettiness (τ 3 /τ 2 ) [50] and the angular correlation function (ACF) [55] were added later. Here, the efficiencies include both the effect of a mass cut as well as a cut on C jets have mass in the fixed window of 160 < m J < 240 GeV, and second, we apply a sliding cut on C (β) 3 . In addition, it was noted in Ref. [14] that ratio observables such as C 3 can be IR-unsafe without an additional cut. We therefore apply a third cut that C (β) 2 > 0.1, which makes C 3 explicitly IR-safe.
In scanning over the range of 0.5 < β < 2.5, we found that the best discrimination over a wide range of signal efficiencies using C (β) 3 is obtained for β = 1.0. This is the same β value that is optimal for N -subjettiness τ (β) 3,2 [50]. A plot of the distribution of C (1) 3 for top jets and QCD jets in the kinematic and mass windows from above is shown in Fig. 14a. The discrimination curves for the different values of β are shown in Fig. 14b, where the quoted efficiencies only include the effect of a cut on the observable C 3 for jets in the mass window of 160 to 240 GeV.
Finally, we compare the performance of C 3 against several other top tagging procedures in Fig. 15. For this plot, the quoted efficiencies include both the effect of the mass cut as well as the effect from a cut on C 3 . While not as powerful as methods like N -subjettiness, the energy correlation function yields comparable discrimination power to other methods. Of course, the performance may be improved by combining information from different values of β, as well as including additional C 2 and C 1 information.

Conclusions
In this paper, we have introduced arbitrary-point energy correlators that are sensitive to Nprong substructure. These correlators are effective when used as part of an energy correlation double ratio C (β) N , though more general combinations deserve further exploration. Through an NLL calculation, we have seen how C 1 yields excellent quark/gluon discrimination, with β 0.2 being most effective at capturing the differences in color charges. We have also shown the power of C 2 for boosted 2-prong objects like Higgs bosons, and the potential power of C 3 for boosted 3-prong objects like top quarks.
Given the explosion of jet substructure methods over the past few years, it is worth asking whether C N is sufficiently novel to merit further experimental and theoretical studies. Indeed, it is a quite unique variable that combines a number of desirable features. Like Nsubjettiness, C N is a variable which tests for N -prong substructure, but can behave more continuously in situations with soft subsets. Like planar flow and related jet shapes, C N can be calculated directly from the energies and angles of the jet constituents without a separate axes finding step, but it is designed for identifying N -prong substructure instead of just exotic kinematic configurations. Finally, like jet angularities, C N is sensitive to higher-order radiation about LO substructure, but because it is a recoil-free observable, it can better probe the collinear physics that distinguishes a jet's color with 0.2 β 1.0. Because C N has a high computational cost for N > 3, we expect C N will be most useful in practice for 1-, 2-, and 3-prong jet studies.
To gain further confidence in the behavior and performance of these observables, further analytic studies are needed. Of particular need is to calculate C 2 for QCD backgrounds. We already saw that the behavior of C 2 for QCD jets depends strongly on the jet's mass over p T ratio, and it is likely that different theoretical descriptions will be needed for C 2 as a function of m/p T . While C 2 is built as a ratio of IRC safe observables, C 2 itself is only IRC safe with a cut on the jet mass (which acts like a cut on the denominator), and it is an interesting question how to best perform NLL resummation for generic ratio observables. Like all jet shape observables, C 2 is sensitive to underlying event, initial state radiation, and pileup, which must be accounted for in determining the optimal β value. Ideally, theoretical progress on understanding C 2 and other jet shapes will match the rapid experimental progress in implementing them, such that jet substructure observables can truly be a robust tool for LHC physics.

A Fixed-Order Calculation
In this appendix, we present the details of the fixed-order calculation of C 1 and matching to the NLL result from Sec. 3.2. The calculation is valid for any jet algorithm that, for configurations involving exactly two partons in some neighborhood, combines two partons into a single jet if they are separated by an angle less than R 0 and otherwise places them in separate jets. At this order we define a quark jet to be a jet that contains a single quark, or a quark and a gluon. A gluon jet is a jet that contains a single gluon, a gluon pair, or a quark-antiquark pair (of identical flavor). 24 24 Algorithms that satisfy the condition for when they pair partons into a single jet include all members of the generalized-kT family, notably the anti-kT algorithm [85]. One subtlety is that the flavor of jets from such algorithms is not infrared safe for configurations with three or more particles in a common neighborhood. There exist algorithms designed specifically to guarantee a safe jet flavor to all orders, the "flavor-kT " algorithms [76]. These, however, have the property that a quark-antiquark pair can be combined into a jet even for angular separations larger than R0, and so they do not yield the same jets at O (αs) as the generalized-kT family and as we assume in the calculation here. We have reason to believe that it is possible to design an algorithm that is both equivalent to generalized-kT at O (αs) and flavor safe to all orders, but leave an investigation of this question to future work.
In the limit where the jet radius R 0 is small, the O(α s ) cumulative distribution of the observable C and we have approximated the full matrix element by the appropriate splitting function, P (z), as is legitimate for R 0 1. The splitting functions are for quarks and To match the fixed-order cumulative distribution to the NLL cumulative distribution, we use the "Log-R" matching scheme [120]. In this matching scheme, the fixed-order corrections are exponentiated with the NLL distribution. The logarithms that appear in the fixedorder expression must be properly subtracted so as to eliminate a double counting with the logarithms that were resummed. Also, at large values of C 1 , we want the distribution to agree with the fixed-order result. This requires "turning off" the logarithms of the resummation properly.
Matching O(α s ) fixed-order to NLL is straightforward. The matching scheme can be written as Σ match = Σ (L) resum e − αs π (R1−G0−G1L−G2L 2 ) . (A.8) Here, R 1 is defined from the fixed-order cumulative distribution as and G 0 , G 1 L and G 2 L 2 are placeholders representing the constant terms, single logarithms and double logarithms that have been resummed, respectively. For quarks, the logarithms are and for gluons the logarithms are The choice of G 0 in Σ(L) resum is arbitrary because these terms are subleading to the NLL resummation. Subtracting these logarithms from R 1 , in addition to the constant terms G 0 , eliminates double counting. Also, to verify that the distribution agrees with the fixed-order result at large values of C 1 , we can shift the argument of the logarithms appropriately to vanish when C 1 = R β 0 4 , which is the maximum value of C 1 . That is, we replace the logarithms in the resummation and subtraction to be We use this expression to determine the quark versus gluon discrimination in Sec. 3.2.

B Breakdown of Perturbative Calculation
In this appendix, we provide some simple quantitative arguments for the breakdown of the perturbative calculation from Sec. 3.2 for values of β less than about 0.2. There are two effects that we will consider: the QCD Landau pole and the breakdown of the independent emission approximation. Of course, there may be other effects that become important at small values of β, but these nevertheless suggest that our perturbative calculation of the quark versus gluon discrimination power ceases to make sense at very small values of β. First, the QCD Landau pole. At small β, the smallest scale Q 0 that the running coupling is sensitive to is Q 0 = p T R 0 e −L/β . (B.1) The perturbative calculation can be trusted when Q 0 Λ NP , where Λ NP is a scale at which α s becomes non-perturbative. We can estimate the value of β at which the non-perturbative effects become important as follows. The logarithm of the observable C 1 can be roughly estimated from the LL, fixed-coupling expression for the cumulative distribution Σ for quarks, where Σ e Demanding that Q 0 > Λ NP and using the expression for L from the above equation we find that β min π ln 1/Σ α s C F suggesting that non-perturbative effects become critical for β 0.2-0.3. One can also perform such an analysis numerically using the full NLL expressions for Σ q , including all runningcoupling effects, and one reaches a similar conclusion. Second, the NLL calculation assumed that emissions could be treated as independent, but multiple emissions cannot be regarded as independent when each emission can take an O(1) fraction of the energy of the jet. That is, if the logarithm of C 1 is not large then our analysis (appropriate for soft-collinear emissions) breaks down. Assuming as above that the cumulative distribution can be written as in Eq. (B.2), the minimal value of β is The precise value of L min at which the soft-collinear analysis is deemed to break down will change this value. Nevertheless, multiple hard, collinear emissions become important and result in a breakdown of the analysis when β is too small. To include the leading effect of energy conservation among emissions, one must match the NLL resummation to NLO (O(α 2 s )) splitting functions.
It should be noted that the fact that the non-perturbative analysis and the multiple emissions analysis give the same ballpark of β min is a coincidence due to the choice of parameters that were made.