Toward Multi-Differential Cross Sections: Measuring Two Angularities on a Single Jet

The analytic study of differential cross sections in QCD has typically focused on individual observables, such as mass or thrust, to great success. Here, we present a first study of double differential jet cross sections considering two recoil-free angularities measured on a single jet. By analyzing the phase space defined by the two angularities and using methods from soft-collinear effective theory, we prove that the double differential cross section factorizes at the boundaries of the phase space. We also show that the cross section in the bulk of the phase space cannot be factorized using only soft and collinear modes, excluding the possibility of a global factorization theorem in soft-collinear effective theory. Nevertheless, we are able to define a simple interpolation procedure that smoothly connects the factorization theorem at one boundary to the other. We present an explicit example of this at next-to-leading logarithmic accuracy and show that the interpolation is unique up to $\alpha_s^4$ order in the exponent of the cross section, under reasonable assumptions. This is evidence that the interpolation is sufficiently robust to account for all logarithms in the bulk of phase space to the accuracy of the boundary factorization theorem. We compare our analytic calculation of the double differential cross section to Monte Carlo simulation and find qualitative agreement. Because our arguments rely on general structures of the phase space, we expect that much of our analysis would be relevant for the study of phenomenologically well-motivated observables, such as $N$-subjettiness, energy correlation functions, and planar flow.

Abstract: The analytic study of differential cross sections in QCD has typically focused on individual observables, such as mass or thrust, to great success. Here, we present a first study of double differential jet cross sections considering two recoil-free angularities measured on a single jet. By analyzing the phase space defined by the two angularities and using methods from soft-collinear effective theory, we prove that the double differential cross section factorizes at the boundaries of the phase space. We also show that the cross section in the bulk of the phase space cannot be factorized using only soft and collinear modes, excluding the possibility of a global factorization theorem in soft-collinear effective theory. Nevertheless, we are able to define a simple interpolation procedure that smoothly connects the factorization theorem at one boundary to the other. We present an explicit example of this at next-to-leading logarithmic accuracy and show that the interpolation is unique up to α 4 s order in the exponent of the cross section, under reasonable assumptions. This is evidence that the interpolation is sufficiently robust to account for all logarithms in the bulk of phase space to the accuracy of the boundary factorization theorem. We compare our analytic calculation of the double differential cross section to Monte Carlo simulation and find qualitative agreement. Because our arguments rely on general structures of the phase space, we expect that much of our analysis would be relevant for the study of phenomenologically well-motivated observables, such as N -subjettiness, energy correlation functions, and planar flow.

Introduction
Historically, there has been significant effort devoted to understanding and computing the all-orders distributions of jet observables in QCD [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. This has led to incredibly precise predictions for the differential cross sections of these observables which have been used, for example, to determine the strong coupling α s to high precision [9,11,21]. For all their successes, though, this program can only answer questions about individual observables. In this paper, we move beyond this paradigm of single differential cross sections, to exploring the full phase space of multi-differential cross sections analyzed on a single jet. 1 Using angularities as a case study, we find that factorization methods are confined to the boundaries on the physical phase space regions, and we propose an interpolation method to connect all the various forms of factorization possible.
There exist multiple motivations for why one might want to study such multi-differential cross sections. Aside from purely formal interest in connecting different effective field theory regimes, we focus on two motivations here: for studying the correlations between different observables and for understanding the properties of observables formed from the ratio of two infrared and collinear (IRC) safe observables. Phenomenologically one would want to know the correlations between different observables so as to determine the extent to which they probe identical physics. However, this cannot be done by studying the differential cross sections of individual observables alone. Correlations are encoded in the multi-differential cross section of the observables and so to understand the correlations between two observables we must study their double differential cross section.
Studying the correlations between two observables is not necessary to make highly precise predictions for QCD. However, with the advent and boom of the jet substructure program [23][24][25] increasingly detailed questions about the dynamics of QCD jets are being asked and probed by experiment . In particular, one of the goals of jet substructure is to design highly efficient observables and procedures for discriminating QCD jets from boosted heavy particle decays. Many of the proposed procedures for doing so involve the measurement of several observables on the jet and making appropriate cuts. Thus, to determine if a QCD jet can fake looking like a boosted W , Z, H or top quark requires a thorough analysis of the correlations of the observables that go into the discrimination procedure.
Several of the most powerful discrimination observables are formed from the ratio of two IRC safe observables. This includes N -subjettiness [53,54], energy correlation functions [7,55], planar flow [56,57] and angular structure functions [58,59]. While it might seem like the ratio of two IRC safe observables is still IRC safe and so calculable order-by-order in perturbation theory, it was shown in Ref. [60] that ratio-type observables are actually IRC unsafe, if the denominator observable can become arbitrarily small. Naïvely, this is an insurmountable barrier to understanding these observables in perturbative QCD. Indeed, this is true with the standard procedure of computing single differential cross sections which require IRC safety to be well-defined in perturbation theory.
However, in Ref. [61] it was shown that ratio observables can actually be made welldefined, if all-orders effects are taken into account. There, the simple observable formed from the ratio of two angularities [6,10,57] measured on a single jet was studied, where the angularity e α is 2 (1.1) E J is the jet energy, R 0 is the jet radius, θ i is the angle between particle i and an appropriately defined jet axis, and α > 0 for IRC safety. The approximation is accurate for R 0 1, which we assume throughout this paper.
The differential cross section of the ratio r of two angularities e α and e β can be found by marginalizing the double differential cross section of the two angularities: (1.2) Ref. [61] showed that, while the ratio observable is not IRC safe and so cannot be computed order-by-order in α s , by resumming the large logarithms present in the double differential cross section to all orders, the differential cross section for the ratio r is well-defined and calculable. This property was called "Sudakov safety" because the calculability of the cross section of r relied on the fact that small values of the angularities e α and e β are exponentially suppressed by the Sudakov factor. The calculation of the double differential cross section of angularities was done to leading logarithmic (LL) accuracy in Ref. [61] with no robust predictions about what happens at higher logarithmic orders. In particular, Sudakov safety was only exhibited to LL accuracy, and some important and subtle physics might arise at higher orders that could change the story. Given these motivations, the double differential cross section of two angularities measured on a jet provides a laboratory for understanding multi-differential cross sections. To have adequate control over large logarithmic corrections, we need to prove a factorization theorem which would provide an order-by-order recipe for resumming to arbitrary accuracy. We will find that establishing such a factorization theorem for all of phase space in the double differential cross section is not possible with identified soft and collinear modes. In particular, a subtlety in the resummation of double differential cross sections is that the two measured observables do not define a unique set of scales for soft and collinear radiation in the jet.
Nevertheless, we will show that there do exist factorization theorems on the boundaries of phase space for the double differential cross section of two angularities using soft-collinear effective theory (SCET) [62][63][64][65][66]. Single differential cross sections factorize when the observable is sufficiently small, when one can say that soft and collinear dynamics dominate the structure of the jet. For the case of the double differential cross section of angularities e α and e β , small values of the angularities does mean that soft and collinear dynamics dominate the jet. However, the physical phase space for the double differential cross section is two dimensional, and the precise scaling of e α and e β with respect to one another emphasizes soft over collinear physics, or vise-versa. Strictly speaking, only on the boundaries of the phase space are the soft and collinear modes on-shell, where the factorization theorems hold. 3 The boundaries of phase space are defined by the requirements of energy conservation and clustering of emissions into the jet of radius R 0 . Energy conservation corresponds to the boundary 4 where e β α = e α β and the jet radius requirement is the boundary line e α = e β . The physical phase space lies in between. We will show that, at these boundaries, the double differential cross section of the angularities e α and e β reduces to the single differential cross section for one of the angularities times a δ-function for the other angularity, depending on the boundary, plus terms that integrate to 0. For example, near the boundary e β α = e α β , the double differential cross section reduces to where denotes the relationship that follows from the factorization theorem and f α + is a function that integrates to zero on e β ∈ [0, e β/α α ]. H represents the hard function, J(e α , e β ) is the double differential jet function and S(e α ) is the soft function for e α alone.
Importantly, this relationship captures the effect of canonical resummation on this boundary as predicted by the factorization theorem and the only non-trivial dependence on e β exists in the ratio e β /e β/α α . The fact that the soft function is independent of e α implies that the ultraviolet singular structure of the cross section exists on the line e β = 0, as enforced by δ(e β ). This is a non-trivial statement of the factorization theorem on this boundary that to all orders the RG evolution does not generate a non-zero value for e β .On the other boundary of phase space, where e α = e β , we find a similar relationship for the singular terms, with the single differential cross section of e β times δ(e α ): Non-trivial dependence on e α only exists in the ratio e α /e β . Because the factorization theorem only applies near the boundaries of the phase space, we cannot formally claim any logarithmic accuracy of the double differential cross section in  Figure 1: Summary of the results of the factorization theorem of the double differential cross section of angularities. The factorization theorems exists near the boundaries of the allowed phase space where the double differential cross section reduces to the appropriate single differential cross section plus terms that integrate to 0. The bulk of the phase space is described by an interpolating function.
the bulk of the phase space. However, we are able to determine a function that interpolates into the bulk of the phase space between the boundaries; crucial to this is the existence of factorization theorems at the boundaries. The interpolation between the boundary factorization theorems can be determined most simply by appropriately setting scales in the logarithms and by adding terms that are subleading at the boundaries. This conjectured double differential cross section must satisfy several consistency conditions, such as correctly reproducing the single differential cross section of one of the angularities. Thus, while we are unable to fully demonstrate formal logarithmic accuracy in all of the phase space, we will present a conjecture for the double differential cross section to next-to-leading-logarithmic accuracy (NLL) which satisfies all consistency conditions. 5 The summary of this factorization theorem discussion is illustrated in Fig. 1. The structure of the cross section as found from the interpolation procedure is fascinating and manifests the barrier to proving a factorization theorem in the bulk of the phase space. The interpolation procedure defines a double cumulative cross section containing the following logarithms: (1.5) log e α and log e 1/β β can naturally be understood as arising from soft and collinear divergences, respectively, and so correspond to the modes that are identified in SCET. The other logarithms, which we refer to as "k T ", 6 are novel, arising neither from soft nor collinear modes over all of the phase space of e α and e β . Indeed, the fact that there are three logarithmic structures in the bulk of the phase space suggests that there must be three modes in a factorization theorem of the double differential cross section that would be valid everywhere. 7 At the boundaries of phase space, the k T logarithms degenerate to soft or collinear logarithms, which is why SCET factorization applies there. This situation is very different than, for example, the recoil convolution in the broadening factorization theorem [12,69]. In that case, the relevant modes were still only soft and collinear. Any factorization theorem of the double differential cross section must be super-SCET.
A possible complaint with the interpolation procedure 8 is that it is not unique and therefore there is no control over the logarithms that appear in the bulk of the phase space in the double differential cross section. This is an especially valid point because there is no factorization theorem in the bulk of the phase space and therefore no formal accuracy of the interpolation conjecture in this region is guaranteed. However, we will show that (under reasonable assumptions on the double differential cross section) to NLL accuracy, the boundary conditions are sufficiently robust to forbid all logarithms that are not generated by our procedure for interpolation up to O(α 4 s ) in the exponent of the double cumulative cross section. This is strong evidence that our interpolation procedure of setting scales can capture all logarithms that exist in the double differential cross section of two angularities to NLL accuracy over all of the phase space.
The outline of this paper is as follows. In Sec. 2 we discuss the relevant phase space for the double differential cross section in the two angularities e α and e β . This will also necessitate a discussion of the definition of the axis about which the angularities are measured. To remove sensitivity to recoil from soft wide angle emissions, we measure angularities about the broadening axis of a jet [20,55,61]. In Sec. 3 we compute the double differential cross section at fixed-order. This will illustrate some of the subtleties of resummation of the double differential cross section. In Sec. 4 we present the factorization theorem of the double differential cross section. We first discuss what can be learned simply from the phase space, then turn to the relevant SCET modes that contribute to the two angularities, and finally explicitly show that the double differential cross section factorizes near the boundaries of the phase space. Because the factorization theorem contains unfamiliar double differential jet and soft functions, we discuss the structure of these objects in Sec. 5 from constraints of power 6 For a single emission, this new logarithm reduces to the relative transverse momentum of the emitted parton. 7 It might seem that the case when α = 1 or β = 1 is special, where the logarithms degenerate, which may suggest that the number of modes that contribute to these cases is reduced. However, as in the case of traditional broadening, just because the contributions from different modes to the observable degenerate does not mean that the number of modes that contribute changes. As was observed with recoil-free angularities in Ref. [20], we expect that there is smooth behavior through α = 1 and β = 1. 8 We thank Jesse Thaler for extensive discussions of this point.
counting and consistency of the factorization. In Sec. 6 we suggest a simple procedure for interpolating the double differential cross section from the boundaries into the bulk of the phase space. We will show that this interpolating conjecture for the double differential cross section satisfies non-trivial consistency conditions and provide evidence that it captures all logarithms to NLL accuracy. In Sec. 7 we compare our expression for the double differential cross section to Monte Carlo simulation and find good qualitative agreement. Finally we close in Sec. 8 and suggest several future directions and applications for studying double differential cross sections.

Angularities Phase Space
We begin with a discussion of the phase space of the differential cross section of two angularities. From the introduction, we define the angularity e α measured on a narrow jet as where E J is the jet energy, R 0 1 is the jet radius, and the sum runs over all constituents in the jet. For IRC safety, α > 0. θ i is the angle between particle i and an appropriately chosen axis. Historically, this has been chosen to be the jet axis, defined as the sum of three-momenta of all particles in the jet. However, recently [55] it has been emphasized that this choice of axis is sensitive to recoil effects from the emission of soft, wide angle particles. At small values of the angular exponent α, the effect of recoil dominates the value of the angularity, significantly reducing its power for quark versus gluon jet discrimination, for example.
Instead, one can define an axis that is insensitive to these recoil effects and one example of this is the broadening axis [55,61]. 9 The broadening axis is defined as the axis in the jet that minimizes the β = 1 measure of N -subjettiness [53,54]; equivalently, the broadening axis is defined as the axis that minimizes the jet broadening [3,71,72]. That is, the broadening axisb corresponds to the axis that minimizes the scalar sum of momentum transverse to it: For a jet with two constituents, the broadening axis aligns with the hardest particle. In general, the broadening axis typically aligns with the direction of the hard core of energy in the jet. We also define the broadening axis to be the center of the jet so that all particles in the jet are closer than the jet radius R 0 to the broadening axis. With this set-up, now consider the allowed phase space of the double differential cross section of two angularities e α and e β . We will assume that α > β and so, because all angles 9 It should be noted that the broadening axis is one definition that results in recoil-free observables. Other recoil-free examples include energy correlation function observables [7,55,58] and the axis defined by the winner-take-all jet algorithm recombination scheme [20,70]. To the accuracy that we work in this paper, all of the recoil-free choices are identical. between particles and the broadening axis are less than R 0 , e α < e β . This implies that as e β → 0, then e α → 0. Also, because the angularities e α and e β are first non-zero at the same order in perturbation theory, then e α → 0 implies that e β → 0. Therefore, in addition to the upper bound on the phase space, there must also be a lower bound on the phase space for two angularities e α and e β . This lower bound of the phase space follows from energy conservation.
These properties can be seen explicitly in a jet with two constituents. The phase space can be described by the splitting angle θ and the energy fraction z of the emission. For the emission to be in the jet, θ < R 0 and for energy to be conserved z < 1. The matrix element then necessarily contains the restrictions In these phase space coordinates, in the soft emission limit, the recoil-free angularity e α is 10 which ranges from 0 to 1. The phase space coordinates z and θ can be rewritten in terms of the two angularities e α and e β as (2.5) The phase space restrictions written in terms of e α and e β are then where the first Θ-function follows from energy conservation and the second Θ-function follows from demanding that the emission is in the jet. The allowed phase space in the (e α , e β ) plane is illustrated in Fig. 2, setting α = 2 and scanning over a range of values for β. 11

Fixed-Order Cross Section
In this section, we will explicitly compute the double differential cross section of two jet angularities at O(α s ). The process we will consider is gluon emission from a quark and we will use the soft emission form of the angularities from Eq. (2.4). For simplicity, we will just use the QCD splitting function as representative of the matrix element, but this only differs 10 Strictly speaking, the recoil-free angularities to this order in αs are Throughout this paper, we will only consider logarithmically-enhanced contributions to the angularities. Therefore, to the accuracy that we consider, the definition of the angularities in Eq. (2.4) is sufficient. 11 This phase space has been discussed previously in Ref. [61]. The allowed phase space of the double differential cross section of angularities e α and e β . The angular exponent α is fixed to be 2 and β is varied. For a given value of β, the phase space consists of the respective shaded region and all shaded regions above.
from the full QCD matrix element at O(α s ) by non-singular terms. To the accuracy that we consider, the (normalized) cumulative distribution of two angularities can be computed from Σ(e α , e β ) = 1 + α s π where P q (z) is the quark splitting function given by The −1 in the first line is the subtraction of the virtual contribution which, by unitarity, we can assume is defined on the same phase space as the real contribution. On the physical phase space defined by e β > e α and e β α > e α β , we find 12 From the double cumulative cross section, the double differential cross section is found by differentiating with respect to e α and e β . Away from the boundaries of the phase space, we find The structures of the cumulative cross section and the differential cross section have some surprising distinctions. In the cumulative distribution, there are several terms which appear power-suppressed in the physical phase space region. For example, consider the term Because e α β < e β α in the physical phase space, this term is suppressed by powers of the angularities. Specifically, it is constant on the curve e β α = e α β , but otherwise vanishes in the physical phase space as e α , e β → 0. However, in the double differential cross section, this term produces Because e α < e β on the physical phase space, this term actually diverges as e α , e β → 0. Clearly, this term is integrable so one would not necessarily think that it needs to be resummed. This term, however, is actually vital to reproduce the single differential cross section of one angularity to single logarithmic accuracy. By marginalizing over one of the angularities, we have Note that the first term is evaluated at the upper limit of the phase space. This means that in this term, e α has been integrated over its entire physical range and so by itself, this term must be the differential cross section of e β . That is, The second term on the second line of Eq. (3.5) therefore must be zero to reproduce the correct cross section. This can be checked explicitly. The derivative of the cumulative distribution with respect to e β is For e α = e β , this produces which is correct to this accuracy. For e α = e α/β β , it indeed vanishes. However, there is a delicate cancelation of terms that is necessary for this term to vanish. Note that if the naïvely power-suppressed terms in the double cumulative cross section are removed, the derivative becomes (3.9) At the boundary where e α = e α/β β , we find which is clearly non-zero. Therefore, to guarantee that the double differential cross section is accurate and consistent to single logarithmic accuracy requires that the cumulative cross section contains terms that are naïvely power-suppressed with respect to logarithmic terms. This is unfamiliar from the calculation of resummed single differential cross sections because there is no analogous consistency condition and will be important in the following sections. Note that the double logarithms are correct, even when all power-suppressed terms are removed.

Factorization Theorem
Having discussed the phase space and fixed-order calculation of the double differential cross section, we now turn to studying its all-orders properties. We present the factorization theorem for the double differential cross section of angularities measured on a single jet. This section consists of three parts ordered in increasing technical detail, but only the first part is necessary to understand the remainder of this paper. First, we return to discussing the phase space of the double differential cross section. Nearly all of the conclusions from this section follow from simple geometric arguments about the limiting behavior of the double cumulative distribution at the boundaries of the phase space. We then discuss the relevant on-shell SCET modes which contribute to the angularities. We will show that there are two soft modes with different invariant mass which are relevant in the bulk of the phase space of the angularities. This will be an obstacle to factorization on the full phase space but, by an appropriate partition, we are able to prove factorization of the cross section at the boundaries of the phase space. The form of the factorization theorem will result in non-trivial identities between the cross sections at the two boundaries and we will use this in the following section to interpolate the cross section from the two boundaries into the bulk region.

A Study of the Phase Space
Consider again the allowed phase space for the two angularities e α and e β , with α > β. The double cumulative distribution Σ(e α , e β ) is the integral of the double differential cross section over a rectangle that includes the origin of the phase space. In particular, the double cumulative distribution can be evaluated at one of the boundaries of the phase space, illustrated in Fig. 3. For example, if the double cumulative distribution is evaluated at the boundary where e β α = e α β , note that e β has been integrated over its entire allowed range: from e β = e α to e β = e β/α α . Therefore, on this boundary, the cumulative distribution can only depend on e α : Σ(e α , e β = e β/α α ) = Σ(e α ) ,  To determine the differential cross section, we differentiate the double cumulative distribution. The boundary behavior of the cumulative distribution implies that which is the single differential cross section for e α . This can be related to the double differential cross section by integration: which holds for all values of e β/α α > 0. For this relationship to be true at this boundary, the double differential cross section should be expressable as 13 where the f + function integrates to zero on e β ∈ [0, e β/α α ]. A similar relationship holds for the other boundary, where e α = e β .
Then, the statement of the boundary factorization theorem is: the double differential cross section simplifies at the boundaries: (4.6) 13 All that is necessary is that the function that multiplies the differential cross section of eα integrates to 1 and the remainder function integrates to 0. Using distributions, a function that integrates to 1 can always be expressed as an appropriate δ-function plus a distribution that integrates to 0. Therefore, the expression in Eq. (4.5) is not unique, but will be justisfied with the factorization theorem in the following sections.
We have made use of the fact that the arguments of the non-trivial functions assume a very specific form dictated by the factorization, as discussed in detail in Sec. 5. Factorization theorems for differential cross sections of individual angularities are well-known [6,10,12,20]; therefore, the double differential cross section factorizes at the boundaries. In the following sections, we argue for factorization by studying the on-shell modes of the double differential cross section in detail. This relationship between the single and double differential cross section is quite remarkable, and can be understood as a precise statement of the UV structure of the effective theory. The e α = 0 and e β = 0 lines are where all UV divergences are localized. This is consistent with the fact that these lines are parametrically far away from the boundary where the factorization theorem is valid, since this is a statement about the UV structure of the theory.

Modes of the Double Differential Cross Section
We now turn to studying the on-shell SCET modes that contribute to the double differential cross section. For small values of a jet angularity e α , the dominant contributions to e α come from collinear and soft radiation in the jet. In general, the contribution to e α from collinear modes scales like θ α , where θ is the characteristic angular size of the collinear splittings. Soft modes, by contrast, contribute an amount that scales like their energy. Therefore, for the soft and collinear modes to be on-shell and contribute comparably to the angularity e α , their momenta must scale like 14 in the −, + and ⊥ lightcone coordinates, respectively. λ is a small parameter which sets the size of the angularity; here λ ∼ e 1/α α . For the double differential cross section, this analysis can be extended to the two angularities, e α and e β . We will only consider on-shell modes, which for small values of e α and e β are only soft and collinear radiation. For on-shell collinear modes, the scaling of their momenta must be the same as for a single observable, from Eq. (4.7). Because the angular scaling of the angularities e α and e β is different, the collinear modes contribute an amount of order λ α to e α and λ β to e β . Soft modes are more subtle. Now, because there are two angularities, there are two possible scalings of the soft modes. Either the momenta of the soft modes scale like λ α or they scale like λ β . Any other scaling would either be off-shell or would not be consistent with the collinear modes. Therefore, while there is a single collinear mode that contributes to the double differential cross section of e α and e β , there are two soft modes whose scalings are set by the angular exponents of the angularities. Table 1: Scaling of the on-shell collinear (C) and soft (S α , S β ) modes of the double differential cross section. z is the energy fraction of the mode and θ is the angle of the mode from the jet axis. This is shown in Table 1 where the scaling in the small parameter λ of the collinear and soft modes is given in terms of their energy fraction z and their splitting angle θ. Also, we show the contribution to the two angularities from each mode. The collinear mode contributes a different amount to each angularity, depending on the angular exponent. By contrast, each soft mode contributes the same amount to the two angularities, because the angularities are linear in the energy of the modes. Thus, in addition to having to deal with two soft modes on a single jet, the scaling of the angularities will be unfamiliar from the single differential cross section.
We will prove in the following section that with this scaling of the modes, the double differential cross section does factorize. Here, we will give a heuristic argument for the factorization of the cross section. If we assume that α > β, we can determine the dominant modes that contribute at leading power in λ to the cross section. For now, we will assume that the cross section can be written in the factorized form: where σ 0 is the Born-level cross section, H is the hard function, J(e α , e β ) is the jet function describing the collinear modes' contribution and S(e α , e β ) is the soft function describing the soft modes' contribution. The ⊗ symbol denotes convolution. This form of the cross section is suggestive, but must be expanded in powers of λ to ensure that the divergences in the hard, jet and soft functions cancel consistently at leading power in λ. This expansion can be done depending on the chosen scaling of the soft modes. By choosing a particular scaling of the soft modes, we restrict ourselves to a small region of the full angularities phase space, described in Sec. 2, where those soft modes are on-shell. If we first choose the S α soft modes, then the soft and collinear contributions to the angularity e α both scale like λ α . Therefore, they both will appear in the leading-power cross section. However, for this choice of soft mode scaling, the contribution from soft modes to e β , which scale like λ α , is power-suppressed with respect to the contribution from collinear modes, which scale like λ β . Explicitly, this is the limit in which e α e β , corresponding to a region of phase space far from the boundary e α = e β . Therefore with this choice of scaling of the soft modes, the leading-power factorized cross section has the form: where the superscript α denotes the scaling of the soft modes. Note that both angularities appear in the jet function and so the angle of the splitting is dominating the double differential cross section. Thus, this form of the cross section is valid in the region of phase space controlled by collinear emissions, near the boundary where e β α = e α β . By similar arguments, choosing the other scaling of the soft modes produces the factorized cross section which corresponds to an expansion with e β e β/α α , which is far from the boundary where e α β = e β α . Therefore, the factorization theorem of Eq. (4.10) is valid near the boundary dominated by soft emissions, where e α = e β .
Therefore, the double differential cross section factorizes near the boundaries of the phase space. The form of the factorization is quite interesting. Near the α boundary, consistency of the renormalization group implies that where γ F denotes the anomalous dimensions of the appropriate function F . The hard function has no dependence on the observable and the soft anomalous dimension only depends on the angularity e α . Thus, near this boundary of phase space the anomalous dimension of the jet function can only have non-trivial dependence on e α . Now we are in a position to see that the analytic forms of the double differential factorization theorems given in Eq. (4.6) capture the UV structure of each factorization. The δ-function term (which multiplies the single differential cross-section) contains all of the divergences of the factorization, and hence dictates its canonical resummation. It must be the case that the UV-divergence structure is localized by these δ-functions, since each factorization contains a single differential function, and between all sectors divergences cancel. Put simply, on a boundary the UV structure of the factorization reduces to that of the single differential cross-section.

Proof of Boundary Factorization Theorem
We now present a proof in SCET that the double differential cross section factorizes near the boundaries of the phase space. In this proof, we will implicitly use many results from Ref. [10] which discussed the factorization of jet observables for the first time, and so, here, will only focus on the novel aspects of the factorization theorem of the double differential cross section. Also, our analysis will focus on jets in e + e − → qq events, but using Ref. [10], this can be extended to jets in e + e − collision events with any number of well-separated jets.
We begin with the double differential cross section in QCD for e + e − → qq: where L µν is the leptonic tensor and J µ (x) is the QCD current at position x. The two δ-functions enforce the measured values of the angularities and the operator O J is the jet algorithm restriction. Here, we will mostly be agnostic to the form and function of this operator. It is defined to return a jet in the event on which the angularities are measured.
To the order to which we work, the jet algorithm in e + e − → qq events can be enforced by integrating over one hemisphere of the event and boosting to constrain radiation of the other hemisphere to exist in a cone of radius R 0 . This setup will be sufficient for our discussion here, with a more detailed discussion of jet algorithm factorization left to Ref. [10]. As mentioned in the introduction, we will only discuss the resummation of global logarithms. Factorization-violating non-global and clustering logarithms will be left to future work.
To be able to factorize the cross section, we match the QCD current with the corresponding current written in terms of fields in SCET as is a light-like Wilson line, and C is the matching coefficient matrix. Spinor indices have been suppressed for simplicity. However, matching SCET to QCD is more subtle than for single differential cross sections. For the case of the double differential cross section of two angularities, the form of the factorization depends not only on the fact that the angularities are small, but also the way in which they scale with respect to one another. This is important because the current matching does not set the virtuality of the soft radiation. For a single differential cross section, the measurement of the observable sets the virtuality of the soft emission, but for the double differential cross section, the soft radiation does not have a unique, well-defined virtuality, as discussed in Sec. 4.2. Only once the relative scaling of the angularities e α and e β is specified does the soft radiation have a well-defined virtuality.
To enforce a virtuality of the soft modes, we can restrict the measurement operator to only have support in the region of phase space where the relative scaling of the angularities e α and e β produce a unique soft mode. From Sec. 4.2 we found two on-shell soft modes, and so to do this, we can partition the phase space into two regions: one in which the soft modes have virtuality λ 2α and the other in which the soft modes have the virtuality λ 2β . That is, the measurement operator can be written as where we have inserted the identity. κ controls the relative scaling of e α with respect to e β . On the physical phase space, κ ∈ [1, α/β]. At this level, Eq. (4.14) is an operator identity, however, each Θ-function constrains the angularities in a region of phase space with a unique, on-shell soft mode. This partitioning is illustrated in Fig. 4 where boundary region α corresponds to e κ β > e α and boundary region β corresponds to e α > e κ β . Inserting Eq. (4.14) into the expression for the full QCD cross section, we have Note that the Θ-functions commute with all operators because they are functions of pure numbers (i.e., the value of the angularities). Because each term after the second equal sign in Eq. (4.15) only has a single on-shell soft mode, each term separately can be factorized by matching currents as defined in Eq. (4.13). For single angularities measured with respect to the jet thrust axis, this was done in Ref. [10] and for angularities measured with respect to the jet broadening axis this was done in Ref. [20]. Because it is a relatively standard and familiar procedure, we do not present the details of the factorization of the QCD cross section into SCET operators.
Performing the factorization and expanding to leading power, we find the following form for the cross section for angularities e α and e β measured on a single jet where the hard function H is the absolute square of the matching coefficient matrix C, H = C † C. The subscript on the symbol ⊗ denotes the appropriate convolution. For example, in the first term, because the soft function is independent of e β , the jet and soft functions are only convolved in e α . The single differential functions are where the jet's − component of momentum is Q and we have assumed that the jet is in thē n direction. The double differential jet and soft functions are The calculation of the double differential jet and soft functions will be discussed in Sec. 5. The definitions of the various operators appearing in these functions can be found Ref. [12] and references therein. Thus, the double differential cross section factorizes in the boundary regions of phase space. By the arguments of the previous section the double differential cross section reduces to a single differential cross section of the appropriate angularity, depending on the boundary. Because the double differential cross section must be independent on the choice of partitioning defined by κ, this will provide powerful constraints on the double differential cross section and will allow us to define an interpolating function from the boundaries into the bulk of the phase space. This will be studied in detail in Sec. 6.

Limit of Soft-Collinear Factorization
We are now in a position to understand why only two factorization theorems can be written down for the double differential cross section, and using the traditional ingredients of softcollinear factorization, no universal factorization formula could be presented. 15 Two separate arguments apply, leading to this conclusion. First, the fact that there exist two distinct soft modes as defined in Table 1 implies that there is no unique singular fixed-order cross section. Rather, there are two different singular cross sections that depend on the scaling of the soft mode. No one soft mode covers all of phase space.
Alternatively, one can be wholly ignorant of the power-counting and still come to the same conclusion. Formally, the SCET Lagrangians for the soft and collinear sectors at leading power are equivalent to full QCD [73]. Thus one can forget about the relative power counting of the low-scale components of the momenta between the soft and jet modes, and simply write down all possible jet and soft functions that could contribute. As long as the number of jets in the process is fixed, and hence also both the hard function and the number of eikonal lines in the soft function, the set of jet and soft functions is finite, and is controlled only by how many angularity measurements are imposed on a sector. So for a two-jet process, the only on-shell functions that can be written down are the single and double differential jet and soft functions from Eqs. (4.17) and (4.18). Fixing the phase space fixes the form of the divergences, regardless of how one power-counts the modes in the sector relative to each other. 16 Given these functions, a simple one-loop calculation is sufficient to show which combinations are RG consistent with each other.
As can be seen from the results of App. A, the only RG consistent combinations are those in Eq. (4.16). In particular, using only soft or collinear modes, there is no sense to the factorization theorem Eq. (4.8), independent of any argument about power counting. This constitutes a remarkable test of the consistency and power of the effective theory approach: fixing the scaling of the soft modes and appropriately expanding the phase space according to the power counting automatically generates RG-consistent combinations of on-shell functions. The precise power counting must be taken seriously to have consistent factorization. Of course one must eventually consider the power counting to know where in the phase space a given factorization formula holds, and this then shows that there is no universal factorization formula using traditional ingredients of soft-collinear factorization.
The non-uniqueness of the low-scale theory has an important consequence new to multidifferential cross sections. Namely, there is no operator product expansion (OPE) from one region to the other that allows a tower of effective theories that one could construct that covers all of phase space. Thus no RG scheme can connect the different regions of phase space. This is in distinction to the single differential cross section, where the singular distribution is unique. Indeed, this uniqueness of the singular terms is what allows the various regions of the differential cross section to be connected by controlling the RG evolution of the sectors. Even when one is in the tail of the distribution, the factorization theorem is correctly reproducing a unique set of terms in the fixed-order cross section, so one only needs to add the non-singular terms in the cross section to achieve the full result.

Double Differential Jet and Soft Functions
The factorized form of the double differential cross section from Eq. (4.16) contains single as well as double differential jet and soft functions. Soft and jet functions for individual angularities measured on a jet have been computed in Refs. [10,20], but the double differential objects are new. As discussed in the previous section, the divergences of the double differential jet and soft functions can only have non-trivial dependence on one of the angularities, for consistency of the factorization. However, the finite terms will have non-trivial dependence on both angularities and these contributions are necessary for improved accuracy of the double differential cross section. Here, we use general arguments to determine the form of the double differential jet and soft functions to all orders. The explicit calculation of the jet and soft functions is presented in App. A.

Jet Function
Much of the structure of the double differential jet function can be determined by power counting and the form of the factorization theorem. From the power counting of the factorization theorem, the jet function must scale as where the angularities scale as e α ∼ λ α and e β ∼ λ β . In addition, for consistency of the factorization theorem, the divergences in the double differential jet function can only have non-trivial dependence on e α , of exactly the same form as the single differential jet function: where div denotes the divergent parts of the jet functions. These two observations imply that the jet function has the following general form to all orders: J(e α , e β ) = C(α s ) δ(e α )δ(e β ) + e α /e β . This last quality is critical so that all of the divergences can be localized at e β = 0 (as required by the factorization theorem) with the +-prescription [74]. For a function f with support on [0, b], where b > 0, the function can be expressed as The b superscript denotes that the +-distribution is defined on (0, b]. It has the property that it integrates to zero:

Soft Function
We apply similar arguments to the the double differential soft function, S(e α , e β ). The scaling of the soft function is different than the jet function, because it exists near the boundary where e α = e β ∼ λ β . Then, from the factorization theorem the soft function scales like For consistency of the factorization theorem, the divergences in the double differential soft function can only have non-trivial dependence on e β and must be of the same form as the single differential soft function: [S(e α , e β )] div = [S(e β )] div δ(e α ) . (5.8) As with the jet function, these observations imply that the soft function has the following general form to all orders: 17 to be consistent with the anomalous dimension. F L is a function that depends on the loop order but scales like λ 0 to all orders. For the scaling of the double differential soft function, the only such combination of e α and e β that scales like λ 0 is e α /e β .
The singularities of the soft function can be localized at e α = 0 by the +-prescription. As discussed with the jet function, because e α is defined on [0, e β ] in the soft function, the function F L can be written as a +-distribution: (5.10) The calculation of the double differential soft function at one-loop is presented in App. A.2.

Interpolating between boundary regions
With the boundary factorization theorem, we would like to determine the double differential cross section throughout the allowed phase space for the two angularities. Because the factorization theorem only holds near the boundaries, we cannot claim any formal accuracy in the bulk of the phase space. Nevertheless, because the double differential cross section must satisfy several non-trivial constraints, these can be used to determine an interpolation from one boundary of the phase space to the other. In this section, we will present the interpolation to NLL accuracy in the boundary factorization theorem.
First, we will define what we mean by "NLL accuracy" for the double differential cross section. Typically, for a single observable e, NLL is defined to capture the leading terms in the exponent of the cumulative distribution with the scaling that α s log e ∼ 1. That is, NLL accuracy is log Σ NLL (e) ⊃ α n s log n+1 e, α n s log n e , (6.1) for all n > 0. For the double cumulative distribution of angularities e α and e β , we define NLL similarly, but include all possible logarithms of e α and e β : log Σ NLL (e α , e β ) ⊃ α n s log n+1−m e α log m e β , α n s log n−l e α log l e β , for all n > 0, 0 ≤ m ≤ n + 1 and 0 ≤ l ≤ n. This definition assumes that the logarithms of the double cumulative distribution exponentiates, which we believe is a reasonable expectation. 18 Also, as we measure the angularities on a jet, there will be non-global logarithms that arise at NLL; however, we will ignore them here. Now, we collect the constraints that were discussed in Sec. 4 on the double differential cross section of two angularities and its factorization theorem. With Σ(e α , e β ) the double cumulative distribution of the angularities e α and e β , it must reduce on the boundaries to: Σ(e α , e β = e β/α α ) = Σ(e α ) , Σ(e α = e β , e β ) = Σ(e β ) . (6. 3) The derivatives of the cumulative distribution are also constrained: The fact that these constraints are satisfied only for the total cross section implies that the factorization into soft and collinear modes cannot occur throughout the allowed phase space. The form of the boundary factorization theorem from Eq. (4.16) is This must be independent of κ for the factorization theorems at the two boundaries to be consistent with one another.
To determine a conjecture for the double differential cross section that interpolates between the boundaries of phase space subject to the above constraints, we will do the simplest thing possible. We will set the scales in the logarithms that appear in the boundary factorization theorem appropriately so that the total cross section constraints in Eqs. (6.3) and (6.4) are satisfied and the two boundary factorization theorems match onto one another continuously. However, because the RG evolution in the double differential cross section can only ever generate a non-zero value for one of the two angularities at the boundaries, this interpolation must be done at the level of the double cumulative cross section.
With this approach, we can then set scales in the logarithms of the double cumulative distribution on the boundary where e β = e β/α α (where it reduces to the cumulative distribution for e α alone) so that when continued to the boundary where e α = e β , it reduces to the cumulative distribution for e β and satisfies the other constraints. As we observed in the fixed-order calculation of Sec. 3, to satisfy the derivative constraints on the double cumulative distribution to single logarithmic accuracy required including naïvely power-suppressed terms in the cumulative distribution. Similar power-suppressed terms will need to be included in the resummed double cumulative distribution, in addition to setting scales, to satisfy all constraints. Because these power-suppressed terms are not exponentiated in the boundary factorization theorem, they are otherwise arbitrary and correspond to an uncertainty in the calculation.
To illustrate our procedure for interpolation, we will study in detail the double cumulative distribution to NLL accuracy. This will allow us to use known results for the cumulative distributions for individual recoil-free angularities at NLL. 19 Higher accuracy can be achieved by matching to fixed-order double differential cross sections, profiling the jet and soft scales in the factorization theorem [11] or resumming the individual angularities to higher logarithmic 19 At NLL, recoil-free angularities are identical to two-point energy correlation functions for the same value of the angular exponent β.
order. Here, we will only consider NLL order and will address improved accuracy in future work.

NLL Interpolation
At the boundaries of the phase space, the double cumulative distribution of two angularities e α and e β must reduce to the cumulative distribution of a single angularity, so we start by considering the form of the cumulative distribution for a single recoil-free angularity. To NLL accuracy, 20 the normalized cumulative distribution of a single recoil-free angularity e β measured on a jet can be expressed as [7,10] R(e β ) is often referred to as the radiator and consists of the cusp pieces of the anomalous dimensions of the jet and soft function and to NLL accuracy, is evaluated at two-loop order. The second term in this exponent, γ i T (e β ), is the non-cusp piece of the anomalous dimensions which result from hard collinear splittings. For NLL accuracy, it is evaluated at one-loop order. The prefactor accounts for the effects of multiple emissions adding together to produce a given value of the angularity e β . R (e β ) is the logarithmic derivative of the radiator: γ E is the Euler-Mascheroni constant. The explicit expression for the cumulative distribution at NLL is given in App. B. Because our strategy for achieving the interpolation is to only change the scale of the logarithms appearing in the single cumulative distribution, the normalized double cumulative distribution must be of the same functional form: Σ(e α , e β ) = e −γ ER (eα,e β ) Γ(1 +R(e α , e β )) e −R(eα,e β )−γ i T (eα,e β ) , (6.8) for some functions R(e α , e β ), T (e α , e β ) andR(e α , e β ). This then enforces the boundary conditions on the double cumulative distribution onto its constituent functions: R(e α , e β = e β/α α ) = R(e α ) , R(e α = e β , e β ) = R(e β ) , T (e α , e β = e β/α α ) = T (e α ) , T (e α = e β , e β ) = T (e β ) , R(e α , e β = e β/α α ) =R(e α ) ,R(e α = e β , e β ) =R(e β ) , up to terms suppressed by positive powers of e α , e β . In addition, the derivatives of each constituent function must satisfy the boundary conditions so as to correctly reproduce the differential cross sections of individual angularities at the boundaries. For example, for the derivative with respect to e α , we have the following derivative boundary conditions: Similar constraints exist for derivatives with respect to e β . With these results, we can consider each function separately and determine how it can be defined so as to interpolate between the boundary regions. As illustration of the interpolation, we will analyze the one-loop cusp component of the radiator R(e α , e β ) and the non-cusp function T (e α , e β ). The complete expression for the double cumulative distribution that satisfies all constraints is given in App. C. An important point to note is that, because they are defined by one-gluon emission, R(e α , e β ) and T (e α , e β ) can be directly computed in QCD. Here, we choose to compute them via the interpolation to illustrate the procedure. Also, we expect that the logarithmic structures generated by the interpolation are generic, and could be tested by computing anomalous dimensions at higher orders directly. On the other hand, the multiple emissions factorR(e α , e β ) can not be interpreted as the logarithmic derivative of the radiator R(e α , e β ) and so the method for computing it directly is not clear. However, its logarithmic structure can be determined by matching to the boundary conditions. This is an illustration of the power of the interpolation.

One-Loop Cusp/Radiator Interpolation
Consider first the one-loop radiator function for the angularity e α : where β 0 is the coefficient of the one-loop β-function. Equivalently, this can be written as an integral over the jet and soft function cusp anomalous dimensions: , (6.10) where β[α s ] is the β-function and the cusp anomalous dimensions of the jet and soft function to one-loop are µ is the renormalization scale and µ J and µ S are the jet and soft scales, which we take to be their canonical values: µ J = e 1/α α Q , µ S = e α Q , (6.12) where Q is the energy of the jet. Making this identification, we will refer to the terms in Eq. (6.9) with log e α as soft logarithms and those with log e 1/α α as collinear. 21 In this section, we will use the expression for the radiator in Eq. (6.9) because we are only working to oneloop order. However, the expression Eq. (6.10) is true to all orders, and so the interpolation obtained in this section could be tested at higher orders, given the cusp anomalous dimensions of the jet and soft functions at higher orders.
This expression in Eq. (6.9) is the one-loop component of the radiator R(e α , e β ) on the boundary e β = e β/α α . To this accuracy, we are free to change the argument of the soft and collinear logarithms by an order-1 number near this boundary of the phase space. The natural such number is e α β /e β α , which will enable the radiator to be continued into the bulk of the phase space, away from the boundary e β = e Note that soft logarithms on one boundary can mix and become soft or collinear logarithms on the other boundary of phase space (and similarly for collinear logarithms). Therefore, there are four possible terms that we must consider: log e α → log e β , log e α → log e β β , log e α α → log e β , log e α α → log e β β , (6.14) where the arrow indicates the interpolation from boundary e β = e β/α α to the boundary e α = e β . For example, consider the interpolation log e α → log e β . We multiply the argument of the soft logarithm on the e β = e The exponent c that satisfies this equation is c = 0. The three other logarithmic interpolations can be determined similarly. 21 Note that the anomalous dimensions are singular at α = 1. However, as shown explicitly in Ref. [20], the cross section is continuous through α = 1 and at α = 1 the relevant divergences transform into ultraviolet and rapidity divergences.
With this prescription for scale setting, the one-loop radiator is for some constants x, y. We have used the short-hand When e β = e β/α α , this reduces correctly to R (1) (e α ), and when e α = e β , only soft and collinear logarithms of e β are produced.
We now enforce the boundary conditions on R (1) (e α , e β ) to determine the constants x and y. When e α = e β , Eq. (6.16) becomes For this to reproduce R (1) (e β ), we must have To fix the remaining coefficient, we consider the derivative boundary conditions. Taking the derivative of the radiator with respect to e α and evaluating it on the boundary e β = e α it must vanish: (6.20) which then requires The other derivative boundary conditions produce the same constraints on x and y. It then follows that , (6.22) and so the radiator function at one-loop is which satisfies all boundary conditions. The form of this expression is interesting and we will discuss it in more detail in Sec. 6.2. For the radiator of a single angularity, there were only two logarithmic structures corresponding to soft or collinear logarithms. However, the interpolating radiator for two angularities has three logarithmic structures: soft (log e α ), collinear (log e 1/β β ) and what we will call "k T " logarithms: (6.24) We use the term k T because this combination of e α and e β reduces to k T = zθ/R 0 for one emission: Near the boundaries of the phase space the k T logarithms appropriately reduce to either soft or collinear logarithms. This would seem to suggest that to fully describe the bulk of the phase space requires introducing an additional mode into the effective theory. However, because there are only two types of singularities in QCD, we do not know how this would be done. The existence of a possible meta-effective theory that is well-defined over the entire phase space would be intriguing and deserves further study.

Non-Cusp Interpolation
As a second example of interpolation from the boundary into the bulk of the phase space, we will study the non-cusp piece, T (e α , e β ). To NLL accuracy, the non-cusp piece for a single angularity e β is T (e β ) = 1 πβ 0 log 1 + 2α s β 0 log e β β . Nevertheless, this expression does not satisfy the derivative boundary conditions. For example, the derivative with respect to e α vanishes, which satisfies the boundary condition when e β = e α . However, it clearly does not reproduce the correct term when e β = e β/α α so as to reproduce the differential cross section of e α . Other terms will need to be added to T (e α , e β ) to accomplish this. The terms that must be added cannot spoil the logarithmic accuracy of T (e α , e β ) and must produce the correct single logarithmic expressions when differentiated. Therefore, we must add a term to T (e α , e β ) that is power suppressed, but when differentiated produces singular terms. This was anticipated in Sec. 3 where it was observed that naïvely powersuppressed terms in the cumulative cross section were necessary to reproduce the correct single logarithms of the differential cross section. Motivated by the expressions there, we add to T (e α , e β ) a term that is suppressed by powers of e α and e β : T (e α , e β ) = 1 πβ 0 log 1 + 2α s β 0 log e β β + 2 α s π c e a α e b β β + 2α s β 0 log e β , (6.28) for exponents a, b and coefficient c. For the added term to be truly power suppressed in the small e α , e β limit, we require a + b > 0. The derivative boundary conditions will constrain a, b, c further. Taking the derivative with respect to e α , we find There are no constraints beyond these from taking the derivative with respect to e β . The constraints on T (e α , e β ) do not fully specify the parameters a, b, c so we must impose an additional, arbitrary condition. This should be interpreted as an uncertainty in the calculation that can be formally corrected by matching to the fixed-order cross section. Note that in the cumulative distribution these power-suppressed terms are beyond NLL accuracy anyway, so are only required to satisfy the boundary conditions and not to obtain formal NLL accuracy. Here, we will fix the parameters by considering a + b = 1 , (6.32) but any positive value for the sum of a and b would work. One could also consider adding several power-suppressed terms to satisfy the boundary conditions. With our choice on the sum of a and b, the non-cusp piece becomes 33) which satisfies all of the boundary conditions to leading power at NLL accuracy. Using the procedures developed here, all other pieces of the double cumulative cross section can be determined that satisfy the boundary conditions. We present the full expression in App. C.

Mixing Structure of Collinear and Soft Logarithms
We now discuss in more detail the mixing of the collinear and soft logarithms found in the radiator interpolation. When performing the interpolation we allowed for the mixing of the soft and collinear logarithms on one boundary into either of the soft or collinear logarithms on the other boundary. This allows for the presence of four possible logarithmic structures in the radiator for the double cumulative distribution. However, the consistency conditions at the phase space boundaries enforced that only three appear: the soft, collinear and k T logarithms, discussed in the previous sections. The fourth possible logarithmic structure arising from the interpolation between a soft logarithm on the e β = e which combines collinear energy scaling and soft angle scaling. The form of the factorization theorems on the boundary guarantees that such a structure does not appear in the bulk of the phase space. The non-appearance of this logarithmic structure implies that as we transition from one boundary to the other, either the collinear or soft logarithms (but not both) split to become a sum of collinear and soft logarithms on the other boundary. In particular, consider the transition from the e β = e β/α α boundary through the bulk to the e α = e β boundary. From Eq. (6.23) for the one loop radiator, we see that the logarithms which reduce to soft logarithms near the e β = e β/α α boundary remain as soft logarithms in the bulk, and on the e α = e β boundary. However, the collinear logarithms near the e β = e β/α α boundary split into a sum of collinear and k T logarithms in the bulk, and the k T logarithms then reduce to soft logarithms near the e α = e β boundary. This can also be phrased in terms of the mixing of the anomalous dimensions of the jet and soft functions appearing in the factorization theorems on the different boundaries. Recall from Eq. (6.10) that the one loop radiator for a single angularity can be written in terms of integrals of the jet and soft function cusp anomalous dimensions: The mixing of the logarithms described in the previous paragraph is equivalent to the following mixing of the jet and soft function anomalous dimensions: In this expression, K β S (µ S,β ), K β J (µ J,β ) are the integrals of the cusp anomalous dimensions for the factorization theorem near the e α = e β boundary, evaluated at their canonical scales, and the K α are the similar integrals of the anomalous dimensions for the factorization theorem near the e β = e β/α α boundary, but evaluated at the appropriately modified scales, and z = β−α α(β−1) . From this expression, we can see that as we move from the region where e β ∼ e β/α α to e α ∼ e β , the anomalous dimension of the jet function near the e β = e β/α α boundary splits into two pieces, one of which contributes to the anomalous dimension of the soft function near the e α = e β boundary and the other to the anomalous dimension of the jet function. This mixing structure is illustrated in Fig. 5. Fig. 5 also makes it transparent how the form of the factorization theorems on the two boundaries dictates that there can only be mixing between the jet functions near the e β = e β/α α boundary and the soft functions near the e α = e β boundary, as only these depend on both angularities. This clarifies why the only new logarithms in the bulk are the k T logarithms, and not the logarithms of Eq. (6.35).
To summarize, the structure of the radiator logarithms found via the interpolation procedure gives rise to a single new logarithmic structure in the bulk of the phase space, the k T logarithms. Unlike the soft (log e α ) and collinear (log e 1 β β ) logarithms, which reduce to, respectively, the soft and collinear logarithms of the two different factorization theorems near the boundaries, the k T logarithms reduce to the soft logarithms near the e α = e β boundary and the collinear logarithms near the e β = e β/α α boundary. This further clarifies the impediment to writing down a factorization theorem valid in the entire bulk region.

Evidence for Uniqueness of Interpolation
While the interpolation between the boundary regions presented above satisfies all constraints on the cross section from Eqs. captures the logarithms to any formal accuracy in the bulk of the phase space and so matching the resummed double differential cross section to the fixed-order cross section would be meaningless. However, making some reasonable assumptions about the structure of the logarithms in the bulk of the phase space, we will argue that the boundary conditions on the cumulative cross section are sufficiently strong to enforce the uniqueness of the interpolation up to O(α 4 s ). To prove this, we assume that the logarithms in the bulk of the phase space exponentiate. Then, the true double cumulative cross section to logarithmic accuracy can be written as log Σ(e α , e β ) = log Σ int (e α , e β ) + where d 1n mj and d 2n mj are coefficients, independent of α s . We assume that the logarithms in f n cannot be rewritten in such a way that factors of log e α e β , log e α β e β α exist. That is, all dependence on these logarithms has been explicitly factored out in Eq. (6.40). Because we assume that the interpolation cross section Σ int (e α , e β ) satisfies all boundary conditions, the second term, corresponding to interpolation-violating contributions, must vanish when either e α = e β or e β α = e α β so that the double cumulative cross section reduces appropriately at the boundaries. In addition, the derivatives of this term must also vanish when evaluated at the boundaries of the phase space to correctly reproduce the single differential cross sections. These boundary conditions are automatically satisfied for the the sum over n in Eq. (6.40) to start at n = 4 and for the sum over i to range from i = 2 to i = n − 2. From Eq. (6.41), this shows that the lowest order at which interpolation-violating terms can arise is α 3 s in the exponent. However, the only possible interpolation-violating (IV) contribution at O(α 3 s ) is leading logarithmic, which has the form log Σ The existence of these logarithms in the double cumulative cross section could be checked explicitly. However, it seems very unlikely that such a term could exist because it is double logarithmic and so should be totally captured by the resummation of the double cumulative cross section presented in Ref. [61]. While the lack of existence of this term would not necessarily prove exponentiation, it would demonstrate that the interpolation for the double differential cross section is significantly robust. Therefore, assuming exponentiation of logarithms in the double cumulative cross section and one-loop running capturing all leading logarithms, the lowest order at which the interpolation-violating contributions can exist is O(α 4 s ) in the exponent. If exponentiation of the logarithms of the double cumulative cross section does not occur, then the bulk of the phase space would only be described at fixed-order. Thus, this is strong evidence that, at least to NLL accuracy, the interpolation captures the dominant logarithmic structure of the double cumulative cross section. We further conjecture that the interpolation presented in Sec. 6.1 correctly resums all logarithms to the accuracy of the single cumulative cross sections at the boundaries. Testing this requires at least an O(α 3 s ) calculation, which is well beyond the state-of-the-art for fixed-order distributions of jet observables.

Comparison to Monte Carlo
With analytic results for the double differential cross section to NLL accuracy as defined by interpolation between the boundaries of phase space, we present a numerical analysis and compare to Monte Carlo simulation. Because we are interested in comparing the logarithmic  Table 2: Table comparing location and height of peaks of the double differential cross sections from the analytic NLL interpolation and Pythia.
structure of the analytic and Monte Carlo double differential cross section of angularities e α and e β , we will plot it in the plane (log e β , log e α ), with α > β. In this plane, the upper and lower boundaries of phase space are straight lines with slopes equal to 1 and α/β, respectively. We will plot the double differential cross section weighted by the two angularities: The Sudakov double logarithms manifest themselves as a concave down paraboloid in the (log e β , log e α ) plane.
We generate e + e − → qq events simulated with Pythia 8.165 [75,76] at a center-of-mass energy of 1 TeV with hadronization turned off, two-loop running of α s , and α s (m Z ) = 0.118. 22 To analyze the jets, we cluster jets with the e + e − anti-k T algorithm [77] with FastJet 3.0.3 [78] with a fat jet radius R 0 = 1.5. We analyze only the hardest jet in the event, requiring that the cosine of the angle between the jet momentum axis and the initial hard parton be greater than 0.9. We only include particles that lie within an angle R 0 = 0.4 from the broadening axis of the hardest jet. The energy of the jets is required to be in the range of Q ∈ [450, 550] GeV. We then measure the recoil-free angularities for various values of the angular exponents α and β of the jets in the sample.
In Figs. 6 and 7 we plot the distributions, fixing one angularity to be thrust, e 2 , and scanning over the other angularity: β = 1.5, 1, 0.5, 0.2. In the NLL interpolation plots, Fig. 6, the double differential distribution has been set to zero at very small values corresponding to scales near the Landau pole of α s . While the scale of the contours in the corresponding NLL interpolation and Pythia plots differ by up to a factor of 2, there is good qualitative agreement between the distributions. Both exhibit a peak in the distribution in the bulk of phase space at approximately the same location. This suggests that the correlations between angularities with different angular exponents are well-modeled in Monte Carlo.
The comparison between the analytic and Monte Carlo results can be made more quantitative by comparing the location and height of the peak of the distribution. In Table 2, 22 The quarks that are produced are only u, d, or s, so mass effects should be minimal.  Figure 6: Plots of the double differential cross section defined from the analytic NLL interpolation of Sec. 6 measured on quark jets with one angularity fixed to be thrust (α = 2) and scanning over the other angularity: β = 1.5, 1, 0.5, 0.2. The energy of the jets is Q = 500 GeV and the jet radius is R 0 = 0.4. The dashed lines on the plot correspond to the expected phase space boundary.
we list the location and height ("Peak") of the peak in (log e α , log e β ) space for α = 2, β = 1.5, 1, 0.5, 0.2. There are several features that illustrate qualitative agreement including: • The location of the peak in log e 2 generally becomes more negative as β decreases.
• The location of the peak in log e β moves to less negative values as β decreases.
• The height of the peak is relatively large for β near α = 2 and β near 0 and smaller for intermediate values of β.
While this qualitative agreement is encouraging, an honest quantitative comparison between Monte Carlo and analytic results would require going to at least NLL accuracy. That is, we would include the contributions from low scale matrix elements convolved with the NLL resummation kernel.
One apparent distinction between the analytic result and Pythia is that the double differential cross section in Pythia vanishes in the region near the line e α = e β , while it does not in the NLL cross section. We attribute this difference to angular ordering/veto imposed in the Monte Carlo. The line e α = e β requires that all emissions contributing there are located at θ = R 0 , the edge of the jet. Such a configuration is exponentially suppressed in a Monte Carlo, but is allowed in our NLL expression for the double differential cross section.
However, one can show that at NLL in each individual factorization, the double differential cross section vanishes at the boundaries of the phase space. The argument is as follows. As noted above Eq. (4.5), the general form of the fixed order singular cross section for the e α ∼ e α β β factorization boundary is: where we have explicitly indicated we are taking the cross sections at fixed order (fo). Note that f α + encodes all non-trivial e β dependence, and is solely fixed by the jet function matrix element for this factorization, and we have made explicit the boundary Θ-function enforced by the phase space of the jet function. 23 If we canonically resum this distribution, the resulting cross section is: where U (e α ) is the resummation kernel for e α . For non-zero e β , only the second term contributes, thus as e α approaches the boundary of phase space c e α β β the cross section vanishes, since the limits of integration become squeezed to zero. Note that this argument does not depend on the particular order to which one has calculated the cross section, and thus is a robust prediction of the factorization theorem for the double-differential cross section. We leave a detailed analysis at NLL , including interpolation into the bulk of the phase space, to future work.

Conclusions
In this paper, we have used the double differential cross section of two angularities measured on a single jet as a case study for understanding the factorization properties of double differential cross sections. We have explicitly shown that the double differential cross section for two angularities factorizes near the boundaries of the phase space, where it reduces to the single differential cross section of one of the angularities. Indeed, we have also shown the impossibility of a factorization theorem valid in the entire phase space region using only soft and collinear modes.
We presented a conjecture for the NLL double differential cross section using an interpolation procedure, based on scale setting and the addition of subleading terms, between the two factorization theorems defined on the boundaries of phase space. This interpolation procedure has the interesting property of introducing what we termed k T logarithms in the bulk region of phase space. These logarithms reduce to soft or collinear logarithms on the boundaries of the phase space where the factorization theorem applies, but are required in the bulk to interpolate between the soft logarithms on one boundary and the collinear logarithms on the other. The conjectured double differential cross section is subject to numerous consistency constraints from the boundary factorization theorems, which guarantee that it is unique to logarithmic accuracy up to at least O(α 4 s ). The interpolation scale choices that we found at NLL could be tested at higher accuracy by computing the anomalous dimensions of the jet and soft functions to higher loop order. We compared our calculation for the double differential cross section of angularities with a parton shower Monte Carlo, and found qualitative agreement, evidence that Monte Carlos model the correlations between angularities well.
While we have only discussed the perturbative aspects of the double differential cross section, the effect of non-perturbative physics is also important. Because the recoil-free angularities are additive and there exists a factorization theorem, this suggests that nonperturbative corrections can be incorporated by some kind of shape function [79,80]. In Ref. [61], a shape function was assumed to exist for the double differential cross section of angularities and it qualitatively agreed with the hadronization corrections in Pythia Monte Carlo. Nevertheless, a rigorous definition of the non-perturbative corrections to the double differential cross section is vital for determining the effect of low energy physics on the correlations of angularities. Angularities are additive observables and so the shape function for the double differential cross section should be similar in form to the shape function for a single differential cross section. However, the phase space constraints can be deformed by the non-perturbative corrections, which could result in subtle, but important, effects on the differential cross section.

Future Directions
This paper presents the first step in a wider program with the goal of understanding the factorization and resummation properties of double differential cross sections of IRC safe ob-servables. We therefore conclude with a number of future directions and possible applications of these techniques.

Extension to Other Observables
Although this paper has focused specifically on the example of angularities, the conclusions and techniques should be applicable to the double differential cross sections of more phenomenologically relevant observables. In particular, the argument for the reduction of the double cumulative distribution to a single cumulative distribution on the boundary of phase space presented in Sec. 4.1 is geometric in nature and does not rely on the detailed form of the boundaries. The only requirement is that the boundary is described by a monotonically increasing function. We therefore believe this reduction to be generic. This simplifies the problem of proving factorization theorems for double differential cross sections to that of proving factorization theorems for single differential cross sections, several of which are already known.
Armed with factorization theorems near the boundaries of phase space one can attempt an interpolation procedure by shifting scales and adding subleading terms in the cumulative distribution, as was done explicitly for the case of angularities in Sec. 6. For the relevant case when the two observables define boundary factorization theorems with different scalings for the soft modes, this interpolation procedure necessarily introduces a new logarithmic structure in the bulk of the phase space, which reduces appropriately on the boundaries to either a soft or collinear logarithm. For the case of two angularities, this was the k T logarithm discussed in Sec. 6.1.1. Furthermore, with certain assumptions on the logarithmic structure, a proof similar to that given in Sec. 6.3 could be used to argue for the uniqueness of the interpolation procedure. We believe that if this procedure is indeed possible for a particular pair of observables, then it gives a strong conjecture for the NLL resummed double differential cross section. This further allows for the computation of the ratio observable through marginalization. It is also interesting to speculate on the existence of a super-SCET formalism allowing for the incorporation of the additional modes required in the bulk, however, we leave this to future study.
One observable of particular phenomenological interest is N -subjettiness, which merits a more detailed discussion due to the interesting structure of its phase space. The Nsubjettiness observable τ (β) N is defined as [53,54] where R 0 is the jet radius and the sums run over all particles in the jet. R n,i is the angle between axis n and particle i and β > 0 for IRC safety. The axes in the jet can be chosen in several ways; the most elegant being to choose the axes so as to minimize the value of τ (β) N . The ratio of τ 2 /τ 1 has proven very powerful for discriminating boosted W jets from massive QCD jets [24,42,52]. 24 Progress has been made in computing the distribution of τ 2 /τ 1 for signal jets by relating it to the event-wide thrust distribution in e + e − collisions [13]. Understanding the background distribution is a formidable challenge that has not been studied for arbitrary values of the ratio of jet mass to jet energy. Nevertheless, we suspect that a boundary factorization theorem exists for the double differential cross section of N -subjettiness observables τ 2 and τ 1 . Because τ 2 is defined about two axes in the jet while τ 1 is only defined about one axis, τ 2 < τ 1 , with no non-trivial lower bound on the phase space. That is, τ 2 can be zero and τ 1 be non-zero, which is different from the angularities considered in this paper. However, as we illustrate in Fig. 8, the double cumulative distribution should still reduce to a single cumulative distribution on the appropriate boundaries. For example, evaluating the double cumulative distribution on the boundary τ 2 = τ 1 should reduce to the cumulative distribution for τ 1 alone as τ 2 has been integrated over its entire range. From the arguments in Sec. 4.1 this then implies that the double differential cross section reduces near this boundary: up to terms that integrate to zero on τ 2 ∈ [0, τ 1 ]. A similar relationship holds near the boundary τ 1 = 1 where the cross section reduces as again, up to terms that integrate to zero on τ 1 ∈ [0, 1]. 25 Therefore, to prove that the double differential cross section of τ 2 and τ 1 factorizes on the boundaries of phase space only requires proving that the single differential cross sections factorize. τ 1 is just the jet angularities, for which there exists a factorization theorem and so the double differential cross section should factorize on the τ 2 = τ 1 boundary. While there is up to now no factorization theorem for the differential cross section of τ 2 for QCD jets, we expect that there is a factorization theorem for τ 2 when τ 1 = 1. On this boundary, when τ 1 = 1, the structure that dominates τ 1 is a hard, perturbative emission. This essentially defines the jet to have two, well-separated subjets. In this configuration, τ 2 is dominated by soft and collinear radiation about those subjets, suggesting that its differential cross section is factorizable. 26 A proof of the factorization of τ 2 and subsequently the calculation of the double differential cross section of τ 1 and τ 2 would provide deep insights into the power of N -subjettiness as a discrimination observable as well as substantial information about the structure of QCD jets. 25 We have assumed that the maximum value of τ1 is 1. More generally, for τ1 near its maximum value, the double differential cross section should reduce to the single differential cross section for τ2. 26 Because the jet has two hard subjets, the factorization theorem would probably follow from SCET+ of Ref. [81]. Right: Evaluated on the boundary τ 2 = τ 1 which reduces the double cumulative distribution to Σ(τ 1 ).

Exclusive PDFs
Our results may also have consequences for the recent program of fully unintegrated or fully exclusive parton distribution functions (PDFs) [82][83][84][85][86][87][88][89][90] that depend on all components of colliding parton momentum, not only the longitudinal component. As we have shown in the context of angularities, the measurement of multiple observables on a single parton defines a region of the allowed phase space that is determined by the scaling of the observables with respect to one another. Typically, analyses that resum the logarithms of the unintegrated PDFs are constrained to a particular region of phase space; in Ref. [90], they essentially take the beam broadening comparable to the square-root of the beam thrust. 27 Though the PDF is more exclusive and should contain more information about the colliding parton, much of that information is lost because one is forced into a small region of the phase space. By studying the boundaries of the PDF phase space, it may be possible to interpolate between the boundary regions, producing a description of the unintegrated PDFs throughout the phase space.

Monte Carlos
Beyond its purely theoretic applications, the double differential cross section of two angularities can be used to tune Monte Carlos. Typically, tuning involves adjusting the arbitrary parameters in a Monte Carlo so as to match the measured differential cross section of several observables that are sensitive to the parameters. Tuning is not a precise science and involves significant art to choose parameters consistently so as to match many different distributions. However, if instead Monte Carlos were tuned to joint differential cross sections of observables, correlations would be naturally incorporated. With theoretical input for the double differential cross sections of angularities, parameters in the Monte Carlo could be adjusted appropriately to correctly model the higher order perturbative effects and separately, nonperturbative physics. This tuning program would also require the measurement of the double differential cross sections from the experiments, something that has not yet been published in studies at the Large Hadron Collider. 28 Because of theoretical progress, its potential application for Monte Carlo tuning, and the information it contains regarding correlations of observables, we strongly advise the ATLAS and CMS experiments to provide the measurements of double differential cross sections of jet observables.
for a total jet momentum of (Q, l + , 0) in light-cone coordinates. In this frame, the double differential jet function of a quark jet is J (1) (e α , e β ) = g 2 µ 2 C F dl + 2π We have assumed that the jet radius R 0 is O(1) and so the jet algorithm constraint is trivial to leading power in λ 1. Evaluating this in d = 4 − 2 dimensions, we find  All x dependence (or equivalently e β ) can be treated with the +-prescription as defined in Eq. (5.4). To the lowest orders in the dimensional regularization parameter , all e β dependence can be expressed as 29 where the divergent term is These expressions are consistent with the general arguments of Sec. 5.1. 29 This expansion has been performed with the Mathematica package HypExp [92,93].

A.2 Soft Function
At one-loop, the double differential soft function of a jet in the process e + e − → qq has the following form: (A.10) The Θ-function is the constraint of the jet algorithm for the definition of the angularities from Eq. (1.1). In d = 4 − 2 dimensions, the result is With this expansion of e α , we can now expand the remaining dependence in +distributions of e β . With the MS prescription, the soft function is S (1) (e α , e β ) = S (1) (e α , e β ) div + S (1) (e α , e β ) fin , (A.14) where the divergent term is S (1) (e α , e β ) div = α s 2π 2 Θ(e β ) e β + (A. 15) and the finite term is S (1) (e α , e β ) fin = α s 2π where higher-order terms in have been ignored. This form of the soft function is in agreement with the general arguments made in Sec. 5.2.

B The Cumulative Distribution of a Single Angularity
To NLL order, the cumulative distribution of a recoil-free angularity e β can be expressed as Σ(e β ) = e −γ E R (e β ) Γ(1 + R (e β )) e −R(e β )−γ i T (e β ) . (B.1) R(e β ) is the radiator and consists of the cusp pieces of the jet and soft function anomalous dimensions. To NLL accuracy, the cusp anomalous dimensions are evaluated at two loop order and the radiator is Here, C i is the color of the jet, λ = 2α s β 0 log e β , β 0 and β 1 are the one-and two-loop β-functions: For NLL accuracy, this only needs to be evaluated at one-loop: R (e β ) NLL = C i πβ 0 1 β − 1 log 1 + 2α s β 0 log e β β − log (1 + 2α s β 0 log e β ) . (B.8)

C The Double Cumulative Distribution for Two Angularities
From Sec. 6.1, the ansatz of the form of the double cumulative cross section for angularities e α and e β to NLL is Σ(e α , e β ) = e −γ ER (eα,e β ) Γ(1 +R(e α , e β )) e −R(eα,e β )−γ i T (eα,e β ) . (C.1) The functions R(e α , e β ), T (e α , e β ) andR(e α , e β ) can be found by setting scales in the logarithms so as to satisfy the boundary conditions. We find, for the radiator R(e α , e β ) R(e α , e β ) = C i 2πα s β 2 (C.4) The power suppressed terms have been chosen so that the sum of the exponents of e α and e β is 1.