Jet Shapes and Jet Algorithms in SCET

Jet shapes are weighted sums over the four-momenta of the constituents of a jet and reveal details of its internal structure, potentially allowing discrimination of its partonic origin. In this work we make predictions for quark and gluon jet shape distributions in N-jet final states in e+e- collisions, defined with a cone or recombination algorithm, where we measure some jet shape observable on a subset of these jets. Using the framework of Soft-Collinear Effective Theory, we prove a factorization theorem for jet shape distributions and demonstrate the consistent renormalization-group running of the functions in the factorization theorem for any number of measured and unmeasured jets, any number of quark and gluon jets, and any angular size R of the jets, as long as R is much smaller than the angular separation between jets. We calculate the jet and soft functions for angularity jet shapes \tau_a to one-loop order (O(alpha_s)) and resum a subset of the large logarithms of \tau_a needed for next-to-leading logarithmic (NLL) accuracy for both cone and kT-type jets. We compare our predictions for the resummed \tau_a distribution of a quark or a gluon jet produced in a 3-jet final state in e+e- annihilation to the output of a Monte Carlo event generator and find that the dependence on a and R is very similar.


Motivation and Objectives
Jets provide troves of information about physics within and beyond the Standard Model of particle physics. On the one hand, jets display the behavior of Quantum Chromodynamics (QCD) over a wide range of energy scales, from the energy of the hard scattering, through intermediate scales of branching and showering, to the lowest scale of hadronization. On the other hand, jets contain signatures of exotic physics when produced by the decays of heavy, strongly-interacting particles such as top quarks or particles beyond the Standard Model.
Recently, several groups have explored strategies to probe jet substructure to distinguish jets produced by light partons in QCD from those produced by heavier particles [1,2,3,4,5,6,7,8], and methods to "clean" jets of soft radiation to more easily identify their origin, such as "filtering" or "pruning" for jets from heavy particles [5,9,10] or "trimming" for jets from light partons [11]. Another type of strategy, explored in [12], to probe jet substructure is the use of jet shapes, which are modifications of event shapes [13] such as thrust. Jet shapes are continuous variables constructed by taking a weighted sum over the four-momenta of all particles constituting a jet. Different choices of weighting functions produce different jet shapes, and can be designed to probe regions closer to or further from the jet axis with greater sensitivity. 1 While such jet shapes may integrate over some of the detailed substructure for which some other methods search, they are better suited to analytical calculation and understanding from the underlying theory of QCD.
In this paper, we consider measuring the shape of one or more jets in an e + e − collision at center-of-mass energy Q producing N jets with an angular size R according to a cone or recombination jet algorithm, with an energy cut Λ on the radiation allowed outside of jets. We use this exclusive characterization of an N -jet final state looking forward to extension of our results to a hadron collider environment, where such a final state definition is more typical. For the jet shape observable we choose the angularity τ a of a jet, defined by (cf. [12,17]), where a is a parameter taking values −∞ < a < 2 (for IR safety, although factorizability will require a < 1), the sum is over all particles in the jet, E J is the jet energy, p T is the transverse momentum relative to the jet direction, and η = − ln tan(θ/2) is the (pseudo)rapidity measured from the jet direction. The jet is defined by a jet algorithm, such as a cone algorithm, the details of which we will discuss below. We complete the calculation for the jet shape τ a for jets defined by cone or recombination algorithms, but our logic and methods could be applied to a wider spectrum of jet shapes and jet algorithms. We have organized our results in such a way that the pieces independent of the choice of jet shape and dependent only on the jet algorithm are easily identifiable, requiring recalculation only of the observable-dependent pieces to extend our results to other choices of jet shapes. Reliable theoretical prediction of jet observables in the presence of jet algorithms is made challenging by the presence of many scales. Logarithms of ratios of these scales can become large and spoil the behavior of perturbative expansions predicting these quantities. These scales are determined by the jet energy ω, the cut on the angular size of a jet R, the measured value of the jet shape such as τ a , and any other cut or selection parameters introduced by the jet algorithm. 1 The original "jet shape," to which the name properly belongs, is the quantity Ψ(r/R), the fraction of the total energy of a jet of radius R that is contained in a subjet of radius r [14,15,16] . This observable falls into the larger class of jet shapes we have described here and for which we have hijacked the name.
Precisely this separation of scales, however, allows us to take advantage of the powerful tools of factorization and effective field theory. Factorization separates the calculation of a hard scattering cross section into hard, jet and soft functions each depending only on physics at a single scale [18,19]. Renormalization group (RG) evolution of these functions between scales resums logarithms of these scales to all orders in α s , with the logarithmic accuracy determined by the order to which the anomalous dimensions in the running are calculated [20]. Effective field theory organizes these concepts and tools into a conceptually simple framework unifying many ingredients going into traditional methods, such as power counting, gauge invariance, and resummation through RG evolution. The rules of effective theory facilitate proofs of factorization and achievement of logarithmic resummation at leading order in the power counting and make straightforward the improvement of results order-by-order in power counting and logarithmic accuracy of resummation.

Soft-Collinear Effective Theory and Factorization
Soft-Collinear Effective Theory (SCET) [21,22,23,24] has been successfully applied to the analysis of many hard scattering cross sections [25] including the production of jets. SCET is constructed by integrating out of QCD all degrees of freedom except those collinear to a lightlike direction n and those which are soft, that is, have much lower energy than the energy of the hard scattering or of the jets. Using this formalism, the factorization and calculation of two-jet cross sections and event shape distributions in SCET were developed in [26,27,28,29]. Later, these techniques were extended to the factorization of jet cross sections and observables using jet algorithms in [30]. Calculations in SCET of two-jet rates using jet algorithms have been performed in [27,31], and more recently in [32]. Calculations of cross sections with more than two jet directions have been given in [33,34,35].
Building on many of the ideas in these previous studies, in this paper, we will demonstrate a factorization theorem for jet shape distributions in e + e − → N jet events, dσ(P 1 , . . . , P N ) dτ 1 · · · dτ M = σ (0) (P 1 , . . . , P N )H(n 1 , ω 1 ; · · · n N , ω N ; µ) (1.2) × J n 1 ,ω 1 (τ 1 ; µ) · · · J n M ,ω M (τ M ; µ) ⊗ S n 1 ···n N (τ 1 , . . . , τ M ; R, Λ; µ) × J n M +1 ,ω M +1 (R; µ) · · · J n N ,ω N (R; µ) , where the N jets have three-momenta P i , and M ≤ N of the jets' shapes τ 1 , . . . , τ M are measured. σ (0) is the Born cross-section, H is a hard function dependent on the directions n i and energies ω i of the N jets, J n,ω (τ ) is the jet function for a jet whose shape is measured to be τ , J n,ω (R) is the jet function for a jet with size R whose shape is not measured, and S is the soft function connecting all N jets, dependent on all jets' shapes τ i , sizes R, and total energy Λ that is left outside of all jets. The symbol "⊗" stands for a set of convolution integrals in the variables τ i between the measured jet functions and the soft function. All terms in the factorization theorem depend on the factorization scale µ. SCET is typically constructed as a power expansion in a small parameter λ formed by the ratio of soft to collinear or collinear to hard scales, determined by the kinematics of the process under study. λ is roughly the typical transverse momentum p T of the constituent of a jet (relative to the jet direction) divided by the jet energy E J . This is set either by the measured value of the jet shape τ a for a measured jet or the algorithm measure R for an unmeasured jet. Thus we encounter in this work the new twist that the size of λ may be different for different jets. We will comment on further implications of this in subsequent sections. Still, in each separate collinear sector, the momentum p n of collinear modes in the light-cone direction n in SCET is separated into a large "label" momentump n containing O(E J ) and O(λE J ) components and a "residual" component of O(λ 2 E J ), the same size as soft momenta. Effective theory fields have dynamical momenta only of this soft or residual scale. This fact, along with the fact that soft quarks and soft gluons can be shown to decouple from collinear modes at the level of the Lagrangian [24], makes possible the factorization of a jet shape distribution into hard, jet, and soft functions depending only on the dynamics at those respective scales.
In using SCET for jets in multiple directions and using jet algorithms to define the jets, we will encounter the need for several additional criteria to ensure the validity of the N -jet factorization theorem.
• First, to ensure that the algorithm does not group final-state particles into fewer than N jets, the jets must be "well separated." This allows us to use as the effective theory Lagrangian a sum of N copies of the collinear part of the SCET Lagrangian for a single direction n and a soft part, and to construct a basis of N -jet operators built from fields from each of these sectors to produce the final state. Our calculations will reveal the precise quantitative condition that jets must satisfy to be "well separated".
• Second, to ensure that the jet algorithm does not find more than N jets, we place an energy cut Λ on the total energy outside of the observed jets. We will take this energy Λ to scale as a soft momentum so that we will be able to identify the total energy of each jet with the "label" momentum on the SCET collinear jet field producing the jet. Corrections to this identification are subleading in the SCET power counting.
• Third (and related to the above two), we will assume that the N -jet restriction on the final state can itself be factorized into a product of N 1-jet restrictions, one in each collinear sector, and a 0-jet restriction in the soft sector. We represent the energetic particles in the ith jet by collinear fields in the SCET Lagrangian in the n i collinear sector and soft particles everywhere with fields in the soft part of the Lagrangian. We then stipulate that the jet algorithm acting on states in the n i collinear sector find exactly one jet in that sector, and when acting on the soft final state find no additional jet in that sector.
• Fourth, the way in which a jet algorithm combines particles in the process of finding a jet must respect the order of steps envisioned by factorization. In particular, factorization requires that the jet directions and energies be determined by the collinear particles alone, so that the soft function knows only about the directions and colors of the jets, not the details of any collinear recombinations. Ideally, all energetic collinear particles should be recombined first, with soft particles within a radius R of the jet axis being recombined into the jet only afterwards. Jet algorithms in use at experiments do not have this precise behavior, but we will discuss in Sec. 3.4 the extent to which common algorithms meet this requirement and estimate the size of the power corrections due to their failure to do so. In general, we will find that for sufficiently large R, infrared-safe cone algorithms and k T -type recombination algorithms satisfy the requirements of factorizability, with anti-k T allowing smaller values of R than k T .
After enforcing the above requirements, a key test of the consistency of Eq. (1.2) will be the independence of the physical cross section on the factorization scale µ. This requires the anomalous dimensions of the hard, jet, and soft functions to sum to zero, 0 = [γ H (µ) + γ J M +1 (R; µ) + · · · + γ J N (R; µ)]δ(τ 1 ) · · · δ(τ M ) + γ J 1 (τ 1 ; µ)δ(τ 2 ) · · · δ(τ N ) + · · · + δ(τ 1 ) · · · δ(τ M −1 )γ J M (τ M ; µ) + γ S (τ 1 , . . . , τ M , R; µ) . (1.3) It seems highly nontrivial that this condition would be satisfied for any number, size, and flavors of jets (and that the soft anomalous dimension be independent of Λ), but we will demonstrate that it does hold at O(α s ), up to corrections of O(1/t 2 ) which violate Eq. (1.3), where t is a measure of the separation between jets. In particular, for a pair of jets, i, j, with 3-vector directions separated by a polar angle ψ ij , the separation t ij is given by t k,l = tan(ψ k,l /2) tan(R/2) . (1.4) Now define t (no indices) as the minimum of t ij over all jet pairs. This quantifies the qualitative condition of jets being well-separated, t 1, that is required to justify the factorization theorem Eq. (1.2). The factorization theorem is valid up to corrections of O(λ) in the SCET power expansion parameter and corrections of O(1/t 2 ) in the separation parameter. As an example of the magnitude of t, for three jets in a Mercedes-Benz configuration (ψ = 2π/3 for all pairs of jets), 1/t 2 = 0.04 for R = 0.7 and 1/t 2 = 0.1 for R = 1, so these corrections are indeed small. More generally, for non-overlapping jets, ψ > 2R, we have 1/t 2 < 1/4.
Notice that for back-to-back jets (ψ = π), t → ∞. Thus, for all cases previously considered in the literature, the jets are infinitely separated according to this measure, and no additional criterion regarding jet separation is required for consistency of the factorization and running. A key insight of our work is that for an N -jet cross-section described by Eq. (1.2), the factorization theorem receives corrections not only in the usual SCET power counting parameter λ, but also corrections due to jet separation beginning at O(1/t 2 ).

Power Corrections to Factorized Jet Shape Distributions
As always, there are power corrections to the factorization theorem which we must ensure are small. One class of power corrections arises from approximating the jet axis of the measured jet with the collinear direction n i , which labels that jet in the SCET Lagrangian. This direction n i is the direction of the parent parton initiating the jet. The jet observable must be such that the difference between the parent parton direction and the jet axis identified by the algorithm makes a subleading correction to the calculated value of the jet observable. In the context of angularity event shapes, such corrections were estimated in [17,29] and found to be negligible for a < 1, and we find the same condition for jet shapes.
In the presence of algorithms, however, there are additional power corrections due to the difference in the soft particles that are included or excluded in a jet by the actual algorithm and in its approximated form in the factorization theorem. We study the effect of this difference on the measurement of jet shapes, and find that for sufficiently large R the power corrections due to the action of the algorithm on soft particles remain small enough not to spoil the factorization for infrared-safe cone and k T -type recombination algorithms. Algorithm-related power corrections to jet momenta were studied more quantitatively in [36], and their estimated R dependence is consistent with our observations. We do not address in this work the issue of power corrections to jet shapes due to hadronization. Event shape distributions are known to receive power corrections of the order 1/(τ a Q), enhanced in the endpoint region but suppressed by large energy. The endpoints of our jet shape distribution near τ a → 0, therefore, will have to be corrected by a nonperturbative shape function. Such functions have been constructed for event shapes in [37,38]. The shift in the first moment of event shape distributions induced by these shape functions was postulated to take a universal form in [39,40] based on the behavior of single soft gluon emission, and the universality was proven to all orders in soft gluon emission at leading order in the SCET power counting in [41,42]. This universality relied on the boost invariance of the soft function describing soft gluon radiation from two backto-back collinear jets. The extent to which such universality may survive for jet shapes with multiple jets in arbitrary directions is an open question that must be addressed in order to construct appropriate soft shape function models to deal adequately with the power corrections to jet shapes from hadronization. Nonperturbative power corrections to jet observables from hadronization and the underlying event in hadron collisions were also studied in [36], and hadronization corrections were found to scale like 1/R. In this work, we focus only on the perturbative calculation and resummation of large logarithms of jet shapes, and leave inclusion of nonperturbative power corrections for future work.

Resummation and Logarithmic Accuracy
Knowing the anomalous dimensions of the hard, jet, and soft functions in the factorization theorem allows us to resum logarithms of ratios of the hard, jet, and soft scales. We take this opportunity to explain the order of accuracy to which we are able to resum these logarithms. For an event shape distribution dσ/dτ (i.e. Eq. (1.2) with two jets and integrated against δ(τ − τ 1 − τ 2 )), the accuracy of logarithmic resummation [43] is typically characterized by counting logs in the exponent ln R(τ ) of the "radiator," where they appear in the form α n s ln m τ with m ≤ n + 1. At leading-logarithmic (LL) accuracy all the terms with m = n + 1 are summed; next-to-leading-logarithmic (NLL) accuracy sums also the m = n terms, and so on. In more traditional methods in QCD, event shapes that have been resummed include NLL resummation of thrust in [43,44], jet masses in [43,45,46,47], jet broadening in [48,49], the C-parameter in [50], and angularities in [17]. Resummation of an event shape distribution using the modern SCET method was first illustrated with the thrust distribution to LL accuracy in [51]. Heavy quark jet mass distributions were resummed in SCET to NLL accuracy as part of a proposed method to extract the top quark mass in [52]. The N 3 LL resummed thrust distribution in SCET was compared to LEP data to extract a value for the strong coupling α s to high precision in [53]. Angularities were resummed to NLL accuracy in SCET in [54] directly in τ a -space instead of in moment space as in [17].
Summation of logarithms in effective field theory is achieved by RG evolution. In the factorized radiator of the thrust distribution Eq. (1.5), one finds that the hard function contains logarithms of µ/Q, the jet functions contain logarithms of µ/(Q √ τ ), and the soft function contains logarithms of µ/(Qτ ). Thus, evaluating these functions respectively at the hard scale µ H = Q, jet scale µ J = Q √ τ and soft scale µ S = Qτ eliminates large logarithms in each function. They can then be RG-evolved to the common factorization scale µ after calculating their anomalous dimensions. The solutions of the RG evolution equations are of the form that logarithms of τ are resummed to all orders in α s to a logarithmic accuracy determined by the order in α s to which the anomalous dimensions and hard/jet/soft functions are known. This underlying hierarchy of scales is illustrated Fig. 1 [in this case, with only one (measured) jet scale and soft scale and ω = Q] along with a table that lists the order in α s to which various quantities must be known in order to achieve a given N k LL accuracy in the exponent of the radiator Eq. (1.5). The power of the EFT framework is to organize of the logs of τ arising in Eq. (1.5) into those that arise from ratios of the jet to the hard scale and those that arise from ratios of the soft to the hard scale, which then allows RG evolution to resum them. For the multijet shape distribution in Eq. (1.2), the strategy to sum logs is the same, but is complicated by the presence of additional scales. This also makes trickier the classification of logarithmic accuracy into the standard N k LL scheme. Our aim will be to sum as many logs of the jet shapes τ a as possible, while not worrying about any others. For instance, phase space cuts induce logs of R and Λ/ω (where ω is a typical hard jet energy), and the presence of multiple jets induces logs of jet separations n i · n j or ratios of jet energies ω i /ω j . We will not aim to sum these types logs systematically in this paper, only those of τ a (though we sum subsets of the others incidentally). In particular, resumming the phase space logs of R or Λ/ω is complicated by how the phase space cuts act order-by-order in perturbation theory 2 , and the fact that a simple angular cut R is less restrictive than a small jet mass or angularity on how collimated a jet must be. That is, an angular cut allows particles in a jet to be anywhere within an angle R of the jet axis regardless of their energy, while a small jet mass or angularity forces harder particles to be closer to the jet axis. The former allows hard particles to lie along the edges of a jet, and hard scale "unmeasured" jet scale soft scales Figure 1: An illustration of generic scales along with a table of log-accuracy versus perturbative order. A cross section with jets of energy ∼ ω, radius R, and energy Λ outside the jets, with some jets' shapes τ a being measured and others' shapes left unmeasured, induces measured and unmeasured jet scales at µ meas J and µ unmeas J . Dynamics at these scales are described by separate collinear modes in SCET. Soft dynamics occur at several soft scales, µ Λ and µ ΛR induced by the soft out-of-jet energy cut Λ and jet radius R, and µ meas S induced by the measured jet shape τ a . RG evolution in SCET resums logs of ratios of jet scales to the hard scale µ H individually, and logs of the ratio of a "common" soft scale µ S to the hard scale. Remaining logs of ratios of soft scales to one another are not resummed in current formulations of SCET. The accuracy of logarithmic resummation of these ratios of scales is determined by the order to which anomalous dimensions and matching coefficients or matrix elements (i.e. hard/jet/soft functions) are calculated in perturbation theory. In this paper we perform the RG evolution indicated by the arrows to NLL accuracy.
soft radiation from such configurations that escapes the jets can lead to logs of Λ/ω that are not captured in our treatment. These are not issues we solve in this paper, in which we focus on resumming logs of jet shapes τ a . (Some exploration of phase space logarithms in SCET was carried out in [31,32].) completely to a desired N k LL order. The complication is in the soft function. The soft function is sensitive to soft radiation into measured and unmeasured jets and outside of all jets. As we will see by explicit calculation, this induces logs of µ tan 1−a (R/2)/(ωτ a ) from radiation into measured jets, and logs of µ/(2Λ) and µ/(2Λ tan R 2 ) from radiation from unmeasured jets cut off by the energy Λ. In addition, though not illustrated in Fig. 1, there can be logs of multiple jet shapes to one another, τ i a /τ j a . No single choice of a soft scale µ S will minimize all of these logs.
In the present work, we will start with the simple strategy of choosing a single soft scale µ S ∼ ωτ a / tan 1−a (R/2) for a jet whose shape τ a we are measuring and logs of which we aim to resum. We will calculate hard/jet/soft functions and anomalous dimensions corresponding to "NLL" accuracy listed in Fig. 1. By this we do not mean all potentially large logs in Eq. (1.2) are resummed to NLL, but only those logs of ratios of a jet scale to the hard scale or of the (common) soft scale to the hard scale. We do not attempt to sum logs of ratios of soft scales to one another completely to NLL accuracy (which can contain τ a ). In the case that all jets' shapes are measured and are similar to one another, τ i a ∼ τ j a , our resummation of large logs of τ i a would be complete to NLL accuracy. We will nevertheless venture to propose a framework to "refactorize" the soft function into further pieces dependent on only a single soft scale at a time and perform additional RG running between these scales to resum the additional logarithms, and will implement it at the level of the O(α s ) soft functions we calculate. However, one cannot really address mixed logarithms such as log(τ i a /τ j a ) that arise for multiple jets until O(α 2 s ), the first order at which two soft gluons can probe two different physical regions. This lies beyond the scope of the present work. (We note, however, that our implementation of refactorization using the one-loop soft function does already seem to tame logarithmic dependence on Λ in our numerical studies of jet shape distributions.) These issues are related to some types of "non-global" logarithms described by [46,56,57,58] that spoil the simple characterization of NLL accuracy. In [59] these were identified as next-to-leading logs of ΛR 2 /(ω i τ a ) and Λ/Q (when R 1) that appear at O(α 2 s ) in jet shape distributions. These authors organized the radiator for a single jet shape distribution into a "global" and "non-global" part [58,59], In this language, the calculations we undertake in this paper resum logs in the global part to NLL accuracy but not in the non-global part. The first argument in R ng is related to ratios of soft scales illustrated in Fig. 1, and the second argument arises when there are unmeasured jets. In the case that all jets are measured, R ∼ 1, and Λ ∼ ω i τ i a , the non-global logs vanish.
While summing all global and non-global logs to at least NLL accuracy will be important for precision jet phenomenology, what we contribute in this paper are key developments and calculations necessary to resum even global logs of jet shapes for jets defined with algorithms. We also believe the effective theory approach and the idea of refactorizing the soft function will help us understand and resum many types of non-global logarithms.

Detailed Outline of This Work
In this paper, we will formulate and prove a factorization theorem for distributions in the jet shape variables we introduced above, calculate the jet and soft functions appearing in the factorization theorem to O(α s ) in SCET, and use the renormalization group evolution of these functions to sum global logs of τ a to NLL accuracy. We consider N jets (defined with a cone or k T algorithm) produced in an e + e − collision, with M of the jets' shapes (angularities) being measured. The key formal result is our demonstration of Eq. (1.3), the consistency of the anomalous dimensions of hard, jet and soft functions to O(α s ) for any number of total jets, any numbers of quark and gluon jets, any number of these jets whose shapes are measured, and any value of the distance measure R in cone or k Ttype algorithms (as long as t 1). These results lead to accurate predictions for the shape of the τ a distribution near the peak, but not necessarily the endpoints for very small τ a (where hadronization corrections dominate) and very large τ a (where fixed-order NLO QCD corrections take over, which are not yet calculated and not captured by NLL resummation). 3 In Sec. 2 we describe in detail the jet shapes and jet algorithms that we use. We describe features of an "ideal" jet algorithm that would respect exactly the order of operations envisioned in the factorization theorem derived in SCET, and the extent to which cone and recombination algorithms actually in use resemble this idealization.
In Sec. 3, using the tools of SCET, we will derive in detail a factorization theorem for exclusive 3-jet production where we measure the angularity jet shape of one of the jets, and then perform the straightforward extension to N -jet production with M ≤ N measured jets. We will give a review of the necessary technical details of SCET in Sec. 3.1. In the process of justifying the factorization theorem, we identify the new requirements listed above on N -jet final states and jet algorithms that must be satisfied for factorization to hold. In Sec. 3.4 we explore in some detail the power corrections to the factorization theorem due to soft radiation and the action of jet algorithms that cause tension with these requirements, and argue that for sufficiently large R in infrared-safe cone and recombination algorithms, these corrections are sufficiently small.
Next we calculate to O(α s ) the jet and soft functions corresponding to N cone or k T -type jets, with M jets' shapes measured.
In Sec. 4 we calculate the jet functions for measured quark jets, J q ω (τ a ), unmeasured quark jets, J q ω , measured gluon jets, J g ω (τ a ), and unmeasured gluon jets, J g ω , where ω = 2E J is the label momentum of the collinear jet field in each jet function. We find that in collinear sectors for measured jets, the collinear scale (and thus the SCET power counting parameter in that sector λ i ) is given by ω i τ 1/(2−a) a shape τ a with R ∼ λ 0 i , while in unmeasured jet sectors, λ i ∼ tan(R/2). Thus one should understand tan(R/2) to be significantly less than 1 but much larger than any jet shape τ a .
In Sec. 5 we calculate the soft function. To do this, we split the soft function into several contributions from different parts of phase space in order to facilitate the calculation and elucidate its intuitive structure. We find it most convenient to split the soft function into an observable-independent part that arises from soft emission out of the jets, S unmeas , and a part that depends on our choice of angularities as the observable that arises from soft emission into measured jet i, S meas (τ i a ). S unmeas is hence sensitive to the scale Λ while S meas (τ i a ) is sensitive to the scale ω i τ i a . In Sec. 6, having calculated all the jet and soft function contributions to O(α s ), we extract the anomalous dimensions and perform renormalization-group (RG) evolution. We find the hard anomalous dimension from existing results in the literature. The hard, jet, and soft anomalous dimensions have to satisfy the consistency condition Eq. (1.3) in order for the physical cross section to be independent of the arbitrary factorization scale µ. Our calculations reveal that, as long as the jet separation parameter t Eq. (1.4) between all pairs of jets is much larger than 1, the condition is satisfied.
Even after requiring t 1, the satisfaction of the consistency condition is non-trivial. The hard function knows only about the direction of each jet and the jet function knows only the jet size R; the soft function knows about both. Furthermore, it is not sufficient only to include regions of phase space where radiation enters the measured jets. We learn from our results in this Section that it is crucial to include soft radiation outside of all jets with an upper energy cutoff of Λ. Only after including all of these contributions from the various parts of phase space do the jet, hard, and soft anomalous dimensions cancel and we arrive at a consistent factorization theorem.
We conclude Sec. 6 by proposing in Sec. 6.4 a strategy to sum logs due to a hierarchy of scales in the soft function, by "refactorizing" it into multiple pieces, each sensitive to a single scale, as suggested by the discussion surrounding Fig. 1. Our current implementation of this procedure does tame the logarithmic dependence of jet shape distributions on the ratio Λ/ω in our numerical studies, but we leave for further development the resummation of all "non-global" logs of ratios of multiple soft scales that begin at NLL and O(α 2 s ). To help the reader find the results of the calculations in Sec. 4 through Sec. 6 just outlined, Table 1 provides a summary with equation numbers.
In Sec. 7 we compare our resummed perturbative predictions for jet shape distributions to the output of a Monte Carlo event generator. We test both the accuracy of these predictions and assess the extent to which hadronization corrections affect jet shapes. We will illustrate our results in the case of e + e − → 3 jets, with the jets constrained to be in a configuration where each has equal energy and are maximally separated. In both the effective theory and Monte Carlo, we can take the jets to have been produced by an underlying hard process e + e − → qqg. After placing cuts on jets to ensure each parton corresponds to a nearby jet, we measure the angularity jet shape of one of the jets. We compare our resummed theoretical predictions with the Monte Carlo output for quark and gluon jet shapes with various values of a and R. We find that the dependencies on a and R of the shapes of the distribution and the peak value of τ a agree well between the theory and  Table 1: Directory of main results: the fixed-order NLO quark and gluon jet functions for jets whose shapes τ a are measured or not; the fixed-order NLO contributions to the soft functions from parts of phase space where a soft gluon enters a measured jet, S meas (τ a ), or not, S unmeas ; their anomalous dimensions; the contributions the jet and soft functions make to the finite part of the NLL resummed distributions; and the full NLL resummed jet shape distribution itself.
Monte Carlo, with small but noticeable corrections due to hadronization. We can estimate these corrections by comparing output with hadronization turned on or off in Monte Carlo. In Sec. 8, we give our conclusions and outlook. We also collect a number of technical details and results for O(α s ) finite pieces of jet and soft functions in the Appendices.
Our work is, to our knowledge, the first achieving factorization and resummation of a jet observable distribution in an exclusive N -jet final state defined by a non-hemisphere jet algorithm. 4 Having demonstrated the consistency of this factorization for any number of quark and gluon jets, measured and unmeasured jets, and phase space cuts in cone and k Ttype algorithms, and having constructed a framework to resum logarithms of jet shapes in the presence of these phase space cuts, we hope to have provided a starting point for future precision calculations of many jet observables both in e + e − and hadron-hadron collisions. The case of pp collisions will require a number of modifications, including turning two of our outgoing jet functions into incoming "beam functions" introduced in [62]. We leave this generalization for future work.
The reader wishing to follow the general structure of our ideas and logic and understand the basis of the final results of the paper without working through all the technical details may read Secs. 1 and 2, and then skip to Sec. 7. Some short less technical discussion also appears in Sec. 3.4.

Jet Shapes
Event shapes, such as thrust, characterize events based on the distribution of energy in the final state by assigning differing weights to events with differing energy distributions. Events that are two-jet like, with two very collimated back-to-back jets, produce values of the observable at one end of the distribution, while spherical events with a broad energy distribution produce values of the observable at the other end of the distribution. While event shapes can quantify the global geometry of events, they are not sensitive to the detailed structure of jets in the event. Two classes of events may have similar values of an event shape but characteristically different structure in terms of number of jets and the energy distribution within those jets.
Jet shapes, which are event shape-like observables applied to single jets, are an effective tool to measure the structure of individual jets. These observables can be used to not only quantify QCD-like events, but study more complex, non-QCD topologies, as illustrated for light quark vs. top quark and Z jets in [12,60]. Broad jets, with wide-angle energy depositions, and very collimated jets, with a narrow energy profile, take on distinct values for jet shape observables. In this work, we consider the example of the class of jet shapes called angularities, defined in Eq. (1.1) and denoted τ a . Every value of a corresponds to a different jet shape. As a decreases, the angularity weights particles at the periphery of the jet more, and is therefore more sensitive to wide-angle radiation. Simultaneous measurements of the angularity of a jet for different values of a can be an additional probe of the structure of the jet.

Jet Algorithms
A key component of the distribution of jet shapes is the jet algorithm, which builds jets from the final state particles in an event. (We are using the term "particle" generically here to refer to actual individual tracks, to cells/towers in a calorimeter, to partons in a perturbative calculation, and to combinations of these objects within a jet.) Since the underlying jet is not intrinsically well defined, there is no unique jet algorithm and a wide variety of jet algorithms have been proposed and implemented in experiments. The details of each algorithm are motivated by particular properties desired of jets, and different algorithms have different strengths and weaknesses. In this work we will calculate angularity distributions for jets coming from a variety of algorithms. Because we calculate (only) at next-to-leading order, there are at most 2 particles in a jet, and jet algorithms that implement the same phase space cuts at NLO simplify to the same algorithm. At this order the two standard classes of algorithms, cone algorithms and recombination algorithms, each simplify to a generic jet algorithm at NLO. At NLO the cone algorithms place a constraint on the separation between each particle and the jet axis, while the recombination algorithms place a constraint on the separation between the two particles.
Cone algorithms build jets by grouping particles within a fixed geometric shape, the cone, and finding "stable" cones. A cone contains all of the particles within an angle R of the cone axis, and the angular parameter R sets the size of the jet. In found jets (stable cones), the direction of the total three-momentum of particles in the cone equals the cone (jet) axis. Different cone algorithms employ different methods to find stable cones and deal with differently the "split/merge" problem of overlapping stable cones. The SISCone algorithm [63] is a modern implementation of the cone algorithm that finds all stable cones and is free of infrared unsafety issues. In the next-to-leading order calculation we perform, there are at most two particles in a jet, and we only consider configurations where all jets are well-separated. Therefore, it is straightforward to find all stable cones, there are no issues with overlapping stable cones, and the phase space cuts of all cone algorithms are equivalent. This simplifies all standard cone algorithms to a generic cone-type algorithm, in which each particle is constrained to be within an angle R of the jet axis. For a two-particle jet, if we label the particles 1 and 2 and the jet axis n, then the cone-like constraints for the two particles to be in a jet are cone jet: θ 1n < R and θ 2n < R . (2.1) This defines our cone-type algorithm. Recombination algorithms build jets by recursively merging pairs of particles. Two distance metrics, defined by the algorithm, determine when particles are merged and when jets are formed. A pairwise metric ρ pair (the recombination metric) defines a distance between pairs of particles, and a single particle metric ρ jet (the beam, or promotion, metric) defines a distance for each single particle. Using these metrics, a recombination algorithm builds jets with the following procedure: 5 0. Begin with a list L of particles.
1. Find the smallest distance for all pairs of particles (using ρ pair ) and all single particles (using ρ jet ).
2a. If the smallest distance is from a pair, merge those particles together by adding their four momenta. Replace the pair in L with the new particle.
2b. If the smallest distance is from a single particle, promote that particle to a jet and remove it from L.
3. Loop back to step 1 until all particles have been merged into jets.
The k T , Cambridge-Aachen, and anti-k T algorithms are common recombination algorithms, and their distance metrics are part of a general class of recombination algorithms. For e + e − colliders, a class of recombination algorithms can be defined by the parameter α: where α = 1 for k T , 0 for Cambridge-Aachen, and −1 for the anti-k T algorithm. The parameter R sets the maximum angle between two particles for a single recombination. 6 In the multijet configurations we consider the jets are separated by an angle larger than R, so that only the pairwise metric is relevant for describing the phase space constraints for particles in each jet. For a two-particle (NLO) jet, the only phase space constraint imposed by this class of recombination algorithms is that the two particles be separated by an angle less than R: This defines a generic recombination algorithm suitable for our calculation. We will denote this as the k T -type algorithm.
The configurations with two widely separated energetic particles best distinguish conetype jets from k T -type jets at NLO. For instance, the case where the two energetic particles are at opposite edges of a cone jet (at an angle 2R apart) is not a single k T jet. However, it is important to note that these configurations will not be accurately described in this SCET calculation for R λ, as such configurations are power suppressed in our description of jets. Our concern is in accurately describing the configurations with narrow jets (small τ a ), and not the wide angle configurations above.
Because jets are reliable degrees of freedom and provide a useful description of an event when they have large energy, in the description of an event we impose a cut Λ on the minimum energy of jets. An N -jet event, therefore, is one where N jets have energy larger than the cutoff Λ, with any number of jets having energy less than the cutoff. In our calculation, we impose the same constraint: any jet with energy less than Λ is not considered when we count the number of jets in the final state. This imposes phase space cuts: for a gluon radiated outside of all jets in the event, that gluon is required to have energy E g < Λ to maintain the same number of jets in the event. The proper division of phase space in calculating the jet and soft functions is a key part of the discussion below, and careful treatment of the phase space cuts is needed.

Do Jet Algorithms Respect Factorization?
The factorization theorem places specific requirements on the structure of jet algorithms used to describe events. As in Eq. (1.2), the factorization theorem divides the cross section into separately calculable hard, jet, and soft functions. The hard function depends only on the configuration of jets, while the jet and soft functions describe the degrees of freedom in each jet in terms of the observable τ . While the soft function is global, the jet function depends only on the collinear degrees of freedom in a single jet. The limited dependence of the hard and jet functions implies constraints on the jet algorithm.
Because jets are built from the long distance degrees of freedom arising from evolution of energetic partons to lower energies, the configuration of jets in an event depends on dynamics across all energy scales. This naively breaks factorization in SCET, since the configuration of jets in the hard function would depend on dynamics at low energy in the soft function. However, we can describe a jet algorithm that respects factorization, and in Sec. 3.4 we will parameterize the power corrections that arise from various algorithms.
The primary constraint on the jet algorithm in order to satisfy the factorization theorem is that the phase space cuts on the collinear particles in the jet are determined only by the collinear degrees of freedom. This essentially ensures that the jet functions are independent of dynamics in the soft function. Correspondingly the soft function can only know about the jet directions and their color representations. The direction of the jet is naturally set by the collinear particles, as soft particles have energy parametrically lower than the collinear ones and change the jet direction by a power suppressed amount. The further restriction that the phase space cuts on the collinear degrees of freedom are independent of the soft degrees of freedom places a strong constraint on the action of the jet algorithm. Cone algorithms already implement this constraint: the jet boundary (the cone) is determined by the location of the jet axis, which is the direction of the sum of all collinear particles up to a power correction. Recombination algorithms, however, are constrained by the factorization theorem to operate in a specific way: all collinear particles must be recombined before soft particles. As discussed in Sec. 3.4, commonly used algorithms obey this constraint up to power corrections in the observable for measured jets. Of particular note is the anti-kT algorithm, which exhibits behavior very close to what is required by the factorization theorem (similarly to the way cone algorithms behave).

Factorization of Jet Shape Distributions in e + e − to N Jets
In this Section we formulate a factorization theorem for jet shape distributions in e + e − annihilation to N jets. All the formal aspects we need to describe an N -jet cross section appear already in the 3-jet cross section, so we will give explicit details only for that case. We will use the framework of Soft-Collinear Effective Theory (SCET), developed in [21,22,23,24], to formulate the factorization theorem. We begin with a basic review of the relevant aspects of the effective theory.

Overview of SCET
SCET is the effective field theory for QCD with all degrees of freedom integrated out, other than those traveling with large energy but small virtuality along a light-like trajectory n, and those with small, or soft, momenta in all components. A particularly useful set of coordinates is light-cone coordinates, which uses light-like directions n andn, with n 2 =n 2 = 0 and n ·n = 2. In Minkowski coordinates, we take n = (1, 0, 0, 1) and n = (1, 0, 0, −1), corresponding to collinear particles moving in the +z direction. A generic four-vector p µ can be decomposed into components In terms of these components, p = (n · p, n · p, p ⊥ ), collinear and soft momenta scale with some small parameter λ as where E is a large energy scale, for example, the center-of-mass energy in an e + e − collision. λ is then the ratio of the typical transverse momentum of the constituents of the jet to the total jet energy. Quark and gluon fields in QCD are divided into collinear and soft effective theory fields with these respective momentum scalings: We factor out a phase containing the largest components of the collinear momentum from the fields q n , A n . Defining the "label" momentump µ n =n ·p n n µ 2 +p µ ⊥ , wheren ·p n contains the O(1) part of the large light-cone component of the collinear momentum p n , andp ⊥ the O(λ) transverse component, we can partition the collinear fields q n , A n into their labeled components, The sums are over a discrete set of O(1, λ) label momenta into which momentum space is partitioned. The binp = 0 is omitted to avoid double-counting the soft mode in Eq. (3.3) [65]. The labeled fields q n,p , A n,p now have spacetime fluctuations in x which are conjugate to "residual" momenta k of the order Eλ 2 , describing remaining fluctuations within each labeled momentum partition [23,65]. It will be convenient to define label operators P µ = n · Pn µ /2 + P µ ⊥ which pick out just the label components of momentum of a collinear field: Ordinary derivatives ∂ µ acting on effective theory fields φ n,p (x) are of order Eλ 2 . The final step to construct the effective theory fields is to isolate the two large components of the Dirac spinor q n,p for a fermion with lightlike momentum along n. The large components ξ n,p and the small Ξ n,p can be separated by the projections ξ n,p = n /n / 4 q n,p , Ξ n,p =n /n / 4 q n,p , (3.6) and we have q n,p = ξ n,p + Ξ n,p . One can show, substituting these definitions into the QCD Lagrangian, that the fields Ξ n,p have an effective mass of order E and can be integrated out of the theory. The effective theory Lagrangian at leading order in λ is [22,23,24] where the collinear quark Lagrangian L ξ is where W n is the Wilson line of collinear gluons, the collinear gluon Lagrangian L An is where c n is the collinear ghost field and α the gauge-fixing parameter; and the soft Lagrangian L s is which is identical to the form of the full QCD Lagrangian (the usual gauge-fixing terms are implicit). In the collinear Lagrangians, we have defined several covariant derivative operators, In addition, there is an implicit sum over the label momenta of each collinear field and the requirement that the total label momentum of each term in the Lagrangian be zero. Note the soft quarks do not couple to collinear particles at leading order in λ. Meanwhile, the coupling of the soft gluon field to a collinear field is in the component n·A s only, according to Eqs. (3.8) and (3.10), which makes possible the decoupling of such interactions through a field redefinition of the soft gluon field given in [24]. We will utilize this soft-collinear decoupling to simplify the proof of factorization below.
The SCET Lagrangian Eq. (3.7) may be extended to include collinear particles in more than one direction [25]. One adds multiple copies of the collinear quark and gluon Lagrangians Eqs. (3.8) and (3.10) together. The collinear fields in each direction n i constitute their own independent set of quark and gluon fields, and are governed in principle by different expansion parameters λ associated with the transverse momentum of each jet, set either by the angular cut R in the jet algorithm or by the measured value of the jet shape τ a . Each collinear sector may be paired with its own associated soft field A s with momentum of order Eλ 2 with the appropriate λ. For the purposes of keeping the notation tractable while proving the factorization theorem in this section, we will for simplicity take all λ's to be the same, with a single soft gluon field A s coupling to collinear modes in all sectors. In Sec. 6.4, we will discuss how to "refactorize" the soft function further into separate soft functions each depending only on one of the various possible soft scales.
The effective theory containing N collinear sectors and the soft sector is appropriate to describe QCD processes with strongly-interacting particles collimated in N well-separated directions. Thus, in addition to the power counting in the small parameter λ within each sector, guaranteeing that the particles in each direction are well collimated, we will find in calculating an N -jet cross section the need for another parameter that guarantees that the different directions n i are well separated. This latter condition requires t 1, where t is defined in Eq. (1.4). 7

Jet Shape Distribution in e + e − → 3 Jets
Consider a 3-jet cross section differential in the jet 3-momenta P 1,2,3 , where we measure the shape τ 1 a of one of the jets, which we will call jet 1. The full theory cross section for e + e − → γ * → 3 jets at center-of-mass energy Q is where the J (X) is the jet algorithm acting on final state X, and N (J (X)) is the number of jets identified by the algorithm [30]. P(jet j) is the 3-momentum of jet j, and is also a function of the output of the jet algorithm J (X). L µ is the leptonic part of the amplitude for e + e − → γ * → qqg. The current j µ is summing over colors a and flavors f . When the three jet directions are well separated, we can match the QCD current j µ (x) onto a basis of three-jet operators in SCET [34,66]. We build these operators from quark jet fields χ n , related to collinear quark fields ξ n by χ n = W † n ξ n , where W n is given by Eq. (3.9), and a gluon jet field B ⊥ n related to gluons A n by The matching relation is (3.16) with sums over Dirac spinor indices α, β and Lorentz index ν, and over label directions n 1,2,3 and label momentap 1,2,3 . Sums over colors and flavors are implied. We have chosen to produce a quark in direction n 1 , antiquark in n 2 , and gluon in n 3 . The matching coefficients C µ αβν are found by equating QCD matrix elements of j µ to SCET matrix elements of the right-hand side of Eq. (3.16). These coefficients have been found at tree level in [66]. The number of independent Dirac and Lorentz structures that can actually appear with nonzero coefficients is considerably smaller than suggested by Eq. (3.16) due to symmetries. We will keep the form of these coefficients general to make extension to N jets transparent, which would require the introduction of a basis of N jet fields in Eq. (3.16), with specified the final state. We could move away from the large-t limit and account for corrections to it by using a basis of operators with arbitrary numbers of jets and properly accounting for the regions of overlap between an N jet operator and (N ± 1)-jet operators. This is outside the scope of the present work, where we limit ourselves to kinematics well described by an N -jet operator, and thus, limit ourselves to the large-t limit.
numbers of quark, antiquark, and gluon fields. We will not write the details for an N -jet cross section here, but the procedures are obvious extensions of all the steps in factorizing the 3-jet cross section below.
As a final step before factorization, we redefine the collinear fields to decouple collinearsoft interactions in the Lagrangian [24]: where Y n is a Wilson line of soft gluons along the light-cone direction n, with A s in the fundamental representation. 8 Y n is similar but in the adjoint representation.
The new fields χ n do not have interactions with soft fields in the SCET Lagrangian at leading order in λ. Henceforth, we use only these redefined fields, but for simplicity drop the (0) superscripts.
The cross section in SCET can now be written, To proceed to factorize this cross section, it is convenient to rewrite the remaining delta functions that depend on the final state X in terms of operators acting on X. Those quantities depending on the jet algorithm J can be rewritten in terms of an operator containing the momentum flow operator, where T µν is the energy-momentum tensor, evaluated at time t and position Rn. The operator E µ (n) measures the flow of four-momentum P µ in the direction n (cf. [29,68,69]), and the jet algorithm J can be written as an operatorĴ acting on the momentum flow in all directions [30]. Correspondingly we can define an operator for the 3-momentum of the jet,P(J j (Ĵ )). In addition, the event shape τ a (jet 1) can also be expressed as an operator τ a (J 1 (Ĵ )), built from the momentum flow operator, acting on the state |X (cf. [29]): The operator is constructed to count only particles actually entering the jet in direction n 1 determined by the action of the jet algorithm (for simplicity we will suppress the argument J 1 (Ĵ ) ofτ a in the following, but add a superscript for the jet number). Using these operators, we can eliminate the X dependence in the delta functions in Eq. (3.19) and perform the sum over states X, obtaining The matrix element can be calculated as the sum over cuts of time-ordered Feynman graphs, with the delta function operators inside the matrix element enforcing phase space restrictions from the jet algorithm and jet shape measurement on the final state created by the cut. The operatorsτ a andĴ depend linearly on the energy-momentum tensor, which itself splits linearly in SCET into separate collinear and soft pieces, which will aid us to factorize the full matrix element in Eq. (3.22) into separate collinear and soft matrix elements. To achieve this factorization, however, we must make some more approximations: 1. The contribution of soft particles and residual collinear momenta to the momentum P(jet j) of each jet can be neglected, and the jet momentum is just given by the label momentump j of the collinear state |X j . Thus the energy and jet axis of each jet is approximated to be that of the parent collinear parton initiating the jet. In particular, the squared mass of the jet is order λ 2 compared to its energy. So in this approximation we take the jet energy to be equal to the magnitude of its 3momentum. On the other hand, we keep the leading non-zero contribution to the angularity even though it is also of order λ 2 . These approximations also require that we treat the energy of any particles outside all of the jets, and thus the cutoff Λ, as a soft or residual energy.
2. The Kronecker delta restricting the total number of jets to 3 can be factored into three separate Kronecker deltas restricting the number of jets in each collinear direction n i to 1, and one Kronecker delta restricting the soft particles not to create an additional jet. This approximation requires the separation between jets to be much larger than the size of any individual jet so that different jets do not overlap. Factoring the restriction on the number of jets in this way is one reason that the parameter t ij in Eq. (1.4) is required to be large.
We describe to what extent the algorithms we consider actually satisfy these two approximations in Sec. 3.4. For now we assume these approximations and facts hold, which allows us to factor the cross section Eq. (3.19), We have rewritten the cross section to be differential in E i (the magnitude of P i ) and Ω i (the direction of P i ). In the sum over label directions, n i can be chosen to align with P i such thatp ⊥ i = 0. In Eq. (3.24) we have written the label momentum as ω i ≡n i ·p i . In Eq. (1.1) we approximate the jet axis by this n i and the jet energy byn i ·p i /2, so that they do not depend on soft momenta at all. The operatorsτ n 1 a andτ s a are defined to count only particles inside the measured jet identified by the algorithm.
In the soft matrix element in Eq. (3.24), we have introduced the soft Wilson line Y n in the antifundamental representation to remove the time-and anti-time-ordering operators T,T in Eq. (3.19) [27], and related Wilson lines Y n in the adjoint representation to those in the fundamental representation by [24] Defining the jet functions by the relations 26c) and the soft function by we can express the cross section Eq. (3.24) as The quark and antiquark jet functions are now for a single flavor, and we have summed over N F flavors to obtain the factor in front. The jet functions depend only on the smallest component of momentum n i · k i in each collinear direction. The residual and soft momenta appearing in the exponentials can be reabsorbed into the label momenta by a series of reparameterizations of the label momenta and directions, under which the SCET Lagrangian is invariant [70]. The three-jet operator Eq. (3.16) will receive corrections of order λ 2 (which we can drop) under the reparameterizations we perform below, or can be constructed from the outset to be explicitly reparameterization invariant (RPI) [66]. First, collect the residual and soft momenta together by defining K = k 1 + k 2 + k 3 + r. We can decompose x in n 1 light-cone coordinates, so Performing a type-A transformation (in the language of [70]) on the label momentum and a type-IB transformation on the vector n 1 itself, shifts the label momentum on the jet function 1 by ω 1 n 1 /2 → (ω 1 +n 1 · K)n 1 /2 + K ⊥1 . The summation variables n 1 , ω 1 can then be shifted to eliminaten 1 · K and K ⊥1 from the exponentials entirely. We drop all corrections suppressed by λ 2 due to these shifts. It remains to absorb the n 1 · K component of residual and soft momentum, appearing in the exponential factor e −in 1 ·Kn 1 ·x/2 . This cannot be achieved by RPI transformations in the n 1 sector since this momentum is purely residual-there is no label momentum in this direction. However, in a multijet cross section,n 1 can be written as a linear combination of n 2 , n 3 , and, say, n ⊥2 (a unit vector transverse to n 2 ,n 2 ), so that RPI transformations on ω 2 , ω 3 and n 2 similar to those above can absorb n 1 · K into the label momenta also. Then, the x-dependent residual and soft exponentials all disappear, and we can combine the nine superfluousn i · k i and k ⊥i integrals with the nine discrete sums over label directions and momenta to give continuous integrals over total momenta. Performing these with the remaining energy and direction delta functions in Eq. (3.28) and the x integral with the remaining exponential gives the momentum conservation delta function The resulting cross section Eq. (3.28) takes the form where we used that the matching coefficients C µ αβν (n i ;p i ) are such that, by construction, the right-hand side at tree-level is simply the Born cross section (denoted by σ (0) ) for e + e − → qqg times δ(τ 1 a ). The hard function H = 1 + O(α s ) is determined by calculating the matching coefficients C order-by-order in perturbation theory.
The above may be easily modified to describe the antiquark or gluon jet angularities, by moving the appropriate delta function δ(τ i a − τ a (jet i)) into the J 2 or J 3 jet functions. In addition, we may choose from among various jet algorithms. The choice determines what Θ-function restrictions must be inserted into the final state phase space integrations created by cutting the jet and soft diagrams to determine which particles make it into the jet. As long as the algorithm is such that the approximations enumerated above hold, it will not violate factorization of the jet shape cross section. We will discuss factorization and its potential breakdown in the context of particular jet algorithms in more detail in Sec. 3.4.
Another check of the validity of the factorization theorem is that the factorized jet and soft functions be separately IR safe, which is a stronger condition than the full cross section being IR safe. If the observable [54,71] or algorithm [32] too strongly weights final states with narrow jets whose invariant masses are the same as the virtuality of soft particles, then the jet and soft functions for the observable in standard SCET with dimensional regularization become IR divergent. When this occurs it does not necessarily mean that factorization is not possible; but at least not in the standard form derived from the version of SCET we utilized above. It could be, for example, that a scheme to further separate modes by defining the theory with an explicit cutoff [32] or by factorizing modes by rapidity instead of virtuality [65] can restore a version of the factorization theorem. We leave an explicit study of which algorithms and observables give IR safe jet and soft functions in SCET in dimensional regularization for a separate publication. However, we note here that the algorithms and observables (τ a for a < 1) that we consider in this paper, at least at NLO, do give rise to IR safe jet and soft functions.

Jet Shapes in e + e − → N jets
To generalize the result Eq. (3.32) to an arbitary number N of jets, we simply add more quark and gluon jet fields to the operator matching in Eq. (3.16), resulting in the corresponding number of additional quark and gluon jet functions in Eq. (3.32), along with a hard coefficient and a soft function for N jet directions. Consider an event with 2N q quark and antiquark jets and N g gluon jets, where 2N q + N g = N . Furthermore, we can choose to measure the angularity shape for any number of these jets. Achieving a factorization theorem that remains consistent for any of these combinations is a nontrivial task and thus a powerful test of the validity of the effective theory.
For an N = 2N q + N g jet event, we generalize the matching of the QCD current in Eq. (3.16) to: with implicit sums over Dirac spinor indices α i , β j , Lorentz indices ν k , (anti-)fundamental color indices a i , b j , and adjoint color indices A k . The N jet cross section differential in M jet shapes, with M < N , factors in SCET into the form where σ (0) is the Born cross section for e + e − to N q quarks, N q antiquarks, and N g gluons; the color indices on the hard and soft functions H and S allow for color mixing; and f i is the flavor of each jet function (quark, antiquark, or gluon). H itself is determined by calculating the matching coefficients C in Eq.

Do Jet Algorithms Induce Large Power Corrections to Factorization?
In this section we explore when power corrections to the factorization theorem above become large, in particular those that are due to the action of jet algorithms. We will argue that power corrections to jet angularities induced by the commonly-used cone and recombination algorithms remain suppressed as long as R is sufficiently large. In particular, we need in general that R satisfies R λ and, for the case of the k T algorithm, we require R λ. 9 These power corrections are associated with assumptions we made about the action of the jet algorithm on final states in deriving Eq. (3.34). In general, the size of these power corrections depends both on the algorithm and the observable.
Power corrections to the p T of a jet arising from perturbative emissions (as well as from hadronization and the underlying event in pp collisions) for various jet algorithms were explored in [36]. These power corrections arise for similar reasons as those we discuss below, namely, perturbative emissions changing which partons get combined into the jet. Ref. [36] finds that such power corrections scale like ln R for small R. This result is consistent with our qualitative discussion below, where we argue that power corrections to jet angularities arising from the jet algorithms we use are minimized when R is sufficiently large. For us, R should be at least O(λ) (and for the case of the k T algorithm, R λ). One set of power corrections that is independent of the choice of algorithm arises from the approximation of setting the jet axis equal to the label direction n. Since this neglects the effects of soft particles, it is valid up to O(λ 2 ) corrections. It was argued in Refs. [17,29,41] for the case of hemisphere jet algorithms that these corrections in turn induce corrections to the angularity τ a of order λ 2(2−a) , which, for τ a ∼ λ 2 , are subleading for a < 1. Essentially the same arguments can be applied to all of the algorithms we consider.
Jet algorithm-dependent power corrections arise from the difference in soft particles included in a jet by a given algorithm and those included in the soft function in the factorization theorem. The algorithms also differ amongst themselves in which soft particles they include in a jet. For observables that scale as O(1), such as the jet energies and 3momenta, the contribution of soft momenta can be neglected since they scale as O(λ 2 ). Clearly then, these observables are not dependent on our choice of jet algorithm and so the assumptions we made about factorization of the algorithm in deriving Eq. (3.34) are trivially satisfied.
However, for observables that scale as O(λ 2 ) such as angularities, soft contributions become important and so the details of the algorithms we consider become relevant. We now estimate how closely the phase space region included in the soft function in the SCET factorization theorem approximates the region of soft particles actually included by the jet algorithm. We will argue that unless R O(λ) for the anti-k T and infrared-safe cone algorithms, and R λ for the k T algorithm, the mismatch in areas is sufficiently large 9 We noted earlier that there may be different λ's for the SCET modes describing different jets. For measured jets, λ 2 ∼ τa, while for unmeasured jets, λ ∼ tan(R/2), and for soft gluons outside jets, λ 2 ∼ Λ/Q. In this section, we mean by λ the expansion parameter associated with measured jets, and ensuring R is much larger than this λ. But if R is too large, the separation parameter t ∝ 1/ tan(R/2) becomes too small. We will consider R ∼ 0.4 to 1 to be safe.
Parent Location (before collinear splitting)

Additional Soft Region of Algorithm
Soft Region Common to Both Algorithm and SCET Factorization Thm.
Daughter Location (after collinear splitting) so as to cause a leading-order power correction to the measured jet angularities. For R larger than these bounds, the corrections are negligible. This miscounting arises due to the fact that factorization requires that collinear particles be combined first, and that the soft function only knows about the parent collinear direction. When algorithms do not obey this ordering, factorization may be violated.
To determine the size of the soft particle phase space region for each jet algorithm that is not included by the factorization theorem, we consider the situation depicted in Fig. 2. A parent collinear particle splits into two daughter collinear particles. In the factorization theorem, since collinear particles are combined first, the region of phase space where soft particles are combined into the jet is a circle of radius R about the parent particle direction. However, in a jet algorithm, soft particles in an additional region outside of this may also be combined into the jet (the hatched regions in Fig. 2). If the area of this region is of the same order as the area included by the factorization theorem, then the power corrections to jet angularities induced by the jet algorithm will be leading-order.
Because soft particles have momenta that are parametrically smaller than collinear particle momenta, we determine the omitted region of soft particle phase space by considering the dominant action of the jet algorithm. The k T algorithm serves as a useful example. The k T metric between a pair of soft particles is O(λ 2 )θ ss , the metric between a soft particle and collinear particle is O(λ 2 )θ cs , and the metric between a pair of collinear particles is O(λ 0 )θ cc . Therefore, collinear-collinear recombinations only occur if the angle θ cc between the collinear particles is smaller than the separation between any soft particle and its nearest neighbor by a factor of O(λ 2 ). Given that the typical angle between collinear particles is O(λ), the dominant action of the jet algorithm is to first merge all soft particles with their nearest neighbors, while collinear-collinear recombinations occur late in the operation of the algorithm. This description will suffice to accurately determine the area of the omitted soft phase space for the k T algorithm. Since collinear particles are combined last, on average soft particles within circles of radius R about the collinear daughter particles are included in the k T algorithm jet, 10 as shown in Fig. 2A. The area that the k T algorithm includes that the soft function does not is represented by the hatched region, which is an area of O(θ ij R). This area must be parametrically smaller than the area included by the soft function (of O(R 2 )) for the associated power corrections to be small. We thus require that R θ ij ∼ λ in the SCET power counting.
The anti-k T algorithm combines particles in a manner that is closer to respecting factorization. It finds the hardest particle first and merges particles at successively larger distances from this particle. For the example of two collinear daughters, it will merge all soft particles with the hardest daughter that are closer than the distance to the softer daughter before merging the two daughters and then merging all soft particles a distance R from the merged daughters (i.e., the parent particle), as shown in Fig. 2B. As the Figure  illustrates, the hatched area of the anti-k T jet tends to be smaller than that of the k T jet. In fact, for R > 3θ ij /2, this region vanishes completely (and this case of having only two collinear daughters is a worst-case scenario). This leads us to expect that, for any number of collinear splittings, for R λ (i.e., not necessarily parametrically larger), power corrections due the action of the anti-k T algorithm vanish.
Cone algorithms such as the SISCone algorithm can also include regions that differ from the lowest-order region at higher orders in perturbation theory. We now argue that an arithmetic bound R λ is sufficient to minimize the power correction from these differences, as for the anti-k T algorithm. This situation arises due to the fact that stable solutions may exist with overlapping cones when collinear splittings are larger than the cone radius, i.e., R < θ ij . In these cases, the amount of radiation that falls into the overlapping region is used to decide whether the cones are split or merged. In either case, the boundary of the resulting jet(s) has roughly the appearance of Fig. 2A and the difference between the region of soft radiation assumed in SCET and that by the actual algorithm is O(1). However, for R > θ ij ∼ λ for the case of a single collinear splitting, all of the collinear radiation lies within a region of size R and there will always be a stable cone that includes this radiation and thus the algorithm and the SCET soft function will be sensitive to soft particles in the same region of phase space.
In summary, we have argued that for all the algorithms we consider (k T , anti-k T , infrared-safe cone), power corrections are negligible for sufficiently large R. While anti-k T and cone allow simply R λ instead of R λ as for k T , we will in fact always consider R λ (for the λ in a measured jet sector) in the remainder of this paper, guaranteeing small power corrections for all these algorithms. (R still determines the scale λ in an unmeasured jet sector.) Our focus will remain on resumming logs of jet shapes such as angularities τ a in the presence of jet algorithms, without worrying about resumming logs of R themselves. 11 10 Soft particles in this region can also be removed from this region by merging with other soft particles outside of the region and vice-versa, but this average area suffices for our discussion.

Jet Functions at O(α s ) for Jet Shapes
In this section, we calculate the quark and gluon jet functions for jet shapes at next-toleading order in perturbation theory. The jet functions can be divided into two categories: those for measured jets, which are fixed to have a specific angularity τ a , and those for unmeasured jets, which are not. We will denote the quark jet function by J q ω , the gluon jet function by J g ω , where ω is the label momentum, and a jet function J q,g (τ a ) with an argument of τ a denotes a measured jet. We will calculate the jet functions for the two classes of jet algorithms, k T -type and cone-type algorithms.
In the course of these calculations, we will demonstrate the crucial role of zero-bin subtractions [65] from collinear jet functions in obtaining the consistent anomalous dimensions and the correct finite parts. In this case zero-bin subtractions are not merely scaleless integrals converting IR to UV divergences, but in fact contribute part (sometimes the most important part) of the correct nonzero result, as was already pointed out by [32,72]. The relation of zero-bin subtractions in SCET to eikonal jet subtractions from soft functions in traditional methods of QCD factorization was explored in [41,73,74]. In addition, we find that the zero-bin subtraction removes the contribution of collinear emissions that escape a jet, leaving only power-suppressed pieces in Λ/ω i .

Phase Space Cuts
To calculate the jet functions for a particular algorithm, Figure 3: A representative diagram for the NLO quark and gluon jet functions. The incoming momentum is l = n 2 ω+n 2 l + and particles in the loop carry momentum q ("particle 1") and l − q ("particle 2").
we must impose phase space restrictions in the matrix element. From the jet function definitions, Eq. (3.26), these cuts take two forms. One kind, imposed by the operator δ N (Ĵ ),1 in Eq. (3.26), is common to every jet function. It is the set of phase space restrictions related to the jet algorithm, and requires exactly one jet to arise from each collinear sector of SCET. The other, imposed by the operator δ(τ a −τ a ), is implemented only on measured jets and restricts the kinematics of the cut final states to produce a fixed value of the jet shape. In this section we describe these phase space cuts in detail.
The typical form of the NLO diagrams in the jet functions is shown in Fig. 3. As shown in the figure, the momentum flowing through the graph has label momentum l − ≡n · l = ω and residual momentum l + ≡ n · l, and the loop momentum is q. We will label "particle 1" as the particle in the loop with momentum q and "particle 2" as the particle in the loop with momentum l − q. For the quark jet, we take particle 1 as the emitted gluon and particle 2 as the quark.
As usual, the total forward scattering matrix element can be written as a sum over all cuts. Cutting through the loops corresponds to the interference of two real emission diagrams, each with two final state particles, whereas cutting through a lone propagator that is connected to a current corresponds to the interference between a tree-level diagram and a virtual diagram, each with a single final state particle. Thus, the phase space restrictions and measurements we impose act differently depending on where the diagrams are cut. In addition, since we will be working in dimensional regularization (with d = 4−2 ), which sets scaleless integrals to zero, the only diagrams that contribute are the cuts through the loops. This means that we only need to focus on the form of phase-space restrictions and angularities in the case of final states with two particles.
The regions of phase space for two particles created by cutting through a loop in the jet function diagrams can be divided into three contributions: 1. Both particles are inside the jet.
In contributions (2) and (3), the jet has only one particle, which is the remaining particle with E > Λ.
It is well known that collinear integrations of jet functions can be allowed to extend over all values of loop momenta so long as a "zero-bin subtraction" is taken from the result to avoid double counting the soft region already accounted for in the soft function [65]. We will demonstrate that contributions (2) and (3) are power suppressed by O(Λ/ω), which scales as λ 2 , after the zero-bin subtraction.
The phase space cuts that enforce both particles to be in the jet depend on the jet algorithm. There are two classes of jet algorithms that we consider, cone-type algorithms and (inclusive) k T -type algorithms, and all the algorithms in each class yield the same phase space cuts. We label the phase space restrictions as Θ cone and Θ k T , generically Θ alg . For the cone-type algorithms, These Θ functions demand that both particles are within R of the label direction. For the k T -type algorithms, the only restriction is that the relative angle of the particles be less than R: In the second line we took the collinear scaling of q (q + q − ). While this is not strictly needed, it makes the calculations significantly simpler.
For the phase space restrictions of zero-bin subtractions, we take the soft limit of the above restrictions. The zero-bin subtractions are the same for all the algorithms we consider. For the case of particle 1, which has momentum q, the zero-bin phase space cuts are given by The zero bin of particle 2 is given by the replacement q → l − q.
For all the jet algorithms we consider, the zero-bin subtractions of the unmeasured jet functions are scaleless integrals. 12 However, for the measured jet functions, the zero-bin subtractions give nonzero contributions that are needed for the consistency of the effective theory.
In the case of a measured jet, in addition to the phase space restrictions we also demand that the jet contributes to the angularity by an amount τ a with the use of the delta function δ R = δ(τ a −τ a ), which is given in terms of q and l by In the zero-bin subtraction of particle 1, the on-shell conditions can be used to write the corresponding zero-bin δ-function as (and for particle 2 with q → l − q).

Quark Jet Function
The diagrams corresponding to the quark jet function are shown in Fig. 4. The fully inclusive quark jet function is defined as and has been computed to NLO (see, e.g., [75,76]) and to NNLO [77]. Below we compute the quark jet function at NLO with phase space cuts for the jet algorithm for both the measured jet, J q ω (τ a ), and the unmeasured jet, J q ω . As discussed above, we will find that the only nonzero contributions come from cuts through the loop when both cut particles are inside the jet. 12 Note that algorithms do exist that give nonzero zero-bin contributions to unmeasured jet functions [32]. Fig. 4. The sum of these contributions is

The measured quark jet function includes contributions from naive Wilson line graphs (A) and (B) and QCD-like graphs (C) and (D) in
(4.7) The contribution proportional to d − 2 comes from the QCD-like graphs (C) and (D) in Fig. 4. Only the Wilson line graphs have a nonzero zero-bin limit, which comes from taking the scaling limit q ∼ λ 2 of the naive contribution: (4.8) All jet algorithms that we use yield the same zero-bin contribution, since the phase space cuts are the same.
To evaluate these integrals, we can analytically extract the coefficient of δ(τ a ) by integrating over τ a and using the fact that the remainder is a plus distribution, as defined in Eq. (A.2). We find the naive contribution is The only difference between the jet algorithms that we consider resides in the finite distri-butionJ q alg (τ a ), which is a complicated function of τ a that we give in Appendix A. Note that the divergent part of the naive contribution is proportional to δ(τ a ). This is due to the fact that the jet algorithm regulates the distribution for τ a > 0. The divergent plus distributions come entirely from the zero-bin subtraction, which is given by (4.10) Adding the leading-order contribution to all of the NLO graphs and expanding in powers of , adopting the MS scheme, we find the total quark jet function, This agrees with the standard jet function J(k + ) given in [75,76] by setting a = 0 and k + = ωτ a . We have shown the divergent terms explicitly, and collect the finite pieces in J q alg (τ a ), which we give in Eq. (A.14). Note that there is no jet algorithm dependence in the divergent parts of the jet function at this order in perturbation theory.

Gluon Outside Measured Quark Jet
In this section we calculate the contribution to the quark jet function from the region of phase space in which the gluon exits the jet carrying an energy E g < Λ. This cut causes the contribution to be power suppressed by Λ/ω, which scales as λ 2 . However, we elect to evaluate this case explicitly as it provides a clear example of the zero-bin subtraction giving the proper scaling to the total contribution. We only evaluate this contribution for the cone algorithm; the details of the k T algorithm calculation are similar. Note that the contribution when the quark is out of the jet is power suppressed at the level of the Lagrangian given in Sec. 3.1, in which soft quarks do not couple to collinear partons at leading order in λ.
For the cone algorithm, the gluon exits the jet when the angle between the jet axis, n 1 , and the gluon is greater than R. When the gluon is not in the jet, the cone axis is the quark direction, and so it makes no contribution to the angularity. Therefore, this region of phase space contributes only to the δ(τ a ) part of the angularity distribution.
For the naive contributions, requiring the gluon to be outside the jet and have energy less than Λ, we have the integral Note that the theta function requiring q − < 2Λ is more restrictive than q − < ω. Evaluating Eq. (4.12) yields a contribution that scales with Λ only below the leading term in 1/ : The zero-bin subtraction of Eq. 4.12 is Evaluating Eq. (4.14), we find the zero bin will exactly remove the leading term in 1/ : Therefore, the difference is power suppressed only after the zero bin is included. Because other contributions when one particle is outside of the jet are similarly power suppressed, we will drop them in our remaining discussion of the jet functions.

Unmeasured Quark Jet
When the angularity of a jet is not measured, the jet function has no τ a dependence. The naive and zero-bin contributions are the same as Eqs. (4.7) and (4.8) except for the factor of δ R . The zero-bin contribution is alg .

(4.16)
This integral is scaleless and therefore equal to 0 in dimensional regularization. This implies that the NLO part of the quark jet function for an unmeasured jet is just the naive result. We find, making the divergent part explicit, in the MS scheme, where the finite part J q alg is given in Eq. (A.18). 13

Gluon Jet Function
The diagrams needed for the gluon jet function at NLO are shown in Fig. 5. The fully inclusive jet function, defined as (with J g ω (l + ) normalized to 2πδ(l + ) at tree-level) has been calculated to NLO in Feynman gauge in [34,78,79] and was reported to give the same result in both R ξ and light-cone gauges in [35]. Since our phase space restrictions and the observables act differently on cuts through loops and on cuts through external propagators, it is useful to reproduce these results by explicitly cutting the diagrams.
After some algebra, we find that the sum of all cuts through the loops of the naïve collinear graphs gives We also record the zero bin that needs to be subtracted from Eq. (4.19). To leading-power it is given by (4.20) 13 The unmeasured jet function Eq. (4.17) is not simply obtained by integrating the measured jet function Eq. (4.11) over τa. This is due to the different relative scaling of R with the SCET expansion parameter λi in a measured and unmeasured jet sector, as noted earlier. Namely, R ∼ λ 0 i in a measured jet sector (where λ ∼ √ τa) while λ k ∼ tan(R/2) in an unmeasured jet sector. Without inserting any additional constraints, this integral is scaleless and zero in dimensional regularization. Therefore, in the absence of phase-space restrictions, the naïve integral Eq. (4.19) gives the standard (inclusive) gluon jet function in the MS scheme. The measured and unmeasured jet functions are obtained by inserting Θ alg δ R and Θ alg , respectively, into Eqs. (4.19) and (4.20).

Measured Gluon Jet
The naive contribution to the measured gluon jet can be written as where, as for the quark jet function, the finite distributionsJ g alg (τ a ) differ among the algorithms we consider. They are given in Appendix A.
The zero-bin result is .

(4.24)
Subtracting the zero-bin from the naive integral and adding the leading-order contribution, we obtain in MS The finite distribution J g alg (τ a ) is given in Eq. (A.14).

Unmeasured Gluon Jet
As for the quark jet function, for unmeasured jets the naive and zero-bin contributions are given by the measured jet contributions without the δ R constraint. Also, as it was for the quark jet function, the zero-bin contribution to the unmeasured jet function is a scaleless integral. Thus, the final answer is just the naive result, which is given by with the finite part J g alg given in Eq. (A.29) in the Appendix.

Soft Functions at O(α s ) for Jet Shapes
In this section we compute the soft function for a generic N jet event. In Sec. 5.1, we describe the phase space cuts that we impose on soft particles emitted into the final state. We then give an outline of how we disentangle the various contributions to the N -jet soft function in Sec. 5.2 and proceed to calculate these contributions in the remaining subsections.

Phase Space Cuts
Soft particles in the final state must satisfy the phase space cuts required by the operator δ N (Ĵ ),0 in Eq. (3.27), which constrains the soft particles to not create an extra jet. A soft particle is allowed in the final state if it is emitted into one of the jets as defined by the jet algorithm, or if it is not in a jet but has energy less than a cutoff Λ. At NLO, only a single soft gluon can be emitted. Therefore, for both cone-type and k T -type algorithms, the constraint that the soft gluon be in a jet is simply that the angle of the gluon with respect to the jet axis satisfies θ gJ < R. Thus, our requirement on soft gluons is that they obey one of the two following conditions: in jet i: θ gJ i < R out of all jets: E g < Λ and θ gJ i > R for all i .  These conditions can be written in terms of theta functions on the gluon momentum k. We denote the energy restriction for out-of-jet gluons as and we denote the requirement that a gluon be in jet i in terms of the light-cone components k ± about the direction of jet i, n i , as For the case when the soft gluon is in a measured jet, we demand that it contributes an amount τ a to the angularity of a jet with label momentum ω with the use of the delta function

Calculation of contributions to the N -Jet Soft Function
The topology of the various graphs that we need to compute the soft function is shown in Fig. 6. The next-to-leading order part S (1) of the soft function S is the sum of the interference of soft gluon emissions from Wilson lines in directions i and j, S ij , over all lines i and j with i = j (since for i = j, the diagram is proportional to n 2 i = 0), It is conceptually straightforward to see that S ij is just the sum of the following three classifications of the final state cut gluon: • The gluon is in a measured jet and thus contributes to the jet angularity.
• The gluon is outside of all the jets and has energy E g < Λ.
• The gluon is in an unmeasured jet and has any energy.
However, the second contribution is technically difficult to compute due to the complicated shape of the space with all jets cut out of it, like Swiss cheese. A division of phase space leading to a simpler calculation is the following: The gluon is in measured jet k and contributes to the jet's angularity τ k a .
• S k ij : The gluon is in jet k with energy E g > Λ (and does not contribute to τ k a ).
•S k ij : The gluon is in jet k with energy E g < Λ (and does not contribute to τ k a ).
• S incl ij : The gluon is anywhere with E g < Λ (and does not contribute to any angularity).
In terms of these pieces, the NLO soft function with M measured jets and N − M unmeasured jets is given by From the definitions above, it is easy to see that the term in large parentheses on the second line is equivalent to the sum of the last two contributions on the original list above, i.e., the contributions from a gluon not in any jet with E g < Λ and from a gluon in an unmeasured jet with any energy. We can simplify this expression by noting that the contribution from a gluon in jet k with no energy restriction involves a scaleless integral over the energy that vanishes in dimensional regularization and thus Using this relation, the soft function simplifies to where the first term in Eq. (5.8) is a universal contribution that is needed for every N -jet observable, defined as The second term, defined as, depends on our choice of angularities as the observable.
We now proceed to set up the one-loop expressions for the contributions in Eq. (5.8). The integrals involved are highly nontrivial and so in this section we simply quote the result of each integral, referring the reader to Appendix B for details. Most of these integrals are most easily written in terms of the variable t ij , defined in Eq. (1.4) as t ij ≡ tan ψ ij 2 / tan R 2 , where ψ ij is the angle between jet directions i and j. (That is, n i · n j = 1 − cos ψ ij .) In Table 2, we summarize the divergent parts of the soft function.
The Feynman rules for the emission of a soft gluon from fundamental-and adjointrepresentation Wilson lines (corresponding to quark and gluon jets, respectively) have the same kinematic structure. The difference is entirely encoded in the color charge of the Wilson lines which, using the color space formalism of [80,81], we denote as T i for emission from Wilson line i. Thus, we can consider the N -jet soft function without specifying the color representation of the final-state jets.

Inclusive Contribution: S incl ij
The contribution to the soft function from a gluon going in any direction with energy E g < Λ is given by the integral We evaulate this integral in Sec. B.1 of the Appendix and find Note that this calculation is related to the inclusive, timelike soft function that has applications elsewhere (see, e.g., [82,83,84]), generalized for non back-to-back jets: Using Eq. (5.7), the contribution S k ij from a gluon emitted into jet k from lines i and j is given by the integral Much like for the S meas ij contribution, if k = i, j, there is an additional divergence (arising from n k · k → 0) relative to the case k = i, j, and so we evaluate these two cases separately below.
Case 1: k = i or j The integrals for this case are performed in Sec. B.2 of the Appendix, with the result that S j ij is .

(5.15)
Case 2: k = i, j These contributions are at most 1/ divergent since the matrix element does not have the n k · k → 0 singularity. We show in Appendix B.3.2 that the result takes the form where β ij is the angle between the i-k and j-k planes and the finite function F is given in Eq. (B.33) and is O(1/t 2 ).

Soft gluon inside measured jet k: S meas
ij (τ k a ) The contribution of a gluon emitted into jet 1 (the measured jet) from lines i and j is given by the integral The singularity structure of this integral depends on whether or not k = i or j. Thus, we evaluate the case k = i or j and the case k = i, j separately below.
where G is O(1/t 2 ) and is given in Eq. (B.36) and, again, β ij is the angle between the i-k and j-k planes. contribution divergent terms Table 2: Summary of the divergent parts of the soft function at NLO. The first block contains the the observable-independent contributions: soft gluons emitted by jets i and j in any direction with energy E g < Λ in row 1; soft gluons entering jet k with E g > Λ (with k = i or j in the second row and k = i, j in the third). In the last row of the first block, we summed over i and j and took the large-t limit to get the total S unmeas (1) . Similarly, in the second block we give the contributions to the angularities τ k a (with k = i or j in the first row and k = i, j in the second) and summed over i and j and took the large-t limit to get S meas (1) in the third row.

Total N -Jet Soft Function in the large-t Limit
In this section, we add together the necessary ingredients calculated above to obtain the N -jet soft function from Eq. (5.8). The results for the divergent pieces are summarized in Table 2. In Sec. 6 we use Table 2 to show that the consistency of factorization is explicitly violated by terms of order 1/t 2 , and so in this section we give the full soft function (including the finite terms) to O(1/t 2 ).
Using color-conservation ( i T i = 0), we find that the observable-independent part, S unmeas (1) , defined in Eq. (5.9), is given for large t by S unmeas We see that the finite parts of this contribution are sensitive to two scales, 2Λ and 2Λ tan R 2 . For simplicity, in this paper, since we take tan(R/2) ∼ O(1), we will choose only a single scale µ Λ S to attempt to minimize logs in Eq. (5.20), where chosen as the geometric mean of the two. The remaining part of the soft function that is dependent on angularities as our choice of jet observable is the sum over S meas (1) (τ i a ) (defined in Eq. (5.10)) for each jet i, where S meas (1) (τ i a ) is given by The finite part of this contribution is sensitive to the scale µ i S , where which, in principle, differs for each jet i and from the scale µ Λ S in the unmeasured part of the soft function Eq. (5.20). After discussing resummation of large logarithms through RG evolution, we will describe in Sec. 6.4 a framework to "refactorize" the soft function into pieces depending on multiple separated soft scales. This framework will enable us to resum logarithms of all of these potentially disparate scales.

Resummation and Consistency Relations at NLL
The factorized cross section Eq. (3.34) is written in terms of hard, jet, and soft functions evaluated at a factorization scale µ. Since the physical cross section is independent of µ, the anomalous dimensions of these functions are closely related. We derive and verify this relation in Sec. 6.3. In Sec. 6.1 and Sec. 6.2, we work out the form of the renormalizationgroup equations (RGEs) satisfied by the hard, jet, and soft functions, and obtain their solutions so that we can express each function at the scale µ in terms of its value at a scale µ 0 where logarithms in it are minimized. In Sec. 6.4, we present a framework to refactorize the soft function and give the total resummed distribution in Sec. 6.5.

General Form of Renormalization Group Equations and Solutions
We will build solutions for the hard, jet, and soft functions from two forms of RGEs. The first form is for a function which does not depend on the observable τ a and is multiplicatively renormalized, and satisfies the RGE, where the anomalous dimension γ F is found from the Z-factor by the relation and takes the general form, We call Γ F [α] the "cusp part" of the anomalous dimension and γ F [α] the "non-cusp part". They have the perturbative expansions The RGE Eq. (6.2) has the solution where the evolution kernel U F is given by where K F and ω F will be defined below in Eq. (6.15). The second form of RGE is for a function dependent on the jet shape τ a and is renormalized through a convolution, F bare (τ a ) = dτ a Z F (τ a − τ a ; µ)F (τ a , µ) , (6.9) and satisfying the RGE with an anomalous dimension calculated from (6.11) and taking the general form The solution of an RGE of the form Eq. (6.10) has the solution [82,85,86,87,52] F (τ a ; µ) = dτ U F (τ a − τ a ; µ, µ 0 )F (τ a ; µ 0 ) , (6.13) where the evolution kernel U F is given to all orders in α s by the expression where γ E is the Euler constant. In Eqs. (6.8) and (6.14), the exponents ω F and K F are given in terms of the cusp and non-cusp parts of the anomalous dimensions by the expressions In the case of Eq. (6.8) or if Γ F [α] happens to be zero, we define j F to be 1. To achieve NLL accuracy in the evolution kernels U F , we need the cusp part of the anomalous dimension to two loops and the non-cusp part to one loop, in which case the parameters ω F , K F in Eq. (6.15) are given explicitly by Here r = αs(µ) αs(µ 0 ) , and β 0 , β 1 are the one-loop and two-loop coefficients of the beta function, where (with T R = 1/2) The two-loop running coupling α s (µ) at any scale µ in terms of its value at another scale Q is given by The ratio of the one-loop and two-loop coefficients of Γ cusp is [88] Γ 1 (6.21) Γ 1 cusp and β 1 are needed in the expressions of ω F and K F for complete NLL resummation since we formally take α 2 s ln τ a ∼ O(α s ).

Hard Function
The hard function is related to the matching coefficient C N of the N -jet operator in Eq. (3.33). If there are multiple operators with different color structures, C N is a vector of coefficients. The hard function is then a matrix H ab = C † b C a . The hard function is contracted in the cross section Eq. (3.34) with a matrix of soft functions.
The anomalous dimensions of the matching coefficients C a have been calculated in the literature (for example, Table III of Ref. [89]). For an operator with N legs with color charges T 2 i , the anomalous dimension is 22) where T i is a matrix of color charges in the space of operators, and γ i is given for quarks and gluons by The coefficient Γ(α s ) is the cusp anomalous dimension and is given to one-loop order by Γ(α s ) = α s /π. The anomalous dimension of the hard function itself is given by γ H = γ † C N + γ C N , and takes the form of Eq. (6.4), with cusp and non-cusp parts Γ H [α s ] and γ H [α s ] given to one loop in Table 3 , with the result i is the sum of all the Casimirs, and the effective hard scaleω H appearing as the scale ω in the logarithm in Eq. (6.4) is given by the color-weighted average of the jet energies,ω ]. Γ is the cusp anomalous dimension itself, equal to α s /π at one loop. γ i is given for quark and gluon jets in Eq. (6.23). The third column gives the value of j F appearing in Eq. (6.12) or Eq. (6.15). The last column gives the values of ω appearing in the logarithm ln µ 2 /ω 2 multiplying the cusp part of the anomalous dimension in Eqs. (6.4) and (6.12). The scaleω H is the color-weighted averages of all jet energies defined in Eq. (6.25). All rows except for the last are given in the large-t limit and in the last row we give the remaining terms that are present for arbitrary t. This last row explictly violates consistency at O(1/t 2 ). The first group of rows are needed for measured jets and the second group for unmeasured jets. In the large-t limit, for any number of measured and unmeasured jets, the consistency relation Eq. (6.33) is satisfied.

Jet Functions
There are two forms of jet functions, those for measured and those for unmeasured jets. Unmeasured jet functions J q,g ω (µ) satisfy multiplicative RGEs, with anomalous dimensions given by the infinite parts of Eqs. (4.17) and (4.26), which falls into the general form Eq. (6.4), with cusp and non-cusp parts of the anomalous dimension given in Table 3, and the scale ω in Eq. (6.4) being ω i tan R 2 . The part γ i is given by Eq. (6.23).
Measured jet functions satisfy RGEs of the form Eq. (6.10), with anomalous dimensions given by the infinite parts of Eqs. (4.11) and (4.25), which takes the form Eq. (6.12) with cusp and non-cusp parts of the anomalous dimension split up as in Table 3, and the scale ω in Eq. (6.12) being ω i .

Soft Function
The total N -jet soft function is given by Eq. (5.20) for unmeasured jets added to the sum over measured jets of Eq. (5.22). This soft function depends on the M jet shapes τ 1 a , . . . , τ M a , and satisfies the RGE From the infinite parts of the soft function given in Table 2, we find the anomalous dimension γ S (τ 1 , . . . , τ M ; µ) decomposes, as required by the consistency condition Eq. (6.33) given below, into a sum of terms, where the pieces γ unmeas S (µ) and γ meas S (τ k ; µ) are given in terms of their cusp and non-cusp parts in Table 3, with the result which takes the form of Eq. (6.4) with no cusp part, and which takes the form of Eq. (6.12) with no non-cusp part, and the scale ω in Eq. (6.12) being ω k / tan 1−a R 2 . The solution of the RGE Eq. (6.28) is given by where U unmeas S is given by the form of Eq. (6.8) and U k S (τ k a ) by the form of Eq. (6.14). The solution Eq. (6.32) is appropriate if all the scales appearing in the soft function are similar, and thus all large logarithms in the finite part can be minimized at a single scale µ 0 . As we noted in Sec. 5.3, however, the potentially disparate scales ω i τ i a , set by the jet shapes of the measured jets, and Λ, set by the cutoff on particles outside jets, appear together in the soft function, and logarithms of ratios of these scales may be large. In this case, the soft function should be "refactorized" into pieces depending only on one of these scales at a time. We describe a framework for doing so below in Sec. 6.4.
But first, we verify the consistency of the anomalous dimensions for the hard, jet, and soft functions to the order we have calculated them.

Consistency Relation among Anomalous Dimensions
We summarize the anomalous dimensions of the hard, jet, and soft functions in Table 3. We separate contributions to the jet and soft anomalous dimensions that arise from measured jets, from unmeasured jets, and those that are universally present. In all rows except the last row, we take the large-t limit and give the additional terms that arise (from the soft function) for arbitrary t.
Consistency of the effective theory requires that the anomalous dimensions satisfy From the results tabulated in Table 3, up to corrections of O(1/t 2 ), we see that, remarkably, this relation is indeed satisfied! This is highly nontrivial, as jet and soft anomalous dimensions depend on the jet radius R and the jet shape τ a , while the hard function does not; in addition, the hard and soft functions know the directions n i of all N jets, while the jet functions do not. These dependencies cancel precisely between the R-dependent pieces of unmeasured jet contributions to the jet and soft functions, between τ a -dependent pieces of the measured jet contributions, and between n i · n j -dependent pieces of the hard and soft functions. The sum of all jet and soft anomalous dimensions then precisely matches the hard anomalous dimensions, satisfying Eq. (6.33). Making the satisfaction of consistency even more nontrivial, individual contributions to the infinite part of the soft function, and therefore its anomalous dimension, given by Table 2 depend on the energy cut parameter Λ as well. However, these terms cancel in the sum over the contributions S incl ij and S i ij in the first two rows of Table 2. The double poles of the form 1 ln Λ arise from regions of phase space where a soft gluon can become both collinear to a jet direction (giving a 1/ ) and soft (giving a ln Λ). These regions exist in the integral over all directions giving S incl ij but are subtracted back out in the contributions S i ij .
In the surviving "Swiss cheese" region the regions giving these double poles are cut out.
The O(1/t 2 ) terms that violate consistency arise whenever there are unmeasured jets. While this limit is not required for the contribution of measured jets to the anomalous dimension to satisfy the consistency condition Eq. (6.33), the finite parts of measured jet contributions to the soft function contain large logarithms of ω/Λ that can not be minimized with a scale choice but are suppressed by O(1/t 2 ) (cf. Eq. (B.37) of Appendix B). This is the manifestation of the fact that jets need to be well-separated, as explained in Sec. 3. For the remainder of the paper, we work strictly in the large-t limit.
It may seem mysterious that the calculations of the hard, jet, and soft functions themselves and requiring their consistency lead to the condition of a large separation parameter t. Although we already specified qualitatively in the proof of factorization the requirement of well-separated jets, it may not be immediately apparent where it is implemented in the actual calculations. It enters in the definition of the collinear jet functions. In the large-t limit, the N jets are infinitely separated from one another according to the measure given by Eq. (1.4). And indeed, when N -jet operators are constructed in SCET, each collinear jet field contains a Wilson line W n , defined below in Eq. (3.9), of collinear gluons in the direction n that were emitted from the back-to-back directionn, for which the separation measure t → ∞. (This is similar to QCD factorization proofs of hard scattering cross sections, e.g. in [17], in which this directionn is chosen to be along some arbitrary path ξ that is separated by an O(1) amount from the jet direction n.) Furthermore, the n icollinear jet function knows only its own color representation, and not those of the other jets. Meanwhile, the hard and soft functions we calculate "know" about all N jets and their precise directions and color representations. Therefore it is no surprise that, when we actually calculate the anomalous dimensions of these functions, we find that they are consistent with one another only in the limit that the separation parameter t → ∞.

Refactorization of the Soft Function
Our results for the soft function in Sec. 5.3 make clear that in general there are multiple scales that appear in the soft function: the µ 1 S , . . . , µ M S associated with the M measured jets and the scale µ Λ S associated with the out-of-jet cutoff Λ (see Eq. (5.21)). When these scales are all comparable, we can RG evolve the soft function from the single scale µ S . However, when any of them differ considerably from the others, we need to "refactorize" the soft function into multiple contributions, each of which is sensitive to a single energy scale. For illustration, take the scales µ i S to be such that µ 1 S µ 2 S · · · µ M S as in Fig. 7. We also take µ l−1 S µ Λ S µ l S for our discussion, but the result is independent of whether µ Λ S lies in the range µ 1 S < µ Λ S < µ M S or not. We can express the soft function appearing in Eq. (3.34) as where the operator τ i a returns the contribution to τ a of final-state soft particles entering jet i, andΛ returns the energy of final-state soft particles outside of all N jets. We have kept the dependence on the scales µ i S and on Λ implicit on the left-hand side.
There are N Wilson lines appearing in the operator O S , corresponding to the M measured jets and N − M unmeasured jets. The scales associated with soft gluons entering the M measured jets whose shapes are measured to be τ 1 , . . . , τ M are µ 1 S , . . . , µ M S , given by Eq. (5.23). The scale associated with soft gluons outside of measured jets is µ Λ S given by Eq. (5.21). We have illustrated the ladder of these scales in Fig. 7. Each of these soft scales can be associated with different soft fields A i s whose momenta scale as λ 2 i ω i where λ i is associated with the typical transverse momentum λ i ω i of the collinear modes for the ith jet. For measured jets, λ i is determined by τ i a , while for unmeasured jets λ i ∼ tan(R/2). For soft gluons outside jets, however, the soft momentum is set by the cutoff scale Λ, which is why µ Λ S appears in the ladder of Fig. 7.
At a scale µ larger than all µ i S and µ Λ S , the soft function is calculated as we presented in Sec. 5. In particular, we take the quantities τ i a and Λ to be nonzero and finite. At a scale µ below µ 1 S , we integrate out soft gluons of virtuality µ 1 S and match onto a theory with soft gluons of virtuality µ 2 S . The scale µ 1 S associated with τ 1 a is taken to infinity, and phase space integrals for soft gluons entering the measured jet 1 become zero (see, e.g., Eq. (B.17)). Therefore, the matching coefficient from the theory above µ 1 S to the theory below is just the measured jet 1 contribution S meas (τ 1 a ) to the soft function given by Eq. (5.22). The same occurs when matching from the theory above each scale µ i S for soft gluons entering measured jet i to the scale below µ i S , giving a matching coefficient S meas (τ i a ). When going through the scale µ Λ S , in the theory above this scale, the calculation of unmeasured contributions to the soft function gives the result Eq. (5.20), by treating Λ as a nonzero, finite cutoff. In the theory below µ Λ S , we take Λ to infinity, making all phase space integrals originally cutoff by Λ to be scaleless and thus zero. So the matching coefficient between the theory above and below µ Λ S is just S unmeas . After performing the above matchings all the way down to the lowest soft scale in Fig. 7, we find that the original soft function S(τ 1 a , . . . , τ M a ; µ) can be expressed to all orders as where to next-to-leading order S meas and S unmeas are given by is given by Eq. (5.22). Note that no operators restricting the jet shape or the phase space appear in the final matrix element of soft fields living at the lowest scale on the ladder in Fig. 7. If all the scales on the ladder are at a perturbative scale, we can now just use O † S O S = 1 to eliminate the final matrix element. If any scale is nonperturbative, we should have stopped the matching procedure before that scale, and defined the surviving soft matrix element still containing additional delta function operators as a nonperturbative shape function.
Since the factors S unmeas (µ) and S meas (τ i a , µ) are now matching coefficients between two theories above and below the respective scales µ Λ S and µ i S , we can run each of the individual factors separately from their natural scale, instead of from a single soft scale µ 0 as in Eq. (6.32). The result for the RG-evolved soft function is then Eq. (6.36) where each factor at NLO is given by the solution of its RGE, These solutions allow us now to resum logarithms of all of the scales appearing in the ladder in Fig. 7 when these scales are widely disparate. However, the result we obtained in Eq. (6.28) when we took all scales to be of the same order and had a single soft scale has the form Eq. (6.39) at NLL accuracy. We will use equation Eq. (6.39) in all cases to interpolate between these two extremes.

Total Resummed Distribution
Collecting together the above results for the running of hard, jet, and soft functions in the factorized cross section Eq. (3.34), we obtain the RG-improved N -jet cross section differential in M jet shapes, Using results from Appendix C, we obtain the functions f i J,S generated by the finite pieces of the measured jet and soft functions, where d i J is the part of the unmeasured jet function containing no large logarithms (given in Eqs. (A.19) and (A.30) in the Appendix).
The finite parts of the measured and unmeasured jet and soft functions given in Eqs. (6.42) and (6.43) show that to minimize large logarithms in the O(α s ) finite parts in the resummed distribution Eq. (6.40), we should choose initial scales for the running to be These choices eliminate all large logarithms in the O(α s ) hard, jet, and soft functions. They still leave logs of tan R 2 and n i · n j in the unmeasured part of the soft function, and logs of tan R 2 in the measured jet function, but we already take R numerically of O(1) 14 to minimize power corrections from our implementation of the jet algorithm as discussed in Sec. 3.4, and n i · n j ≈ 1 since the jet separation parameter t ij is large compared to 1. All logs of R, Λ, and τ i a are eliminated in the unmeasured jet function and measured part of the soft function. 14 We still consider tan(R/2) to be of order λ k in the collinear sectors describing unmeasured jets, as implied by Eq. (6.44). This means λ k is effectively much larger than the parameter λi in a measured jet sector. In fact, note that Eq. (6.44) tells us that tan R 2 must be parametrically larger than (τ i a )

Plots of Distributions and Comparisons to Monte Carlo
Having resummed the jet shape distributions in τ a to NLL accuracy, in this section we plot the distributions for various values of a and R, compare to Monte Carlo simulated events, and perform scale variation on the resummed distribution. We use the process e + e − → 3 jets to study our predictions of jet shapes, where the jets arise from partons in the "Mercedes-Benz" configuration, with each parton having equal energy. In these configurations, the partons lie in a plane and are equally separated with a pairwise angle of 2π/3. This allows us to study event shape distributions of well-separated jets where t is reasonably large. We choose three values of R to study, R = 1.0, 0.7, and 0.4. With these values of R, the 1/t 2 suppression factor for corrections to the large-t limit are 0.10, 0.044, and 0.014 respectively. We will measure the angularity of only one of the three jets; the other two jets will be unmeasured. In general, the T i · T j color correlations in the soft and hard functions lead to operator mixing in color space under RG evolution. This implies that the RG kernels U S and U H are matrices in color space and must be studied on a process-by-process basis (see, e.g., [89,90,91,92,93,94]). For the case of N = 2, 3 jets there is only one color-singlet operator and hence no mixing. This can be seen, for example, by noting that all color correlations reduce to the Casimir invariants (C F and C A ) in this case (cf. Appendix D). We have restricted the example process we use in this work to N = 3 jets, avoiding the additional complication of color-correlations that comes with a larger number of jets.
The NLL resummed distribution for one quark or gluon jet shape (jet 1) in a three-jet final state, written as the derivative of the radiator Eq. (1.5), is where the various evolution parameters ω i J,S , Ω, K are all defined in Eqs. (6.15) and (6.41), andf J,S are given by f J,S in Eq. (6.42) with the d J,S terms set to zero (accurate to NLL). The best scale choices Eq. (6.44) for this case are In Eq. (7.1) we have used tree-level initial conditions for the hard, jet, and soft functions at these scales. Eq. (7.1) evolves these functions to the arbitrary scale µ at NLL accuracy.  With these choices, we plot Eq. (7.1) in Fig. 8 for several values of a and R for a quark or gluon jet shape in a three-jet final state in e + e − annihilation at center-of-mass energy Q = 500 GeV. 15 The jets are chosen to be in a Mercedes-Benz configuration, where all jets have equal energies of 150 GeV. We choose the jet energy cutoff Λ to be 15 GeV. We choose the factorization scale to be µ = µ H .  Figure 9: Quark vs. gluon jet shapes with comparison to Monte Carlo. Solid, straight curves represent the resummed jet shape distribution in Eq. (7.1), and jagged curves are histograms from the Monte Carlo, normalized as described in the text. The solid histogram has no hadronization, while the dashed histogram includes the effects of hadronization. The distributions are plotted for a = − 1 2 , 0, 1 2 with quark (blue) and gluon (red) jets compared on the same plot, for jets of size R = 1.0, 0.7, 0.4. Gluon jet shape distributions peak at larger values of τ a than quark jets, indicative of their broader shape. The plots are for jets in e + e − annihilation at center-of-mass energy Q = 500 GeV, with three jets produced with equal energies E J = 150 GeV, angular separation ψ = 2π/3 between all pairs of jets, and minimum threshold Λ = 15 GeV to produce a jet.
We compare the results of a jet algorithm implemented on Monte Carlo simulated events with our NLL resummed predictions for a variety of a and R values in Fig. 9. Because the resummed NLL distribution we choose to study applies to an exclusive process, threejet events in the Mercedes-Benz configuration, we implement cuts on the simulated events to obtain a sample that matches onto this configuration. We use MadGraph/MadEvent v.4.4.21 [95] to generate parton-level e + e − → qqg events at a center-of-mass energy Q = 500 GeV, with cuts imposed to obtain partons in the Mercedes-Benz configuration. We shower and hadronize the events with Pythia v.6.414 [96] using p T -ordered parton showers. The process of hadronization will induce a shift in the angularity distribution, which we do not model in our resummed distribution. Therefore, we produce two samples: one sample with only QCD final-state showering, no initial-state radiation, and no hadronization, and another sample with complete showering, initial-state radiation, and hadronization. The anti-k T jet algorithm is run on the final state particles from Pythia, and we use FastJet [97] to implement the jet algorithm. The anti-k T algorithm is particularly well suited for this comparison, as very few particles at an angle θ > R to the jet axis are included in the jet. With anti-k T , the phase space cut on an individual particle matches well with the phase space cuts in the next-to-leading order calculation.
To select a sample of events to compare to our resummed distributions, we make cuts on the final state jets, requiring each of the three hard, well-separated partons from MadGraph to be associated with a jet. This involves a cut on the jet direction and momentum: We analyze events passing these cuts, and tag each associated jet as coming from a quark or a gluon based on which parton it matches onto. The angularity value for each jet is computed from the constituent particles of the jet, using the matching parton direction as the jet axis. The jet direction only differs from the parton direction by a power correction (see Sec. 3.2). In Fig. 9, we isolate some of the quark and gluon jet shapes in Fig. 8 and compare to Monte Carlo events. The relative normalization between the distribution of Monte Carlo events and the NLL resummed angularity distribution requires some explanation. Both our calculation and the Monte Carlo simulation are most accurate in the region that includes the peak of the distribution and the larger-τ side of the peak, but both are inaccurate as τ → 0 and in the tail region. Therefore, each will differ in the relative normalization between the peak region and the tail region. An accurate prediction of the tail region requires matching onto a calculation at fixed-order in α s in full QCD as in [43,53,54]. In Fig. 9, we choose to normalize the area of the Monte Carlo distribution to the total area of the NLL  Figure 10: Location of peak of jet shape distribution as a function of a for quark and gluon jets. We plot the value of τ a at the peak of the jet shape distribution for a between -1.0 and 0.8 for quark vs. gluon jets, using R = 1, 0.7, 0.4.  Figure 11: Scale variation of quark and gluon jet shapes. For a = 0 and R = 0.7, we display the variation of the NLL resummed jet shape distributions with the hard scale µ H , the jet cutoff scale µ Λ S , the unmeasured jet scales µ 2,3 J , the measured jet scale µ 1 J (τ a ), and the measured soft scale µ S (τ a ). In each case we vary the scale between 1/2 and 2 times the natural choices in Eq. (6.44), except for the measured soft scale, which we varied between 1 and 2 times the choice in Eq. (6.44). We keep the factorization scale fixed at the default hard scale given by Eq. (7.2), µ = ω i . resummed theory distribution. We find the area under the theory curves for quark and gluon jets to be approximately 0.3 for R = 0.4, 0.5 for R = 0.7, and 0.7 for R = 1. A more accurate prediction of the normalizations may require summing remaining unsummed logs of the phase space cuts in the theory and Monte Carlo predictions. These plots should be interpreted as comparisons of the predictions of the shapes in τ a and these shapes' scaling as we vary a and R, rather than the overall normalization.
The shapes of the theory and Monte Carlo distributions are largely similar, though they display noticeable differences at the leftmost endpoint near τ a = 0 and in the "sharpness" of the peak. These may be due to the different ways the two approaches deal with the growth of the strong coupling for small τ a , the different orders of log resummation (LL vs. NLL) and the need to match the tails onto fixed-order QCD predictions. Since neither the Monte Carlo nor theory partonic predictions without inclusion of hadronization effects is yet a prediction of a physically observable quantity, we use this comparison as an intermediate diagnostic rather than a conclusive test of either method. Nevertheless, comparing the way the shapes of the distributions and locations of the peaks vary over the range of values of a and R that we sample, the behavior agrees very well between the theory distributions and the Monte Carlo distributions without hadronization for both quark and gluon jets.
In Fig. 10 we plot the location of the peak of the jet shape distributions as a function of a for three values of R, displaying the different variation of the peak of quark and gluon jet shape distributions. The peak value increases with increasing R and a, as wide angle radiation is included (increasing R) and less suppressed (increasing a). Although the difference in the peak value between the quark and gluon jet angularity distributions is large, the width of each distribution creates substantial overlap in angularity values between quark and gluon jets. Distinguishing between quark and gluon jets using jet angularities is a complex task which we will explore in future work; for now, we note only that the NLL resummed distributions indicate that discrimination between quark and gluon jets using jet angularities is possible.
As a rough estimate of the theoretical uncertainty in our NLL resummed predictions, we show in Fig. 11 the change in the a = 0 quark and gluon τ a jet shape distributions for R = 0.7 when the various scales that appear in the resummed cross section Eq. (7.1) are varied. These are the initial scales at which the hard, jet, and soft functions are evaluated to minimize logarithms in the NLO fixed-order part, from which the evolution kernels run them to the common factorization scale µ. In the top row of Fig. 11, we vary µ between ω H /2 and 2ω H . The tiny variation is a sign of the consistency condition satisfied by the anomalous dimensions in Eq. (6.33). In the next four rows, we vary the hard scale µ H , the soft jet energy cutoff scale µ Λ S , the unmeasured jet scales µ 2,3 J , and the measured jet scale µ 1 J (τ 1 a ) between half and twice the natural values given in Eq. (7.2). In the last row, we vary the measured soft scale µ 1 S (τ 1 a ) between one and two times the value in Eq. (7.2). This is because too low a value of µ 1 S (τ 1 a ) as τ a → 0 brings it into the nonperturbative region where α s (µ 1 S ) blows up, so that the perturbative estimate of uncertainty is not so meaningful. We note that, while the uncertainty in the vertical scale of the distributions is considerable in some cases, the location of the peak is much more stable.
Finally, in Fig. 12 we give a sense for how robust our theoretical predictions are for other  Figure 12: Jet shapes for other kinematic configurations. We compare our theoretical predictions to Monte Carlo simulations for the shape τ 1 0 (a = 0) for a quark or gluon jet found in a three-jet configuration where the two jets with narrowest separation angle ψ near have equal energy. We consider the two cases ψ near = π/2 and π/3. In the first row, we plot shapes of one of the jets in the "near" pair. The blue solid curve is the shape of quark jet found near a gluon jet, the green dotted curve is a quark found near an antiquark, and the red solid curve is a gluon found near a quark. In the second row, we compare shapes of a quark or gluon jet found far from the near pair. kinematic configurations. We consider e + e − → qqg events where the angle ψ near between two partons is either π/2 or π/3, and these partons have equal energy. We find jets using the anti-kT algorithm with R = 0.4, and plot jet shapes for a = 0. The selection cuts to choose events from the Monte Carlo are the same as the Mercedes-Benz configuration. In these events there are five distinct characterizations for a single parton. If the event has the quark (or antiquark) as the "far" (most well separated) parton, then each parton in the event is distinct: there is the far quark, the near quark, and the near gluon. If the event has the gluon as the far parton, then there are only two distinct types of partons: the far gluon and the near quark (antiquark). In Fig. 12, we plot all these configurations for both ψ near = π/2 and ψ near = π/3. The agreement between the theory predictions and the Monte Carlo are as good as in the Mercedes-Benz case, a good indication that our calculation applies to a broad range of kinematic configurations of multijet events. Additionally, we observe features consistent with our intuition about the relative differences between the jet shape distributions between different jets in these configurations. As one would expect, the distribution of near jet shapes is weighted more heavily towards larger τ a than the far jet shapes, due to the enhanced soft radiation in the near jet system. When the near quark is near a gluon, the distribution is weighted more heavily towards larger τ a than when the near quark is near an antiquark, due to the enhanced radiation coming from a gluon rather than a quark. These distributions serve as further evidence that jet shapes may be an effective discriminant between quark and gluon jets.

Conclusions
In this work, we have factorized an N -jet exclusive cross section differential in M ≤ N jet observables and resummed global logarithms of the jet observable τ a to NLL accuracy, leaving summation of non-global logarithms to future work. We demonstrated that the anomalous dimensions of the hard, jet, and soft functions in the factorization theorem satisfy the nontrivial consistency condition Eq. (6.33) to O(α s ), for any number of quark and gluon jets, any number of jets whose shapes are measured, and any size R of the jets, as long as the jets are well-separated, meaning t 1. This condition ensures the validity of an effective theory with N collinear directions that are assumed to be distinct. We identified and estimated important power corrections to the factorized form of the cross section. We also illustrated that zero-bin subtractions give nonzero contributions to the anomalous dimensions crucial for consistency.
Armed with consistent factorization and the fixed-order jet and soft functions, we resummed large logarithms in the jet shape distribution by running each individual function from the scale where logs in it are minimized to the common factorization scale µ. We thereby resummed, to NLL accuracy, global logs of the jet shape τ a and logs of the scale Λ/E J of soft radiation outside of jets, but leaving some non-global logs and logs of the angular cut R (but we took R to be numerically of order 1). This is the first such calculation of a resummed jet shape distribution in an exclusive multijet cross section.
We constructed a framework to deal with all the scales that appear in the multijet soft function which depends on the values τ i a of all M jet shapes and the phase space cuts Λ, R. By refactorizing the full soft function into individual pieces depending on one of these scales at a time, we were able to sum logs of ratios of these scales.
We demonstrated the accuracy of our results by comparing our resummed prediction for quark and gluon jet shapes in e + e − → 3 jets to the output of Monte Carlo event generators, MadGraph/MadEvent and Pythia. We compared our predictions with the Monte Carlo output without hadronization. The changes in shape and location of the peak value as functions of a and R match quite well between the theory and Monte Carlo.
Our results provide a basis for future studies of other jet observables at both e + e − and hadron colliders, requiring recalculation of those parts of our jet and soft functions that depend on the choice of observable. Studying jets at hadron colliders requires constructing observables appropriate for that environment and the switching of two of our outgoing jets to incoming beams, which can be described by beam functions in SCET [62].
Precision calculations of jet shapes will allow improved discrimination of jets of different origins. We are applying the results of our predictions of light quark and gluon jet shapes to distinguish quark and gluon jets with greater efficiency than achieved before. Extensions to the shapes of heavy jets and calculations of other types of jet shapes such as the Ψ(r/R) shape introduced in [14,15,16] can also be performed.
Note added in final preparation: As this paper was being completed, Ref. [98] appeared reporting the calculation of a quark jet function for a jet defined with a Sterman-Weinberg algorithm and whose invariant mass s is measured. This jet function is related to our measured jet function J q ω (τ a ) for a cone jet at a = 0 given in Eq. (4.11), since s = ω 2 τ 0 . We have checked that the corresponding results between the two papers agree.

Acknowledgments
We are grateful to C. Bauer for valuable discussions and review of the draft. The authors at the Berkeley CTP and in the Particle Theory Group at the University of Washington thank one another's groups for hospitality during portions of this work. AH was supported in part by an LHC Theory Initiative Graduate Fellowship, NSF grant number PHY-0705682. The work of AH and CL was supported in part by the U.S. Department of Energy under Contract DE-AC02-05CH11231, and in part by the National Science Foundation under grant PHY-0457315. The work of SDE, CKV, and JRW was supported in part by the U.S. Department of Energy under Grants DE-FG02-96ER40956.

A.1 Finite Pieces of the Quark Jet Function
Measured Quark Jet Function The finite pieces the jet functions, which depend on the jet algorithm, share common features. For cone-type algorithms, the finite piece of the naive part of the quark jet function,J q alg (τ a ), is given bỹ where in this Appendix, plus distributions are defined by [62] [Θ(x)g(x)] + = lim depends implicitly on τ a and R and is given by The parameter x cone = x cone (τ a ) is the lower limit on the x = q − /ω scaled gluon momentum integral imposed by the cone restriction. It is given by the solution of the equation in the range 0 < x < 1/2, which is plotted in Fig. 13A. The limit τ max a is given by the maximum value over x of Eq. (A.5). Similarly, for k T -type algorithms,J q k T (τ a ) is given bỹ (A.6) I q k T is given by where R is the region in x where the constraint is satisfied. We plot this region in Fig. 13B and C for the cases a > −1 and a < −1, repsectively. The boundaries of this region are the points x 1,2 illustrated in the figure, and are given by the equation f k T (x 1,2 ) = τ a tan 2−a R where we take x 2 > x 1 if x 2 exists. The upper limit τ max a is given by the maximum value over x of the right-hand side of Eq. (A.8). In general, the constraint Eq. (A.8) is symmetric about x = 1 2 , and so the region R is symmetric about the same point. In general, if a > −1 or τ a < 2 a−2 tan (2−a) R 2 , then R is a single range in x. Otherwise, R is two disjoint ranges in x. Since τ a ≥ 2 a−2 tan (2−a) R 2 can only occur for a < −1, we can write I q k T as Note that I q cone and I q k T involve the same integrand, but for each algorithm the integral is over different ranges. In addition, both x cone and x 1 approach the same limiting value for small τ a , where x = x cone or x 1 for the cone and k T algorithms, respectively. Defining r q (x) = 3x + 2 ln 1 − x x , (A. 13) using Eq. (A.12), and including the zero-bin subtraction in Eq. (4.10), we find that the finite distributions of the full measured quark jet functions are Going through similar steps as for the quark jet function, defining and using Eq. (A.12) to make all logarithmic dependence on τ a explicit, we find for the cone and k T -type jet function finite distributions where z = (1−n 2 )(1−u 2 ) (1−un) 2 . The integration over u = cos θ has singularities at the points u = 1 and u = n which correspond to z = 1 and z = 0, respectively. To isolate these singularities, we split the integration over u into the ranges (−1, δ) and (δ, 1) where n < δ < 1, I incl (n i · n j ) = I incl 1 (n i · n j ) + I incl 2 (n i · n j ) , (B.3) where I incl 1 (n i · n j ) ≡ Over the range of integration of u in I incl 1 , z ∈ [0, 1) for δ < 1. For I incl 2 , z ∈ (0, 1]. Furthermore, the singularity at u = n in I incl 1 is made more explicit through the use of the identity f a (z) gives an O(1/ ) contribution and we proceed by using the following trick that we exploit multiple times throughout the Appendix.
To integrate a product of functions f (x, )g(x, ) where f diverges at the point x 0 as (x − x 0 ) −1+O( ) , we write the integation as dx f (x, )g(x, ) = dx f (x, )g(x 0 , ) + dx f (x, ) g(x, ) − g(x 0 , ) . (B.6) The first integral has relatively simple x dependence since g(x 0 , ) does not depend on x.
The term in parenthesis in the second integral vanishes as x − x 0 for regular functions g and so the entire integrand can be expanded in .
We can now evaluate f a (z) by adding and subtracting the non-singular part of the integrand (which is the hypergeometric function) evaluated at u = n as in Eq. (B.6), whereas f b (z) is O( ) and so we can simply expand about = 0. Adding these contributions, we find that B.2.2 S meas ij (τ i a ) To evaluate Eq. (5.17) for the case that k = i, we use light cone coordinates in the frame of jet i, k + = n i · k and k − =n i · k. In terms of these variables, the on-shell condition can be used to give n j · k = k + cos 2 ψ ij 2 + k − sin 2 ψ ij 2 − √ k + k − sin ψ ij cos φ, (B. 16) with cos ψ ij = 1 − n i · n j , and φ the angle in k ⊥ -space (the azimuthal angle about n i ). We can do the k ⊥ and k + integrals using the on-shell and τ a delta functions respectively. The resulting S meas ij (τ i a ) has non-trivial integrals over k − and φ:

B.2.3 S i ij
The Θ-functions in Eq. (5.14) are easiest to deal with if we shift to variables where each Θ-function is in a different variable. The simplest choices are just the arguments of the Θ functions Θ Λ and Θ i R , k 0 and u = cot R 2 k + /k − , respectively, where k ± are defined with respect to direction n i . This gives a form similar to the integral in S meas ij (τ i a ), where I(α, β, t) is defined in Eq. (B.10) and evaluates to Eq. (B.15). This gives Eq. (5.15).

B.3 S meas
ij (τ k a ) and S k ij for k = i, j We again use light cone coordinates centered on jet k. The integrations involved in S meas ij (τ k a ) and S k ij only give rise to a 1/ pole as explained in the text, but integrating the eikonal factor 1/(n i · k)(n j · k) is more complicated than for the other cases since there is a third direction, n k , involved.
For unmeasured jets when there are n ≥ 3 total final state jets, S k ij is needed. However, as we explain in the text, measured jets violate consistency at O(1/t 2 ) even for n = 2 (non back-to-back) jets and the contribution of S k ij does not ameliorate this fact when n ≥ 3. To show this, we need to evaluate the divergent contribution of S k ij . In addition, we give the form of the finite pieces which are O(1/t 2 ).
For each measured jet when there are n ≥ 3, the sum S meas ij (τ k a ) + S k ij δ(τ k a ) is needed. However, in this case the 1/ pole cancels in this sum and we are left with only a single, finite integral to evaluate. This is clear from the expressions for S meas ij (τ k a ) and S k ij which we derive in Sec. B.3.2 and Sec. B.3.3, respectively. We evaluate the sum explicitly in Sec. B.3.4.

B.3.1 Common Integrals
We find the following integral arising in both S meas ij (τ k a ) and S k ij : We can evaluate I (0) straightforwardly. For the range of our interest, t a,b > 1 and 0 < u < 1, it gives .

C. Convolutions and Finite Terms in the Resummed Distribution
In evaluating the final resummed distribution Eq. (6.40), each measured jet function must be convolved against a corresponding soft function piece S meas . These convolutions take the form.
(C.1) For the class of functions of the form x −1−ω with ω = 0 and ω < 1, we define the plus distribution by where the plus functions on the second line are given by Eq. (A.2), From these definitions, we can derive the identities (see, e.g., Appendix B of [54]) Using the above identities, we find that the final result can be written in the form Eqs. (6.40) and (6.42) with the functions d J (τ i a ) given by

D. Color Algebra for n = 2, 3 Jets
For the two and three jet cases, there are no color correlations since all color generator inner products T i · T j can be expressed in terms of the Casimir invariants C A and C F . For n = 2, there is a quark jet with charge T q and an anti-quark jet with charge Tq that each square to C F . There is only one inner-product in this case and using color conservation ( i T i = 0), we have that For n = 3 jets color conservation gives that, for example, Referring to the quark, anti-quark, and gluon generators as T q , Tq, and T g , respectively, using T 2 q = T 2 q = C F and T 2 g = C A in Eq. (D.2) gives