Factorization and resummation for jet processes

From a detailed analysis of cone-jet cross sections in effective field theory, we obtain novel factorization theorems which separate the physics associated with different energy scales present in such processes. The relevant low-energy physics is encoded in Wilson lines along the directions of the energetic particles inside the jets. This multi-Wilson-line structure is present even for narrow-cone jets due to the relevance of small-angle soft radiation. We discuss the renormalization-group equations satisfied by these operators. Their solution resums all logarithmically enhanced contributions to such processes, including non-global logarithms. Such logarithms arise in many observables, in particular whenever hard phase-space constraints are imposed, and are not captured with standard resummation techniques. Our formalism provides the basis for higher-order logarithmic resummations of jet and other non-global observables. As a nontrivial consistency check, we use it to obtain explicit two-loop results for all logarithmically enhanced terms in cone-jet cross sections and verify those against numerical fixed-order computations.


Introduction
A crucial ingredient for the factorization of high-energy processes is the simple form of soft emissions. Their eikonal structure forms the basis of the resummation of soft-photon effects in QED [1][2][3] and proved to be equally relevant when analyzing factorization in JHEP11(2016)019 QCD [4][5][6]. In position space, soft emissions are described by Wilson lines along the classical trajectories of the energetic particles [7,8] and the renormalization properties of the corresponding Wilson-line matrix elements [9][10][11][12][13][14] play an important role when resumming soft-gluon effects. An important property of soft radiation is that wide-angle radiation off a jet of collinear particles does not resolve the individual energetic partons. Instead of probing the charge distribution inside the jet, the long-wavelength soft radiation is sensitive only to the total charge, because the eikonal factors associated with the collinear momenta p µ i can be expanded around a common reference vector n µ pointing along the jet direction, where k and ǫ are the momentum and polarization of the emitted soft gluon. This approximation fails for a soft gluon which is collinear to the jet, but typically this small region of phase space does not give an unsuppressed contribution to the cross section. As long as this is true, the soft function for a dijet cross section involves only two Wilson lines, accounting for soft radiation emitted from the two jets. The cross section then factorizes into a product (in the appropriate space; more generally a convolution) of this soft function S, jet functions J and J for the collinear radiation inside the two jets, and a hard function H encoding the virtual corrections to the production of the two leading partons: Computations based on this structure were successfully used to resum large logarithms in many event-shape variables, in particular thrust, heavy jet mass or the Cparameter [15,16]. In the meantime, using renormalization-group (RG) methods in Soft-Collinear Effective Theory (SCET) [17][18][19] (see [20] for a review), the resummations for selected observables has been carried out up to next-to-next-to-next-to-leading logarithmic accuracy [21][22][23]. However, it was found that certain observables contain single-logarithmic terms which are not captured by resummation techniques based on (1.2). The classic example is the light-jet mass event shape in e + e − collisions, or equivalently, the jet mass in a single hemisphere. In [24], Dasgupta and Salam traced the problem back to the fact that these observables are insensitive to radiation in certain regions of phase space (the right hemisphere in the example of the left-hemisphere mass). The additional logarithms first appear at two-loop order and are referred to as non-global logarithms (NGLs). From an analysis of soft emissions, Dasgupta and Salam extracted the leading two-loop logarithm analytically and gave an algorithm to numerically compute the leading higher-order logarithms in the large-N c limit. Based on the simple form of amplitudes in the strongly energy-ordered limit (see e.g. [6,25]), Banfi, Marchesini and Smye (BMS) derived a non-linear integral equation which can be used to perform the resummation of the leading NGLs in the large-N c limit [26]. However, since strong energy ordering is a crucial ingredient for the BMS equation, its logarithmic accuracy cannot easily be increased. What has been achieved is leading-logarithmic resummation of these logarithms for N c = 3 [27,28]. It turns out that the corrections to the large-N c limit are numerically quite small [29].

JHEP11(2016)019
Since the vast majority of collider observables include hard phase-space cuts or, more generally, regions of phase space in which radiation is not restricted, the presence of NGLs severely limits the applicability of higher-order resummation techniques. For this reason, a lot of effort was put into trying to get a better understanding of these types of logarithms. For example, several groups computed hemisphere soft functions up to two-loops to obtain the full result for the non-global structure at this order [30][31][32][33]. Also, using the BMS equation, the analytic result for the leading-logarithmic terms up to five-loop order was extracted [25]. This analysis was also extended beyond the large-N c limit by computing the higher-order terms directly from strongly-ordered soft amplitudes [34]. While these fixedorder computations provide important insights into the form of NGLs, ultimately one is interested in their all-order structure. Steps towards a resummation of such terms were recently taken in [35,36]. The authors claim that the NGLs arise from soft subjets near phase-space boundaries and propose a set of factorization theorems which resum global logarithms in the presence of subjets. They then argue that this resummation will capture a large part of the NGLs in more inclusive cross sections. They propose an expansion in the number of soft subjets, which they call "dressed gluons". Since the dressed gluons include Sudakov factors, the expansion in dressed gluons does not suffer from the same divergence as standard perturbation theory when the logarithms become large. However, an arbitrary number of soft subjets contributes even to the leading NGLs, and it is not clear what expansion parameter governs the expansion in subjets and whether there is any parametric suppression of the higher-multiplicity terms. 1 From a numerical point of view (see figure 8 in [35]), the expansion in dressed gluons appears to provide only a modest improvement over a pure fixed-order treatment for moderately large values of the logarithms. A second interesting proposal to go beyond leading-logarithmic resummation is the functional RG of Caron-Huot [37]. We will comment in more detail on both approaches and their relation to our results below.
An important example of non-global observables are jet cross sections, in particular those involving cone jets, which are insensitive to radiation inside the jet cone. In the present paper we analyze dijet cross sections in e + e − collisions. In addition to the narrowjet case analyzed in our previous work [38], we also treat the case where the opening angle of the jet cone is large. For brevity, we will refer to these as "wide-angle jets". We find that in both situations, the simple factorization theorem (1.2) is incorrect. This is immediately obvious for wide-angle jets since the jet opening angle is as large as the typical angle of the soft radiation. The approximation (1.1) is therefore not appropriate and each hard parton inside the jet produces its own Wilson line. In the next section we show that the relevant factorization theorem for the cross section takes the form (1. 3) The function H m is the squared amplitude for having m particles inside the two jets, integrated over their energies but at fixed angles. The function S m contains soft Wilson

JHEP11(2016)019
lines along the directions of the m hard partons. The symbol ⊗ indicates an integral over the angles specifying these directions. The functions H m and S m are matrices in the color space of the m hard partons and the angle brackets indicate that one has to take the trace of their product to get the cross section. One could naively expect to recover the form (1.2) in the narrow-jet case, but this turns out not to be true because small-angle soft radiation gives a leading-power contribution to the cross section. The relevance of a new mode describing small-angle soft radiation was demonstrated in [38,39]. Since the radiation is simultaneously collinear and soft, we call it "coft" to distinguish it from the former two types of radiation. In the narrow-jet case we find that, for k energetic partons inside the left jet and l inside the right one, the additional factorization is valid, up to corrections suppressed by the small cone angle. The first equation is the usual collinear factorization, the statement that the hard function factorizes into a hard function for the two jets times collinear splitting functions J l and J k for the partons inside the jets. The second equation states that the soft radiation splits up into wide-angle soft radiation consisting of two Wilson lines and coft radiation in the direction of the left or right jet, which resolves the individual energetic partons inside the jets. The tildes in this equation indicate that the product form holds in Laplace space; in momentum space the factorization would involve a convolution. The coft function U l contains l Wilson lines from particles inside the right jet and a single Wilson line along the left jet, and analogously for U k . We use the method of regions [40,41] and consistently expand both the amplitudes and the phase-space constraints in small momentum components to avoid double counting between different momentum modes. The expansion of the phase-space constraints is an important ingredient of our approach, which leads to complete scale separation and avoids the necessity of performing zero-bin subtractions. The renormalization of the operators S m is quite interesting and nontrivial [38]. The factorization theorem (1.3) splits the cross section into a sum of ingredients with a fixed number of hard partons. Scattering amplitudes with a fixed number of partons are not finite, because the virtual corrections to lower multiplicities are needed to cancel the infrared (IR) divergences of the real-emission amplitudes. In the effective field-theory framework we are using, the IR divergences of the hard functions are in one-to-one correspondence to the ultraviolet (UV) divergences of the Wilson-line matrix elements S m , see e.g. [12,13]. The fact that the cancellation of divergences involves different multiplicities then implies that the renormalization matrix Z km which absorbs the divergences of S m is not diagonal, or more specifically, that higher-multiplicity Wilson-line operators mix into lower-multiplicity ones under renormalization. The reason for this mixing are UV divergences arising from soft emissions inside the jets, which are unconstrained. The scale dependence of the different components of the factorization theorem is driven by anomalous dimensions obtained from the respective Z-factors. Solving the relevant RG equations resums all large logarithms in cone-jet processes [38].
The focus of the present paper is on the derivation of the factorization theorems and the RG equations which govern its ingredients. A large fraction of our paper is devoted JHEP11(2016)019 to computing the ingredients for the next-to-next-to-leading (NNLO) jet cross section in both the narrow-jet and wide-angle cases. These computations provide us with nontrivial consistency conditions on our results and show that our formalism correctly predicts all logarithmic terms. Moreover, we establish that RG evolution in the effective theory reproduces the known results at leading logarithmic accuracy in the large-N c limit. A detailed study of the RG equation and its solution beyond this accuracy is left for future work. Given that all jet observables as well as most other experimental measurements are non-global, there are many potential applications for our framework, once methods have been developed to solve the associated RG evolution equations.
The structure of the paper is as follows. In section 2, we analyze the factorization properties of cone-jet cross sections, both in the wide-angle case and for narrow jets. We also discuss the renormalization of the operators which appear in these formulae. As a check of the factorization theorem for narrow jets, we compute in section 3 all ingredients to two-loop accuracy. The cancellation of divergences provides an important consistency check and we verify that our analytical result is compatible with numerical fixed-order results. In section 4, we perform the analogous two-loop analysis for wide-angle jets. We again check our results against numerical computations and then verify their factorization properties in the limit small-angle limit. We give explicit results for the one-loop anomalous dimensions matrices which govern the resummation in section 5. We also verify that for the leading NGLs at large N c our results reduce to the ones obtained from the BMS equation. After a short discussion of methods to solve the RG equations and a comparison to the approaches of [35,37], we conclude. Some lengthy two-loop expressions are relegated into appendices. Appendix A provides expressions relevant for the narrow-cone case, while those for the wide-angle case are given in appendix B.

Factorization of jet cross sections
In this section we will derive a factorization formula for the cross section for e + e − → 2 jets at center-of-mass energy √ s = Q. We use the thrust axis, defined as the unit-vector n in the direction of maximum momentum flow, as the jet axis and define two light-like vectors n µ = (1, n) andn µ = (1, − n) along the jets. Using these vectors, we can rewrite any four-momentum as p µ =n · p n µ 2 + n · pn We use the thrust axis to split each event in two hemispheres and call particles with n · p <n · p right-moving. The definition of the thrust axis implies that the total transverse momentum in each hemisphere vanishes. Particles are considered to be part of the right jet if n · p < δ 2n · p, where the parameter δ is related to the opening angle α of the jet via 2) see figure 1. To define the cross section, we impose that the total energy emitted outside the left and right jet cones, fulfills the condition 2E out < βQ. Except for the choice of the jet axis, our definition is identical to the one in the seminal paper of Sterman and JHEP11(2016)019 n α δ = tan(α/2) 2E out < βQ Figure 1. Definition of the parameters δ and β of the dijet cross section. We use the thrust axis n as the jet axis.
Weinberg [42]. Using the thrust vector as the jet axis leads to a simpler form of the phase-space constraints and enables us to use existing two-loop results for the cone-jet soft function obtained in [32,33].

Wide-angle jets
Let us first consider wide-angle jets with δ ∼ 1. In this case the effective theory contains only two relevant momentum regions, whose components (n · p,n · p, p ⊥ ) scale as follows: The hard mode describes the energetic particles inside the jet. Since we are dealing with wide jets, the energetic radiation inside the jet covers a large angular range. It is thus not collinear to n but has a homogenous scaling of all components. Given their large energy, these particles can never go outside the jet, in contrast to the soft partons which can be emitted inside or outside. Since there are no collinear singularities for large cone size, the cross section is single-logarithmic, i.e. the leading logarithms have the form α n s ln n β. The factorization of an amplitude with m hard partons and an arbitrary number of soft partons is of course well known. Each hard parton gets dressed with a Wilson line along its direction. For an outgoing particle in the color representation T i propagating along the direction n i , the appropriate Wilson line is given by the path-ordered exponential The Wilson line S i is a matrix in color space, which acts on the color index of particle i. The operator for the emission from an amplitude with m hard partons then takes the form where n µ i = p µ i /E i , and we use the compact notation {p} ≡ {p 1 , p 2 , . . . , p m }. This equation is analogous to the factorization for amplitudes with coft particles [38], but while the coft case involves splitting amplitudes, we are now dealing with ordinary amplitudes |M m ({p}) . In writing (2.5) we use the color-space formalism of [43,44], in which amplitudes are treated as n-dimensional vectors in color space. Since they act on different particles, the different generators trivially commute [T a i , T b j ] = 0 for i = j. The same is therefore true for the JHEP11(2016)019 associated Wilson lines. Note that the gluon fields in the product of Wilson lines are timeordered, but for brevity, we do not indicate this explicitly. Since they commute, Wilson lines along common directions immediately combine into single Wilson lines, for example This property ensures that collinear particles only produce a single Wilson line carrying the total color charge. However, since we deal with large-angle jets, the individual Wilson lines do not combine in our example. To derive formula (2.5) in the effective field theory we introduce a separate collinear field for each of the energetic particles in the final state, i.e. we write down the SCET operators for processes with m jets. This is possible since on the amplitude level there is no difference between collinear and hard on-shell particles. The relevant purely collinear SCET Lagrangian consists of m copies of the ordinary QCD Lagrangian. Operators in the effective theory are conveniently expressed in terms of gauge-invariant fields χ i and A µ i⊥ , which are related to the usual quark and gluon fields via [45] The i-collinear Wilson lines in the fundamental representation are defined analogously to the soft Wilson lines in (2.4) as The argument denotes the direction of the Wilson line, which is conjugate to the direction n i of the collinear particle. These Wilson lines ensure that these fields are invariant under collinear gauge transformations in each sector [17,18]. At leading order in power counting, m-jet operators in this effective theory involve exactly one collinear field Φ i ∈ {χ i ,χ i , A µ i⊥ } from each sector i = 1, . . . , m. Performing the usual decoupling transformation with the appropriate color representation T i for each field, yields the Wilson-line structure shown in (2.5). Finally, one evaluates the matrix element of the operator with one collinear particle in each sector, using (2.10) Along with the Wilson coefficient of the m-jet operator this gives the amplitude |M m ({p}) , see [13] for details. Since the particles are on the mass shell, the higher-order corrections to the relations (2.10) are all scaleless and vanish.

JHEP11(2016)019
To get the amplitude for the emission of l soft partons in the final state with momenta k 1 , . . . , k l , one computes the matrix element k 1 , . . . , k l | S 1 (n 1 ) S 2 (n 2 ) . . . S m (n m ) |0 (2.11) of the Wilson-line operator. To obtain the contribution of an arbitrary number of soft partons to the jet cross section, one first defines the squared matrix element for the emissions from m partons as (2.12) This is the same as the coft function which arises for narrow-angle jets [38], up to the fact that the constraint now acts on the out-of-cone energy E out of the soft radiation, as opposed ton · p out , the large component of the total momentum of the coft fields. Since the soft function depends on the outside energy, it depends on the cone size δ. In terms of the matrix element (2.12), the jet cross section takes the form 13) up to terms suppressed by powers of β. The integration is over the m-dimensional phasespace of the hard partons, which are all constrained to lie inside the two jet cones. The function Θ nn in p ensures that the hard partons are either inside the right jet along the direction n or the left jet alongn. In the narrow-cone case, we will encounter constraints which involve only one of the jets. Note that, due to the multipole expansion, the contribution of soft particles must be neglected in the momentum-conservation δ-functions.
In order to write the cross section in a more transparent way, we now define hard functions which are obtained by integrating over the energies of the hard particles subject to the constraint that their sum is equal to the center-of-mass energy Q, while keeping their directions n µ i fixed, (2.14) These hard functions are distribution-valued in the angles of the particles, since they contain additional divergences which arise when particles become collinear. These realemission divergences get cancelled by the divergences associated with the virtual corrections to amplitudes with fewer legs. In contrast, the soft function  In contrast to the standard formula (1.2), the factorization formula (2.15) does not involve jet functions. The reason is that there is no collinear scale in our problem. The collinear matrix elements in (2.10) are scaleless and do not receive higher-order corrections. In dimensional regularization, there is thus a cancellation between IR and UV singularities in these functions. When added to the hard functions, the IR divergences in the collinear matrix elements cancel against IR divergences of the hard functions, so that the net effect is to convert the IR divergences of the hard function into UV divergences. It would be possible to separate IR and UV singularities at each step by introducing parton masses or additional subjet resolution parameters into our analysis, however this would only complicate the problem in an unnecessary way.
We have obtained the formula (2.15) from an analysis of QCD amplitudes in the hard and soft momentum regions. In the context of a low-energy effective theory, the hard momentum modes are integrated out and the functions H m ({n}, Q, δ) are Wilson coefficients of the soft operators S m ({n}, Qβ, δ). Because all components of the soft modes are smaller than their hard counterparts, all hard-soft interactions in the Lagrangian are power suppressed. As in standard SCET, the soft Lagrangian is the same as the QCD Lagrangian. To resum large logarithms of β, one has to solve the RG equations of the Wilson coefficients. This will be discussed in detail below, after we have analyzed the factorization properties for the narrow-jet case.
In our derivation of the factorization theorem we first analyzed the factorization properties on the amplitude level and then computed the cross section with the factorized amplitudes. As we discussed, the relevant factorization of the amplitude shown in (2.5) can easily be obtained from SCET with collinear fields along the directions of the energetic on-shell particles. However, in standard SCET derivations one is usually working directly on the level of the cross section. For our approach to be valid, it was important to know the relevant momentum regions in the cross section, so that the appropriate expansion of the amplitude could be performed. In other words, it was important that the phasespace integrations indeed only involved the regions we considered for the expansion of the amplitude.

JHEP11(2016)019
Qβδ 2n · p n · p Figure 2. Momentum regions relevant for narrow-angle jet production. The plot shows the scaling of the light-cone components n · p andn · p, and we assume that β ≪ δ (we use β ∼ δ 2 in the narrow-jet case to ensure this condition). The meshed gray area shows the veto in the out-of-jet region which forbids the presence of energetic modes. In the wide-angle limit δ ∼ 1, soft and coft modes coincide and the collinear and hard scalings are the same.
What is special about the present case is that every single component of the coft momentum is suppressed compared to its collinear counterpart, while in standard applications the small component of the collinear mode is commensurate with the soft scaling. In our case, the relative scaling is the same as for the soft and hard modes present in the wide-angle case. Because of this fact, the coft Wilson lines associated with different collinear particles cannot be combined and we end up with multi-Wilson-line operators. We emphasize that the coft modes have very low virtuality p 2 t = Λ 2 t = (Qβδ) 2 , much lower than the virtuality of the collinear and soft modes. The presence of this low physical scale might have important implications for the relevance of non-perturbative effects. These

JHEP11(2016)019
are suppressed by the ratio Λ QCD /Λ t , where Λ QCD ∼ 0.5 GeV is a scale associated with strong QCD dynamics. Non-perturbative corrections to jet processes can thus be much larger than the naive expectation Λ QCD /Q. For example, for a jet opening angle α = 10 • (δ ≈ 0.09) and 5% of the collision energy outside the jets (β = 0.1), one obtains Λ t ≈ 1 GeV for Q = 100 GeV. It would be interesting to explore phenomenological consequences of this low-scale physics.
One way to construct the effective theory containing coft modes is to first match QCD onto standard SCET containing collinear fields along the jet directions and soft fields. One then matches the QCD quark vector current onto the vector current in the effective theory and obtainsψ The field χ c ) describes an energetic collinear quark propagating in the n direction (n direction), as defined in (2.7). The Wilson coefficient C V contains matching corrections from hard virtual particles. In the above equation we have already decoupled the soft fields which gives rise to the product S(n) S † (n) of two soft Wilson lines in the fundamental representation. Note that, as in (2.4), we use outgoing soft Wilson lines. After this decoupling the soft and collinear fields no longer interact. From now on we drop the superscript "(0)" on the fields, since we will need to perform a second decoupling transformation below.
In a next step, one splits the collinear field into two submodes 20) and analogously for the quark fields. More explicitly, one matches the purely collinear theory onto a second effective theory, which distinguishes genuine collinear from coft momenta. Note that the relative scaling of the momentum components of the collinear and coft particles is exactly the same as the one of the hard and soft modes analyzed in the wide-angle jet case. Indeed, if we perform a boost n · p → n · p/δ,n · p →n · p δ and identifyQ = Qδ with the hard scale, we end up with exactly the same configuration as in the wide-angle case discussed in the previous subsection: after the boost the collinear fields become hard fields at the scaleQ and the coft fields correspond to soft particles with momentum scaleQβ. We can thus analyze the problem with the same method used in section 2.1. For the case with m particles inside the jet, we introduce collinear fields along their directions with scaling p 2 = (Qβ) 2 and then write down the corresponding leadingpower operators. One important difference is that the "full theory" expression for this second matching step contains collinear quark fields which contain collinear Wilson lines. Specifically, let us consider the second term in (2.19) which contains the field χ c . This field involves a Wilson line W † c (n) along then direction, see (2.7). Since the full theory contains gluons in the direction n 0 ≡n, we also need to introduce a corresponding collinear field along the n 0 direction in the effective theory. We thus define a new collinear building block This building block is the remnant of the anti-quark field (labeled as particle "0") in the original theory and behaves like an anti-quark under the decoupling transformation. It

JHEP11(2016)019
is of O(1) in SCET power counting, while the fields χ 0 and A 0⊥ count as O(β). As a consequence, the general leading-power operator for m partons inside the jet is After performing the decoupling transformations in complete analogy to expression (2.5) obtained in the wide-angle case. The relevant collinear amplitude is We have added the argument p 0 to indicate that these matrix elements are the splitting amplitudes of the collinear quark into the given final state. Since they connect the color of the quark field to the color of the final-state particles, they are strictly speaking matrices in color space, and a more precise but also somewhat cumbersome notation would be to denote them by Sp m (p 0 ; {p}).
In order to check these arguments diagrammatically and to confirm the factorization of soft, collinear and coft fields shown in (2.19) and (2.24), we have explicitly expanded the squared γ * (q) → q(p 1 )q(p 2 ) g(p 3 ) g(p 4 ) tree-level QCD amplitude in all relevant momentum regions and compared the results with (2.19) and (2.24). For the leading contributions in the different regions, we obtain We have indicated the assumed scaling of the gluon momenta as a superscript, while the (anti-)quark momentum has always (anti-)collinear scaling. The functions are squared one-and two-gluon matrix elements of an operator with Wilson lines along the directions {n} = {n 1 , . . . , n m }. The one-gluon matrix elements A (1) are given by the expressions shown in (3.8) below, while the two-gluon matrix elements A (2) can be found e.g. in appendix C of [48]. The splitting functions F (1) and F (2) describe the emission of one or two collinear gluons from the quark with momentum p. The function F (1) (p; p 1 , p c 3 ) is given by the integrand in (3.23) below. We can now compare the results in (2.26) with the factorized structures derived above. The first three relations are quite simple. The matrix elements factorize into the Bornlevel matrix element times the squared two-gluon amplitudes for each sector. The next JHEP11(2016)019 two factorization relations are more interesting. They demonstrate that for two energetic particles with momenta p 1 and p 3 inside the right jet one ends up with a two-Wilson-line matrix element for a soft gluon, as in (2.19), but a three-Wilson-line matrix element for a coft gluon, confirming our relation (2.24). Finally, for the case of M(q; p 1 , p 2 , p s 3 , p t 4 ) one ends up with a product of a soft and a coft matrix element, demonstrating the coft-soft factorization. We see that all relations in (2.26) are fully compatible with the structure of (2.19) and (2.24).
Squaring the factorized amplitude and expanding the phase-space constraints in each region, we obtain a factorization theorem for the narrow-angle jet cross section. However, before presenting its final form, let us discuss an alternative way to derive the mode factorization in the narrow-jet case. To this end, we start with the wide-angle result (2.15) and narrow the jet angle. To take the limit δ → 0, we can group the final-state particles into k left and l right ones, {p} = {p L } ∪ {p R }. For small cone size, the amplitude factors into the amplitude for the production of a quark and an anti-quark, multiplied by splitting amplitudes The hard functions are the square of these amplitudes, and hence we obtain the factorized result While the hard function H k+l depends both on the cone angle δ and the collision energy Q, this dependence is factorized on the right-hand side. The jet functions only depend on the collinear scale Qδ and for l partons they are defined as where p Xc = m i=1 p i is the total momentum of the collinear particles. The theta function Θ n in ensures that every collinear parton is inside the right jet, while the analogous constraint in the definition of the hard function (2.14) was allowing for particles inside both jets. Again, we only integrate over the energies of the particles, keeping their directions fixed. The integrals over the directions will be performed after multiplication with the coft matrix elements. The function J k describes the energetic particles inside the left jet and is defined exactly as (2.29) but with n ↔n.
A similar factorization holds also for the soft functions S m in (2.12). However, because the energy deposited outside the jet cones is shared between the soft and coft modes, we need a Laplace transformation to factorize the corresponding phase-space constraints. We define the Laplace-transformed functions as

JHEP11(2016)019
Note that τ ∼ β for power-counting purposes. The derivative in this definition is present because we introduce the Laplace transform for the differential soft function, defined with a fixed energy, rather than the integrated function S m given in (2.12). We can then factorize the Laplace transformed soft function as where the soft function S(Qτ ) only includes two Wilson lines along the jet directions. To arrive at this factorized expression one starts with the original soft function and splits up the soft field as A s → A s + A t . One can then decouple the coft field using the usual field redefinition. This redefinition yields a new soft field A ′ s , which is split up into A ′ s → A ′ s +At. The anti-coft-soft interactions can then be decoupled in turn. The l soft Wilson lines in the right hemisphere can then be expanded around a common reference vector n and combined using (2.6), and similarly for the soft Wilson lines in the left hemisphere. This yields the same Wilson-line structure as in (2.19). Squaring the soft amplitudes yields (2.32) Because the soft radiation has parametrically large angle, it is always outside the jet and the energy constraint is imposed on the total energy E Xs . The coft function U m ({n R }, Qδβ) with m Wilson lines is given by The right-moving coft particles are always outside the left jet in the sense that the out-ofjet constraint is always fulfilled after the multipole expansion, independent of the angle of the coft particle. The momentum p out is therefore the total momentum outside the right jet. The anti-coft function U k has the Wilson line U 0 along the n instead of then direction and the constraint is placed on n · p out . Putting these ingredients together, the cross section in Laplace space takes the form [38] σ(τ, δ) where we have used the fact that both jets give an identical contribution. In figure 3 we show a pictorial representation of the structure of the factorization formula and the different types of radiation relevant in both the wide-angle and narrow-jet cases.

Renormalization and resummation
The factorization theorems we have obtained involve operators with an arbitrary number of Wilson lines, both in the wide-angle and narrow-jet case. We now discuss the renormalization of these operators. The associated RG equations form the basis for the resummation of logarithmically-enhanced contributions to all orders in perturbation theory. This resummation is achieved by evolving the Wilson coefficients of these operators from the high scale µ ∼ Q down to the scale where the low-energy physics takes place. Let us first discuss the wide-angle cross section for which the factorization theorem has been given in (2.15). In our effective theory, the hard functions H m are the Wilson coefficients of the Wilson-line matrix elements S m and we regularize both quantities in d = 4−2ǫ dimensions. The effective field theory matrix elements contain UV divergences since the short-distance structure of the full theory is not resolved. The corresponding 1/ǫ poles can be removed by renormalizing the hard Wilson coefficients according to In practice, it is easiest to obtain the bare Wilson coefficients from on-shell matching calculations, where the poles arise from IR divergences. However, these IR poles are in one-to-one correspondence to UV divergences since the effective-theory loop-integrals in such matching computations are scaleless, see e.g. [13] for a detailed explanation of this point within SCET. We have discussed this correspondence after (2.15). It implies that we can understand the UV divergences of H m from the structure of the IR divergences in the real and virtual diagrams which contribute to these quantities. Given that the coefficients H m are fixed-multiplicity QCD amplitudes squared, integrated over energy, it is clear that the matrix Z H lm ({n}, Q, δ, ǫ, µ) cannot be diagonal: lower-multiplicity virtual diagrams are needed to cancel the divergences of real-emission diagrams. In order to achieve this cancellation, the renormalization matrix must have the form where we indicate the perturbative order of each element. At each higher order in perturbation theory, more off-diagonal contributions fill in. We have anticipated the upper JHEP11(2016)019 diagonal structure of the matrix in (2.35) by restricting the sum to l ≤ m. Note that Z H lm ({n}, Q, δ, ǫ, µ) has logarithmic Q dependence, because the fixed-multiplicity amplitudes involve both soft and collinear divergences. This dependence is a familiar feature of Sudakov-type processes.
By consistency, the matrix Z H must render the soft functions finite, i.e. we must find that the functions are finite for ǫ → 0. The structure of this result is at first sight quite surprising, since Wilson-line matrix elements can usually be renormalized multiplicatively. However, in the present case additional UV divergences in the real-emission diagrams arise because the soft radiation is not constrained inside the jet. It is precisely those types of divergences which lead to NGLs. Furthermore, the upper triangular form of Z H lm implies that highermultiplicity soft functions are needed to absorb the divergences of matrix elements with fewer Wilson lines. The symbol⊗ indicates that in (2.37) one has to integrate over the (m− l) additional directions of the unresolved partons on which the bare function S m depends.
The scale dependence of the renormalized hard and soft functions is governed by the RG equations which ensure that the cross section (2.15) is scale independent. The anomalous-dimension matrix is obtained from the standard relation and it has linear dependence on ln(Q/µ) as is familiar from Sudakov-type problems. However, the wide-angle cross section we consider contains only a single large logarithm at each order. The Sudakov double logarithms must cancel in the sum over multiplicities in (2.15). A related observation is that the RG equation (2.39) for the soft functions is only consistent if the Q-dependence of the anomalous dimension drops out after the integrals over the unresolved partons have been performed, since the expression on the left-hand side only involves the soft scale Qβ. This implies a set of highly nontrivial consistency relations among the entries of the anomalous-dimension matrix. At one-loop order this will be studied in section 5. Solving the RG equations (2.38) and (2.39) one can resum all large logarithms in the wide-angle jet cross section (2.15). At the soft scale µ s ≈ Qβ the soft functions do not involve large logarithms, and hence they can be calculated in a perturbative series in powers of α s (µ s ). Likewise, at the hard scale µ h ≈ Q the hard functions do not involve JHEP11(2016)019 large logarithms, and hence they can be calculated in a perturbative series in powers of α s (µ h ). The large logarithms of the scale ratio µ h /µ s are resummed by evolving the soft functions up to the hard scale (or vice versa), with an evolution matrix of the form The path-ordering symbol P is necessary since Γ H is a matrix. The resummed cross section is then given by The renormalization procedure in the narrow-jet case is quite similar, except that the underlying factorization formula (2.34) is more complicated. The RG evolution of the functions H and S is well known. Both are renormalized multiplicatively, and explicit expressions for γ V , γ W and the cusp anomalous dimension Γ cusp at threeloop order can be found in the same reference. The jet functions J m for a fixed number of partons cannot be renormalized multiplicatively, as was the case for H m . In analogy to (2.35), we write where Z J is an upper triangular matrix with the same structure as in (2.36). RG invariance of the cross section in (2.34) implies that the product of renormalization factors must render the coft matrix elements finite, i.e.

JHEP11(2016)019
From this one can derive the evolution equation where the anomalous dimension Γ J is defined in analogy to (2.40), and 2γ φ = γ W − γ V governs the non-cusp part of the RG evolution of the quark parton distribution function near x → 1 [50,51]. Similarly to (2.39), the RG equation for the coft functions implies nontrivial consistency conditions on Γ J . The dependence of this function on the jet scale Qδ must contain a universal piece proportional to Γ cusp , which can be factored out and conspires with ln τ term in (2.50) to produce a logarithm of the coft scale Qδτ . The remaining dependence on the jet scale drops out once Z U is applied to the coft functions. The resummation of large logarithms works in an analogous way to the wide-angle case discussed earlier. In practice, it is most convenient to evolve the coft functions from the coft scale µ t ≈ Qδτ to the jet scale µ c ≈ Qδ by means of an equation analogous to (2.41). The evolution of the hard and soft functions to the jet scale is readily obtained by solving the standard evolution equations (2.45). The explicit computations presented below allow us to verify that the renormalization of the operators indeed works as we have described. In the next section, we demonstrate at two-loop order that (2.48) renders the coft function U 1 finite. In section 5, we verify at one-loop order that the matrix Z H lm , which renormalizes the Wilson coefficients H l , indeed renders all soft functions S m in (2.37) finite. In appendix C, we perform the same check also for the coft functions U m and relate the one-loop Z-factors for the two cases. Also in section 5 we show that at large N c the RG equation (2.39) is equivalent to the BMS equation at the leading logarithmic level.

Two-loop result for the narrow-jet cross section
In the following two sections we derive all ingredients entering the factorization theorems (2.15) and (2.34) at NNLO and verify that we reproduce the correct cross sections at this order. While the cross sections are finite after coupling renormalization, the different component functions in the factorization theorems contain divergences, which we regularize dimensionally. The 1/ǫ n poles cancel when the various components are combined, and hence we can test the factorization theorems working with bare functions. We start with the analysis of the narrow-jet cross section, for which our NNLO expression can be compared with numerical results obtained using the event generator Event2 [44]. The wide-angle case will be studied in section 4.

Ingredients of the factorization theorem
where we have used that at lowest order U m = 1 + O(α s ). Recall that the symbol ⊗ is defined as in (2.16) as an integral over the particle directions with measure Π i dΩ(n i )/(4π). For each parton in the jet, we introduce a polar angle θ i with respect to the jet axis n and an azimuthal angle φ i in the plane orthogonal to n. Rotational invariance in the transverse plane allows us to choose one of the azimuthal angles at will, and we choose to set φ 1 = 0 for the first parton. For the small-angle case it will be convenient to rescale the polar angles to new variablesθ i ≡ θ i /(2δ) ∈ [0, 1] (see below). We will thus parameterize the directions {n 1 , . . . , n m } with the angular variables {θ 1 , . . . ,θ m , φ 2 , . . . , φ m } and introduce new jet and soft functions (with the same names but different variables) such that Note that for m > 2 there are additional azimuthal angles in d = 4 dimensions, see (3.12) below.
The function H(Q, ǫ) is the standard hard function for two-jet processes and can be obtained from a matching calculation of the QCD vector current onto the corresponding SCET current operator at time-like momentum transfer. The renormalized two-loop expression can be found, for example, in [49,50]. The bare function is identical to the two-loop on-shell quark form factor, which has the form [52][53][54] The Laurent series in ǫ of the coefficients h F , h 2F , h A and h f to the required order is given in appendix A. Note that we expand the bare functions in the bare coupling constant α 0 = Z α α s , with where α s ≡ α s (µ) is the renormalized coupling constant. We write the d-dimensional coupling as α 0μ 2ǫ to make α 0 dimensionless, where the scaleμ 2 = µ 2 e γ E /(4π) is chosen such that the subtraction of divergences leads to renormalization in the MS scheme. The soft function S(Qτ, ǫ) defined in (2.32) is the same as the soft function for threshold resummation in Drell-Yan production, except that in our case the Wilson lines are outgoing instead of incoming. Since the soft function only depends on the product of the two lightlight vectors, its result is unchanged under the simultaneous change n → −n andn → −n (see [55] for a detailed discussion of relations among time-like and space-like soft functions). We can therefore simply use the Drell-Yan soft function at NNLO obtained in [49,56]. It reads The coefficients W i can again be found in appendix A. We now turn to the new ingredients in the factorization theorem, namely the coft and jet functions. For convenience we write their perturbative expansions in the form Coft functions. In (3.1) we need the coft functions U 1 to two-loop order and U 2 with one-loop accuracy. For the momentum scaling of coft particles in (2.18), the phase-space constraint allows for emission both inside and outside the jet cones. The energy is only constrained for emissions outside of the jet, because the coft momentum inside the jet is negligible compared to the momentum of the collinear particles. It is therefore dropped in the multipole expansion of the energy-conservation δ-function. Because of this fact, coft functions with all particles inside the jet are scaleless (their energy can be arbitrarily large). Also, a coft particle in the right-moving jet does not see the left-moving jet, since the out-of-left-jet condition is always fulfilled once the multipole expansion is performed. According to the definition (2.33) the coft function U 1 contains two Wilson lines, one along the direction n 1 of the particle inside the right jet and a second one along thē n direction, which describes emissions from the left jet. Similarly, the coft function U 2 contains three Wilson lines, two along the direction of the particles inside the right jet and a third one along then direction. We first discuss the calculation of a general coft function U m at one-loop order. The relevant Feynman diagrams contributing for the special case m = 2 are shown in figure 4. Analogous diagrams can be drawn for U 1 and for all higher coft functions U m . The general one-loop expression for (2.33) involves a sum over all pairs of emissions and absorptions from directions i and j, such that (ij)

JHEP11(2016)019
where (ij) with i = j denotes an unordered pair of numbers in the range 0 . . . m, and the scaleμ is defined after (3.4). Here, T 0 and n 0 =n are the color generators and the direction of the single Wilson line in the left hemisphere. Since the contribution from radiation inside the jet cone is scaleless, we have restricted the emission to lie outside the cone. For the special cases m = 1, 2 the sum over pairs yields where we used color conservation (i.e., the fact that 3 i=1 T a i = 0) to express the product of color generators T i · T j ≡ a T a i T a j through Casimir invariants. The loop integrals can be expressed in terms of the invariant scalar products of the external light-like reference vectors n,n, n 1 and n 2 . The result for the coft function U m must be invariant under the reparameterization transformation which leaves the product n ·n = 2 unchanged, as well as under arbitrary rescalings of the light-like vectors n µ i . It follows that the answer can be expressed in terms of the invariant variablesθ i = 1 δ n · n ī n · n i , Φ ij = 2 δ 2 n i · n j n · n in · n j , βδ . (3.10) The product βδ enters the definition of the coft scale. The variablesθ i are closely related to the polar angles θ i of the quark and antiquark with respect to the thrust axis, namelŷ where the last step holds in the small-angle limit. The jet-cone constraint implies that 0 ≤θ i ≤ 1. In four dimensions, the variables Φ ij are related to the azimuthal angle Performing the loop integral in terms of these variables and applying the Laplace transformation (2.30), we obtain

JHEP11(2016)019
The Laurent expansions of the coefficient functions to the required order in ǫ read where (3.15) and in general f i (0) = 0 and g i (0, 0, π) = 0. We also note that From momentum conservation and the definition of the jet axis it follows that the direction of particle 1 (the quark) in the jet function J 1 must be aligned with the jet axis (i.e. n 1 = n), which enforcesθ 1 = 0. Likewise, the azimuthal angles of particles 1 (quark) and 2 (gluon) in the jet function J 2 must differ by 180 • , which enforces φ 2 = π for our choice φ 1 = 0. These kinematic relations are not encoded in the coft functions, but they will be enforced when we evaluate the convolutions of the coft functions with the jet functions. We can thus simplify the calculation of the coft functions by implementing these constraints from the beginning. It is for this reason that we only need the function g 0 (θ 1 ,θ 2 , π) given in (3.15), and in the expression for U 1 we can setθ 1 = 0, in which case f i (0) = 0 for all i and we end up with a simple expression in terms of ζ values. In order to obtain the required two-loop contribution to U 1 settingθ 1 = 0 means a significant simplification. We can then exploit the fact that, after a Lorentz boost along the jet axis, the coft function U 1 with two opposite Wilson lines in directions n andn can be mapped onto the hemisphere soft function with p L Xs and p R Xs the total momenta in the respective hemispheres. In our case the energy in one hemisphere (the one corresponding to the inside of the jet) can be arbitrarily large.

JHEP11(2016)019
The hemisphere soft function is known at two-loop accuracy [25]. In order to obtain the two-loop coft function from this result, we first map the plane used to separate the two hemispheres onto the right jet cone via a Lorentz boost along the thrust axis. Under this boost ω L → ω L /δ and ω R → ω R δ. In other words, computing the hemisphere soft function with ω L < βδQ and arbitrarily large ω R is the same as computing the coft function with energy βQ outside the cone, i.e.
The integration over ω R needs to be performed before renormalization, since it leads to additional singularities. Taking the Laplace transform as defined in (2.30), we then obtain where L t = ln Qδτ µ and the explicit expressions for the coefficients V i are again collected in appendix A. Figure 5. Sample Feynman diagrams contributing to the jet function J 1 at different orders in perturbation theory. Note that only a single propagator is cut.

JHEP11(2016)019
Feynman diagram shown in figure 6, integrated over the appropriate phase space. We find where p 1 = E 1 n 1 and p 2 = E 2 n 2 are the momenta of the quark and gluon, respectively. We now split the integration domain into two regions, such that: Region II : In region I the integration generates double poles, because the gluon can be collinear to the quark and it can become soft (E 2 → 0), while in region II only single poles corresponding to a collinear divergence arise. These divergences are overlapping and must be separated using suitable choices of variables. We find it convenient to keep the larger of the two anglesθ max as one variable and use their ratio y =θ min /θ max as a second one. Both run from 0 to 1. In terms of these variables we obtain with L c = ln Qδ µ and C ǫ = e ǫγ /Γ(1 − ǫ). The 1/ǫ n singularities arise from the regions y → 0 andθ max → 0 and are non-overlapping. One can thus easily expand the jet function a Laurent series using the standard formula We quote the resulting expression up to O(ǫ 0 ), reinterpreted as a distribution on functions of the angular variablesθ 1 ,θ 2 and φ 2 . We find The differentials dy/dθ i in the last two lines change one of the integration variables from dθ i to dy; for example, when applied to a function F (θ 1 ,θ 2 ) the term in the third line gives The contributions of order ǫ and ǫ 2 can be obtained in an analogous way but are too lengthy to be presented here.
Convolutions. We finally consider the convolution J 2 ⊗ U 2 in (3.1), which we need to evaluate to O(α 2 s ). Performing the convolution with the tree-level coft function U (0) 2 = 1, and adding up the contributions from the two regions, we find

JHEP11(2016)019
Next we need to evaluate the convolution with the one-loop coft function. Doing so, we obtain the NNLO collinear-coft mixing contribution This result was presented earlier in the supplemental material to [38], but the finite terms were given only in numerical form.
Higher-order jet functions. The last unknown ingredients in the factorization formula (3.1) involve the one-loop corrections to the 2-particle jet function as well as the 3-particle jet functions with parton content qgg and qq ′ q ′ (summed over flavors q ′ ). Their combined contribution to the cross section can schematically be written as Some sample Feynman diagrams for these contributions are shown in figure 6. We have not computed these contributions individually but have inferred their divergent parts from the finiteness of the cross section. The explicit result is given in appendix A.

NNLO cross section
We now have all the ingredients at hand to obtain the full NNLO result for the cone-jet cross section. The bare functions need to be combined according to the NNLO expansion (3.1) of the factorization formula (2.34). After coupling renormalization all divergences cancel and we get a finite result for the Laplace-transformed cross section σ(τ, δ). Despite the fact that we have not explicitly computed the two-loop divergences of the jet functions, this provides a highly nontrivial check of the factorization formula (2.34), since the individual two-loop ingredients all depend on different scales. After expanding in ǫ, the divergences then involve logarithms of the different scales, which must cancel in the cross section. We stress that one would not obtain a finite result starting from the "standard" factorization formula (1.2) involving only two soft Wilson lines. Beyond one-loop order the nontrivial Wilson-line structure in (2.34) becomes an essential feature.

JHEP11(2016)019
Up to the desired order, the Laplace-transformed cross section is a quadratic polynomial in ln τ . For such a function, the Laplace transformation (2.30) can be inverted by means of the simple substitutions ln τ → ln β , ln 2 τ → ln 2 β − π 2 6 . (3.33) We choose µ = Q for the renormalization scale of the strong coupling and write  and c J,f 2 in expression (A.5) for the higher-order jet functions (3.32). In the upper panels in figure 7 we compare our analytic expressions for the three color structures in B(β, δ) without the contributions from the unknown constants to high-precision numerical results obtained using the event generator Event2 [44]. We choose very small values for δ and β, so that the logarithmic terms become dominant and the power corrections in δ and β can be ignored. At the scale of these plots we find perfect agreement. In the lower panels we show the differences ∆B between the numerical results from Event2 and our predictions, where the vertical scale is now expanded by a factor 1000. For sufficiently small δ and β these differences should be equal to the unknown constants c F,A,f The uncertainty on the last constant is fairly large due to numerical instabilities in the four-fermion channel. The same numerical problems were observed for thrust [57], where the analytic result for the constant is known [30,58].
The above fixed-order expression for the cross section makes it clear that the cross section depends in a most complicated way on ratios of the four relevant and hierarchical physical scales Q ≫ Qδ ≫ Qβ ≫ Qδβ. Our factorization formula separates these scale and factorizes this complicated dependence to all orders in perturbation theory.

Renormalization at NNLO
We now perform the renormalization of the jet and coft functions, which are the new ingredients in our factorization formula (2.34), and verify the properties we have discussed earlier in section 2.3. The jet function J 1 in (3.21) is trivial and does not require any renormalization. It follows that selves. According to the same relation the jet function J 2 in (3.25) requires a subtraction proportional to J 1 . We obtain 38) where the dots refer to higher-order terms in α s . The renormalization factor contains logarithmic dependence on the collinear scale Qδ, as is typical for Sudakov problems. The renormalized jet function is now obtained as (recall that y =θ min /θ max ) As we discussed earlier, the product of the jet, hard and soft renormalization factors must renormalize the coft functions, see (2.48). The factors Z H and Z S can be reconstructed from the RG evolution equations (2.45) or read off from the bare results for the corresponding functions, which are listed in appendix A. To the order we are working, the renormalization of the coft function takes the form where we have used the fact that U 3 = 1 + O(α s ). For simplicity we only indicate the arguments µ and ǫ in order to differentiate between renormalized and bare functions. The off-diagonal contributions depend on the directions of the additional partons, and the symbol⊗ indicates that one has to integrate over these. In the first relation above we do not know the two-loop contributions to the renormalization factors Z J 12 and Z J 13 which enter in (3.40). Similarly, in the second relation we do not know the one-loop contribution to the renormalization factors Z J 23 . The corresponding contribution to (3.41) is a series of 1/ǫ n pole terms, whose coefficients depend on the collinear logarithm L c . It is a nontrivial check of our method that in the remaining pole terms all reference to other logarithms cancels out. After some algebra, we find (atθ 1 = 0) To compute the coefficient u A 1 we have made use of relations (3.16). Note that the renormalized coft function only depends on the single coft scale, via L t = ln Qδτ µ , and the two-loop coefficient u F 1 is one half of the square of the one-loop result, as required by non-abelian exponentiation. The one-loop expression for the renormalized coft function U 2 reads

Two-loop logarithmic structure for wide-angle jet processes
In this section, we calculate all two-loop ingredients of the cross section for producing two wide-angle jets. This allows us to fully predict all logarithmically enhanced pieces in the cross section also for δ ∼ 1. We then verify that the relations (2.31) and (2.28) among the ingredients in the wide-and narrow-jet cases are indeed fulfilled. Up to two-loop order, the factorized wide-angle cross section (2.15) is given by For convenience, we factor out the Born cross section, so that H Inserting these expansions into (4.1), we find for the NLO coefficient in the cross section (3.34) and for the NNLO coefficient we obtain In the following, we first evaluate the hard function H 3 and then describe in detail the calculation of the soft functions S 2 and S 3 . The finiteness of the cross section can be used to infer the higher-order unknown logarithmic terms in the contribution from the hard functions H (2) 3 and H (2) 4 . As in the narrow-angle case, we compare the resulting oneand two-loop coefficients A(β, δ) and B(β, δ) with the numerical results obtained using the event generator Event2. Finally, we study all the two-loop ingredients in the small-δ limit and verify the factorization formulas (2.28) and (2.31) at this order.

Hard function
Following the operator definition in (2.14), the hard function H 3 describing the process γ * → q(p 1 )q(p 2 ) g(p 3 ) starts at O(α s ) and is given by The integrations over the energies are performed keeping the directions of the three particles fixed. The crucial point is that for three-particle final states, the thrust axis always points opposite to the direction of the most energetic particle, and the jet cone therefore centers JHEP11(2016)019 Figure 8. Three-body phase space for the γ * (q) → q(p 1 )q(p 2 ) g(p 3 ) process in different kinematic regions, corresponding to different energy hierarchies among the three partons. The gray region represents the part of phase-space in which the jet-angle constraint is fulfilled.
around it. For this reason, it is natural to decompose the phase-space integration into the following regions: where we parameterize the particle energies as x i = 2E i /Q. The corresponding threebody phase space is shown in figure 8. In figure 9 the kinematical configurations in the different regions are illustrated. Regions I and VI suffer from overlapping soft and collinear divergences, while regions II and V contain collinear divergences only. Regions III and IV are infrared safe. Due to the x 1 ↔ x 2 symmetry, only the contributions of regions I, II and III need to be computed. It is obvious from figure 9 that the phase-space constraints for H 3 take a relatively complicated form. To perform the relevant integrations, it is convenient to use the energy and angle of the least energetic particle as variables for the hard function. In region I, the gluon with momentum p 3 is the particle with the lowest energy. We parameterize its energy fraction x 3 and angle θ 3 with respect to the thrust axis as (4.8) The variables u and v are integrated from 0 to 1 and disentangle overlapping soft and collinear divergences. Using momentum conservation one can express x 3 as a function of the JHEP11(2016)019 angles θ 2 and θ 3 , so that the variables u and v play the role of the angular variablesθ i used in the previous section. Indeed, in the small-angle limit the variable u is directly related to the ratio of the angles introduced for the jet function: y =θ min /θ max = u/(2 − u). Analogous parameterizations are introduced in regions II and III; see appendix B for more details. The nontrivial part of the angular integration is then written as an integral over u and v The one-loop hard functions H 3 in the different sectors are given by  S 2 containing two Wilson lines. In momentum space, the bare soft function is given by

(4.12)
Anticipating the convolution with the hard function we have set the reference vectors n 1 and n 2 of the Wilson lines equal to n andn. If the momentum k is clustered into a jet the resulting integral is scaleless. Therefore, the one-loop function only receives a nonvanishing contribution from the out-of-jet region. After convolution with the LO hard function, the result reads with L s = ln Qβ µ . The full expression for this function, without expanding in ǫ, can be found in [32]. For the two-loop soft function S (2) 2 three color structures arise. The color-averaged expression has the form (4.14) The C 2 F term in this expression can be derived from the one-loop result by means of nonabelian exponentiation, while the C F C A and C F T F n f terms can be extracted from the results given in [32]. In appendix B.2 we describe the necessary steps in detail and list the explicit two-loop expressions.
At NNLO we also need the convolution of the one-loop soft function S 3 with the oneloop hard function H 3 . First we calculate the expression S T i · T j dΩ(n k ) 4π n i · n j n i · n k n k · n j Θ nn out (n k ) , 3 for three different cone sizes α.
with i, j ∈ {1, 2, 3}. The divergence arises from the energy integral. The angular integral is finite, since the gluon is always outside the cone and can therefore never become collinear to a Wilson-line vector whose direction tracks a hard parton inside the jet. In region I the hard function H 3 suffers from double divergences. Therefore, in order to obtain all divergent terms after the convolution, we need to calculate S Performing the convolution of the one-loop hard and soft functions, we obtain We have computed the singly-divergent terms in numerical form. In table 1 their values for some specific cone sizes are listed. Finite terms do not give rise to large logarithms in the cross section and are thus omitted. Now we have all two-loop ingredients at hand, except for the purely hard contributions given by the two-loop hard functions H (2) 3 and H (2) 4 . Finiteness of the cross section implies that we can obtain the divergent terms of the sum of these two-loop hard functions. These terms only depend on the hard scale Q, which provides a nontrivial check of our computation of the various other ingredients. Explicitly, they must be of the form where the two-loop coefficients are functions of δ, but independent of β. Our explicit results, which are presented in appendix B.4, satisfy this condition.

NNLO cross section
After combining all ingredients and performing charge renormalization, we obtain a finite NNLO cross section. Similar to the narrow-angle case, the coefficients A and B in (3.34) can also be computed numerically using the event generator Event2 [44], and a comparison with these results provides a check of the factorization formula at this order.
(4.19) Since we now count δ ∼ 1, this result holds for arbitrary values of δ, up to terms suppressed by powers of β. In figure 11, we compare our analytical result for A(β, δ) (red line) to the numerical results obtained using Event2 (blue dotted). As it must be, the difference ∆A between the logarithmic terms and the full result goes to zero at small values of ln β within the numerical uncertainty of the Monte-Carlo integration.

The small-δ limit
As a final check, we now evaluate all two-loop bare ingredients in the small-δ limit and verify that they fulfill the factorization formulas (2.28) and (2.31). The hard function H 2 is independent of the jet opening angle and the factorization (2.28) becomes trivial. Interesting relations first arise for the hard function H 3 . In the small-δ limit, the contributions from regions III and IV are power suppressed. For the remaining regions, we indeed find    [44]. In the lower panels we show the difference ∆B between Event2 and our result, which should be equal for small values of β. The cone size is chosen as α = π/4, corresponding to δ ≈ 0.414.
This factorization property of the cone-jet soft function was observed also in [39]. Similarly, for S 3 we find at one-loop accuracy S 3, I(II) (Qτ, δ) = S(Qτ ) U 1 (Qδτ ) U 2, I(II) (Qδτ ) (4.24) in each region. Since all relevant relations are fulfilled, we indeed reproduce the corresponding NNLO cross section (3.34) in the small-δ limit. We have also verified that in this limit the one-loop coefficient (4.19) coincides with the one in (3.34), and that all the logarithms ln β in the two-loop coefficient (4.20) are the same as those in the narrow-angle coefficient (3.35).

Renormalization-group evolution
In this section, we focus on the wide-angle case and check that renormalization works at the one-loop level as described in section 2.3. To this end, we determine Z H ij , the Z-factor for the hard function, and verify that the same matrix renders the soft function S m with m Wilson lines finite. We then give results for the one-loop anomalous dimensions and show that the lowest-order RG evolution equation is equivalent to the BMS equation.

Renormalization at one-loop order
Let us write the expansion of the Z-factor defined in (2.35)

in the form
with z (0) i,j = 1. The entries z i,j are matrices in the color space of the partons in the amplitude and its conjugate. We denote the color generators T a i acting on i-th particle in the amplitude on the left-hand-side of H m in (2.14) as T a i,L , and those acting on the conjugate amplitude on the right-hand side as T a i,R . Because of the structure of (2.15), the roles of T a i,L and T a i,R are reversed for the case of the soft function: the generators T a i,L act on the right-hand side of S m .
Let us now verify that Z H , which is introduced to absorb the divergences of the hard function, can indeed be used to renormalize the one-loop soft function. If this is true, we must find that where we have used S m = 1 + O(α s ), so that the Z-factors multiply the identity matrix.
In the second term we integrate over the angle of the additional emission. One can easily obtain the divergent part of the one-loop soft functions, since it is given by a sum of exchanges between two legs. A sample Feynman diagram is shown in figure 10. We get where we have introduced the dipole radiator The function Θ nn out (n k ) = 1 − Θ nn in (n k ) ensures that the gluon is outside the two jet cones around the n andn directions. Note that the angular integral does not suffer from collinear divergences, since the vectors n i and n j lie inside the jet cones, while the direction n k associated with the soft emission points outside the cone. (The soft radiation can also be JHEP11(2016)019 emitted inside the cone, but as mentioned earlier this contribution is scaleless, since it does not have an upper limit on the energy of the emission.) In (5.3), the quantity z m,m represents the divergences of the virtual corrections to the amplitude with m legs, while z m,m+1 gives the divergences from an additional real emission. Let us now consider the real and virtual corrections together, since all collinear divergences drop out and only a single soft divergence remains. The leading divergence can be obtained by using the soft approximation for the emitted (real or virtual) gluon. In the soft approximation, the real-emission contribution factorizes as In this approximation, one can write the virtual correction in the same form as the realemission contribution, because the principal-value part of the propagator of the emission does not contribute. The virtual correction then reads (5.7) The virtual gluon can be either inside or outside the cone. Since the gluon is soft, we can replace Q − E k → Q in the hard function. Adding up the two terms and extracting the divergence from the lower end of the energy integration, the total result for the divergent part becomes Since the color factors are contracted with the trivial tree-level soft function, we do not need to distinguish the left and right color generators. Note that inside the cone the real and virtual corrections have cancelled, so that the net result only gets contributions from out-of-cone radiation and precisely cancels against the divergence of the soft function. We see that the renormalization indeed works at the one-loop level. We have repeated the same exercise also for the narrow-jet case, see appendix C. In this case, we can give explicit expressions for the angular integrals. Again, we find that the divergences cancel as they should.

Renormalization-group evolution at leading logarithmic level
We now discuss the anomalous-dimension matrix Γ H defined in (2.40), which governs the RG evolution of the hard (2.38) and soft functions (2.39), and verify the agreement between the perturbative expansion of the BMS equation and our RG-based resummation method. In order to resum the leading logarithmic terms, the anomalous-dimension matrix is needed up to O(α s ). It can be expressed as

JHEP11(2016)019
where It follows from the discussion in the previous section that, in the soft approximation, the corresponding matrix elements are given by The anomalous dimensions V m and R m depend on the directions {n} = {n 1 , . . . , n m } and colors of the hard partons, and the indices i, j in the sum run from 1 to m. The quantities R m also depend on the additional direction n m+1 of the real emission. The integration over this direction is performed after the multiplication with the soft function. At first sight, the expressions in (5.11) look problematic, since the angular integrals involve collinear divergences. However, we know on general grounds that the collinear divergences must cancel when the anomalous-dimension matrix is applied to the soft functions, see section 2.3. We have observed this cancellation at the one-loop level in the previous section, and we will see explicitly in the following that the same pattern continues at higher orders. The expressions (5.11) are therefore valid for the RG evolution of the soft function. To obtain an expression which is suitable for evolving the hard functions, one would need to regularize the angular integrations and extract the collinear divergences. Furthermore, the soft approximation would then not be appropriate to obtain the anomalous-dimension matrix. At the soft scale µ s ≈ Qβ the soft functions do not involve large logarithms, and the higher-order corrections are thus suppressed by powers of α s . These can be neglected at leading logarithmic accuracy, so that the soft functions reduce to the identity matrix S m ({n}, Qβ, δ, µ s ) = 1 + . . . . (5.12) To get the resummed cross section at leading logarithmic accuracy, we evolve the soft functions to the hard scale µ h ≈ Q using (2.41). At this scale the hard functions do not involve any large logarithms, and hence all higher-order hard functions are suppressed, ). The only non-vanishing contribution arises from the lowest-order function The net effect of the convolution with H 2 is that the reference vectors n 1 and n 2 are set equal to the jet directions n andn, together with a multiplication by the Born cross section. We thus find that at leading logarithmic order the resummed cross section is equal to 14)

JHEP11(2016)019
where the symbol⊗, which was introduced in (2.37), indicates that one has to integrate over the additional directions present in the higher-multiplicity anomalous dimensions R m and V m . It is of course highly nontrivial to obtain an explicit form for the formal expression (2.42), but it is easy to write down explicit expressions order by order. To do so we rewrite the exponent of the evolution matrix (2.42) at leading order in RG-improved perturbation theory in the form which corresponds to leading logarithmic accuracy. By using the prefactor as our expansion parameter, we automatically include running coupling effects. Using the structure of (5.10) and expanding the exponential in (2.42), we find for the first three coefficients of the expansion in t S (1) As explained after (5.14), we have to integrate over the directions of the additional emissions but for brevity, suppress the integrations in the above expressions. Including the integration, the one-loop term reads The structure of the result (5.17) is very simple. To obtain the result at the next order, one takes the existing result and adds an additional real emission plus a virtual correction to each term. This type of iterative structure is reminiscent of a parton shower, and it should therefore be possible to solve the evolution equation numerically, using Monte Carlo methods. Indeed such Monte Carlo methods have been used to perform resummations of NGLs, see e.g. [24,59]. One nontrivial complication is that the anomalous dimensions are matrices in color space and that the color algebra becomes nontrivial for high multiplicities. This difficulty can be avoided by taking the large-N c limit, which is also useful to compare (5.17) to the result obtained using the BMS equation. To take the large-N c limit, it is simplest to adopt the trace basis (see for example [60]), i.e. to write the color structure of the m-particle amplitudes in the form  Figure 13. The action of the operator R m on an amplitude in the large-N c limit, where (i 3 , i 4 , · · · , i m ) is a permutation of {3, 4, · · · , m}. The double and single lines represent gluons and quarks, respectively. The sum on the right-hand side represents all the contributions from the planar diagrams.
with color-ordered amplitudes A σ ({p}). The indices a and b are the color indices of thē qq pair. We only include single trace terms, since contributions with multiple traces are suppressed at large N c . At large N c , emissions arise only between nearest-neighbour legs, since all other attachments would lead to non-planar contributions which are suppressed. Based on the above simplification, the effect of R m in the large-N c limit is shown diagrammatically in figure 13. The action of V m simplifies analogously, as shown in figure 14. The large-N c color factor from squaring the amplitudes is simply a factor of N c for each color loop, and the number of additional color loops is equal to the number of powers of α s , so that the color factor is obtained by switching to the 't Hooft coupling λ = N c α s . We now plug the explicit results (5.11) for the anomalous-dimension coefficients V m and R m into the expressions (5.17). For the coefficients of the expansion in t, we then obtain where Ω 3 Out = dΩ(n 3 ) 4π Θ nn out (n 3 ), and we have used the abbreviation The above expressions include all leading logarithms, i.e. the global and non-global logarithmic terms appear together. Let us now relate the above expressions to the leading logarithmic resummation of NGLs at large N c , which can be obtained by solving the BMS equation [26] ∂LG kl (L) = d Ω(n j ) 4π  with boundary condition G kl (0) = 1. The function G kl (L) depends on two light-like reference vectors n k and n l . After solving the equation, the resummed soft function is obtained as S({n}, Qβ, µ) = G 12 (L) withL = 4N c t. While the non-linear integral equation (5.22) can in general only be solved numerically, it is easy to solve it iteratively, order by order inL. This was done in [25], where the resulting non-global terms were listed up to three-loop order. Our expressions (5.20) agree with these results. This demonstrates the equivalence of the RG equation (2.39) driven by the one-loop anomalous dimensions (5.10) and the BMS equation in the large-N c limit. However, an important advantage of our RG framework is that it is valid for arbitrary N c and logarithmic accuracy. We have verified that the two-loop expression S 2 indeed reproduces the leading logarithmic term in the NNLO result (4.20). In the literature one often distinguishes global from non-global logarithms. The global part is defined as the exponentiated one-loop logarithm, where the parameter t was defined in (5.16). After removing the global piece, we find that the leading non-global terms in (4.20) have the form This agrees with the result of [61], which analyzed non-global logarithms in the presence of rapidity gaps, after one relates the rapidity gap ∆η to our cone parameter using δ = exp (−∆η/2). From our point of view, the separation into global and non-global logarithms is somewhat artificial. Beyond leading logarithmic accuracy it is also not unique: in general, the non-global logarithms are defined by dividing the cross section through the cross section of a global observable, which has the same leading logarithms. Different reference observables have different subleading logarithms, so that the non-global subleading logarithms depend on the choice of the global observable. The RG evolution approach developed in this work can be compared to the functional RG approach by Simon Caron-Huot [37]. In this formalism, unitary matrices U (n) in color space, which are functions of a direction vector n, are used to track the contributions from JHEP11(2016)019 different particle multiplicities. Using these matrices, the author defines the color-density matrix functional σ[U ], which is given by the cross section modified by multiplying finalstate partons along directions n by the matrix U (n). Taking functional derivatives with respect to these matrices then separates out the contributions to the cross section from particles traveling along the corresponding directions. The one-loop expressions for the anomalous dimensions governing the RG evolution of the color-density matrix defined in this reference are in one-to-one correspondence to the anomalous dimensions in (5.11). At leading logarithmic order, the resummed result is obtained from solving the corresponding RG evolution equation for trivial initial conditions and is the same as in our formalism. Beyond leading logarithmic accuracy, the relation between the two formalisms is less obvious. The approach [37] does not distinguish hard and soft partons but multiplies every parton by a matrix U (n), and the hard function contains arbitrary hard radiation outside the jet. Furthermore, the Wilson-line structure, which is an important feature of our approach, is not immediately manifest. A one-to-one mapping should arise after an averaging procedure over U (n), which sets those matrices corresponding to outside radiation to zero, effectively vetoing hard outside radiation, and multiplies the inside partons by soft Wilson lines [62]. The details of this procedure were not discussed in [37], which was mainly concerned with the computation of the two-loop anomalous dimension. Recently the result for the anomalous dimension was even extended to three-loop order for planar N = 4 Super-Yang-Mills theory [63]. Translated into our formalism, these results should provide an important ingredient for the resummation of subleading non-global logarithms.
Let us also compare our method in more detail with the approach put forward by Larkoski, Moult and Neill in [35]. Rather than attempting to directly arrive at a facorization theorem for a non-global observable, these authors imagine performing a set of increasingly differential measurements on the jet. The goal of these measurements is to isolate the regions of phase space where the jet contains soft subjets near the jet boundary, which then give rise to NGLs. The basic strategy is that resolving the substructure of the jet down to lower scales reduces the size of the NGLs. The resummation of the global logarithms associated with the subjet observables then resums part of the NGLs present in the original, inclusive cross section. However, at each stage non-factorized logarithms remain. A complete factorization theorem for a jet cross section with additional measurements would be at least as complicated as the factorization theorem we obtain for the inclusive cross section, and would include multi-Wilson-line operators. By construction, the method of [35] involves increasingly differential measurements, each of which requires a dedicated effective theory containing the relevant modes for the associated factorization theorem. The method thus involves a tower of effective theories with more and more degrees of freedom. This is rather different from our approach, which factorizes the original cross section and allows for the resummation of the NGLs with RG methods. Rather than with a tower of effective theories, we work with a given theory with a fixed number of degrees of freedom. The wide-angle case is especially simple in this regard, since it only involves a single, soft low-energy mode. Instead of soft subjets, we encounter soft Wilson lines along the directions of the energetic particles, and the NGLs get factorized into hard logarithms from vetoing hard emissions outside the jet, and soft logarithms due to the fact JHEP11(2016)019 that soft radiation is not constrained inside the jet. The physical picture emerging from our approach is rather different from that put forward in [35]: we find that the NGLs arise from "subjets" which are generically hard. The leading logarithms can be obtained by using strongly-ordered amplitudes in which all gluons are soft, but this is only a feature of this particular approximation.
Solving the RG evolution equation (2.41) to obtain the evolution between the hard and soft scales is obviously complicated. This is not surprising, given that even in the large-N c limit and at leading logarithmic order our equation is equivalent to the non-linear integral equation derived by BMS. However, from a conceptual point of view our RG resummation method is completely standard, and it is clear which ingredients are required for a given logarithmic accuracy. The expansion in subjets, on the other hand, is not immediately related to logarithmic accuracy. Even for leading logarithmic accuracy, terms with arbitrary many subjets contribute. The authors of [35] argue that the higher-order terms in the expansion are phase-space suppressed, but it is not clear to us if there is any parametric suppression associated with this assertion.

Summary and outlook
In this paper, we have derived all-order factorization theorems for cone-jet cross sections, both in the wide-angle case and in the narrow-jet limit. With a veto on out-of-jet energy, these cross sections suffer from large logarithms of the ratio of the energy outside and inside the jets and, for narrow jets, from large logarithms of the opening angle of the jet cone. Our factorization theorems separate the physics associated with the different energy scales relevant in these processes and provide operator expressions for the various ingredients of the cross section. We believe that this is an important step forward, since it opens the door for higher-logarithmic resummations for such observables. In the context of the effective field-theory framework we are using, this resummation is achieved by solving the RG evolution equations of the relevant operators.
The low-energy structure we encounter in jet cross sections is significantly more complicated than what is present for event-shape variables such as thrust. In the case of global dijet event shapes, soft emissions are obtained from an operator consisting of two Wilson lines along the jet directions. For cone-jet cross sections, we instead encounter a separate Wilson line for each energetic particle inside the jet. The presence of this multi-Wilsonline structure is quite intuitive for wide-angle jets, where the characteristic angle between the energetic particles is parametrically large and of the same size as the typical angle of soft emissions. The emissions thus resolve the individual energetic partons, and the soft Wilson lines for the particles inside a jet cannot be combined into a single Wilson line. Interestingly, we find that the same structure also persists for narrow jets, because it turns out that narrow-angle soft radiation gives an important contribution to the cross section. The effective theory in the narrow-jet limit therefore contains a an additional mode, which is simultaneously soft and collinear. Boosted soft modes have appeared in other contexts in SCET; the new mode we find is different, because every component of its momentum is parametrically smaller than the corresponding one of the energetic, collinear particles.

JHEP11(2016)019
It is this basic property which leads to the multi-Wilson structure not present in earlier SCET applications.
Starting from our factorization theorems, we have computed the full logarithmic structure of the NNLO jet cross sections both in the wide-angle and narrow-jet cases and have verified that they agree with numerical results obtained from the fixed-order event generator Event2. We have also shown how the wide-angle cross section can be refactorized in the narrow-angle limit, and we have verified the resulting relations at NNLO. These results demonstrate that our approach indeed captures all logarithms present in these cross sections.
Since our factorization theorem involves operators with an arbitrarily large number of Wilson lines, the relevant RG evolution equation involves an infinite matrix, and it will be quite challenging to solve it. This complexity is perhaps not that surprising, given that jet cross sections are prime examples of non-global observables, and it is well known that even the leading logarithms in such processes have quite a complicated structure, which is captured by the BMS equation. Indeed, we have shown that our RG evolution equation is equivalent to this non-linear integral equation at the leading logarithmic level. We expect that a natural approach to solving our RG equation will be based on using Monte Carlo methods, since its structure is reminiscent of a parton shower. Also, for many practical applications it could be sufficient to predict the first few NGLs instead of resumming the entire series. In any case, as a next step, it will be important to develop methods to solve the RG equations and to obtain the relevant anomalous dimensions needed to resum subleading NGLs.
Jet observables are the most important observables at colliders and there are countless applications for our framework once the methods for solving the associated RG equations have been developed. In particular, it would be interesting to use RG methods to systematically improve jet substructure predictions, which are currently mostly based on parton-shower programs. In our paper, we have analyzed the simplest example of a jet cross section, but it is clear that the basic multi-Wilson-line structure is present in all nonglobal observables. Nevertheless, there can be additional difficulties when trying to apply our approach to multi-jet cross sections. For example, it was stressed in [64,65] that the complexity of the phase-space constraints for recombination jet algorithms might present a stumbling block when trying to derive all-order statements for such cross sections. While it is interesting to extend the resummation to more complicated observables, it is equally important to find observables whose structure is simple enough to allow for systematic higher-order resummations. For example, in the case we analyzed, it turned out to be advantageous to choose the thrust axis as the jet axis rather than choosing the axis so as to maximize the energy inside the jet.
For simplicity, we have studied dijet cone cross sections in e + e − collisions, but it will obviously be important to apply the same methods to hadronic jet cross sections. In fact, since there are no detectors at very large rapidities, every hadron collider cross section contains two narrow cone jets along the beam direction. Their angular size is given by δ = e −ymax , where y max is the maximum rapidity up to which particles are measured. If a jet veto is imposed in the measured region, then the presence of these beam jets JHEP11(2016)019 will induce logarithms of the same type as the ones we have studied in this paper. So far, these logarithms were only studied using parton showers [66] (a method to mitigate their effects was proposed in [67]), and it will be interesting to analyze them within our framework. An interesting complication in the case of hadronic collisions is the presence of Glauber gluons. An comprehensive effective-theory analysis of these has been presented recently [68], following earlier exploratory studies [69][70][71][72][73]. Glauber effects are interesting in the context of non-global observables, because their presence can lead to enhanced powers of NGLs, the so-called super-leading logarithms [74,75]. We look forward to studying these issues and applying our framework to many phenomenologically relevant problems.
The real-virtual corrections S R−V 2 can be written in terms of the one-loop soft function S The remaining two contributions can be extracted from [32], where the two-loop soft function for jet thrust was computed. This quantity is defined as S R (k L , k R , λ) = 1 N c Xs 0| Sn S n M(k L , k R , λ) |X s X s | S † n S † n |0 , (B.10)

JHEP11(2016)019
where the measurement function is given by M(k L , k R , λ) |X s = δ(k L −n · p L Xs ) δ(k R − n · p R Xs ) δ(λ − E out Xs ) |X s . (B.11) The soft radiation inside the right (left) jet has momentum P R Xs (P L Xs ), and E out Xs represents the total out-of-jet energy. Therefore, after integration over k L , k R and λ our two-loop function is obtained. Explicitly, we find is given in [32]. When performing the integral for S In−Out 2 , we need to split the k L(R) integration at k L(R) = Qβ. Combining all contributions, we obtain for the two-loop coefficients defined in (4.14)